text:
# Day 18: Snailfish
mutable struct SnailfishDigit
    number::Int
    depth::Int
    magnitude::Int # multiplier for calculation of magnitude
end
const SnailfishNumber = Vector{SnailfishDigit}

function input2sfn(input)
    # Luckily, today's snailfish numbers are valid Julia code which can be parsed!
    input = eval(Meta.parse(input))
    out = SnailfishNumber()

    _append!(leaf::Int, depth, mag) = push!(out, SnailfishDigit(leaf, depth, mag))
    function _append!(tuple, depth, mag)
        _append!(first(tuple), depth + 1, 3 * mag)
        _append!(last(tuple), depth + 1, 2 * mag)
        return nothing
    end
    _append!(input, 0, 1)
    return out
end

load(::Day{18}, path) = input2sfn.(readlines(path))

# Part 1 – compute magnitude of sum of snailfish numbers
function check_explosion!(sn::SnailfishNumber)
    len = length(sn)
    for (i, s) in enumerate(sn)
        if s.depth > 4
            # Add left value of pair to left neighbor
            i > 1 && (sn[i - 1].number += s.number)
            # Add right value of pair to right neighbor
            i < len - 1 && (sn[i + 2].number += sn[i + 1].number)
            # Remove and replace by zero
            popat!(sn, i + 1)
            s.number = 0
            s.depth -= 1
            s.magnitude = div(s.magnitude, 3)
            return true
        end
    end
    return false
end

function check_split!(sn::SnailfishNumber)
    len = length(sn)
    for (i, s) in enumerate(sn)
        if s.number >= 10
            n, d, m = s.number, s.depth, s.magnitude
            s.number = div(n, 2, RoundDown)
            s.depth += 1
            s.magnitude *= 3
            insert!(sn, i + 1, SnailfishDigit(div(n, 2, RoundUp), d + 1, m * 2))
            return true
        end
    end
    return false
end

function increase_depth!(sn::SnailfishNumber, mag)
    for s in sn
        s.depth += 1
        s.magnitude *= mag
    end
    return sn
end

function add!(a::SnailfishNumber, b::SnailfishNumber)
    append!(increase_depth!(a, 3), increase_depth!(b, 2))
    while true
        check_explosion!(a) && continue
        check_split!(a) && continue
        break
    end
    return a
end

magnitude(sn::SnailfishNumber) = mapreduce(s -> s.number * s.magnitude, +, sn)

solve1(::Day{18}, data) = magnitude(reduce(add!, deepcopy(data)))

# Part 2 – largest magnitude of any sum of two different snailfish numbers
function solve2(::Day{18}, data)
    largest = 0
    len = length(data)
    for i in 1:(len - 1)
        for j in (i + 1):len
            mag = magnitude(add!(deepcopy(data[i]), deepcopy(data[j])))
            mag > largest && (largest = mag)
            mag = magnitude(add!(deepcopy(data[j]), deepcopy(data[i])))
            mag > largest && (largest = mag)
        end
    end
    return largest
end

testresult1(::Day{18}) = 4140
testresult2(::Day{18}) = 3993


text:
export read_xml_plist_string, read_xml_plist

using Dates
using Base64

function parse_obj(n::ElementNode)
    tag = nodename(n)

    if     "plist" == tag
        # Quick and dirty. We just try to get out data somehow
        # not caring too much about correctness of format.
        children = elements(n)
        if isempty(children)
            return nothing
        else
            parse_obj(first(children))
        end
    elseif "integer" == tag || "real" == tag
        Meta.parse(nodecontent(n))
    elseif "string" == tag
        nodecontent(n)
    elseif "true" == tag
        true
    elseif "false" == tag
        false
    elseif "date" == tag
        df = Dates.DateFormat("yyyy-mm-dd")
        # TODO: include time info
        date_string, time_string = split(nodecontent(n), 'T')
        return Date(date_string, df)
    elseif "data" == tag
        base64decode(nodecontent(n))
    elseif "array" == tag
        parse_array(n)
    elseif "dict" == tag
        parse_dict(n)
    end
end

function parse_array(parent::ElementNode)
    map(elements(parent)) do child
        parse_obj(child)
    end
end

function parse_dict(parent::ElementNode)
    dict =  Dict{Any, Any}()

    children = elements(parent)
    it = iterate(children)
    while it != nothing
        keynode, i = it
        if nodename(keynode) == "key"
            valuenode, i = iterate(children, i)
            dict[nodecontent(keynode)] = parse_obj(valuenode)
        else
            error("Expected XML node to be 'key' not '$(nodename(keynode))'")
        end
        it = iterate(children, i)
    end
    dict
end

function read_xml_plist_string(text::AbstractString)
    doc = parsexml(text)
    if hasroot(doc)
        parse_obj(root(doc))
    else
        nothing
    end
end

function read_xml_plist(stream::IO)
    text = read(stream, String)
    read_xml_plist_string(text)
end

function read_xml_plist(filename::AbstractString)
    open(read_xml_plist, filename)
end

ext: jl


text:
"""
```julia
gradient(O, vars::AbstractVector; simplify = true)
```

A helper function for computing the gradient of an expression with respect to
an array of variable expressions.
"""
function gradient(O, vars::AbstractVector; simplify = true)
    Num[expand_derivatives(Differential(v)(value(O)),simplify) for v in vars]
end

"""
```julia
jacobian(ops::AbstractVector, vars::AbstractVector; simplify = true)
```

A helper function for computing the Jacobian of an array of expressions with respect to
an array of variable expressions.
"""
function jacobian(ops::AbstractVector, vars::AbstractVector; simplify = true)
    Num[expand_derivatives(Differential(value(v))(value(O)),simplify) for O in ops, v in vars]
end

"""
```julia
sparsejacobian(ops::AbstractVector, vars::AbstractVector; simplify = true)
```

A helper function for computing the sparse Jacobian of an array of expressions with respect to
an array of variable expressions.
"""
function sparsejacobian(ops::AbstractVector, vars::AbstractVector; simplify = true)
    I = Int[]
    J = Int[]
    du = Num[]

    sp = jacobian_sparsity(ops, vars)
    I,J,_ = findnz(sp)

    exprs = Num[]

    for (i,j) in zip(I, J)
        push!(exprs, expand_derivatives(Differential(vars[j])(ops[i]), simplify))
    end
    sparse(I,J, exprs, length(ops), length(vars))
end

"""
```julia
jacobian_sparsity(ops::AbstractVector, vars::AbstractVector)
```

Return the sparsity pattern of the Jacobian of an array of expressions with respect to
an array of variable expressions.
"""
function jacobian_sparsity(du, u)
    du = map(value, du)
    u = map(value, u)
    dict = Dict(zip(u, 1:length(u)))

    i = Ref(1)
    I = Int[]
    J = Int[]

    # This rewriter notes down which u's appear in a
    # given du (whose index is stored in the `i` Ref)
    r = [@rule ~x::(x->haskey(dict, x)) => begin
        push!(I, i[])
        push!(J, dict[~x])
        nothing
    end] |> Rewriters.Chain |> Rewriters.Postwalk

    for ii = 1:length(du)
        i[] = ii
        r(du[ii])
    end

    sparse(I, J, true, length(du), length(u))
end

"""
    exprs_occur_in(exprs::Vector, expr)

Return an array of booleans `finds` where `finds[i]` is true if `exprs[i]` occurs in `expr`
false otherwise.
"""
function exprs_occur_in(exprs, expr)
    vec(jacobian_sparsity([expr], exprs))
end

"""
```julia
hessian(O, vars::AbstractVector; simplify = true)
```

A helper function for computing the Hessian of an expression with respect to
an array of variable expressions.
"""
function hessian(O, vars::AbstractVector; simplify = true)
    vars = map(value, vars)
    first_derivs = map(value, vec(jacobian([values(O)], vars, simplify=simplify)))
    n = length(vars)
    H = Array{Num, 2}(undef,(n, n))
    fill!(H, 0)
    for i=1:n
        for j=1:i
            H[j, i] = H[i, j] = expand_derivatives(Differential(vars[i])(first_derivs[j]))
        end
    end
    H
end

isidx(x) = x isa TermCombination

"""
```julia
hessian_sparsity(ops::AbstractVector, vars::AbstractVector)
```

Return the sparsity pattern of the Hessian of an array of expressions with respect to
an array of variable expressions.
"""
function hessian_sparsity end

let
    # we do this in a let block so that Revise works on the list of rules

    _scalar = one(TermCombination)

    linearity_propagator = [
          @rule +(~~xs) => reduce(+, filter(isidx, ~~xs), init=_scalar)
          @rule *(~~xs) => reduce(*, filter(isidx, ~~xs), init=_scalar)
          @rule (~f)(~x::(!isidx)) => _scalar
          @rule (~f)(~x::isidx) => if haslinearity_1(~f)
              combine_terms_1(linearity_1(~f), ~x)
          else
              error("Function of unknown linearity used: ", ~f)
          end
          @rule (^)(~x::isidx, ~y) => ~y isa Number && isone(~y) ? ~x : (~x) * (~x)
          @rule (~f)(~x, ~y) => begin
              if haslinearity_2(~f)
                  a = isidx(~x) ? ~x : _scalar
                  b = isidx(~y) ? ~y : _scalar
                  combine_terms_2(linearity_2(~f), a, b)
              else
                  error("Function of unknown linearity used: ", ~f)
              end
          end] |> Rewriters.Chain |> Rewriters.Postwalk |> Rewriters.Fixpoint

    global hessian_sparsity

    """
    ```julia
    hessian_sparsity(ops::AbstractVector, vars::AbstractVector)
    ```

    Return the sparsity pattern of the Hessian of an array of expressions with respect to
    an array of variable expressions.
    """
    function hessian_sparsity(f, u)
        @assert !(f isa AbstractArray)
        f = value(f)
        u = map(value, u)
        idx(i) = TermCombination(Set([Dict(i=>1)]))
        dict = Dict(u .=> idx.(1:length(u)))
        found = []
        f = Rewriters.Prewalk(x->haskey(dict, x) ? dict[x] : x)(f)
        _sparse(linearity_propagator(f), length(u))
    end
end

"""
```julia
islinear(ex, u)
```
Check if an expression is linear with respect to a list of variable expressions.
"""
function islinear(ex, u)
    isempty(hessian_sparsity(ex, u).nzval)
end

"""
```julia
sparsehessian(O, vars::AbstractVector; simplify = true)
```

A helper function for computing the sparse Hessian of an expression with respect to
an array of variable expressions.
"""
function sparsehessian(O, vars::AbstractVector; simplify = true)
    O = value(O)
    vars = map(value, vars)
    S = hessian_sparsity(O, vars)
    I, J, _ = findnz(S)
    exprs = Array{Num}(undef, length(I))
    fill!(exprs, 0)
    prev_j = 0
    d = nothing
    for (k, (i, j)) in enumerate(zip(I, J))
        j > i && continue
        if j != prev_j
            d = expand_derivatives(Differential(vars[j])(O), false)
        end
        expr = expand_derivatives(Differential(vars[i])(d), simplify)
        exprs[k] = expr
        prev_j = j
    end
    H = sparse(I, J, exprs, length(vars), length(vars))
    for (i, j) in zip(I, J)
        j > i && (H[i, j] = H[j, i])
    end
    return H
end

function toexpr(O::Term)
  if isa(O.op, Differential)
     return :(derivative($(toexpr(O.args[1])),$(toexpr(O.op.x))))
  elseif isa(O.op, Sym)
    isempty(O.args) && return O.op.name
    return Expr(:call, toexpr(O.op), toexpr.(O.args)...)
  end
  if O.op === (^)
      if length(O.args) > 1  && O.args[2] isa Number && O.args[2] < 0
          return Expr(:call, ^, Expr(:call, inv, toexpr(O.args[1])), -(O.args[2]))
      end
  end
  return Expr(:call, O.op, toexpr.(O.args)...)
end
toexpr(s::Sym) = nameof(s)
toexpr(s) = s

function toexpr(eq::Equation)
    Expr(:(=), toexpr(eq.lhs), toexpr(eq.rhs))
end

toexpr(eq::AbstractArray) = toexpr.(eq)
toexpr(x::Integer) = x
toexpr(x::AbstractFloat) = x
toexpr(x::Num) = toexpr(value(x))


text:
#!/usr/bin/env julia

using Pkg
using MPI
using ICI_simulateur: Parallelisation, City, Population, Simulateur, Save
using TimerOutputs

const to = TimerOutput() ;

######################## Calcul du temps reel ##################################
deb = time() ;
################################################################################

####################### Initialisation de la parallélisation ###################
@timeit to "Parallelisation" comm,rank = Parallelisation() ;
################################################################################

######################## Définition de la ville ################################
@timeit to "City" Ville = City() ;
################################################################################

################### Définition de la population ################################
@timeit to "Population" Individu = Population() ;
################################################################################

############################### Simulateur #####################################
@timeit to "Simulateur" Propagation = Simulateur(Ville,Individu,rank) ;
################################################################################

########################### Partie des résultats ###############################
@timeit to "Save" Save(Propagation,rank,comm) ;
################################################################################

if(ICI_simulateur.ParamSimu.MPI)
    MPI.Barrier(comm) ;
end

if(rank == 0)
    elapsed = time() - deb ;
    println("Temps reel : ", elapsed, " secondes.") ;
end

if(ICI_simulateur.ParamSimu.MPI)
    MPI.Finalize() ;
end

print_timer(to)


text:
module XLSXasJSON

using Printf, REPL
using JSON
using XLSX
using DataStructures

include("index.jl")
include("pointer.jl")
include("worksheet.jl")
include("workbook.jl")
include("writer.jl")

export JSONWorkbook, JSONWorksheet,
        hassheet, sheetnames,
        xlsxpath, dropnull

end # module


text:
function num_var(model::Particle)
    nq = model.nq
    nc = model.nc
    nb = 2
    nq + nc + nc + nc + nc * nb + nc + nc * nb
end

function indices_z(model::Particle)
    nq = model.nq
    nc = model.nc
    nb = 2

    q = collect(1:nq)
    γ = collect(nq .+ (1:nc))
    sγ = collect(nq + nc .+ (1:nc))
    ψ = collect(nq + nc + nc .+ (1:nc))
    b = collect(nq + nc + nc + nc .+ (1:(nc * nb)))
    sψ = collect(nq + nc + nc + nc + nc * nb .+ (1:nc))
    sb = collect(nq + nc + nc + nc + nc * nb + nc .+ (1:(nc * nb)))
    IndicesZ(q, γ, sγ, ψ, b, sψ, sb)
end

function indices_optimization(model::Particle)
    nq = model.nq
    nc = model.nc
    nb = 2
    nz = num_var(model)

    IndicesOptimization(
        nz,
        nz,
        [collect(nq .+ (1:nc)), collect(nq + nc .+ (1:nc))],
        [collect(nq .+ (1:nc)), collect(nq + nc .+ (1:nc))],
        [
            [collect([nq + nc + nc + 1, (nq + nc + nc + nc .+ (1:nb))...]), collect([nq + nc + nc + nc + nb + 1, (nq + nc + nc + nc + nb + nc .+ (1:nb))...])],
        ],
        [
            [collect([nq + nc + nc + 1, (nq + nc + nc + nc .+ (1:nb))...]), collect([nq + nc + nc + nc + nb + 1, (nq + nc + nc + nc + nb + nc .+ (1:nb))...])],
        ],
        collect(1:(nq + nc + nc * nb + nc)),
        collect(nq + nc + nc * nb + nc .+ (1:nc)),
        collect(nq + nc + nc * nb + nc + nc .+ (1:3)),
        [collect(nq + nc + nc * nb + nc + nc .+ (1:3))],
        collect(nq + nc + nc * nb + nc .+ (1:(nc + 3))))
end

Trajectory(model::Particle, T) = Trajectory(model, T, nc=model.nc, nb=2)
GradientTrajectory(model::Particle, T) = GradientTrajectory(model, T, nc=model.nc, nb=2)

function residual(model::Particle, mass_matrix, dynamics_bias, kinematics, kinematics_jacobians, z, θ, μ)

    # dimensions
    nq = model.nq
    nu = model.nu
    nw = model.nw
    nc = model.nc
    nf = length(friction_coefficients(model))
    nb = 2

    # unpack data
    q0 = θ[1:nq]
    q1 = θ[nq .+ (1:nq)]
    u1 = θ[2nq .+ (1:nu)]
    w1 = θ[2nq + nu .+ (1:nw)]
    friction_body_world = θ[2nq + nu + nw .+ (1:nf)]
    h = θ[2nq + nu + nw + nf .+ (1:1)]

    # unpack variables
    q2 = z[1:nq]
    γ1 = z[nq .+ (1:nc)]
    sγ1 = z[nq + nc .+ (1:nc)]
    ψ1 = z[nq + nc + nc .+ (1:nc)]
    b1 = z[nq + nc + nc + nc .+ (1:(nc * nb))]
    sψ1 = z[nq + nc + nc + nc + nc * nb .+ (1:nc)]
    sb1 = z[nq + nc + nc + nc + nc * nb + nc .+ (1:(nc * nb))]

    # tangential velocity
    v1 = (q2 - q1) / h[1]
    vT = v1[1:2]

    # contact forces
    J = contact_jacobian(model, q1)
    λ1 = transpose(J) * [b1; γ1]

    [
    dynamics(model, mass_matrix, dynamics_bias, h, q0, q1, u1, w1, λ1, q2);
    signed_distance(model, q2) - sγ1;
    vT - sb1;
    friction_body_world[1] * γ1[1] - ψ1[1];
    γ1 .* sγ1 .- μ[1];
    cone_product([ψ1; b1], [sψ1; sb1]) - [μ[1]; 0.0; 0.0];
    ]
end


text:
# Tests for generate.jl

@testset "Generate" begin
    lower = [1.0]
    upper = [2.0]
    pres_inds = [1 3]
    constraints = Constraints(Atom[], lower, upper, pres_inds)
    coords = [
        0.0 0.0 2.5;
        0.0 1.0 0.0;
        0.0 0.0 0.0;
    ]
    @test errorscore(constraints, coords) == 0.25

    # Test L-amino acids are not inverted
    atom_ca = Atom("CA", "ALA", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_co = Atom("C", "ALA", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_n = Atom("N", "ALA", 'A', 20, [0.0, 0.0, 0.0], "N")
    atom_cb = Atom("CB", "ALA", 'A', 20, [0.0, 0.0, 0.0], "C")
    atoms = [atom_ca, atom_co, atom_n, atom_cb]
    coords_l = [
        0.0 1.0 -1.0 -1.0;
        0.0 -1.0 -1.0 1.0;
        0.0 1.0 -1.0 1.0;
    ]
    coords_fixed = [
        0.0 1.0 -1.0 -1.0;
        0.0 -1.0 -1.0 1.0;
        0.0 1.0 -1.0 1.0;
    ]
    fixchirality!(coords_l, atoms)
    @test coords_l == coords_fixed
    # Test D-amino acids are inverted
    coords_d = [
        0.0 -1.0 1.0 -1.0;
        0.0 -1.0 -1.0 1.0;
        0.0 -1.0 1.0 1.0;
    ]
    coords_fixed = [
        0.0 -1.0 1.0 -1.0;
        0.0 -1.0 -1.0 1.0;
        0.0 1.0 -1.0 -1.0;
    ]
    fixchirality!(coords_d, atoms)
    @test coords_d == coords_fixed

    atom_one = Atom("CA", "ALA", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_two = Atom(mod_atom_info["atom_name"], mod_atom_info["res_name"], mod_atom_info["chain_id"], mod_atom_info["res_n"], [0.0, 0.0, 0.0], mod_atom_info["element"])
    atom_three = Atom("CA", "ALA", 'A', 21, [0.0, 0.0, 0.0], "C")
    atom_four = Atom(mod_atom_info["atom_name"], mod_atom_info["res_name"], mod_atom_info["chain_id"], mod_atom_info["res_n"], [0.0, 0.0, 0.0], mod_atom_info["element"])
    atoms = [atom_one, atom_two, atom_three, atom_four]
    coords_one = [
        1.0 10.0 5.0 11.0;
        9.0 4.0 8.0 2.0;
        7.0 3.0 6.0 12.0;
    ]
    coords_two = [
        9.0 6.0 11.0 3.0;
        1.0 2.0 5.0 12.0;
        7.0 8.0 10.0 4.0;
    ]
    struc_one = ModelledStructure(10.0, coords_one)
    struc_two = ModelledStructure(20.0, coords_two)
    strucs = [struc_one, struc_two]
    ensemble = ModelledEnsemble(atoms, strucs)
    removemodulator!(ensemble)
    test_atoms = ensemble.atoms
    @test length(test_atoms) == 2
    @test test_atoms[1].atom_name == "CA"
    @test test_atoms[2].atom_name == "CA"
    test_strucs = ensemble.strucs
    @test strucs[1].score == 10.0
    real_coords_one = [
        1.0 5.0;
        9.0 8.0;
        7.0 6.0;
    ]
    real_coords_two = [
        9.0 11.0;
        1.0 5.0;
        7.0 10.0;
    ]
    @test strucs[1].coords == real_coords_one
    @test strucs[2].coords == real_coords_two

    struc_one = ModelledStructure(40.0, [0.0 0.0; 0.0 0.0; 0.0 0.0])
    struc_two = ModelledStructure(20.0, [0.0 0.0; 0.0 0.0; 0.0 0.0])
    struc_three = ModelledStructure(10.0, [0.0 0.0; 0.0 0.0; 0.0 0.0])
    struc_four = ModelledStructure(30.0, [0.0 0.0; 0.0 0.0; 0.0 0.0])
    strucs = [struc_one, struc_two, struc_three, struc_four]
    ensemble = ModelledEnsemble([], strucs)
    trimbyscore!(ensemble, 2)
    trim_strucs = ensemble.strucs
    @test length(trim_strucs) == 2
    @test trim_strucs[1].score == 10.0
    @test trim_strucs[2].score == 20.0
end


text:
using Nemo
using EllipticCurves
using ClassPolynomials

using Test

######################################################################
# Testing tools.jl
######################################################################

#This file should be ok if the tests below pass

######################################################################
# Testing weierstrass.jl
######################################################################

print("Testing weierstrass.jl... ")

E1 = Weierstrass(QQ(0), QQ(0), QQ(0), QQ(4), QQ(5))
E2 = ShortWeierstrass(QQ(0), QQ(0), QQ(0), QQ(4), QQ(5))
E3 = SeparatedWeierstrass(QQ(0), QQ(0), QQ(0), QQ(4), QQ(5))
E4 = EllipticCurve(QQ(4), QQ(5))

@test base_ring(E1) == QQ
@test base_ring(E2) == QQ
@test base_ring(E3) == QQ
@test E1 == E2 == E3 == E4

E = EllipticCurve(QQ(0), QQ(17))
@test a_invariants(E) == (QQ(0), QQ(0), QQ(0), QQ(0), QQ(17))
@test b_invariants(E) == (QQ(0), QQ(0), QQ(68), QQ(0))
@test c_invariants(E) == (QQ(0), QQ(-14688))
@test discriminant(E) == QQ(-124848)
@test j_invariant(E) == QQ(0)
@test isvalid(E)

eq = equation(E)
Y = gen(parent(eq))
X = gen(base_ring(eq))
@test eq == Y^2 - X^3 - 17

preq = projective_equation(E)
X, Y, Z = gens(parent(preq))
@test preq == Z * Y^2 - X^3 - 17 * Z^3

f3 = divisionpolynomial(E, 3)
f9 = divisionpolynomial(E, 9)
@test f9 % f3 == 0

print("done\n")

######################################################################
# Testing points.jl
######################################################################

print("Testing points.jl... ")

E = EllipticCurve(QQ(0), QQ(17))
P1 = Point(QQ(-2), QQ(3), QQ(1), E)
P2 = Point(QQ(-4), QQ(6), QQ(2), E)
P3 = Point(QQ(2), QQ(5), QQ(1), E)

@test coordinates(P1) == (-2, 3, 1)
@test base_curve(P1) == E
@test base_ring(P1) == QQ
@test isvalid(P1)
@test isvalid(P3)

@test samefields(normalized(P2), P1)
@test P1 == P2

print("done\n")

######################################################################
# Testing weierstrasspoints.jl
######################################################################

print("Testing weierstrasspoints.jl... ")

#Example is taken from Silverman 1 p.59.

E = EllipticCurve(QQ(0), QQ(17))
P1 = Point(QQ(-2), QQ(3), QQ(1), E)
P2 = Point(QQ(-1), QQ(4), QQ(1), E)
P3 = Point(QQ(2), QQ(5), QQ(1), E)
Pinf = zero(E)

@test P1 + P1 == 2 * P1 == Point(QQ(8), QQ(-23), QQ(1), E)
@test Pinf + Pinf == 2 * Pinf == Pinf
@test P1 + Pinf == Pinf + P1 == P1
@test P2 + P2 == 2 * P2 == Point(QQ(137//64), QQ(-2651//512), QQ(1), E)
@test P2 + P3 == Point(QQ(-8//9), QQ(-109//27), QQ(1), E)
@test -P3 == Point(QQ(2), QQ(-5), QQ(1), E)
@test P3 - P3 == 0 * P3 == Pinf
@test P1 - P3 == 1 * P1 + (-1) * P3 == Point(QQ(4), QQ(9), QQ(1), E)

@test projective_add(P2, P3) == P2 + P3
@test projective_scalar_mul(P2, 7) == 7 * P2

E = EllipticCurve(QQ(2), QQ(1))
P = Point(QQ(1), QQ(2), QQ(1), E)
@assert projective_add(P, P) == 2 * P
@assert projective_scalar_mul(P, 17) == 17 * P
print("done\n")

######################################################################
# Testing montgomery.jl
######################################################################

#Really trivial code, no need for testing

######################################################################
# Testing montgomerypoints.jl and curves
######################################################################

print("Testing montgomerypoints.jl... ")

E = Montgomery(QQ(7), QQ(1))
@test base_ring(E) == QQ
@test j_invariant(E) == 24918016//45
@test isvalid(E)

P1 = Point(QQ(1), QQ(3), QQ(1), E)
P0 = Point(QQ(0), QQ(0), QQ(1), E)
@test isvalid(P1) & isvalid(P0)

xP1 = XZPoint(P1)
xP0 = XZPoint(P0)
xinf = XZzero(E)

@test isfixedtorsion(xP0)
@test !isfixedtorsion(xP1)
@test fixedtorsion(E) == xP0

@test xdouble(xP0) == xinf
@test xdouble(xP1) == xP0
@test xdouble(xinf) == xinf

@test xadd(xP1, xP0, xP1) == xP1

#Scalar multiplications

@test 0 * xP1 == xinf
@test 1 * xP1 == xP1
@test 2 * xP1 == xP0
@test 1 * xP0 == xP0
@test 2 * xP0 == xinf
@test 45 * xinf == xinf
@test 4 * xP1 == xinf

@test 8 * xP1 == xinf

print("done\n")

######################################################################
# Testing edwards.jl
######################################################################

#Do we really want to keep this code ?

######################################################################
# Testing edwardspoints.jl
######################################################################

#Do we really want to keep this code ?

######################################################################
# Testing pairings.jl
######################################################################

#Do we really want to keep this code ?

######################################################################
# Testing isogenies.jl
######################################################################

print("Testing isogenies.jl... ")

F, x = FiniteField(101, 1, "x")
E = EllipticCurve(F(1), F(3))
P = Point(F(3), F(29), F(1), E)
Q = Point(F(4), F(24), F(1), E) #another point

#P is a point of three-torsion
@test isvalid(P)
@test 3 * P == zero(E)
@test 3 * Q != zero(E)

f3 = divisionpolynomial(E, 3)
@test f3(P.X) == 0

phi = Isogeny(E, P, 3)
Eprime = image(phi)
@test domain(phi) == E
@test Eprime == EllipticCurve(F(24), F(24)) #given by Sage
@test phi(P) == zero(Eprime)
@test isvalid(phi(Q))
@test phi(Q) == Point(F(91), F(20), F(1), Eprime)

K = subgroup(P, 3)

@test kernel(phi) == K
@test Isogeny(E, K) == phi
@test Isogeny(E, 3, Eprime) == phi

@test isvalid(phi)
psi = dual(phi)

@test compose(phi, psi) == multiplication_isogeny(E, 3)
@test multiplication_isogeny(E, 3)(P) == 3 * P

id = Isogeny(E, 1, E)
@test kernel(id) == 1
@test id(P) == P

@test isvalid(psi)
@test isvalid(id)

phix, phiy = rational_fractions(phi)
@test Isogeny(E, Eprime, phix, phiy) == phi

jprime = j_invariant(Eprime)
@test Isogeny(E, 3, jprime) == phi

#With modular equations

ell = 5
FX, X = PolynomialRing(F, "X")
bb, j1 = any_root(ClassicalModularPolynomial(ell, j_invariant(E), X))
@test bb
@test isvalid(Isogeny(E, ell, j1))

@test isvalid(Isogeny(E, ell, j1; USE_ATKIN = true))
print("done\n")

######################################################################
# Testing finfields.jl
######################################################################

print("Testing finfields.jl... ")

#Conversions

F, u = FiniteField(233, 1, "u")
G, z = FiniteField(233, 3, "z")
@test convert(G(1), F) == F(1)
@test convert(F(1), G) == G(1)

#Roots of polynomials

R, X = PolynomialRing(F, "X")
(bool, root) = any_root(X^2 - 1)
@test bool
@test root^2 == 1

#Cardinality

E = EllipticCurve(F(3), F(4))
#@test cardinality(E) == 260: no point counting algorithm yet
#@test frobeniustrace(E) == -26
E2 = Montgomery(F(3), F(1))
#@test cardinality(E2) == 240

#Random points

P = rand(E)
@test isvalid(P)
@test 260 * P == zero(E)
P2 = randXZ(E2)
@test isvalid(P2)
@test 240 * P2 == XZzero(E2)

#Torsion points

Q = torsionpoint(E, 13, 260)
@test Q !== zero(E)
@test 13 * Q == zero(E)
Q2 = torsionXZ(E2, 5, 240)
@test Q2 !== XZzero(E2)
@test 5 * Q2 == XZzero(E2)

#Base extensions

P = rand(base_extend(E, G))
@test isvalid(P)

@test card_over_extension(210, 233, 3) == 12652290

#Montgomery models
E = EllipticCurve(F(204), F(56))
(bool, Eprime) = has_montgomery(E)
@test bool
@test j_invariant(E) == j_invariant(Eprime)

print("done\n")

print("\n")


text:
module TimeseriesSurrogates

using Random
using Distributions
using Distances # Will be used by the LombScargle method
using StatsBase
using AbstractFFTs
using DSP
using Interpolations
using Wavelets
using Requires

include("api.jl")

include("utils/testsystems.jl")

# The different surrogate routines
include("methods/randomshuffle.jl")
include("methods/large_shuffle.jl")
include("methods/randomfourier.jl")
include("methods/aaft.jl")
include("methods/iaaft.jl")
include("methods/truncated_fourier.jl")
include("methods/partial_randomization.jl")
include("methods/wavelet_based.jl")
include("methods/pseudoperiodic.jl")
include("methods/pseudoperiodic_twin.jl")
include("methods/multidimensional.jl")
include("methods/ar.jl")
include("methods/trend_based.jl")

# Methods for irregular time series
include("methods/lombscargle.jl")

# Visualization routine for time series + surrogate + periodogram/acf/histogram
using Requires
function __init__()
    @require UncertainData="dcd9ba68-c27b-5cea-ae21-829cd07325bf" begin
        include("utils/uncertaindatasets.jl")
    end
    @require Makie="ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" begin
        include("plotting/surrogate_plot.jl")
    end
end

export surroplot_path
function surroplot_path()
    joinpath(@__DIR__, "plotting", "surrogate_plot.jl")
end

end # module


text:
"""
    x,y = squaredalitz(k,σs,ms)

    Calculate normalized values in square coordinates, -1 ≤ x,y ≤ 1
"""
function squaredalitz(k,σs,ms::ThreeBodyMasses)
    cθ = cosθij(k,σs,ms^2)
    cθ = cθ ≥ 1.0 ? 1.0 : (cθ ≤ -1.0 ? -1.0 : cθ)
    y = acos(cθ) / π
    mlimsk = sqrt.(lims(k, ms))
    xn = 2*(sqrt(σs[k]) - mlimsk[1]) / (mlimsk[2]-mlimsk[1])-1
    x = acos(xn) / π
    return (x, y)
end

squaredalitz1(σs,ms) = squaredalitz(1,σs,ms)
squaredalitz2(σs,ms) = squaredalitz(2,σs,ms)
squaredalitz3(σs,ms) = squaredalitz(3,σs,ms)

"""
    invsquaredalitz(k,x,y,ms)

    Inverse transformation from square coordinates
"""
function invsquaredalitz(k,x,y,ms::ThreeBodyMasses)
    mlimsk = sqrt.(lims(k, ms))
    #
    _mk = mlimsk[1] + (mlimsk[2]-mlimsk[1]) * (cos(x*π) + 1) / 2
    _σk = _mk^2
    _cosθij = cos(y*π)
    #
    k == 1 && return Invariants(ms, σ1 = _σk,  σ3 = σ3of1(_cosθij,_σk,ms^2))
    k == 2 && return Invariants(ms, σ2 = _σk,  σ1 = σ1of2(_cosθij,_σk,ms^2))
    k != 3 && error("k should be 1,2, or 3, while it is $k")
    return Invariants(ms, σ3 = _σk,  σ2 = σ2of3(_cosθij,_σk,ms^2))
end

invsquaredalitz1(x,y,ms) = invsquaredalitz(1,x,y,ms)
invsquaredalitz2(x,y,ms) = invsquaredalitz(2,x,y,ms)
invsquaredalitz3(x,y,ms) = invsquaredalitz(3,x,y,ms)

"""
    jacobian_squaredalitz(k,σs,ms)

calculates jacobian of transformation to square coordinates
"""
function jacobean_squaredalitz(k,σs,ms::ThreeBodyMasses)
    (i,j) = ij_from_k(k)
    #
    cθ = cosθij(k,σs,ms^2)
    cθ = cθ ≥ 1.0 ? 1.0 : (cθ ≤ -1.0 ? -1.0 : cθ)
    dydcθ = 1 / sqrt(1-cθ^2) / π
    #
    m = sqrt(σs[k])
    sqrtlimsk = sqrt.(lims(k, ms))
    xn = 2*(m - sqrtlimsk[1]) / (sqrtlimsk[2]-sqrtlimsk[1]) - 1
    dxdm = 2 / (sqrtlimsk[2]-sqrtlimsk[1]) / sqrt(1-xn^2) / π
    ρρm = sqrt(Kallen(σs[k],ms[i]^2,ms[j]^2)*Kallen(ms[4]^2,σs[k],ms[k]^2)) / (σs[k]*ms[4]^2) * m
    return ρρm / dydcθ / dxdm #
end

jacobean_squaredalitz1(σs,ms) = jacobean_squaredalitz(1,σs,ms)
jacobean_squaredalitz2(σs,ms) = jacobean_squaredalitz(2,σs,ms)
jacobean_squaredalitz3(σs,ms) = jacobean_squaredalitz(3,σs,ms)

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
"""
    abstract type AbstractFinder <: Any
Supertype for various "lookup" functors.
"""
abstract type AbstractFinder end

"""
    struct DeclFinder <: AbstractFinder
Functor for decl lookup.
"""
struct DeclFinder <: AbstractFinder
    spec::CXXScopeSpec
    result::LookupResult
    kind::CXLookupNameKind
end
DeclFinder(r::LookupResult, k::CXLookupNameKind) = DeclFinder(CXXScopeSpec(), r, k)

function dispose(x::DeclFinder)
    dispose(x.spec)
    dispose(x.result)
end

function DeclFinder(ci::CompilerInstance, kind::CXLookupNameKind=CXLookupNameKind_LookupOrdinaryName)
    sema = getSema(ci)
    loc = get_main_file_begin_loc(getSourceManager(ci))  # used as a fake loc
    result = LookupResult(sema, DeclarationName(), loc, kind)
    return DeclFinder(result, kind)
end

function reset(x::DeclFinder)
    clear(x.result, x.kind)
    clear(x.spec)
end

function get_decl(x::DeclFinder)
    isempty(x.result) && error("no lookup result.")
    return getRepresentativeDecl(x.result)
end

function (x::DeclFinder)(ci::CompilerInstance, parser::Parser, decl::String, scope::String="")
    reset(x)
    parse_cxx_scope_spec(ci, parser, scope, x.spec)
    setLookupName(x.result, DeclarationName(get_name(getASTContext(ci), decl)))
    if isValid(x.spec)
        is_found = LookupParsedName(getSema(ci), x.result, getCurScope(parser), x.spec)
        !is_found && error("failed to find decl(`$decl`) in the CXXScopeSpec.")
    else
        is_found = LookupName(getSema(ci), x.result, getCurScope(parser))
        !is_found && error("failed to find decl(`$decl`) in the current scope.")
    end
    return !isempty(x.result)
end


text:
export sinc_interpolate, sinc_interpolate_sum, downsample
export resample

"""
    resample(arr, new_size [, normalize]; take_real=true, boundary_handling=true)

Calculates the `sinc` interpolation of an `arr` on a new array size
`new_size`.
It is a re-evaluation of the Fourier series at new grid points.
`new_size` can be arbitrary. Each dimension is then independently either up or downsampled.

This method is based on FFTs and therefore implicitly assumes periodic
boundaries and a finite frequency support.
`normalize=true` by default multiplies by an appropriate factor so that
the average intensity stays the same.
If `size(new_size)[i] > size(arr)[i]`, we apply zero padding in Fourier space.
If `size(new_size)[i] < size(arr)[i]`, we cut out a centered part of the
Fourier spectrum.

We apply some tricks at the boundary to increase accuracy of highest frequencies.
If you set the keyword argument `boundary_handling=false` you turn this boundary handling off.
This can increase inaccuracies for the highest frequency (not a big problem in real images/signals).
However, `resample` is then faster on CPUs and especially on CUDA GPUs.

# Examples
```jldoctest
julia> resample([1.0, 2.0, 3.0, 4.0], 8)
8-element Array{Float64,1}:
 1.0
 1.085786437626905
 2.0
 2.5
 3.0
 3.914213562373095
 4.0
 2.5

julia> resample([1.0  2.0; 3.0 4.0], (4,4))
4×4 Array{Float64,2}:
 1.0  1.5  2.0  1.5
 2.0  2.5  3.0  2.5
 3.0  3.5  4.0  3.5
 2.0  2.5  3.0  2.5

julia> resample([1.0, 0.0, 1.0, 0.0, 1.0, 0.0], (3))
3-element Array{Float64,1}:
 0.5
 0.5
 0.5

julia> resample([1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0], (6))
6-element Array{Float64,1}:
  1.0
 -0.3333333333333333
  1.0
 -0.3333333333333333
  1.0
 -0.3333333333333333

julia> resample([1 2 3; 4 5 6], (3, 2))
3×2 Array{Float64,2}:
 1.0   3.0
 3.25  5.25
 3.25  5.25
```
"""
function resample(arr::AbstractArray{T, N}, new_size, normalize=true; take_real=true, boundary_handling=true) where {T<:Real, N}

    if typeof(new_size) <: Number
        new_size = Tuple(new_size)
    end

    # go to fourier space
    arr_f = fftshift(fft(arr))
    # create fourier space new array

    # the idea is the following
    # 1) first we handle all the upsampling dimensions. By doing so, we leave the size
    # or increase it. downsampling is done afterwards. Hence we use max to create
    # a new array which is a least as large as the initial one or larger
    # 2) we then handle all the downsampling dimensions
    new_size_interp = Tuple(max(x[1], x[2]) for x in zip(size(arr), new_size))
    out_f = similar(arr_f, new_size_interp)
    fill!(out_f, zero(eltype(arr_f)))

    # change arr_f to be a hermitian array because we want a purely real result
    # after iffting

    # in that case, make_hermitian is useless since we only cut out
    if boundary_handling && ~(all(new_size < size(arr)))
        arr_f = make_hermitian(arr_f)
    end

    # it can happen, that arr_f was now padded with an extra row
    # but there is no extra column in out_f
    # therefore, cut it
    inds = []
    for (a, o) in zip(size(arr_f), size(out_f))
        push!(inds, 1:min(a, o))
    end
    arr_f = arr_f[inds...]

    # set the old array into the new 0-padded array
    center_set!(out_f, arr_f)
    arr_f = out_f

    # part 2). We handle now the downsampling
    # if the new_size[d] is even, we need to add the highest positive frequency
    # of the initial spectrum
    # to the highest negative one. In that way, we get a purely real result
    if boundary_handling
        arr_f = add_high_frequencies(size(arr_f), arr_f, new_size, N)
    end
    # return arr_out_f

    # do the cutting in Fourier space
    arr_f_n = center_extract(arr_f, new_size)
    # back to real space
    arr_out = ifft(ifftshift(arr_f_n))

    if normalize
        arr_out .*= length(arr_out) ./ length(arr)
    end
    if take_real
        return real(arr_out)
    else
        return arr_out
    end
end

 # for a complex signal: split into real and imaginary part and solve for each
 # part individually
function resample(arr::AbstractArray{T}, new_size, normalize=true) where T<:Complex
    # array of same shape but with Complex element type
    arr_r = resample(real(arr), new_size, normalize)
    arr_i = 1im .* resample(imag(arr), new_size, normalize)
    return arr_r .+ arr_i
end

"""
    sinc_interpolate_sum(arr, new_length)

Calculates the `sinc` interpolation of an 1D `arr` on a new array size
`new_size`.
This method is slow, because of an explicit sum evalulation and
not a FFT based evaluation.

# Examples
```jldoctest
julia> sinc_interpolate_sum([1.0, 2.0, 3.0, 4.0], 8)
8-element Array{Float64,1}:
 1.0
 1.7825353626292277
 2.0
 2.1220659078919377
 3.0
 4.1592491794681985
 4.0
 2.0735615442829793
```
"""
function sinc_interpolate_sum(arr, new_size)
    out = zeros(eltype(arr), new_size)
    T = 1 / size(arr)[1]
    t_arr = range(0, 1.0, length=size(out)[1]+1)[1:size(out)[1]]
    for (j, t) in enumerate(t_arr)
        v = zero(eltype(arr))
        for n = 1:size(arr)[1]
            v += arr[n] * sinc((t - (n-1) * T) / T)
        end
        out[j] = v
    end
    return out
end


text:
function f(x, y) end
f(println(1), println(2))

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
function storeinstr_preprend(op::Operation, vloopsym::Symbol)
    # defaultstoreop = :vstore!
    # defaultstoreop = :vnoaliasstore!
    isvectorized(op) && return Symbol("")
    vloopsym ∉ reduceddependencies(op) && return Symbol("")
    # vectorized is not a loopdep, but is a reduced dep
    opp = first(parents(op))
    # while vectorized ∉ loopdependencies(opp)
    while !isvectorized(opp)
        oppold = opp
        for oppp ∈ parents(opp)
            if vloopsym ∈ reduceddependencies(oppp)
                @assert opp !== oppp "More than one parent is a reduction over the vectorized variable."
                opp = oppp
            end
        end
        @assert opp !== oppold "Failed to find any parents "
    end
    reduction_to_scalar(reduction_instruction_class(instruction(opp)))
end

function reduce_expr_u₂(toreduct::Symbol, instr::Instruction, u₂::Int)
    t = Expr(:tuple)
    for u ∈ 0:u₂-1
        push!(t.args, Symbol(toreduct, u))
    end
    Expr(:call, lv(:reduce_tup), reduce_to_onevecunroll(instr), t)
end
function reduce_expr!(q::Expr, toreduct::Symbol, instr::Instruction, u₁::Int, u₂::Int, isu₁unrolled::Bool, isu₂unrolled::Bool)
    if isu₂unrolled# u₂ != -1
        _toreduct = Symbol(toreduct, 0)
        push!(q.args, Expr(:(=), _toreduct, reduce_expr_u₂(toreduct, instr, u₂)))
    else
        _toreduct = Symbol(toreduct, '_', u₁)
    end
    if (u₁ == 1) | (~isu₁unrolled)
        push!(q.args, Expr(:(=), Symbol(toreduct, "##onevec##"), _toreduct))
    else
        push!(q.args, Expr(:(=), Symbol(toreduct, "##onevec##"), Expr(:call, lv(reduction_to_single_vector(instr)), _toreduct)))
        # push!(q.args, :(@show $_toreduct))
        # push!(q.args, Expr(:(=), Symbol(toreduct, "##onevec##"), :(@show $(Expr(:call, lv(reduction_to_single_vector(instr)), _toreduct)))))
    end
    nothing
end

function lower_store_collection!(
    q::Expr, ls::LoopSet, op::Operation, ua::UnrollArgs, mask::Bool, inds_calc_by_ptr_offset::Vector{Bool}
)
    omop = offsetloadcollection(ls)
    batchid, bopind = omop.batchedcollectionmap[identifier(op)]
    collectionid, copind = omop.opidcollectionmap[identifier(op)]
    opidmap = offsetloadcollection(ls).opids[collectionid]
    idsformap = omop.batchedcollections[batchid]

    @unpack u₁, u₁loopsym, u₂loopsym, vloopsym, vloop, u₂max, suffix = ua
    ops = operations(ls)
    # __u₂max = ls.unrollspecification[].u₂
    nouter = length(idsformap)

    t = Expr(:tuple)
    # u = Core.ifelse(isu₁unrolled(op), u₁, 1)

    for (i,(opid,_)) ∈ enumerate(idsformap)
        opp = first(parents(ops[opidmap[opid]]))

        isu₁, isu₂ = isunrolled_sym(opp, u₁loopsym, u₂loopsym, vloopsym)#, __u₂max)
        u = Core.ifelse(isu₁, u₁, 1)
        mvar = Symbol(variable_name(opp, ifelse(isu₂, suffix, -1)), '_', u)
        # mvar = Symbol(variable_name(_op, suffix), '_', u)
        # mvar = Symbol(variable_name(_op, ifelse(isu₂, suffix, -1)), '_', u)
        # @show mvar, isu₂, isu₂unrolled(opp)
        push!(t.args, mvar)
    end

    offset_dummy_loop = Loop(first(getindices(op)), MaybeKnown(1), MaybeKnown(1024), MaybeKnown(1), Symbol(""), Symbol(""))
    unrollcurl₂ = unrolled_curly(op, nouter, offset_dummy_loop, vloop, mask, 1)
    inds = mem_offset_u(op, ua, inds_calc_by_ptr_offset, false)

    falseexpr = Expr(:call, lv(:False));
    trueexpr = Expr(:call, lv(:True));
    rs = staticexpr(reg_size(ls));
    if isu₁unrolled(op) && u₁ > 1 # both unrolled
        if first(getindices(op)) === vloopsym
            interleaveval = -nouter
        else
            interleaveval = 0
        end
        unrollcurl₁ = unrolled_curly(op, u₁, ua.u₁loop, vloop, mask, interleaveval)
        inds = Expr(:call, unrollcurl₁, inds)
    end
    uinds = Expr(:call, unrollcurl₂, inds)
    vp = vptr(op)
    storeexpr = Expr(:call, lv(:_vstore!), vp, Expr(:call, lv(:VecUnroll), t), uinds)
    # not using `add_memory_mask!(storeexpr, op, ua, mask)` because we checked `isconditionalmemop` earlier in `lower_load_collection!`
    mask && push!(storeexpr.args, MASKSYMBOL)
    push!(storeexpr.args, falseexpr, trueexpr, falseexpr, rs)
    push!(q.args, storeexpr)
    nothing
end
function lower_store!(
    q::Expr, ls::LoopSet, op::Operation, ua::UnrollArgs, mask::Bool,
    reductfunc::Symbol = storeinstr_preprend(op, ua.vloop.itersymbol), inds_calc_by_ptr_offset = indices_calculated_by_pointer_offsets(ls, op.ref)
)
    @unpack u₁, u₁loopsym, u₂loopsym, vloopsym, vloop, u₂max, suffix = ua

    omop = offsetloadcollection(ls)
    batchid, opind = omop.batchedcollectionmap[identifier(op)]
    if ((batchid ≠ 0) && isvectorized(op)) && (!rejectinterleave(op))
        (opind == 1) && lower_store_collection!(q, ls, op, ua, mask, inds_calc_by_ptr_offset)
        return
    end

    isunrolled₁ = isu₁unrolled(op) #u₁loopsym ∈ loopdependencies(op)
    # isunrolled₂ = isu₂unrolled(op)
    falseexpr = Expr(:call, lv(:False)); trueexpr = Expr(:call, lv(:True)); rs = staticexpr(reg_size(ls));
    opp = first(parents(op))
    if ((opp.instruction.instr === reductfunc) || (opp.instruction.instr === :identity)) && isone(length(parents(opp)))
        opp = only(parents(opp))
    end
    # __u₂max = ls.unrollspecification[].u₂
    isu₁, isu₂ = isunrolled_sym(opp, u₁loopsym, u₂loopsym, vloopsym)#, __u₂max)
    # @show isu₁, isu₂, opp, u₁loopsym, u₂loopsym, vloopsym
    u = isu₁ ? u₁ : 1
    mvar = Symbol(variable_name(opp, ifelse(isu₂, suffix, -1)), '_', u)
    if all(op.ref.loopedindex)
        inds = unrolledindex(op, ua, mask, inds_calc_by_ptr_offset)
        storeexpr = if reductfunc === Symbol("")
            Expr(:call, lv(:_vstore!), vptr(op), mvar, inds)
        else
            Expr(:call, lv(:_vstore!), lv(reductfunc), vptr(op), mvar, inds)
        end
        add_memory_mask!(storeexpr, op, ua, mask)
        push!(storeexpr.args, falseexpr, trueexpr, falseexpr, rs)
        push!(q.args, storeexpr)
    elseif u₁ > 1
        mvard = Symbol(mvar, "##data##")
        isu₁ && push!(q.args, Expr(:(=), mvard, Expr(:call, lv(:data), mvar)))
        for u ∈ 1:u₁
            mvaru = :(getfield($mvard, $u, false))
            inds = mem_offset_u(op, ua, inds_calc_by_ptr_offset, true, u-1)
            # @show isu₁unrolled(opp), opp
            storeexpr = if isu₁
                if reductfunc === Symbol("")
                    Expr(:call, lv(:_vstore!), vptr(op), mvaru, inds)
                else
                    Expr(:call, lv(:_vstore!), lv(reductfunc), vptr(op), mvaru, inds)
                end
            elseif reductfunc === Symbol("")
                Expr(:call, lv(:_vstore!), vptr(op), mvar, inds)
            else
                Expr(:call, lv(:_vstore!), lv(reductfunc), vptr(op), mvar, inds)
            end
            domask = mask && (isvectorized(op) & ((u == u₁) | (vloopsym !== u₁loopsym)))
            add_memory_mask!(storeexpr, op, ua, domask)# & ((u == u₁) | isvectorized(op)))
            push!(storeexpr.args, falseexpr, trueexpr, falseexpr, rs)
            push!(q.args, storeexpr)
        end
    else
        inds = mem_offset_u(op, ua, inds_calc_by_ptr_offset, true, 0)
        storeexpr = if reductfunc === Symbol("")
            Expr(:call, lv(:_vstore!), vptr(op), mvar, inds)
        else
            Expr(:call, lv(:_vstore!), lv(reductfunc), vptr(op), mvar, inds)
        end
        add_memory_mask!(storeexpr, op, ua, mask)
        push!(storeexpr.args, falseexpr, trueexpr, falseexpr, rs)
        push!(q.args, storeexpr)
    end
    nothing
end

function lower_tiled_store!(
    blockq::Expr, opsv1::Vector{Operation}, opsv2::Vector{Operation}, ls::LoopSet, unrollsyms::UnrollSymbols, u₁::Int, u₂::Int, mask::Bool
)
    ua = UnrollArgs(ls, u₁, unrollsyms, u₂, 0)
    for opsv ∈ (opsv1, opsv2)
        for op ∈ opsv
            lower_tiled_store!(blockq, op, ls, ua, u₁, u₂, mask)
        end
    end
end

function donot_tile_store(ls::LoopSet, op::Operation, vloop::Loop, reductfunc::Symbol, u₂::Int)
    (!((reductfunc === Symbol("")) && all(op.ref.loopedindex))) || (u₂ ≤ 1) || isconditionalmemop(op) && return true
    rejectcurly(op) && return true
    omop = offsetloadcollection(ls)
    batchid, opind = omop.batchedcollectionmap[identifier(op)]
    return ((batchid ≠ 0) && isvectorized(op)) && (!rejectinterleave(op))
end

# VectorizationBase implements optimizations for certain grouped stores
# thus we group stores together here to allow for these possibilities.
# (In particular, it tries to replace scatters with shuffles when there are groups
#   of stores offset from one another.)
function lower_tiled_store!(blockq::Expr, op::Operation, ls::LoopSet, ua::UnrollArgs, u₁::Int, u₂::Int, mask::Bool)
    @unpack u₁loopsym, u₂loopsym, vloopsym, u₁loop, u₂loop, vloop = ua
    reductfunc = storeinstr_preprend(op, vloopsym)
    inds_calc_by_ptr_offset = indices_calculated_by_pointer_offsets(ls, op.ref)

    if donot_tile_store(ls, op, vloop, reductfunc, u₂)
        # If we have a reductfunc, we're using a reducing store instead of a contiuguous or shuffle store anyway
        # so no benefit to being able to handle that case here, vs just calling the default `lower_store!` method
        @unpack u₁, u₂max = ua
        for t ∈ 0:u₂-1
            unrollargs = UnrollArgs(u₁loop, u₂loop, vloop, u₁, u₂max, t)
            lower_store!(blockq, ls, op, unrollargs, mask, reductfunc, inds_calc_by_ptr_offset)
        end
        return
    end
    opp = first(parents(op))
    if (opp.instruction.instr === reductfunc) && isone(length(parents(opp)))
        throw("Operation $opp's instruction is $reductfunc, shouldn't be able to reach here.")
        # opp = only(parents(opp))
    end
    isu₁, isu₂ = isunrolled_sym(opp, u₁loopsym, u₂loopsym, vloopsym)#, u₂)
    @assert isu₂
    # It's reasonable forthis to be `!isu₁`
    u = Core.ifelse(isu₁, u₁, 1)
    tup = Expr(:tuple)
    for t ∈ 0:u₂-1
        mvar = Symbol(variable_name(opp, t), '_', u)
        push!(tup.args, mvar)
    end
    vut = Expr(:call, lv(:VecUnroll), tup) # `VecUnroll` of `VecUnroll`s
    inds = mem_offset_u(op, ua, inds_calc_by_ptr_offset, false)
    unrollcurl₂ = unrolled_curly(op, u₂, u₂loop, vloop, mask)
    falseexpr = Expr(:call, lv(:False)); trueexpr = Expr(:call, lv(:True)); rs = staticexpr(reg_size(ls));
    if isu₁ && u₁ > 1 # both unrolled
        unrollcurl₁ = unrolled_curly(op, u₁, u₁loop, vloop, mask)
        inds = Expr(:call, unrollcurl₁, inds)
    end
    uinds = Expr(:call, unrollcurl₂, inds)
    storeexpr = Expr(:call, lv(:_vstore!), vptr(op), vut, uinds)
    if mask && isvectorized(op)
        # add_memory_mask!(storeexpr, op, ua, mask)
        # we checked for `isconditionalmemop` earlier, so we skip this check
        # and just directly take the branch in `add_memory_mask!`
        push!(storeexpr.args, MASKSYMBOL)
    end
    push!(storeexpr.args, falseexpr, trueexpr, falseexpr, rs)
    push!(blockq.args, storeexpr)
    nothing
end


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 7f583ea2-f14f-11eb-3f1a-df193d312557
begin
	using Pkg
	Pkg.activate()
end

# ╔═╡ 4ba5f862-0e30-4a37-9dba-a47c4426e552
begin
	Pkg.develop("TSAnalysis")
	using Dates, DataFrames, LinearAlgebra, FredData, Optim, Plots, Measures;
	using TSAnalysis;
	plotlyjs();
	fred = Fred("6d4cdbed05fa2d25c171e82d73c77e64");
end

# ╔═╡ 6f90864a-1d27-45f6-8a1a-73df7442b952
"""
    download_fred_vintage(tickers::Array{String,1}, transformations::Array{String,1})

Download multivariate data from FRED2.
"""
function download_fred_vintage(tickers::Array{String,1}, transformations::Array{String,1})

    # Initialise output
    output_data = DataFrame();

    # Loop over tickers
    for i=1:length(tickers)

        # Download from FRED2
        fred_data = get_data(fred, tickers[i], observation_start="1984-01-01", units=transformations[i]).data[:, [:date, :value]];
        rename!(fred_data, Symbol.(["date", tickers[i]]));

        # Store current vintage
        if i == 1
            output_data = copy(fred_data);
        else
            output_data = join(output_data, fred_data, on=:date, kind = :outer);
        end
    end

    # Return output
    return output_data;
end

# ╔═╡ 9487ce4c-7093-4e7b-8371-c8ebc27e620f
begin
	# Download data of interest
	Y_df = download_fred_vintage(["INDPRO"], ["log"]);

	# Convert to JArray{Float64}
	Y = Y_df[:,2:end] |> JArray{Float64};
	Y = permutedims(Y);
end

# ╔═╡ 60c09cf6-fbde-4298-b5f6-a2fafec2e679
names(TSAnalysis)

# ╔═╡ 11904f9f-1004-4635-b6e2-f09ea3d042bb
TSAnalysis

# ╔═╡ 886a4bfb-bd85-42ae-99c5-e8d4e8b49620
begin
	d = 2;
	p = 20;
	q = 7;
	arima_settings = ARIMASettings(Y, d, p, q);

	# Estimation
	arima_out = arima(arima_settings, NelderMead(), Optim.Options(iterations=10000, f_tol=1e-2, x_tol=1e-2, g_tol=1e-2, show_trace=true, show_every=500));
end

# ╔═╡ 81485459-b72d-407f-be69-7edd60874e7b
begin
	max_hz = 100;
	fc = forecast(arima_out, max_hz, arima_settings);
end

# ╔═╡ cf73fad4-197b-4fd3-895b-844fb61ce6d5
begin
	# Extend date vector
	date_ext = Y_df[!,:date] |> Array{Date,1};

	for hz=1:max_hz
	    last_month = month(date_ext[end]);
	    last_year = year(date_ext[end]);

	    if last_month == 12
	        last_month = 1;
	        last_year += 1;
	    else
	        last_month += 1;
	    end

	    push!(date_ext, Date("01/$(last_month)/$(last_year)", "dd/mm/yyyy"))
	end

	my_font = font(8, "Helvetica Neue")

	# Generate plot
	p_arima = plot(date_ext, [Y[1,:]; NaN*ones(max_hz)], label="Data", color=RGB(0,0,200/255),
	               xtickfont=my_font, ytickfont=my_font,
	               title="INDPRO", titlefont=my_font, framestyle=:box,
	               legend=:right, size=(800,250), dpi=300, margin = 5mm);

	plot!(date_ext, [NaN*ones(length(date_ext)-size(fc,2)); fc[1,:]], label="Forecast", color=RGB(0,0,200/255), line=:dot)
end

# ╔═╡ Cell order:
# ╠═7f583ea2-f14f-11eb-3f1a-df193d312557
# ╠═4ba5f862-0e30-4a37-9dba-a47c4426e552
# ╠═6f90864a-1d27-45f6-8a1a-73df7442b952
# ╠═9487ce4c-7093-4e7b-8371-c8ebc27e620f
# ╠═60c09cf6-fbde-4298-b5f6-a2fafec2e679
# ╠═11904f9f-1004-4635-b6e2-f09ea3d042bb
# ╠═886a4bfb-bd85-42ae-99c5-e8d4e8b49620
# ╠═81485459-b72d-407f-be69-7edd60874e7b
# ╠═cf73fad4-197b-4fd3-895b-844fb61ce6d5


text:
using Redux
using CImGui

include("Counter.jl")
using .Counter

include("../Renderer.jl")
using .Renderer

const store = create_store(Counter.counter, Counter.State(0))

function counter_ui(store)
    flag = CImGui.ImGuiWindowFlags_NoTitleBar |
           CImGui.ImGuiWindowFlags_NoResize |
           CImGui.ImGuiWindowFlags_AlwaysAutoResize |
           CImGui.ImGuiWindowFlags_NoSavedSettings |
           CImGui.ImGuiWindowFlags_NoFocusOnAppearing |
           CImGui.ImGuiWindowFlags_NoNav
    CImGui.Begin("Counter", Ref(true), flag)
        spacing = CImGui.GetStyle().ItemInnerSpacing.x
        CImGui.PushButtonRepeat(true)
        CImGui.ArrowButton("##left", CImGui.ImGuiDir_Left) && dispatch!(store, Counter.DECREMENT)
        CImGui.SameLine(0.0, spacing)
        CImGui.ArrowButton("##right", CImGui.ImGuiDir_Right) && dispatch!(store, Counter.INCREMENT)
        CImGui.PopButtonRepeat()
        CImGui.SameLine()
        value = get_state(store).counter
        CImGui.Text("$value")
    CImGui.End()
end

Renderer.render(()->counter_ui(store), width=180, height=50, title="App: Counter")


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/MLJTutorials/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/MLJTutorials/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("MLJTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

using MLJ, PrettyPrinting

@load KNNRegressor
# input
X = (age    = [23, 45, 34, 25, 67],
     gender = categorical(['m', 'm', 'f', 'm', 'f']))
# target
height = [178, 194, 165, 173, 168];

scitype(X.age)

pipe = @pipeline MyPipe(X -> coerce(X, :age=>Continuous),
                       hot = OneHotEncoder(),
                       knn = KNNRegressor(K=3),
                       target = UnivariateStandardizer());

pipe.knn.K = 2
pipe.hot.drop_last = true;

evaluate(pipe, X, height, resampling=Holdout(),
         measure=rms) |> pprint

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
# Postcompose with scaling and sum

"""
  Postcompose(f::ProximableFunction, a::Real=1.0, b::Real=0.0)

Returns the function `g(x) = a*f(x) + b`.
"""

immutable Postcompose{T <: ProximableFunction, S <: Real} <: ProximableFunction
  f::T
  a::S
  b::S
  function Postcompose(f::T, a::S, b::S)
    if a <= 0.0
      error("parameter a must be positive")
    else
      new(f, a, b)
    end
  end
end

is_prox_accurate(f::Postcompose) = is_prox_accurate(f.f)

Postcompose{T <: ProximableFunction, S <: Real}(f::T, a::S=1.0, b::S=0.0) = Postcompose{T, S}(f, a, b)

function (g::Postcompose){T <: RealOrComplex}(x::AbstractArray{T})
  return g.a*g.f(x) + g.b
end

function prox!{T <: RealOrComplex}(g::Postcompose, x::AbstractArray{T}, y::AbstractArray{T}, gamma::Real=1.0)
  v = prox!(g.f, x, y, g.a*gamma)
  return g.a*v + g.b
end

function prox_naive{T <: RealOrComplex}(g::Postcompose, x::AbstractArray{T}, gamma::Real=1.0)
  y, v = prox_naive(g.f, x, g.a*gamma)
  return y, g.a*v + g.b
end

fun_name(f::Postcompose) = string("Postcomposition of ", fun_name(f.f))
fun_dom(f::Postcompose) = fun_dom(f.f)
fun_expr(f::Postcompose) = "x ↦ a*f(x)+b"
fun_params(f::Postcompose) = string("f(x) = ", fun_expr(f.f), ", a = $(f.a), b = $(f.b)")


text:
using Plots
using ColorTypes
using Images
gr()

function anim_viz(x_spikes::Array{Int64}, y_spikes::Array{Int64},
                  img::Array{RGB{Float64},2}, n_neurons::Int64, step::Int64,
                  window::Int64 = 60)
    if step > window
        p1 = scatter(
            x_spikes[step-window:step],
            y_spikes[step-window:step],
            ylim = [0,n_neurons],
            xlabel = "Step",
            ylabel = "Neuron",
            legend = false
        )
        p2 = plot(img, ticks = nothing)
        plot(p1, p2, layout = 2)
    end
end

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
# Check consistency with the original WORLD
# Ground truth data was genearted using `generate_groundtruth.sh` in
# https://github.com/r9y9/WORLD
# Note that all of the results of the original WORLD were dumped as %.16lf
# accuracy on Ubuntu 14.04 64bit machine. g++ v4.9.2 was used to compile the
# WORLD.

@assert WORLD.version >= v"0.2.1-4"

info("Check consistency with the original WORLD")

x = vec(readdlm(joinpath(dirname(@__FILE__), "data", "x.txt")))

fs = 16000
period = 5.0
opt = DioOption(71.0, 800.0, 2, period, 1, 0.1) # 0.1 = 0.02 * 5.0

# Fundamental frequency (f0) estimation by DIO
f0, timeaxis = dio(x, fs, opt)
f0_org = vec(readdlm(joinpath(dirname(@__FILE__), "data", "f0.txt")))

println("Maximum error in DIO is $(maximum(abs(f0-f0_org)))")
@test length(f0) == length(f0_org)
@test_approx_eq_eps f0 f0_org 1.0e-10

# F0 refienment by StoneMask
f0 = stonemask(x, fs, timeaxis, f0)
f0_refined_org = vec(readdlm(joinpath(dirname(@__FILE__), "data", "f0_refined.txt")))

println("Maximum error in StoneMask is $(maximum(abs(f0-f0_refined_org)))")
@test length(f0) == length(f0_refined_org)
@test_approx_eq_eps f0 f0_refined_org 1.0e-10

# Spectral envelope estimation by CheapTrick
spectrogram = cheaptrick(x, fs, timeaxis, f0; opt=CheapTrickOption(-0.15))
spectrogram_org = readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram.txt"))'

println("Maximum error in CheapTrick is $(maximum(abs(spectrogram - spectrogram_org)))")
@test size(spectrogram) == size(spectrogram_org)
@test_approx_eq_eps spectrogram spectrogram_org 1.0e-10

aperiodicity = d4c(x, fs, timeaxis, f0; opt=D4COption())
aperiodicity_org = readdlm(joinpath(dirname(@__FILE__), "data", "aperiodicity.txt"))'
println("Maximum error in D4C is $(maximum(abs(aperiodicity-aperiodicity_org)))")
@test size(aperiodicity) == size(aperiodicity)
@test_approx_eq_eps aperiodicity aperiodicity_org 1.0e-10

# Synthesis
y_length = convert(Int, ((length(f0)-1)*period/1000 * fs) + 1)
y = synthesis(f0, spectrogram, aperiodicity, period, fs, y_length)
y_org = vec(readdlm(joinpath(dirname(@__FILE__), "data", "x_synthesized.txt")))
println("Maximum error in Synthesis is $(maximum(abs(y-y_org)))")
@test length(y) == length(y_org)
@test_approx_eq_eps y y_org 1.0e-10


text:
function ortho_gram_schmidt( vin::Array{Float64,2} )
    Ncol = size(vin)[2]
    v = copy(vin)
    for ii = 1:Ncol
        zz = dot( v[:,ii], v[:,ii] )
        v[:,ii] = v[:,ii]/sqrt(zz)
        for jj = ii+1 : Ncol
            puv = prj( v[:,ii], v[:,jj] )
            v[:,jj] = v[:,jj] - puv*v[:,ii]
        end
    end
    return v
end

function ortho_gram_schmidt( Ncol, vin )
    v = copy(vin)
    for ii = 1:Ncol
        zz = dot( v[:,ii], v[:,ii] )
        v[:,ii] = v[:,ii]/sqrt(zz)
        for jj = ii+1 : Ncol
          puv = prj( v[:,ii], v[:,jj] )
          v[:,jj] = v[:,jj] - puv*v[:,ii]
        end
    end
    return v
end

function ortho_gram_schmidt!( v )
    Ncol = size(v)[2]
    for ii = 1:Ncol
        zz = dot( v[:,ii], v[:,ii] )
        v[:,ii] = v[:,ii]/sqrt(zz)
        for jj = ii+1 : Ncol
            puv = prj( v[:,ii], v[:,jj] )
            v[:,jj] = v[:,jj] - puv*v[:,ii]
        end
    end
end

function ortho_gram_schmidt!( Ncol, v )
    for ii = 1:Ncol
        zz = dot( v[:,ii], v[:,ii] )
        v[:,ii] = v[:,ii]/sqrt(zz)
        for jj = ii+1 : Ncol
          puv = prj( v[:,ii], v[:,jj] )
          v[:,jj] = v[:,jj] - puv*v[:,ii]
        end
    end
end

function prj( u, v )
    return dot( u, v )/dot( u, u )
end


text:
using LightGraphs
using LightGraphs.SimpleGraphs

include("SI\\metaplex.jl")

N = 100
M = 10

g = random_configuration_model(N*M, fill(6,N*M))
g = complete_graph(N*M)
h = random_configuration_model(M,fill(3,M))
h = erdos_renyi(M, 0.5)
h = path_graph(M)
h = complete_graph(M)

using GraphPlot
gplot(h)

mpx = Metaplex(g,h)

Xi0 = falses(N*M)
Xi0[1:10] .= true

Xμ0 = [div(i-1,N)+1 for i in 1:N*M]
Xμ0 = [rand(1:M) for i in 1:N*M]

β = 0.2
D = [0.1, 0.1]
tmax = 3.0

@time t, pcs = metaplex_gillespie(mpx, Xi0, Xμ0, β, D, nmax = 500000, tmax=tmax)

using Plots
gr()
pyplot()

plot(t, pcs[:,2,:], label="", line=:steppre)

plot(t, sum(pcs[:,2,:], dims=2)/(N*M),
    label="",
    line=:steppre,
    xlabel="Time",
    ylabel="Fraction of infected individuals"
    )

plot(t, vec(sum(pcs[:,:,:], dims=[2,3])), label="")

t_mf, s_mf, i_mf, sμ_mf, iμ_mf, sol = metaplex_ode(mpx, Xi0, Xμ0, β, D, tmax = tmax)

plot(t_mf, reshape(sum(i_mf, dims=3), length(sol.t), N*M), label="")

plot(t_mf, iμ_mf, label="")

sum(iμ_mf, dims=2)/(N*M)

plot(t, sum(pcs[:,2,:], dims=2)/(N*M),
    label="Gillespie",
    line=:steppre,
    xlabel="Time",
    ylabel="Fraction of infected individuals")
plot!(t_mf, sum(iμ_mf, dims=2)/(N*M), label="Mean Field")

t_mc, s_mc, i_mc, sd_s_mc, sd_i_mc = metaplex_montecarlo(mpx, Xi0, Xμ0, β, D, tmax=tmax, nmax=500000)

plot(t_mc, i_mc, label="")
plot!(t_mf, iμ_mf, label="", linestyle=:dash)

plot(t_mc, sum(i_mc, dims=2)/(N*M), label="Average", legend=:bottomright)
plot!(t_mf, sum(iμ_mf, dims=2)/(N*M), label="Mean Field")

sum(i_mf, dims=2)/(N*M)
plot(t_mc[1:20], sum(i_mc[1:20,:], dims=2)/(N*M), label="Average", legend=:topleft)
plot!(t_mf[1:50], sum(iμ_mf, dims=2)[1:50]/(N*M), label="Mean Field")

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
example_ticket_str = ".--------------------------------------------------------.
| ????: 101    ?????: 102   ??????????: 103     ???: 104 |
|                                                        |
| ??: 301  ??: 302             ???????: 303      ??????? |
| ??: 401  ??: 402           ???? ????: 403    ????????? |
'--------------------------------------------------------'"

test_str = "class: 1-3 or 5-7
row: 6-11 or 33-44
seat: 13-40 or 45-50

your ticket:
7,1,14

nearby tickets:
7,3,47
40,4,50
55,2,20
38,6,12"

test_str2 = "class: 0-1 or 4-19
row: 0-5 or 8-19
seat: 0-13 or 16-19

your ticket:
11,12,13

nearby tickets:
3,9,18
15,1,5
5,14,9"

test_list = split(test_str, '\n')

test_list2 = split(test_str2, '\n')

list = readlines("day16_input.txt")

function count_invalid(list)
    rules_end = findfirst(==(""), list) - 1
    rules = map(l -> split(l, r": |-| or ") |> e -> (e[1], parse(Int, e[2]), parse(Int, e[3]), parse(Int, e[4]), parse(Int, e[5])), list[1:rules_end])
    total = 0
    for i in rules_end+6:length(list)
        nums = parse.(Int, split(list[i], ','))
        for n in nums
            if !any(x -> x[2] ≤ n ≤ x[3] || x[4] ≤ n ≤ x[5], rules)
                total += n
            end
        end
    end
    return total
end

@show count_invalid(test_list) == 71
@show count_invalid(list) == 27802

function count_ticket(list)
    num_rules = findfirst(==(""), list) - 1
    rules = map(l -> split(l, r": |-| or ") |> e -> Tuple([e[1], [parse(Int, e[j]) for j in 2:5]]), list[1:num_rules])
    rule_decoder = Dict(map(x -> 2^(x[1]-1) => x[2][1], enumerate(rules)))

    # get valid tickets
    other_tickets = fill(0, length(list) - num_rules - 5, num_rules)
    num_others = 0
    for i in num_rules+6:length(list)
        m = parse.(Int, split(list[i], ','))
        ok = true
        j = 1
        while ok == true && j ≤ num_rules
            if !any(x -> x[2][1] ≤ m[j] ≤ x[2][2] || x[2][3] ≤ m[j] ≤ x[2][4], rules)
                ok = false
            end
            j += 1
        end
        if ok
            num_others += 1
            other_tickets[num_others,:] = m
        end
    end
    other_tickets = other_tickets[1:num_others, :]

    # build array with possible rules (row) for each position (column)
    ruling_off = fill(true, num_rules, num_rules)
    for k in 1:num_rules
        for j in 1:num_rules
            x = rules[j][2]
            for i in 1:num_others
                n = other_tickets[i,k]
                if !(x[1] ≤ n ≤ x[2] || x[3] ≤ n ≤ x[4])
                    ruling_off[j,k] = false
                end
            end
        end
    end

    # find order from array of possible rules
    order = fill(0, num_rules)
    for k in 1:num_rules
        sums = sum(ruling_off, dims=1)
        argmin_sum = argmax(sum(ruling_off, dims=1))[2]
        for j in 1:num_rules
            if 0 < sums[j] ≤ sums[argmin_sum]
                argmin_sum = j
            end
        end
        aux = argmax(ruling_off[:,argmin_sum])
        order[argmin_sum] = 2^(aux-1)
        ruling_off[aux,:] = fill(0, num_rules)
    end

    # get product of values with rules containing "departure"
    product = 1
    ordered_rules = Tuple(rule_decoder[k] for k in order)
    myticket = parse.(Int, split(list[num_rules+3], ','))
    for j in 1:num_rules
        if occursin("departure", ordered_rules[j])
            product *= myticket[j]
        end
    end

    return product, ordered_rules
end

@show count_ticket(test_list2)[2] == ("row", "class", "seat")
@show count_ticket(list)[1] == 279139880759

nothing

#=
```julia
julia> @btime count_invalid(list);
  907.626 μs (1491 allocations: 440.66 KiB)

julia> @btime count_ticket(list);
  15.494 ms (214734 allocations: 3.86 MiB)
```
 =#


text:
# * Norms

_norm(B::FiniteDifferencesOrRestricted, c::AbstractArray, p::Real=2) =
    norm(c, p)*(step(B)^(inv(p)))

LinearAlgebra.norm(v::FDVecOrMat, p::Real=2) = _norm(v.args..., p)
LinearAlgebra.norm(v::Mul{<:Any, <:Tuple{<:FiniteDifferencesOrRestricted, <:AbstractArray}},
                   p::Real=2) = _norm(v.args..., p)

function LinearAlgebra.normalize!(v::FDVecOrMat, p::Real=2)
    v.args[2][:] /= norm(v,p)
    v
end

function LinearAlgebra.normalize!(v::Mul{<:Any, <:Tuple{<:FiniteDifferencesOrRestricted, <:AbstractArray}},
                                  p::Real=2)
    v.args[2][:] /= norm(v,p)
    v
end

# * Inner products

function _inner_product(a::Adjoint{<:Any,<:AbstractVector}, A::QuasiAdjoint{<:Any,<:FD},
                        B::FD, b::AbstractVector) where {FD<:FiniteDifferencesOrRestricted}
    axes(A.parent) == axes(B) || throw(ArgumentError("Incompatible axes"))
    a*b*step(B)
end

LazyArrays.materialize(inner_product::FDInnerProduct{T,U,FD}) where {T,U,FD<:FiniteDifferencesOrRestricted{U}} =
    _inner_product(inner_product.args...)

function LazyArrays.materialize(inner_product::LazyFDInnerProduct{FD}) where {FD<:FiniteDifferencesOrRestricted}
    aA,Bb = inner_product.args
    _inner_product(aA.args..., Bb.args...)
end

function LazyArrays.materialize(s::Mul{<:Any, <:Tuple{
    <:Mul{<:Any, <:Tuple{
        <:Adjoint{<:Any,<:AbstractVector},
        <:ContinuumArrays.QuasiArrays.QuasiAdjoint{<:Any, <:FD}}},
    <:Mul{<:Any, <:Tuple{
        <:FD,
        <:Diagonal,
        <:ContinuumArrays.QuasiArrays.QuasiAdjoint{<:Any, <:FD}}},
    <:Mul{<:Any, <:Tuple{
        <:FD,
        <:AbstractVector}}}}) where {FD<:FiniteDifferencesOrRestricted}
    a,o,b = s.args
    axes(a.args[2].parent) == axes(o.args[1]) &&
        axes(o.args[3].parent) == axes(b.args[1]) ||
        throw(ArgumentError("Incompatible axes"))
    av = first(a.args)
    ov = o.args[2].diag
    bv = last(b.args)
    v = zero(promote_type(eltype(av),eltype(ov),eltype(bv)))
    @inbounds for i in eachindex(bv)
        v += av[i]*ov[i]*bv[i]
    end
    v*step(first(b.args))
end


text:
# DEFAULT 2019/20 system
function load()
    sys = TaxBenefitSystem( name="Scotland Default 2019/20")
    sys.it.earnings_rates = [19.0,20.0,21.0,41.0,46.0]./100.0
    sys.it.earnings_bands = [2_049.0, 12_444.0, 30_930.0, 150_000.0]
    sys.it.savings_rates  = [0.0, 20.0, 40.0, 45.0]./100.0
    sys.it.savings_bands  = [5_000.0, 37_500.0, 150_000.0]
    sys.it.dividend_rates = [0.0, 7.5,32.5,38.1]./100.0
    sys.it.dividend_bands = [2_000.0, 37_500.0, 150_000.0]
    sys.it.personal_allowance          = 12_500.00
    sys.it.blind_persons_allowance     = 2_450.00
    sys.it.married_couples_allowance   = 8_915.00
    sys.it.mca_minimum                 = 3_450.00
    sys.it.marriage_allowance          = 1_250.00
    sys.it.personal_savings_allowance  = 1_000.00
    sys
end


text:
# ellipse to multiple forward simulated shapes
# only final landmark positions observed

using BridgeLandmarks
using Random
using Distributions
using StaticArrays
using LinearAlgebra
using JLD

workdir = @__DIR__
cd(workdir)

n = 20
nshapes = 1
T = 1.0

xobs0 = [PointF(2.0cos(t), 2.0sin(t))  for t in collect(0:(2pi/n):2pi)[2:end]]
# peach_xcoord(s) = (2.0 + sin(s)^3) * cos(s)
# peach_ycoord(s) = (2.0 + sin(s)^3) * sin(s)

heart_xcoord(s) = 0.2*(13cos(s)-5cos(2s)-2cos(3s)-cos(4s))
heart_ycoord(s) = 0.2*16(sin(s)^3)
xobsT = [PointF(heart_xcoord(t), heart_ycoord(t))  for t in (0:(2pi/n):2pi)[1:n]]

JLD2.@save "data_exp-ckl.jld2" xobs0 xobsT n nshapes


text:
function normalize_select(idx, @nospecialize(cols...))
    selected_cols = [Int[], Int[]]
    for i in 1:length(cols)
        normalize_select!(selected_cols, idx, cols[i])
    end
    unique!(selected_cols)
end
function normalize_select!(selected_cols, idx, cols::ColumnIndex)
    push!(selected_cols[1], idx[cols])
end

function normalize_select!(selected_cols, idx, cols::MultiColumnIndex)
    colsidx = idx[cols]
    for i in 1:length(colsidx)
        push!(selected_cols[1], colsidx[i])
    end
    if cols isa Not
        colsidxnot = idx[cols.skip]
        for i in 1:length(colsidxnot)
            push!(selected_cols[2], colsidxnot[i])
        end
    end
end

# sell and sell! will replace select and select!
# Dataset shouldn't support copycols since it causes modifying a data set without telling other alias data sets
function select(ds::Dataset, @nospecialize(args...))
    selected_cols_all = normalize_select(index(ds), args...)
    selected_cols = setdiff!(selected_cols_all[1], selected_cols_all[2])
    res = AbstractVector[]

    for j in 1:length(selected_cols)
        push!(res, copy(_columns(ds)[selected_cols[j]]))
    end
    newds = Dataset(res, _names(ds)[selected_cols], copycols = false)
    # else
    #     for j in 1:length(selected_cols)
    #         push!(res, _columns(ds)[selected_cols[j]])
    #     end
    #     newds = Dataset(res, _names(ds)[selected_cols], copycols = false)
    # end
    for j in 1:length(selected_cols)
        setformat!(newds, j => getformat(ds, selected_cols[j]))
    end
    if all(index(ds).sortedcols .∈ Ref(selected_cols))
        scols = Int[]
        revs = copy(index(ds).rev)
        scols_ds = index(ds).sortedcols
        for k in 1:length(scols_ds)
            push!(scols, findfirst(isequal(scols_ds[k]), selected_cols))
        end
        _copy_grouping_info!(newds, ds)
        # need to reorder things
        empty!(index(newds).sortedcols)
        empty!(index(newds).rev)
        append!(index(newds).sortedcols, scols)
        append!(index(newds).rev, revs)
    # else
    #     _reset_grouping_info!(ds)
    end
    setinfo!(newds, _attributes(ds).meta.info[])
    return newds
end

function select!(ds, @nospecialize(args...))
    selected_cols_all = normalize_select(index(ds), args...)
    selected_cols = setdiff!(selected_cols_all[1], selected_cols_all[2])
    unwanted_cols = setdiff(1:ncol(ds), selected_cols)
    sort!(unwanted_cols, rev = true)

    newnames = _names(ds)[selected_cols]
    newlookup = newnames .=> 1:length(selected_cols)
    res = AbstractVector[]

    for j in 1:length(selected_cols)
        push!(res, _columns(ds)[selected_cols[j]])
    end
    newformat = Dict{Int, Function}()
    for j in 1:length(selected_cols)
        push!(newformat, j => getformat(ds, selected_cols[j]))
    end
    if all(index(ds).sortedcols .∈ Ref(selected_cols))
        scols = Int[]
        # revs = copy(index(ds).rev)
        scols_ds = index(ds).sortedcols
        for k in 1:length(scols_ds)
            push!(scols, findfirst(isequal(scols_ds[k]), selected_cols))
        end
        # _copy_grouping_info!(ds, ds)
        # need to reorder things
        empty!(index(ds).sortedcols)
        # empty!(index(ds).rev)
        append!(index(ds).sortedcols, scols)
        # append!(index(ds).rev, revs)
    # else
    #     _reset_grouping_info!(ds)
    end
    for i in unwanted_cols
      popat!(_columns(ds), i)
      delete!(index(ds), i)
    end
    for j in 1:length(selected_cols)
      _columns(ds)[j] = res[j]
    end
    empty!(index(ds).lookup)
    empty!(index(ds).names)
    empty!(index(ds).format)
    for j in newlookup
      push!(index(ds).lookup, j)
    end
    for j in newnames
      push!(index(ds).names, j)
    end
    for j in newformat
      push!(index(ds).format, j)
    end
    _modified(_attributes(ds))
    ds
end


text:
push!(LOAD_PATH,"../src/")
using Base.Filesystem
using Revise
using DickeModel
using Documenter,DocumenterCitations
while true
print("Press enter to compile docs")
readline()
try
	rm("build", recursive=true,force=true)
catch
	print("folder locked, close tabs")
	continue
end
revise(DickeModel)
sleep(3)
bib = CitationBibliography("refs.bib")

try
makedocs(bib,
         sitename = "DickeModel.jl",
         modules  = [DickeModel],
         pages=[
                "The DickeModel.jl package" => "index.md",
                "Documentation" => [
                    "ClassicalDicke" => "classicaldicke.md",
					"ClassicalSystems" => "classicalsystems.md"
                ],
				"References" => "references.md"

               ])
catch e
        showerror(stdout, e)
end
end


text:
using F16Model

d2r = pi/180;
npos = 0;
epos = 0;
alt = 10000;
phi = 0;
theta = 0;
psi = 0;
Vt = 300;
alp = 15*d2r;
bet = 0;
p = 0;
q = 0;
r = 0;

x0 = [npos,epos,alt,phi,theta,psi,Vt,alp,bet,p,q,r];
u0 = [9000,0,0,0,0];

# Example 1: Determine xdot
# -------------------------
# There are two ways. See below.

xdot1 = zeros(12);
F16Model.Dynamics!(xdot1,x0,u0); # This is inplace implementation -- use this with DifferentialEquations package.

xdot2 = F16Model.Dynamics(x0,u0); # Use this for linearization of dynamics, etc.

# Example 2: Least restricted trim
# --------------------------------
# Trim the aircraft for SteadyLevel at h0,V0
# xbar = trim state
# ubar = trim control
# status = status of the optimization, status = 0 means optimization found solution and (xbar, ubar) defines a valid trim  point/
# prob = data structure from IpOpt.

h0 = 10000; # ft
Vt0 = 500;   # ft/s
xbar, ubar, status, prob = F16Model.Trim(h0,Vt0); # Default is steady-level

# Example 3: Restricted trim
# ---------------------------
# Trim for steady-climb and coordinated turn can be achieved by calling it with optional arguments
# γ = flight path angle
# ψdot = turn rate
# Defaults for these are zero, as shown in example 2.
# For example, to achieve trim at 5 degrees climb without turning, we can call Trim(...) as
xbar, ubar, status, prob = F16Model.Trim(h0, Vt0, γ=5*pi/180, ψdot=0);

# Example 4: Most restricted trim
# -------------------------------
# Values for ϕ, ψ, θ, α, β, p, q, r are defined as pair (initial_guess, isFixed).
# The nonlinear optimization uses initial_guess to warm start the iterations.
# If isFixed = 1, the optimization fixes the value of the state to initial_guess.
# For example, for the most restricted steady-level flight, we can call Trim(...) as follows.
# Note α and θ are left as free, since they must satisfy γ = θ - α. With γ = 0, we get θ = α.
xbar, ubar, status, prob = F16Model.Trim(h0, Vt0, γ=0, ψdot=0, ϕ=(0,1), ψ=(0,1), β=(0,1), p=(0,1), q=(0,1), r=(0,1));

# Example 5: Linearization about a given trim (xbar,ubar)
# -------------------------------------------------------
A, B = F16Model.Linearize(xbar,ubar); # Linearization is done using ForwardDiff


text:
"""
    HR*_PAT

Pattern to match horizontal rule indicators.
"""
const HR1_PAT = r"\-{3}[ \t\-]*$"
const HR2_PAT = r"\_{3}[ \t\_]*$"
const HR3_PAT = r"\*{3}[ \t\*]*$"

"""
    *_WHITESPACE_PAT

Pattern to match the whitespaces (tabs or spaces) at the start of a line, see
[`dedent`](@ref).
"""
const LEADING_WHITESPACE_PAT = r"^([ \t]*)\S"
const NEWLINE_WHITESPACE_PAT = r"\n([ \t]*)\S"

"""
    LX_COMMAND_PAT

Allowed latex command name.
Underscore are allowed inside the command but not at extremities.
The star `*` is not allowed anywhere.

## Examples:
* \\com
* \\ab1_cd*
"""
const LX_COMMAND_PAT = r"^\\[a-zA-Z](?:[_a-zA-Z0-9]*[a-zA-Z0-9])?$"

"""
    CODE_LANG*_PAT
"""
const CODE_LANG3_PAT = r"^`{3}[a-zA-Z][a-zA-Z0-9-]*$"
const CODE_LANG4_PAT = r"^`{4}[a-zA-Z][a-zA-Z0-9-]*$"
const CODE_LANG5_PAT = r"^`{5}[a-zA-Z][a-zA-Z0-9-]*$"

"""
    HTML_ENTITY_PAT

Pattern for an html entity.
Ref: https://dev.w3.org/html5/html-author/charref.

## Examples:
* ⊓
* ⊓
* ⊓
* ⊓

Note: longest entity is ∳ so capping the max
number of characters to 32.
"""
const HTML_ENTITY_PAT = r"^&(?:[a-zA-Z]{1,32}[0-9]{0,2}|#[0-9]{1,6}|#x[0-9a-f]{1,6});$"

const EMOJI_PAT = r"^\:[a-zA-Z0-9+-_]+\:$"

const FOOTNOTE_PAT = r"^\[\^[\p{L}0-9][\p{L}0-9_]*\](:)?$"

const OL_ITEM_PAT = r"^[0-9]{1,9}(?:[\.\)])[ \t]"


text:
@testset "Regularizer" begin

using ACE1pack

@info("Test generating laplacian regularizer")
# construct basis
ACE_basis = ACE1pack.generate_ace_basis(ace_basis_params(species = :Si, N = 2, maxdeg = 10))
pair_basis = ACE1pack.generate_pair_basis(pair_basis_params(species = :Si, maxdeg = 4))
basis = JuLIP.MLIPs.IPSuperBasis([pair_basis, ACE_basis])

rlap_scal = 2.0
regularizer = regularizer_params(type = "laplacian", rlap_scal = rlap_scal)
@test regularizer["rlap_scal"] == rlap_scal

expected_size = sum(length(b) for b in basis.BB)
P = generate_regularizer(basis, regularizer)
@test size(P) == (expected_size, expected_size)

end


text:
using OrthographicProjection
using Test

workdir = dirname(@__FILE__)
filename = "files/directory.txt" # point format 0

@testset "util" begin
    @testset "get_potree_dirs" begin
        potreedirs = OrthographicProjection.get_potree_dirs(joinpath(workdir,filename))
        @test potreedirs[1] == "C:\\Users\\Documents\\PROJ1"
        @test potreedirs[2] == "C:\\Users\\Documents\\PROJ2"
        @test potreedirs[3] == "C:\\Users\\Documents\\PROJ3"
        @test potreedirs[4] == "C:\\Users\\Documents\\PROJ4"
    end

    @testset "PO" begin
        PO = "XY+"
        @test OrthographicProjection.PO2matrix(PO) == [ 1.0  0.0  0.0; 0.0  1.0  0.0; 0.0  0.0  1.0]

        PO = "XY-"
        @test OrthographicProjection.PO2matrix(PO) == [ -1.0  0.0   0.0; 0.0  1.0   0.0; 0.0  0.0  -1.0]

        PO = "XZ+"
        @test OrthographicProjection.PO2matrix(PO) == [  -1.0  0.0  0.0; 0.0  0.0  1.0; 0.0  1.0  0.0]

        PO = "XZ-"
        @test OrthographicProjection.PO2matrix(PO) == [  1.0   0.0  0.0; 0.0   0.0  1.0; 0.0  -1.0  0.0]

        PO = "YZ+"
        @test OrthographicProjection.PO2matrix(PO) == [ 0.0  1.0  0.0; 0.0  0.0  1.0; 1.0  0.0  0.0]

        PO = "YZ-"
        @test OrthographicProjection.PO2matrix(PO) == [ 0.0  -1.0  0.0; 0.0   0.0  1.0; -1.0   0.0  0.0]
    end

    @testset "Raster array" begin
        PO = "XY+"
        coordsystemmatrix = OrthographicProjection.PO2matrix(PO)
        GSD = 0.1
        V = [ 0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0;
              0.0  0.0  1.0  1.0  0.0  0.0  1.0  1.0;
              0.0  1.0  0.0  1.0  0.0  1.0  0.0  1.0]
        EV = [[1, 2], [3, 4], [5, 6], [7, 8], [1, 3], [2, 4], [5, 7],
              [6, 8], [1, 5], [2, 6], [3, 7], [4, 8]]
        FV = [[1, 2, 3, 4], [5, 6, 7, 8], [1, 2, 5, 6], [3, 4, 7, 8],
              [1, 3, 5, 7], [2, 4, 6, 8]]
        model = V, EV, FV
        BGcolor = rand(3)

        RGBtensor, rasterquote, ref_X, ref_Y = OrthographicProjection.init_raster_array(coordsystemmatrix, GSD, model, BGcolor)
        @test RGBtensor[1,:,:] == fill(BGcolor[1],size(RGBtensor[1,:,:]))
        @test RGBtensor[2,:,:] == fill(BGcolor[2],size(RGBtensor[2,:,:]))
        @test RGBtensor[3,:,:] == fill(BGcolor[3],size(RGBtensor[3,:,:]))
        @test rasterquote == fill(-Inf,size(RGBtensor[1,:,:]))
        @test ref_X == 0.
        @test ref_Y == 1.
    end
end


text:
# this is a script to convert the representations to an easily readable format
using PdfFraud
using JLD2
using FileIO
using ProgressMeter

# first gather data and their labels
data_path = joinpath(PdfFraud.proto_path, "protos")
protos = joinpath.(data_path, readdir(data_path))

function load_representation(file::String)
	data = PdfFraud.read_pb(file)
	y = PdfFraud.visual_representation(data)
	(length(y) == 0) ? nothing : PdfFraud.reshape_representation(y)
end

outfile = joinpath(PdfFraud.proto_path, "extracted_representations/data.jld2")

sha256 = []
representations = []
N = length(protos)
p = Progress(N,1)
for (i, proto) in enumerate(protos)
	x = load_representation(proto)
	if x != nothing
		push!(representations, x)
		push!(sha256, split(basename(proto), ".")[1])
	end
	next!(p)
end

# also get labels
indexf = joinpath(PdfFraud.proto_path, "index.json")
labelf = "/home/vit/vyzkum/bulletproof/pdf_experiment/gt-clustering-pdfs.txt"
labels, labeled_sha256, filenames = PdfFraud.get_labels(sha256, indexf, labelf);

@save outfile sha256 representations labels


text:
import Term: remove_markup,
        has_markup,
        remove_ansi,
        has_ansi,
        get_last_ANSI_code,
        cleantext,
        textlen,
        escape_brackets,
        unescape_brackets,
        replace_text,
        nospaces,
        remove_brackets,
        unspace_commas,
        chars,
        join_lines,
        split_lines,
        textwidth,
        fillin,
        truncate,
        reshape_text
import Term.style: apply_style

@testset "TU_markup" begin
        strings = [
            ("this is [red] some [blue] text [/blue] that I like[/red]",
                "this is  some  text  that I like",
            ), (
                "[bold underline] this is [red on_green] text [/red on_green] I like [/bold underline]",
                " this is  text  I like ",
            )
        ]

        for (s1, s2) in strings
            @test has_markup(s1) == true
            @test remove_markup(s1) == s2
            @test cleantext(s1) == s2
            @test textlen(s1) == textwidth(s2)
        end

        @test remove_markup("text with [[double]] squares") == "text with [[double]] squares"
        @test has_markup("text with [[double]] squares") == false
end

@testset "TU_ansi" begin
    apply_style("test[(.2, .5, .6)]coloooor[/(.2, .5, .6)]")
    strings = [
        (
            "this is \e[31m some \e[34m text \e[39m\e[31m that I like\e[39m",
            "this is  some  text  that I like",
            "\e[39m"
        ), (
            "\e[1m\e[4m this is \e[31m\e[42m text \e[39m\e[49m\e[4m I like \e[22m\e[24m",
            " this is  text  I like ",
            "\e[24m"
        ), (
            "test\e[38;2;51;128;153mcoloooor\e[39m and white",
            "testcoloooor and white",
            "\e[39m"
        )
    ]

    for (s1, s2, ltag) in strings
        @test has_ansi(s1) == true
        @test remove_ansi(s1) == s2
        @test get_last_ANSI_code(s1) == ltag
    end
end

@testset "TU_brackets" begin
    strings = [
        ("test [vec] nn", "test [[vec]] nn"),
        ("[1, 2, 3]", "[[1, 2, 3]]")
    ]

    for (s1, s2) in strings
        escaped = escape_brackets(s1)
        @test escaped == s2
        @test unescape_brackets(escaped) == s1
        @test unescape_brackets(s1) == s1
    end
end

@testset "TU_replace_text" begin
    text = "abcdefghilmnopqrstuvz"

    @test replace_text(text, 0, 5, "aaa") == "aaafghilmnopqrstuvz"
    @test replace_text(text, 0, 5, ',') ==  ",,,,,fghilmnopqrstuvz"

    @test replace_text(text, 18, 21, "aaa") == "abcdefghilmnopqrstaaa"
    @test replace_text(text, 18, 21, ',') ==  "abcdefghilmnopqrst,,,"

    @test replace_text(text, 10, 15, "aaa") == "abcdefghilaaarstuvz"
    @test replace_text(text, 10, 15, ',') ==  "abcdefghil,,,,,rstuvz"

    @test nospaces("a (1, 2, 3) 4") == "a(1,2,3)4"
    @test remove_brackets("aaa (asdsd) BB") == "aaa asdsd BB"

    @test unspace_commas("a, 2, 3") == "a,2,3"
end

@testset "TU_misc" begin
    @test chars("abcd") == ['a', 'b', 'c', 'd']

    strings = [
        "aaa\nadasda\nasdasda",
        """
        asdasd
adsada
asddsa"""
    ]
    for str in strings
        @test join_lines(split_lines(str)) == str
    end
end

@testset "TU_reshape" begin
    str =     """
    Lorem ipsum [bold]dolor sit[/bold] amet, consectetur adipiscing elit,
    ed do e[red]iusmod tempor incididunt[/red] ut [bold]labore et [underline]dolore[/underline] magna aliqua.[/bold] Ut enim ad minim
    veniam, quis[green] nostrud exercitation [on_black]ullamco laboris nisi ut aliquip ex [/on_black]
    ea commodo consequat.[blue] Duis aute irure dolor in[/blue] reprehenderit
    in voluptate velit[/green] esse [italic]cillum dolore[/italic][red] eu[/red][italic green] fugiat [/italic green]nulla
    pariatur. Excepteur[red] sint[/red][blue] occaecat cupidatat [/blue]non proident,
    sunt in culpa qui [italic]officia[/italic] deserunt mollit anim
    id est laborum."""

    str_reshaped = "Lorem ipsum \e[1mdolor sit\e[22m amet, conse\e[0m\nctetur adipiscing elit,\e[0m\ned do e\e[31miusmod tempor incididunt\e[39m\e[0m\nut \e[1mlabore et \e[4mdolore\e[24m\e[1m magna aliqua\e[0m\n\e[1m\e[1m. Ut enim ad minim\e[0m\nveniam, quis\e[32m nostrud exercitation\e[0m\n\e[32m\e[40mullamco laboris nisi ut aliquip\e[0m\n\e[32m\e[40mex \e[49m\e[0m\nea commodo consequat.\e[34m Duis aute\e[0m\n\e[34mirure dolor in\e[39m reprehenderit\e[0m\nin voluptate velit[/green] esse \e[3mc\e[0m\n\e[3millum dolore\e[23m\e[31m eu\e[39m\e[23m\e[3m\e[32m fugiat nulla\e[0m\npariatur. Excepteur\e[31m sint\e[39m\e[34m occaecat\e[0m\n\e[34mcupidatat \e[39mnon proident,\e[0m\nsunt in culpa qui \e[3mofficia\e[23m deserun\e[0m\nt mollit anim\e[0m\nid est laborum."

    strings = [
        (
            "Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua.",
            "Lorem ipsum dolor sit amet, conse\nctetur adipiscing elit, sed do\neiusmod tempor incididunt ut\nlabore et dolore magna aliqua."
        ), (
            "Lorem [red]ipsum dolor sit [underline]amet, consectetur[/underline] adipiscing elit, [/red][blue]sed do eiusmod tempor incididunt[/blue] ut labore et dolore magna aliqua.",
            "Lorem \e[31mipsum dolor sit \e[4mamet, conse\e[0m\n\e[31m\e[4mctetur\e[24m\e[31m adipiscing elit, \e[39m\e[34msed do\e[0m\n\e[31m\e[34meiusmod tempor incididunt\e[39m ut\e[0m\n\e[31mlabore et dolore magna aliqua.\e[0m"
        ), (
            "Lorem[red]ipsumdolorsit[underline]amet, consectetur[/underline] adipiscing elit, [/red]seddoeiusmo[blue]dtemporincididunt[/blue]ut labore et dolore magna aliqua.",
            "Lorem\e[31mipsumdolorsit\e[4mamet, consectet\e[0m\n\e[31m\e[4mur\e[24m\e[31m adipiscing elit, \e[39mseddoeiusmo\e[34mdt\e[0m\n\e[31m\e[34memporincididunt\e[39mut labore et dolor\e[0m\n\e[31me magna aliqua.\e[0m"
        ), (
            "ต้าอ่วยวาทกรรมอาว์เซี้ยว กระดี๊กระด๊า ช็อปซาดิสต์โมจิดีพาร์ตเมนต์ อินดอร์วิว สี่แยกมาร์กจ๊อกกี้ โซนี่บัตเตอร์ฮันนีมูน ยาวีแพลนหงวนสคริปต์ แจ็กพ็อตต่อรองโทรโข่งยากูซ่ารุมบ้า บอมบ์เบอร์รีวีเจดีพาร์ทเมนท์ บอยคอตต์เฟอร์รี่บึมมาราธอน ",
            "ต้าอ่วยวาทกรรมอาว์เซี้ยว กระดี๊กระด๊า\nช็อปซาดิสต์โมจิดีพาร์ตเมนต์ อินดอร์วิว สี่แย\nกมาร์กจ๊อกกี้ โซนี่บัตเตอร์ฮันนีมูน ยาวีแพลนห\nงวนสคริปต์ แจ็กพ็อตต่อรองโทรโข่งยากูซ่ารุม\nบ้า บอมบ์เบอร์รีวีเจดีพาร์ทเมนท์ บอยคอตต์เ\nฟอร์รี่บึมมาราธอน"
        ), (
            "ต้าอ่วยวาท[red]กรรมอาว์เซี้ยว กระดี๊กระด๊า [/red]ช็อปซาดิสต์โมจิดีพาร์ตเม[blue underline]นต์ อินดอร์วิว สี่แยกมาร์กจ๊อกกี้ โซนี่บัตเต[/blue underline]อร์ฮันนีมูน ยาวีแพลนหงวนสคริปต์ แจ็กพ็อตต่อรองโทรโข่งยากูซ่ารุมบ้า บอมบ์เบอร์รีวีเจดีพาร์ทเมนท์ บอยคอตต์เฟอร์รี่บึมมาราธอน ",
            "ต้าอ่วยวาท\e[31mกรรมอาว์เซี้ยว กระดี๊กระด๊า \e[39mช็อป\e[0m\nซาดิสต์โมจิดีพาร์ตเม\e[4m\e[34mนต์ อินดอร์วิว สี่แยกมา\e[0m\n\e[4m\e[34mร์กจ๊อกกี้ โซนี่บัตเตอร์ฮันนีมูน ยาวีแพลนหงวน\e[0m\n\e[4m\e[34mสคริปต์ แจ็กพ็อตต่อรองโทรโข่งยากูซ่ารุมบ้า\e[0m\n\e[4m\e[34mบอมบ์เบอร์รีวีเจดีพาร์ทเมนท์ บอยคอตต์เฟอร์\e[0m\n\e[4m\e[34mรี่บึมมาราธอน\e[0m"
        ), (
            "국가유공자·상이군경 및 전몰군경의 유가족은 법률이 정하는 바에 의하여",
            "국가유공자·상이군경 및 전몰군경의\n유가족은 법률이 정하는 바에\n의하여"
        ), (
            "국[red]가유공자·상이군[bold]경 및 전[/bold]몰군경의 유[/red]가족은 법률이 정하는 바에 의하여",
            "국\e[31m가유공자·상이군\e[1m경 및 전\e[22m\e[31m몰군경의\e[0m\n\e[31m\e[31m유\e[39m가족은 법률이 정하는 바에\e[0m\n\e[31m의하여\e[0m"
        ), (
            "┌────────────────┬────────────────┬────────────────┬────────────────┬──────────────",
            "┌────────────────┬───────────────\n─┬────────────────┬──────────────\n──┬──────────────"
        ), (
            "┌────────────abcde────┬──────────── ────┬────────abcde fghi────────┬────────────────┬──────────────",
            "┌────────────abcde────┬──────────\n── ────┬────────abcde fghi───────\n─┬────────────────┬──────────────"
        ), (
            "┌─────────[red]───ab[/red]cde────┬──────[green]────── ────┬────────abcde fghi[/green]────────┬────────────────┬──────────────",
            "┌─────────\e[31m───ab\e[39mcde────┬──────\e[32m────\e[0m\n\e[32m── ────┬────────abcde fghi\e[39m───────\e[0m\n─┬────────────────┬──────────────\e[0m\n\e[0m"
        ), (
            "┌──────────[red]────[/red]──┬[blue bold]────────────────┬──[/blue bold]──────────────┬────────────────┬──────────────end",
            "┌──────────\e[31m────\e[39m──┬\e[1m\e[34m───────────────\e[0m\n\e[1m\e[34m─┬────────────────┬──────────────\e[0m\n\e[1m\e[34m──┬──────────────end\e[0m"
        ), (
            "."^100  ,
            ".................................\n.................................\n.................................\n."
        ), (
            ".[red]|||[/red]...."^10,
            ".\e[31m|||\e[39m.....\e[31m|||\e[39m.....\e[31m|||\e[39m.....\e[31m|||\e[39m....\e[0m\n.\e[31m|||\e[39m.....\e[31m|||\e[39m.....\e[31m|||\e[39m.....\e[31m|||\e[39m.....\e[31m\e[0m\n\e[31m|||\e[39m.....\e[31m|||\e[39m....\e[0m"
        ), (
            ".|||...."^10,
            ".|||.....|||.....|||.....|||.....\n|||.....|||.....|||.....|||.....|\n||.....|||...."
        ), (
            str,
            str_reshaped
        )
    ]

    for (s1, s2) in strings
        reshaped = reshape_text(s1, 33)
        @test reshaped == s2
    end

    for width in 33:100
        rh = reshape_text(str, width)
        @test any(textlen.(split(rh, "\n")) .> width) == false
    end
end