text:
module DivGradCurls

export Gradient
export Divergence
export Curl

using ..Common
using LinearAlgebra
using SparseArrays

struct Gradient{N}
    components::NTuple{N,SparseMatrixCSC{Float64,Int}}
end

Gradient(components::Vararg{T,N}) where {T<:AbstractArray,N} = Gradient{N}(components)

function Gradient(sim::Simulation{1,Common.Symmetric})
    dx = sim.dx
    n = length(sim)
    m = n+1
    ∇ = spdiagm(m,n,0=>vcat(2/dx,fill(1/dx,n-1)),-1=>vcat(fill(-1/dx,n-1),-2/dx))
    return Gradient(∇)
end

function Gradient(sim::Simulation{2,Common.Symmetric})
    dx = sim.dx
    dy = sim.dy
    nx, ny = sim.latticesize
    mx = nx+1; my = ny+1
    ∇ₓ = spdiagm(mx,nx,0=>vcat(2/dx,fill(1/dx,nx-1)),-1=>vcat(fill(-1/dx,nx-1),-2/dx))
    ∇ᵤ = spdiagm(my,ny,0=>vcat(2/dy,fill(1/dy,ny-1)),-1=>vcat(fill(-1/dy,ny-1),-2/dy))
    return Gradient(kron(sparse(I,ny,ny),∇ₓ),kron(∇ᵤ,sparse(I,nx,nx)))
end

struct Divergence{N}
    gradient::Gradient{N}
end

function Divergence(sim::Simulation{1,Common.Symmetric})
    dx = sim.dx
    m = length(sim)
    n = length(sim)+1
    ∇ = spdiagm(m,n,0=>vcat(fill(-1/dx,m)),1=>vcat(fill(1/dx,m)))
    return Divergence(Gradient(∇))
end

function Divergence(sim::Simulation{2,Common.Symmetric})
    dx = sim.dx
    dy = sim.dy
    mx, my = sim.latticesize
    nx = mx+1
    ny = my+1
    ∇ₓ = spdiagm(mx,nx,0=>vcat(fill(-1/dx,mx)),1=>vcat(fill(1/dx,mx)))
    ∇ᵤ = spdiagm(my,ny,0=>vcat(fill(-1/dy,my)),1=>vcat(fill(1/dy,my)))
    return Divergence(Gradient(kron(sparse(I,my,my),∇ₓ),kron(∇ᵤ,sparse(I,mx,mx))))
end

struct Curl{N}

end

function Curl(sim::Simulation{N}) where N

end

################################################################################
# extend mul!

function LinearAlgebra.mul!(Y::NTuple{N,ScalarField{N}}, A::Gradient{N}, B::ScalarField{N}) where N
    for i ∈ 1:N
        mul!(Y[i], A.components[i], B.values)
    end
    return Y
end

function LinearAlgebra.mul!(Y::ScalarField{N}, A::Divergence{N}, B::NTuple{N,ScalarField{N}}) where N
    if N==1
        mul!(Y.values, A.gradient.components[1], B[1].values)
    else
        fill!(Y.values,0)
        for i ∈ 1:N
            mul!(Y.values, A.gradient.components[i], B[i],1,1)
        end
    end
    return Y
end

end # module

using .DivGradCurls


text:
module DataFramesMeta

using Reexport

@reexport using DataFrames

# Basics:
export @with, @where, @orderby, @transform, @by, @combine, @select,
       @transform!, @select!,
       @eachrow, @eachrow!,
       @byrow, @byrow!, @based_on # deprecated

global const DATAFRAMES_GEQ_22 = isdefined(DataFrames, :pretty_table) ? true : false

include("parsing.jl")
include("macros.jl")
include("linqmacro.jl")
include("eachrow.jl")

end # module

ext: jl


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb04ts(datfile, io=stdout)
    ZERO = 0.0e0
    ONE = 1.0e0
    NIN = 5
    NOUT = 6
    NMAX = 200
    LDA = NMAX
    LDB = NMAX
    LDG = NMAX
    LDQ = NMAX
    LDRES = NMAX
    LDU1 = NMAX
    LDU2 = NMAX
    LDV1 = NMAX
    LDV2 = NMAX
    LDWORK = NMAX
    A = Array{Float64,2}(undef, LDA,NMAX)
    B = Array{Float64,2}(undef, LDB,NMAX)
    G = Array{Float64,2}(undef, LDG,NMAX)
    Q = Array{Float64,2}(undef, LDQ,NMAX)
    CSL = Array{Float64,1}(undef, 2*NMAX)
    CSR = Array{Float64,1}(undef, 2*NMAX)
    TAUL = Array{Float64,1}(undef, NMAX)
    TAUR = Array{Float64,1}(undef, NMAX)
    U1 = Array{Float64,2}(undef, LDU1,NMAX)
    U2 = Array{Float64,2}(undef, LDU2,NMAX)
    V1 = Array{Float64,2}(undef, LDV1,NMAX)
    V2 = Array{Float64,2}(undef, LDV2,NMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    RES = Array{Float64,2}(undef, LDRES,5*NMAX)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    TRANA =  vs[2][1]
    TRANB =  vs[3][1]
    if ( N<=0 || N>NMAX )
        @error "illegal N=$N"
    end

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
    RES[1:N,1:N] .= A[1:N,1:N]

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       B[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
# interp call 1

    #FOREIGN.dlacpy!( 'All', N, N, B, LDB, RES(1,N+1), LDRES )
    RES[1:N,N+1:2N] .= B[1:N,1:N]

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       G[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
# interp call 2

    #FOREIGN.dlacpy!( 'All', N, N, G, LDG, RES(1,2*N+1), LDRES )
    RES[1:N,2N+1:3N] .= G[1:N,1:N]

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       Q[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
# interp call 3

    #FOREIGN.dlacpy!( 'All', N, N, Q, LDQ, RES(1,3*N+1), LDRES )
    RES[1:N,3N+1:4N] .= Q[1:N,1:N]
    close(f)

# interp call 4

    INFO = SLICOT.mb04ts!(TRANA, TRANB, N, 1, A, B, G, Q, CSL, CSR, TAUL, TAUR, LDWORK)
    @test INFO == 0
    if ( INFO!=0 )
        @warn "mb04ts returned info=$INFO"
        return
    end
    U1[1:N,1:N] .= A[1:N,1:N]
    U2[1:N,1:N] .= Q[1:N,1:N]
# interp call 5

    INFO = SLICOT.mb04wr!('U', TRANA, N, 1, U1, U2, CSL, TAUL)
    @test INFO == 0
    if ( INFO!=0 )
        @warn "mb04wr returned info=$INFO"
        return
    end

    V2[1:N,1:N] .= Q[1:N,1:N]
    V1[1:N,1:N] .= B[1:N,1:N]
# interp call 6

    INFO = SLICOT.mb04wr!('V', TRANB, N, 1, V1, V2, CSR, TAUR)
    @test INFO == 0
    if ( INFO!=0 )
        @warn "mb04wr returned info=$INFO"
        return
    end

# interp output 1
    _nc = N
    _nr = N
    println(io, "U1:")
    show(io, "text/plain", U1[1:_nr,1:_nc])
    println(io)
    println(io, "U2:")
    show(io, "text/plain", U2[1:_nr,1:_nc])
    println(io)
    if ( LSAME( TRANA, 'N' ) )
        resnorm = SLICOT.ma02jd!(false,false,N,U1,U2)
        println(io, "U orth. residual norm = $resnorm")
        @test resnorm < 1e-12
        # unable to translate write statement:
        #  write MA02JD( .FALSE., .FALSE., N, U1, LDU1, U2, LDU2, RES(1,4*N+1), LDRES )
    else
        resnorm = SLICOT.ma02jd!(true,false,N,U1,U2)
        println(io, "U orth. residual norm = $resnorm")
        @test resnorm < 1e-12
        # unable to translate write statement:
        #  write MA02JD( .TRUE., .FALSE., N, U1, LDU1, U2, LDU2, RES(1,4*N+1), LDRES )
    end # if
# interp call 7

    #FOREIGN.dlaset!( 'All', N, N, ZERO, ZERO, Q, LDQ )
    Q[1:N,1:N] .= ZERO

    if ( LSAME( TRANA, 'N' ) )
# interp call 8

        #FOREIGN.dlaset!( 'Lower', N-1, N-1, ZERO, ZERO, A(2,1), LDA )
        triu!(view(A,1:N,1:N-1))

        # interp output 3
        _nc = N
        _nr = N
        println(io, "A:")
        show(io, "text/plain", A[1:_nr,1:_nc])
        println(io)
        println(io, "G:")
        show(io, "text/plain", G[1:_nr,1:_nc])
        println(io)

    else
        # interp call 9

        #FOREIGN.dlaset!( 'Upper', N-1, N-1, ZERO, ZERO, A(1,2), LDA )
        tril!(view(A,1:N,1:N-1))

        # interp output 4
        _nc = N
        _nr = N
        println(io, "A:")
        show(io, "text/plain", A[1:_nr,1:_nc])
        println(io)
        println(io, "G:")
        show(io, "text/plain", G[1:_nr,1:_nc])
        println(io)
    end # if
    if ( LSAME( TRANB, 'N' ) )
        if ( N>1 )
            # interp call 10

            #FOREIGN.dlaset!( 'Upper', N-2, N-2, ZERO, ZERO, B(1,3), LDB )
            tril!(view(B,1:N-1,2:N))

        end # if
# interp output 5
        _nc = N
        _nr = N
        println(io, "Q:")
        show(io, "text/plain", Q[1:_nr,1:_nc])
        println(io)
        println(io, "B:")
        show(io, "text/plain", B[1:_nr,1:_nc])
        println(io)

    else
        if ( N>1 )
            # interp call 11

            #FOREIGN.dlaset!( 'Lower', N-2, N-2, ZERO, ZERO, B(3,1), LDB )
            triu!(view(B,2:N-1,1:N-2))

        end # if
        # interp output 6
        _nc = N
        _nr = N
        println(io, "Q:")
        show(io, "text/plain", Q[1:_nr,1:_nc])
        println(io)
        println(io, "B:")
        show(io, "text/plain", B[1:_nr,1:_nc])
        println(io)

    end # if
    if ( LSAME( TRANB, 'N' ) )
        TRANV1 = 'T'
    else
        TRANV1 = 'N'
    end # if
    # interp call 12

    #FOREIGN.dgemm!( TRANA, TRANV1, N, N, N, ONE, RES, LDRES, V1, LDV1, ZERO, RES(1,4*N+1), LDRES )
    BLAS.gemm!(TRANA,TRANV1,ONE,RES[1:N,1:N],V1[1:N,1:N],ZERO,view(RES,1:N,4N+1:5N))

# interp call 13

    #FOREIGN.dgemm!( 'No Transpose', 'Transpose', N, N, N, -ONE, RES(1,2*N+1), LDRES, V2, LDV2, ONE, RES(1,4*N+1), LDRES )
    BLAS.gemm!('N','T',-ONE,RES[1:N,2N+1:3N],V2[1:N,1:N],ONE,view(RES,1:N,4N+1:5N))

# interp call 14

    #FOREIGN.dgemm!( TRANA, TRANA, N, N, N, -ONE, U1, LDU1, A, LDA, ONE, RES(1,4*N+1), LDRES )
    BLAS.gemm!(TRANA,TRANA,-ONE,U1[1:N,1:N],A[1:N,1:N],ONE,view(RES,1:N,4N+1:5N))

    #TEMP = DLANGE( 'Frobenius', N, N, RES(1,4*N+1), LDRES, DWORK )
    TEMP = norm(RES[1:N,4N+1:5N])
# interp call 15

    #FOREIGN.dgemm!( TRANA, 'Transpose', N, N, N, ONE, RES, LDRES, V2, LDV2, ZERO, RES(1,4*N+1), LDRES )
    BLAS.gemm!(TRANA,'T',ONE,RES[1:N,1:N],V2[1:N,1:N],ZERO,view(RES,1:N,4N+1:5N))

# interp call 16

    #FOREIGN.dgemm!( 'No Transpose', TRANV1, N, N, N, ONE, RES(1,2*N+1), LDRES, V1, LDV1, ONE, RES(1,4*N+1), LDRES )
    BLAS.gemm!('N',TRANV1,ONE,RES[1:N,2N+1:3N],V1[1:N,1:N],ONE,view(RES,1:N,4N+1:5N))

# interp call 17

    #FOREIGN.dgemm!( TRANA, 'No Transpose', N, N, N, -ONE, U1, LDU1, G, LDG, ONE, RES(1,4*N+1), LDRES )
    BLAS.gemm!(TRANA,'N',-ONE,U1[1:N,1:N],G[1:N,1:N],ONE,view(RES,1:N,4N+1:5N))

# interp call 18

    #FOREIGN.dgemm!( 'No Transpose', TRANB, N, N, N, -ONE, U2, LDU2, B, LDB, ONE, RES(1,4*N+1), LDRES )
    BLAS.gemm!('N',TRANB,-ONE,U2[1:N,1:N],B[1:N,1:N],ONE,view(RES,1:N,4N+1:5N))

    #TEMP = DLAPY2( TEMP, DLANGE( 'Frobenius', N, N, RES(1,4*N+1), LDRES, DWORK ) )
    TEMP = hypot(TEMP,norm(RES[1:N,4N+1:5N]))
# interp call 19

    #FOREIGN.dgemm!( 'No Transpose', TRANV1, N, N, N, ONE, RES(1,3*N+1), LDRES, V1, LDV1, ZERO, RES(1,4*N+1), LDRES )
    BLAS.gemm!('N',TRANV1,ONE,RES[1:N,3N+1:4N],V1[1:N,1:N],ZERO,view(RES,1:N,4N+1:5N))

# interp call 20

    #FOREIGN.dgemm!( TRANB, 'Transpose', N, N, N, -ONE, RES(1,N+1), LDRES, V2, LDV2, ONE, RES(1,4*N+1), LDRES )
    BLAS.gemm!(TRANB,'T',-ONE,RES[1:N,N+1:2N],V2[1:N,1:N],ONE,view(RES,1:N,4N+1:5N))

# interp call 21

    #FOREIGN.dgemm!( 'No Transpose', TRANA, N, N, N, ONE, U2, LDU2, A, LDA, ONE, RES(1,4*N+1), LDRES )
    BLAS.gemm!('N',TRANA,ONE,U2[1:N,1:N],A[1:N,1:N],ONE,view(RES,1:N,4N+1:5N))

    #TEMP = DLAPY2( TEMP, DLANGE( 'Frobenius', N, N, RES(1,4*N+1), LDRES, DWORK ) )
    TEMP = hypot(TEMP, norm(RES[1:N,4N+1:5N]))
# interp call 22

    #FOREIGN.dgemm!( 'No Transpose', 'Transpose', N, N, N, ONE, RES(1,3*N+1), LDRES, V2, LDV2, ZERO, RES(1,4*N+1), LDRES )
    BLAS.gemm!('N','T',ONE,RES[1:N,3N+1:4N],V2[1:N,1:N],ZERO,view(RES,1:N,4N+1:5N))

# interp call 23

    #FOREIGN.dgemm!( TRANB, TRANV1, N, N, N, ONE, RES(1,N+1), LDRES, V1, LDV1, ONE, RES(1,4*N+1), LDRES )
    BLAS.gemm!(TRANB,TRANV1,ONE,RES[1:N,N+1:2N],V1[1:N,1:N],ONE,view(RES,1:N,4N+1:5N))

# interp call 24

    #FOREIGN.dgemm!( 'No Transpose', 'No Transpose', N, N, N, ONE, U2, LDU2, G, LDG, ONE, RES(1,4*N+1), LDRES )
    BLAS.gemm!('N','N',ONE,U2[1:N,1:N],G[1:N,1:N],ONE,view(RES,1:N,4N+1:5N))

# interp call 25

    #FOREIGN.dgemm!( TRANA, TRANB, N, N, N, -ONE, U1, LDU1, B, LDB, ONE, RES(1,4*N+1), LDRES )
    BLAS.gemm!(TRANA,TRANB,-ONE,U1[1:N,1:N],B[1:N,1:N],ONE,view(RES,1:N,4N+1:5N))

    #TEMP = DLAPY2( TEMP, DLANGE( 'Frobenius', N, N, RES(1,4*N+1), LDRES, DWORK ) )
    TEMP = hypot(TEMP, norm(RES[1:N,4N+1:5N]))
    println(io, "residual H V - U R norm = $TEMP")
    @test TEMP < 1e-12

        # interp output 7
    _nc = N
    _nr = N
    println(io, "V1:")
    show(io, "text/plain", V1[1:_nr,1:_nc])
    println(io)
    println(io, "V2:")
    show(io, "text/plain", V2[1:_nr,1:_nc])
    println(io)

    if ( LSAME( TRANB, 'N' ) )
        # unable to translate write statement:
        #  write MA02JD( .TRUE., .TRUE., N, V1, LDV1, V2, LDV2, RES(1,4*N+1), LDRES )
        resnorm = SLICOT.ma02jd!(true,true,N,V1,V2)
        println(io, "V orth. residual norm = $resnorm")
        @test resnorm < 1e-12
    else

        # unable to translate write statement:
        #  write MA02JD( .FALSE., .TRUE., N, V1, LDV1, V2, LDV2, RES(1,4*N+1), LDRES )
        resnorm = SLICOT.ma02jd!(false,true,N,V1,V2)
        println(io, "V orth. residual norm = $resnorm")
        @test resnorm < 1e-12
    end # if
end # run_mb04ts()


text:
## reader

function Base.parse(::Type{Module}, ir::String, ctx::Context=GlobalContext())
    data = unsafe_wrap(Vector{UInt8}, ir)
    membuf = MemoryBuffer(data, "", false)

    out_ref = Ref{API.LLVMModuleRef}()
    out_error = Ref{Cstring}()
    status = convert(Core.Bool, API.LLVMParseIRInContext(ctx, membuf, out_ref, out_error))

    if status
        error = unsafe_message(out_error[])
        throw(LLVMException(error))
    end

    Module(out_ref[])
end

## writer

Base.string(mod::Module) = unsafe_message(API.LLVMPrintModuleToString(mod))


text:
#=
   Sudoku in Julia + ConstraintSolver

   This is based on the Sudoku example in
   https://wikunia.github.io/ConstraintSolver.jl/stable/tutorial/

   Model created by Hakan Kjellerstrand, hakank@gmail.com
   See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")
include("sudoku_problems.jl")

function sudoku(grid,print_solution=true,timeout=Inf,all_solutions=false)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    # Solver options:
    # https://wikunia.github.io/ConstraintSolver.jl/stable/options/
    # m = Model(optimizer_with_attributes(CS.Optimizer, "all_solutions"=>true,"logging"=>[]))
    m = Model(optimizer_with_attributes(CS.Optimizer,
                                        "all_solutions"=>all_solutions,
                                        "logging"=>[], # :true,

                                        "traverse_strategy"=>:BFS,
                                        # "traverse_strategy"=>:DFS,
                                        # "traverse_strategy"=>:DBFS,

                                        "branch_split"=>:Smallest,
                                        # "branch_split"=>:Biggest,
                                        # "branch_split"=>:InHalf,

                                        # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                        "branch_strategy" => :IMPS, # default
                                        # "branch_strategy" => :ABS, # Activity Based Search
                                        # "activity.decay" => 0.999, # default 0.999
                                        # "activity.max_probes" => 100, # default, 10
                                        # "activity.max_confidence_deviation" => 20, # default 20

                                        # "simplify"=>false,
                                        # "simplify"=>true,
                                        "time_limit"=>timeout,

                                        # "lp_optimizer" => cbc_optimizer,
                                        # "lp_optimizer" => glpk_optimizer,
                                        # "lp_optimizer" => ipopt_optimizer,

                                        ))

    n,_ = size(grid)  # 9
    ns = round(Int, sqrt(n)) # 3
    n2 = ns-1 # 2

    # define the variables
    @variable(m, 1 <= x[1:n,1:n] <= n, Int)

    # hints
    for r=1:n, c=1:n
        if grid[r,c] != 0
            @constraint(m, x[r,c] == grid[r,c])
        end
    end

    # Row and column constraints
    for rc = 1:n
        @constraint(m, x[rc,:] in CS.AllDifferent())
        @constraint(m, x[:,rc] in CS.AllDifferent())
    end

    # Alterntive using transpose (just checking)
    # xt = x' # Transposed version
    # for rc = 1:n
    #    @constraint(m, x[rc,:] in CS.AllDifferent())
    #    @constraint(m, xt[rc,:] in CS.AllDifferent())
    # end

    # Cell constraints
    for i in 1:ns:n, j in 1:ns:n
        @constraint(m, vec([x[i+k,j+l] for k in 0:n2, l in 0:n2]) in CS.AllDifferent())
    end

    # Solve the problem
    optimize!(m)

    status = JuMP.termination_status(m)
    println("status:$status")
    if status == MOI.OPTIMAL
        x = convert.(Integer,JuMP.value.(x))
        if print_solution
            # display(@show x)
            print_grid(x)
        end
        # num_sols = MOI.get(m, MOI.ResultCount())
        # println("num_sols:$num_sols")
        return x
    end
end

#
# Solve the problems in list a
# The problems are in sudoku_problems.jl
#
function check_sudoku(problems, a::Array,print_solutions=true, timeout=Inf)
    println("a:$a")
    for p in a
        println("p:$p")
        p2 = resize_matrix(problems[p])
        @time sudoku(p2,print_solutions,timeout)
    end
end

#
# Solve all problems of size n x n
# The problems are in sudoku_problems.jl
#
function check_sudoku(problems, n::Int,print_solutions=true,timeout=Inf)
    num_problems = 0
    total_time = 0
    for p in sort(collect(keys(problems)),by=k->Symbol(k))
        pp = problems[p]
        len = length(pp)
        if len == n
            num_problems += 1
            println("\n\nproblem $p")
            p2 = resize_matrix(problems[p])
            t = @timed @time sudoku(p2,print_solutions,timeout)
            total_time += t.time
        end
    end
    println("\nnum_problems:$num_problems total_time:$total_time")
end

#
# Solve Norvig's 95 Sudoku's from http://norvig.com/top95.txt
#
# 0.47s
function norvig_problems(print_solutions=false)
   file = "top95.txt"
   # Try to download the file if not exist
   if !isfile(file)
       println("$file does not exist. Trying to download it from http://norvig.com/top95.txt")
       tmp = download("http://norvig.com/top95.txt")
       mv(tmp,file)
       println("downloaded $file")
   end
   if !isfile(file)
       println("$file does not exist and download() didn't work. Please download from http://norvig.com/top95.txt")
       return
   end
   lines = readlines(file).|>line->replace(line,"."=>"0")
   c = 0
   for line in lines
       c += 1
       p = reshape(split(line,"").|>i->parse(Int,i),(9,9))
       @time sudoku(p,print_solutions)
   end
   println("Solve $c Sudoku instances.")

end

timeout = 600.0 # seconds. Only needed for the 25x25 instances
println("Timeout: $timeout")
print_solutions=false

#
# Run the 9x9 instances
#
@time check_sudoku(all_sudoku_problems, 9, print_solutions,timeout) # 40 problems, 0.33s

#
# Run the 16x16 instances
#
@time check_sudoku(all_sudoku_problems, 16,print_solutions,timeout) # 55 problems: 0.59s

#
# Run the 25x25 instances
# Note The 25x25 takes a lots of time so we have to use a timeout when running them.
#
#=
   With a timeout of 600s (10 min)
   - #18: Timeout -> 118s  --> 115s (v0.6.7)
   - #19: 146.5s -> 52.9 (v0.6.7)
   - #20: 515.9s -> 108s (v0.6.7)
   - #21-#29: Timeout
   - #30: 582s
   - #31-#33: Timeout
   - #89: 3.8s -> 1.6s
   - #90: 7.3s -> 8.4s (v0.6.7)
=#
# @time check_sudoku(all_sudoku_problems, 25, true,timeout) # test all instances
# @time check_sudoku(all_sudoku_problems, [19,89,90],true, timeout) # just test the simplest problems
# @time check_sudoku(all_sudoku_problems, [18,19,20,21,89,90],true, timeout) # just test the simplest problems
# @time check_sudoku(all_sudoku_problems, [30],true, timeout) # just test the simplest problems

#
# Run Norvig's 95 Sudoku problems: 0.47s
#
# @time norvig_problems(false)
#

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
"""
Extensible Framework for Quantum Algorithm Design for Humans.
"""
module Yao

export 幺

"""
Extensible Framework for Quantum Algorithm Design for Humans.

简单易用可扩展的量子算法设计框架。

幺 means normalized but not orthogonal in Chinese.
"""
const 幺 = Yao

using Reexport
@reexport using YaoBase, YaoArrayRegister, YaoBlocks, YaoSym

include("deprecations.jl")

end # module


text:
println(filesize("input.txt"))
println(filesize("/input.txt"))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    MicrostructureNoise.Prior(; N, α1, β1, αη, βη, Πα, μ0, C0)
    MicrostructureNoise.Prior(; kwargs...)

Struct holding prior distribution parameters.
`N` is the number of bins,
`InverseGamma(α1, β1)` is the prior of `θ[1]` on the first bin,
the prior on the noise variance `η` is `InverseGamma(αη, βη)`,
the hidden state ``X_0`` at start time is `Normal(μ0, C0)`,
and `Πα` is a prior `Distribution` for `α`,
for example `Πα = LogNormal(1., 0.5)`.

Note: All keyword arguments `N, α1, β1, αη, βη, Πα, μ0, C0`
are mandatory.

Example:
```
prior = MicrostructureNoise.Prior(
N = 40, # number of bins

α1 = 0.0, # prior for the first bin
β1 = 0.0,

αη = 0.3, # noise variance prior InverseGamma(αη, βη)
βη = 0.3,

Πα = LogNormal(1., 0.5),
μ0 = 0.0,
C0 = 5.0
)
```
"""
struct Prior
    N

    α1
    β1

    αη
    βη

    Πα

    μ0
    C0

    function Prior(;Π_...)
        Π = Dict(Π_)
        return new(
            Π[:N],
            Π[:α1],
            Π[:β1],
            Π[:αη],
            Π[:βη],
            Π[:Πα],
            Π[:μ0],
            Π[:C0]
            )
    end
end

"""
    MCMC(Π::Union{Prior,Dict}, t, y, α0::Float64, σα, iterations;
        subinds = 1:1:iterations, η0::Float64 = 0.0, printiter = 100,
        fixalpha = false, fixeta = false, skipfirst = false) -> td, θ, ηs, αs, pacc

Run the Markov Chain Monte Carlo procedure for `iterations` iterations,
on data `(t, y)`, where `t` are observation times and `y` are observations.
`α0` is the initial guess for the smoothing parameter `α` (necessary),
`η0` is the initial guess for the noise variance (optional),
and `σα` is the stepsize for the random walk proposal for `α`.

Prints verbose output every `printiter` iteration.

Returns `td, θs, ηs, αs, pacc`,
`td` is the time grid of the bin boundaries,
`ηs`, `αs` are vectors of iterates,
possible subsampled at indices `subinds`,
`θs` is a Matrix with iterates of `θ` rows.
`paccα` is the acceptance probability for the update step of `α`.

`y[i]` is the observation at `t[i]`.

If `skipfirst = true` and `t` and `y` are of equal length,
the observation `y[1]` (corresponding to `t[1]`) is ignored.

If `skipfirst = true` and `length(t) = length(y) + 1`,
`y[i]` is the observation at `t[i + 1]`.

Keyword args `fixalpha`, `fixeta` when set to `true` allow fixing
`α` and `η` at their initial values.
"""
function MCMC(Π::Union{Prior,Dict}, t, y, α0::Float64, σα, iterations; subinds = 1:1:iterations, η0::Float64 = 0.0, printiter = 100,
    fixalpha = false, fixeta = false, skipfirst = false)

    N = Π.N

    α1 = Π.α1
    β1 = Π.β1
    αη = Π.αη
    βη = Π.βη
    Πα = Π.Πα
    μ0 = Π.μ0
    C0 = Π.C0

    Πx0 = Normal(μ0, sqrt(C0))

    if skipfirst
        if length(t) == length(y)
            shift = 0
            @info "skip observation y[1] at t[1] (skipfirst == true)"
        elseif length(t) == length(y) + 1
            shift = 1
        else
            throw(DimensionMismatch("expected length(t) or length(t) - 1 observations (skipfirst == true)"))
        end
    else
        length(t) != length(y) && throw(DimensionMismatch("expected length(t) observations"))
        shift = 0
    end

    n = length(t) - 1 # number of increments

    N = Π.N
    m = n ÷ N

    Πη = InverseGamma(Π.αη, Π.βη)

    # Initialization
    if shift == 1
        x = [y[1]; y]
    else
        x = copy(y)
    end
    η = η0
    α = α0

    θ = zeros(N)
    ζ = zeros(N - 1)

    Z = zeros(N)
    ii = Vector(undef, N) # vector of start indices of increments
    td = zeros(N+1)
    for k in 1:N
        if k == N
            ii[k] = 1+(k-1)*m:n # sic!
        else
            ii[k] = 1+(k-1)*m:(k)*m
        end

        tk = t[ii[k]]
        td[k] = tk[1]
        Z[k] = sum((x[i+1] - x[i]).^2 ./ (t[i+1]-t[i]) for i in ii[k])
        θ[k] = mean(InverseGamma(α1 + length(ii[k])/2, β1 + Z[k]/2))
    end
    td[end] = t[end]

    acc = 0
    αs = Float64[]
    si = 1

    μ = zeros(n+1)
    C = zeros(n+1)
    θs = Any[]
    ηs = Float64[]

    samples = zeros(N, length(subinds))

    for iter in 1:iterations
        # update Zk (necessary because x changes)
        if !(η == 0.0 && fixeta)
            for k in 1:N
                Z[k] = sum((x[i+1] - x[i]).^2 ./ (t[i+1]-t[i]) for i in ii[k])
            end
        end

        # sample chain
        for k in 1:N-1
            ζ[k] = rand(InverseGamma(α + α, (α/θ[k] + α/θ[k+1])))
        end
        for k in 2:N-1
            θ[k] = rand(InverseGamma(α + α + m/2, (α/ζ[k-1] + α/ζ[k] + Z[k]/2)))
        end
        θ[1] = rand(InverseGamma(α1 + α + m/2, β1 + α/ζ[1] + Z[1]/2))
        θ[N] = rand(InverseGamma(α + length(ii[N])/2, α/ζ[N-1] + Z[N]/2))

        if !fixalpha
            # update parameter alpha using Wilkinson II
            α˚ = α + σα*randn()
            while α˚ < eps()
                α˚ = α + σα*randn()
            end

            lq = logpdf(Πα, α)
            lq += (2*(N-1))*(α*log(α) - lgamma(α))
            s = sum(log(θ[k-1]*θ[k]*ζ[k-1]*ζ[k-1]) + (1/θ[k-1] + 1/θ[k])/ζ[k-1]  for k = 2:N)
            lq += -α*s

            lq˚ = logpdf(Πα, α˚)
            lq˚ += (2*(N-1))*(α˚*log(α˚) - lgamma(α˚))
            lq˚ += -α˚*s

            mod(iter, printiter) == 0 && print("$iter \t α ", α˚)
            if rand() < exp(lq˚ - lq)*cdf(Normal(0, σα), α)/cdf(Normal(0, σα), α˚) # correct for support
                acc = acc + 1
                α = α˚
                mod(iter, printiter) == 0 && print("✓")
            end

        end
        if !fixeta
            # update eta
            z = sum((x[i] - y[i - shift])^2 for i in (1 + skipfirst):(n+1))
            η = rand(InverseGamma(αη + n/2, βη + z/2))

            mod(iter, printiter) == 0 && print("\t √η", √(η))

        end
        mod(iter, printiter) == 0 && println()

        # Sample x from Kalman smoothing distribution

        # Forward pass
        if η == 0.0 && fixeta
            # do nothing
        else
            if skipfirst # ignore observation y[1]
                C[1] = C0
                μ[1] = μ0
                μi = μ0
                Ci = C0
            else
                wi = 0.0
                Ki = C0/(C0 + η)
                μi =  μ0 + Ki*(y[1] - μ0) # shift = 0

                Ci = Ki*η
                C[1] = Ci
                μ[1] = μi
            end
            for k in 1:N
                iik = ii[k]
                for i in iik # from 1 to n
                    wi = θ[k]*(t[i+1] - t[i])
                    Ki = (Ci + wi)/(Ci + wi + η) # Ci is still C(i-1)
                    μi =  μi + Ki*(y[i + 1 - shift] - μi)

                    Ci = Ki*η
                    C[i+1] = Ci # C0 is at C[1], Cn is at C[n+1] etc.
                    μ[i+1] = μi
                end
            end
            hi = μi
            Hi = Ci

            # Backward pass
            x[end] = rand(Normal(hi, sqrt(Hi)))

            for k in N:-1:1
                iik = ii[k]
                for i in iik[end]:-1:iik[1] # n to 1
                    wi1 = θ[k]*(t[i+1] - t[i])
                    Ci = C[i]
                    μi = μ[i]

                    Hi = (Ci*wi1)/(Ci + wi1)
                    hi = μi + Ci/(Ci + wi1)*(x[i+1] - μi) # x[i+1] was sampled in previous step
                    x[i] = rand(Normal(hi, sqrt(Hi)))
                end
            end
        end
        if iter in subinds
            push!(αs, α)
            push!(ηs, η)
            samples[:, si] = θ
            si += 1
        end
    end

    td, samples, ηs, αs, round(acc/iterations, digits=3)
end

"""
```
struct Posterior
    post_t # Time grid of the bins
    post_qlow # Lower boundary of marginal credible band
    post_median # Posterior median
    post_qup # Upper boundary of marginal credible band
    post_mean # Posterior mean of `s^2`
    post_mean_root # Posterior mean of `s`
    qu # `qu*100`-% marginal credible band
end
```

Struct holding posterior information for squared volatility `s^2`.
"""
struct Posterior
    post_t
    post_qlow
    post_median
    post_qup
    post_mean
    post_mean_root
    qu
end

"""
    posterior_volatility(td, θs; burnin = size(θs, 2)÷3, qu = 0.90)

Computes the `qu*100`-% marginal credible band for squared volatility `s^2` from `θ`.

Returns `Posterior` object with boundaries of the marginal credible band,
posterior median and mean of `s^2`, as well as posterior mean of `s`.
"""
function posterior_volatility(td, samples; burnin = size(samples, 2)÷3, qu = 0.90)
    p = 1.0 - qu
    A = view(samples, :, burnin:size(samples, 2))
    post_qup = mapslices(v-> quantile(v, 1 - p/2), A, dims=2)
    post_mean = mean(A, dims=2)
    post_mean_root = mean(sqrt.(A), dims=2)
    post_median = median(A, dims=2)
    post_qlow = mapslices(v-> quantile(v,  p/2), A, dims=2)
    Posterior(
        td,
        post_qlow[:],
        post_median[:],
        post_qup[:],
        post_mean[:],
        post_mean_root[:],
        qu
    )
end

"""
    piecewise(t, y, [endtime]) -> t, xx

If `(t, y)` is a jump process with piecewise constant paths and jumps
of size `y[i]-y[i-1]` at `t[i]`, piecewise returns coordinates path
for plotting purposes. The second argument
allows to choose the right endtime of the last interval.
"""
function piecewise(t_, y, tend = t_[end])
    t = [t_[1]]
    n = length(y)
    append!(t, repeat(t_[2:n], inner=2))
    push!(t, tend)
    t, repeat(y, inner=2)
end


text:
function generate_data_path()
    root_path = "/" * joinpath(split(pathof(Garamond), "/")[1:end-2]...)
    data_path = joinpath(root_path, "test", "data", "generated_data_100_samples.csv")
end

function generate_embeddings_path()
    root_path = "/" * joinpath(split(pathof(Garamond), "/")[1:end-2]...)
    embeddings_path = joinpath(root_path, "test", "embeddings", "word2vec", "sample_model.bin")
end

# The names of the configuration functions below
const CONFIG_FUNCTIONS = [:generate_sample_config_1, :generate_sample_config_2]

function generate_sample_config_1()
    text = """
	{
        "data_loader_name": "juliadb_loader",
        "data_loader_arguments": ["$(generate_data_path())"],
        "data_sampler_name": "identity_sampler",
        "id_key": "id",
        "vectors_eltype": "Float32",
        "embedders": [
            {
                "id": "embedder_1",
                "description": "BM25+RP embedder",
                "embeddable_fields": ["RandString", "StringField", "StringField2", "IntField"],
                "stem_words": false,
                "language": "english",
                "vectors": "bm25",
                "vectors_transform": "rp",
                "vectors_dimension": 50,
                "oov_policy" : "large_vector"
            }
        ],
        "searchers": [
            {
                "id": "searcher_1",
                "id_aggregation": "aggid",
                "description": "A searcher using BM25+RP embeddings and naive indexing",
                "enabled": true,
                "indexable_fields": ["RandString", "StringField", "StringField2", "IntField"],
                "data_embedder": "embedder_1",
                "search_index": "naive",
                "score_alpha": 0.4,
                "score_weight": 0.8
            }
        ]
    }
    """
end

function generate_sample_config_2()
    text = """
    {
        "data_loader_name": "juliadb_loader",
        "data_loader_arguments": ["$(generate_data_path())"],
        "data_sampler_name": "identity_sampler",
        "id_key": "id",
        "vectors_eltype": "Float32",
        "embedders": [
            {
                "id": "embedder_1",
                "description": "BM25+RP embedder",
                "embeddable_fields": ["RandString", "StringField", "StringField2", "IntField"],
                "stem_words": false,
                "language": "english",
                "vectors": "bm25",
                "vectors_transform": "rp",
                "vectors_dimension": 50,
                "oov_policy" : "large_vector"
            },
            {
                "id": "embedder_2",
                "description": "Word2Vec BOE embedder",
                "stem_words": false,
                "language": "english",
                "vectors": "word2vec",
                "vectors_transform": "none",
                "embeddings_path": "$(generate_embeddings_path())",
                "embeddings_kind": "binary",
                "doc2vec_method": "boe",
                "embedder_kwarguments": {},
                "oov_policy" : "large_vector"
            }
        ],
        "searchers": [
            {
                "id": "searcher_1",
                "id_aggregation": "aggid",
                "description": "A searcher using BM25+RP embeddings and naive indexing",
                "enabled": true,
                "indexable_fields": ["RandString", "StringField", "StringField2", "IntField"],
                "data_embedder": "embedder_1",
                "search_index": "naive",
                "score_alpha": 0.4,
                "score_weight": 0.8
            },
            {
                "id": "searcher_2",
                "id_aggregation": "aggid",
                "description": "A searcher using Word2Vec embeddings and naive indexing",
                "enabled": true,
                "indexable_fields": ["StringField2"],
                "data_embedder": "embedder_2",
                "search_index": "ivfadc",
                "search_index_arguments": [],
                "search_index_kwarguments": {"kc":4, "m":4},
                "score_alpha": 0.4,
                "score_weight": 0.8
            }
        ]
    }
"""
end


text:
# keep this separate because it may change between IDL versions

# some constants from idl_export.h
const IDL_TRUE = convert(Int32, 1)
const IDL_FALSE = convert(Int32, 0)
const IDL_MAX_ARRAY_DIM = 8

# IDL types from idl_export.h
const IDL_TYP_UNDEF =      0
const IDL_TYP_BYTE =       1
const IDL_TYP_INT =        2
const IDL_TYP_LONG =       3
const IDL_TYP_FLOAT =      4
const IDL_TYP_DOUBLE =     5
const IDL_TYP_COMPLEX =    6
const IDL_TYP_STRING =     7
const IDL_TYP_STRUCT =     8
const IDL_TYP_DCOMPLEX =   9
const IDL_TYP_PTR =        10
const IDL_TYP_OBJREF =     11
const IDL_TYP_UINT =       12
const IDL_TYP_ULONG =      13
const IDL_TYP_LONG64 =     14
const IDL_TYP_ULONG64 =    15

# translating IDL/C types to julia
typealias IDL_MEMINT Int
typealias IDL_UMEMINT UInt
typealias UCHAR Cuchar
# NOTE: IDL_ARRAY_DIM is fixed length array IDL_MEMINT[IDL_MAX_ARRAY_DIM] (i.e, Int[8])
typealias IDL_ARRAY_DIM Ptr{IDL_MEMINT}
typealias IDL_ARRAY_FREE_CB Ptr{Void}
typealias IDL_FILEINT Int     # possibly different on Windows
typealias IDL_STRING_SLEN_T Cint
const IDL_STRING_MAX_SLEN = 2147483647   # should you check this?

# /***** IDL_VARIABLE flag values ********/
const IDL_V_CONST =       1
const IDL_V_TEMP =        2
const IDL_V_ARR =         4
const IDL_V_FILE =        8
const IDL_V_DYNAMIC =     16
const IDL_V_STRUCT =      32
const IDL_V_NULL =        64

function idl_type(jl_t)
    # IDL type index from julia type
    t = typeof(jl_t)
    if t <: AbstractArray
        t = eltype(jl_t)
    end
    idl_t = -1
    if t == UInt8
        idl_t = IDL_TYP_BYTE
    elseif t == Int16
        idl_t = IDL_TYP_INT
    elseif t == Int32
        idl_t = IDL_TYP_LONG
    elseif t == Float32
        idl_t = IDL_TYP_FLOAT
    elseif t == Float64
        idl_t = IDL_TYP_DOUBLE
    elseif t == Complex64
        idl_t = IDL_TYP_COMPLEX
    elseif t <: AbstractString
        idl_t = IDL_TYP_STRING
    elseif t == Complex128
        idl_t = IDL_TYP_DCOMPLEX
    elseif t == UInt16
        idl_t = IDL_TYP_UINT
    elseif t == UInt32
        idl_t = IDL_TYP_ULONG
    elseif t == Int64
        idl_t = IDL_TYP_LONG64
    elseif t == UInt64
        idl_t = IDL_TYP_ULONG64
    end
    if idl_t < 0 error("IDL.idl_type: type not found: " * string(t)) end
    return idl_t
end

function jl_type(idl_t)
    # julia type from IDL type index
    jl_t = Any
    if idl_t == IDL_TYP_BYTE
        jl_t = UInt8
    elseif idl_t == IDL_TYP_INT
        jl_t = Int16
    elseif idl_t == IDL_TYP_LONG
        jl_t = Int32
    elseif idl_t == IDL_TYP_FLOAT
        jl_t = Float32
    elseif idl_t == IDL_TYP_DOUBLE
        jl_t = Float64
    elseif idl_t == IDL_TYP_COMPLEX
        jl_t = Complex64
    elseif idl_t == IDL_TYP_STRING
        jl_t = Compat.String
    elseif idl_t == IDL_TYP_DCOMPLEX
        jl_t = Complex128
    elseif idl_t == IDL_TYP_UINT
        jl_t = UInt16
    elseif idl_t == IDL_TYP_ULONG
        jl_t = UInt32
    elseif idl_t == IDL_TYP_LONG64
        jl_t = Int64
    elseif idl_t == IDL_TYP_ULONG64
        jl_t = UInt64
    end
    if jl_t == Any
        error("IDL.jl_type: type not found: " * string(idl_t))
    end
    return jl_t
end

#*************************************************************************************************#
# some IDL types from extern.jl
# sizeof(buf) is max size of IDL_ALLTYPES Union (64x2=128 bits or 16 bytes on all platforms)
typealias IDL_ALLTYPES UInt128
immutable IDL_Variable
    vtype::UCHAR
    flags::UCHAR
    flags2::UCHAR
    buf::IDL_ALLTYPES
end

# works as a fixed length array
immutable IDL_DIMS
    d1::IDL_MEMINT
    d2::IDL_MEMINT
    d3::IDL_MEMINT
    d4::IDL_MEMINT
    d5::IDL_MEMINT
    d6::IDL_MEMINT
    d7::IDL_MEMINT
    d8::IDL_MEMINT
end

dims(d::IDL_DIMS) = (d.d1,d.d2,d.d3,d.d4,d.d5,d.d6,d.d7,d.d8)
dims(d::IDL_DIMS, ndims::Integer) = (d.d1,d.d2,d.d3,d.d4,d.d5,d.d6,d.d7,d.d8)[1:ndims]

immutable IDL_Array
    elt_len::IDL_MEMINT                 # Length of element in char units
    arr_len::IDL_MEMINT                 # Length of entire array (char)
    n_elts::IDL_MEMINT                  # total # of elements
    data::Ptr{UCHAR}                    # ^ to beginning of array data
    n_dim::UCHAR                        # # of dimensions used by array
    flags::UCHAR                        # Array block flags
    file_unit::Cshort                   # # of assoc file if file var
    dim::IDL_DIMS                       # dimensions
    free_cb::IDL_ARRAY_FREE_CB          # Free callback
    offset::IDL_FILEINT                 # Offset to base of data for file var
    data_guard::IDL_MEMINT              # Guard longword
end

immutable IDL_String
    slen::IDL_STRING_SLEN_T             # Length of string, 0 for null
    stype::Cshort                       # type of string, static or dynamic
    s::Ptr{Cchar}                       # Addr of string
    IDL_String() = new(0, 0, Base.unsafe_convert(Ptr{Cchar}, Array(Cchar, IDL_RPC_MAX_STRLEN)))
end

# From idl_rpc.h
const IDL_RPC_MAX_STRLEN = 512		# max string length
immutable IDL_RPC_LINE_S
    flags::Cint
    buf::Ptr{Cchar}
    IDL_RPC_LINE_S() = new(0, Base.unsafe_convert(Ptr{Cchar}, Array(Cchar, IDL_RPC_MAX_STRLEN)))
end

const IDL_TOUT_F_STDERR = 1	# Output to stderr instead of stdout
const IDL_TOUT_F_NLPOST = 4	# Output a newline at end of line


text:
#just a function, no struct
function UMRPR(components::Vector{String}; idealmodel=BasicIdeal,
    alpha = MTAlpha,
    mixing = UMRRule,
    activity = UNIFAC,
    translation=MTTranslation,
    userlocations=String[],
    ideal_userlocations=String[],
    alpha_userlocations = String[],
    mixing_userlocations = String[],
     verbose=false)

     return PR(components;
     idealmodel = idealmodel,
     alpha = alpha,
     mixing=mixing,
     activity = activity,
     translation=translation,
     ideal_userlocations = ideal_userlocations,
     alpha_userlocations = alpha_userlocations,
     mixing_userlocations = mixing_userlocations,
     verbose = verbose)
end
export UMRPR


text:
##############################################################################
##
## Finite Difference
##
##############################################################################

type BansalYaronFDMethod
    byp::BansalYaronProblem

    # algorithm parameters
    invdμ::Float64
    invdσ::Float64

    # grid
    μs::Vector{Float64}
    σs::Vector{Float64}
end

function BansalYaronFDMethod(byp::BansalYaronProblem; μn = 30, σn = 50)

    # create volatility grid. square root process has stationary Gamma distribution
    σmin = 0.0
    α =  2 * byp.κσ / byp.νσ^2
    σmax = quantile(Gamma(α, α), 0.999)
    σs = collect(linspace(σmin, σmax, σn))
    invdσ = (σn - 1)/(σmax - σmin)

    # create drift grid
    σ = sqrt(byp.νμ^2 / (2 * byp.κμ))
    μmin = quantile(Normal(byp.μ, σ), 0.001)
    μmax = quantile(Normal(byp.μ, σ), 0.999)
    μs = collect(linspace(μmin, μmax, μn))
    invdμ = (μn - 1)/(μmax - μmin)

    BansalYaronFDMethod(byp, invdμ, invdσ, μs, σs)
end

function solve(x::Union{Type{Val{:ode}}, Type{Val{:nl}}}, byp::BansalYaronProblem; μn = 30, σn = 50,  kwargs...)
    byfd = BansalYaronFDMethod(byp, μn = μn, σn = σn)
    solution = solve(x, byfd, kwargs...)
    return byfd.μs, byfd.σs, solution
end

##############################################################################
##
## Solve
##
##############################################################################
type Grid{T}
    x::Vector{T}
    I::Int
    J::Int
end
Base.getindex(g::Grid, i, j) = g.x[i + g.I * (j - 1)]

function F!{T}(byfd::BansalYaronFDMethod, y::Vector{T}, ydot::Vector{T})
    byp = byfd.byp ; μs = byfd.μs ; σs = byfd.σs
    μn = length(μs)
    σn = length(σs)
    fy = Grid(y, μn, σn)
    ij = 0
    for σi in 1:σn, μi in 1:μn
        ij += 1
        V = fy[μi, σi]
        V∂σ = zero(T)
        V∂2σ = zero(T)
        if (byp.κσ * (1.0 - σs[σi])) >= 0
            V∂σ = (fy[μi, σi+1] - fy[μi, σi]) * byfd.invdσ
        else
            V∂σ = (fy[μi, σi] - fy[μi, σi-1]) * byfd.invdσ
        end
        if σi == 1
            # does not matter since volatility == 0
            V∂2σ = zero(T)
        elseif σi == σn
            V∂2σ = (fy[μi, σi-1] - fy[μi, σi]) * byfd.invdσ^2
        else
            V∂2σ = (fy[μi, σi+1] + fy[μi, σi-1] - 2 * fy[μi, σi]) * byfd.invdσ^2
        end
        V∂μ = zero(T)
        V∂2μ = zero(T)
        if (byp.κμ * (byp.μ - μs[μi])) >= 0
            V∂μ = (fy[μi+1, σi] - fy[μi, σi]) * byfd.invdμ
        else
            V∂μ = (fy[μi, σi] - fy[μi-1, σi]) * byfd.invdμ
        end
        if μi == 1
            V∂2μ = (fy[μi+1, σi] - fy[μi, σi]) * byfd.invdμ^2
        elseif μi == μn
            V∂2μ = (fy[μi-1, σi] - fy[μi, σi]) * byfd.invdμ^2
        else
            V∂2μ = (fy[μi+1, σi] + fy[μi-1, σi] - 2 * fy[μi, σi]) * byfd.invdμ^2
        end
        ydot[ij] = (byp.ρ * byp.θ * max(V, zero(T))^(1-1/byp.θ)
                    + (- byp.ρ * byp.θ + (1-byp.γ) * μs[μi] - 0.5 * (1-byp.γ) * byp.γ * byp.νD^2 * σs[σi]) * V
                    + (byp.κμ * (byp.μ - μs[μi])) * V∂μ
                    + (byp.κσ * (1.0 - σs[σi])) * V∂σ
                    + (0.5 * byp.νμ^2 * σs[σi]) * V∂2μ
                    + (0.5 * byp.νσ^2 * σs[σi]) * V∂2σ
                    )
    end
    return ydot
end

function solve(::Type{Val{:ode}}, byfd::BansalYaronFDMethod; iterations = 1000, kwargs...)
    oldV = fill(byfd.byp.V, length(byfd.μs) * length(byfd.σs))
    newV = oldV
    distance = Inf
    i = 0
    ydot = deepcopy(oldV)
    while i < iterations
        i += 1
        xout, yout = ODE.ode45((t, y) -> F!(byfd, y, ydot), oldV, [0.0 ; 100.0]; kwargs...)
        newV = yout[end]
        distance = chebyshev(newV, oldV)
        @show i, distance
        if distance < 1e-5
            break
        else
            oldV = newV
        end
    end
    return newV
end

function solve(::Type{Val{:nl}}, byfd::BansalYaronFDMethod; kwargs...)
    oldV = fill(byfd.byp.V, length(byfd.μs) * length(byfd.σs))
    out = nlsolve((y, ydot) -> F!(byfd, y, ydot), oldV, method = :trust_region, show_trace = true, autodiff = true, kwargs...)
    return out.zero
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related to lighting condition analysis.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/lighting_condition.jl
# ==============================================================================

# Function: lighting_condition
# ------------------------------------------------------------------------------

@testset "Function lighting_condition" begin
    jd = date_to_jd(2021, 12, 8, 0, 0, 25)
    s_i = sun_position_i(jd)

    r_i = SVector(1.06746e6, 3.22132e6, -6.2841e6)
    lc = lighting_condition(r_i, s_i)
    @test lc === :sunlight

    r_i = SVector(2.52449e6, 4.89401e6, -4.54648e6)
    lc = lighting_condition(r_i, s_i)
    @test lc === :penumbra

    r_i = SVector(3.9575e6, 5.91988e6, -489577.0)
    lc = lighting_condition(r_i, s_i)
    @test lc === :umbra

    lc = lighting_condition(s_i / norm(s_i) * 7000e3, s_i)
    @test lc === :sunlight
end


text:
export SeisMWNI,
SeisPOCS,
SeisSpitzFX
include("SeisMWNI.jl")
include("SeisPOCS.jl")
include("SeisSpitzFX.jl")


text:
function dummycoding(df)
    T = Union{Number,Missing}
    tf = mapreduce(vcat, names(df)) do name
        col = df[!, name]
        vals = unique(col)
        if !(eltype(col) <: T)
            @. name => ByRow(isequal(vals)) .=> Symbol(name * "_", vals)
        else
            Pair[]
        end
    end
    removed = filter(names(df)) do name
        col = df[!, name]
        !(eltype(col) <: T)
    end
    select(transform(df, tf...), Not(removed))
end


text:
# Precompose with diagonal scaling and translation

export PrecomposeDiagonal

"""
**Precomposition with diagonal scaling/translation**

    PrecomposeDiagonal(f, a, b)

Returns the function
```math
g(x) = f(\\mathrm{diag}(a)x + b)
```
where ``f`` is a convex function. Furthermore, ``f`` must be separable, or `a` must be a scalar, for the `prox` of ``g`` to be computable. Parametes `a` and `b` can be arrays of multiple dimensions, according to the shape/size of the input `x` that will be provided to the function: the way the above expression for ``g`` should be thought of, is `g(x) = f(a.*x + b)`.
"""
struct PrecomposeDiagonal{T <: ProximableFunction, R <: Union{Real, AbstractArray}, S <: Union{Real, AbstractArray}} <: ProximableFunction
    f::T
    a::R
    b::S
    function PrecomposeDiagonal{T,R,S}(f::T, a::R, b::S) where {T <: ProximableFunction, R <: Union{Real, AbstractArray}, S <: Union{Real, AbstractArray}}
        if !is_convex(f)
            error("`f` must be convex")
        end
        if !(eltype(a) <: Real)
            error("`a` must have real elements")
        end
        if any(a == 0.0)
            error("elements of `a` must be nonzero")
        else
            new(f, a, b)
        end
    end
end

is_separable(f::PrecomposeDiagonal) = is_separable(f.f)
is_prox_accurate(f::PrecomposeDiagonal) = is_prox_accurate(f.f)
is_convex(f::PrecomposeDiagonal) = is_convex(f.f)
is_set(f::PrecomposeDiagonal) = is_set(f.f)
is_singleton(f::PrecomposeDiagonal) = is_singleton(f.f)
is_cone(f::PrecomposeDiagonal) = is_cone(f.f)
is_affine(f::PrecomposeDiagonal) = is_affine(f.f)
is_smooth(f::PrecomposeDiagonal) = is_smooth(f.f)
is_quadratic(f::PrecomposeDiagonal) = is_quadratic(f.f)
is_generalized_quadratic(f::PrecomposeDiagonal) = is_generalized_quadratic(f.f)
is_strongly_convex(f::PrecomposeDiagonal) = is_strongly_convex(f.f)

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

PrecomposeDiagonal(f::T, a::R, b::S=0.0) where {T <: ProximableFunction, R <: AbstractArray, S <: Real} = PrecomposeDiagonal{T, R, S}(f, a, b)

PrecomposeDiagonal(f::T, a::R, b::S) where {T <: ProximableFunction, R <: Union{AbstractArray, Real}, S <: AbstractArray} = PrecomposeDiagonal{T, R, S}(f, a, b)

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

function gradient!(y::AbstractArray{T}, g::PrecomposeDiagonal, x::AbstractArray{T}) where T <: RealOrComplex
    z = g.a .* x .+ g.b
    v = gradient!(y, g.f, z)
    y .*= g.a
    return v
end

function prox!(y::AbstractArray{T}, g::PrecomposeDiagonal, x::AbstractArray{T}, gamma::Union{R, AbstractArray{R}}=R(1)) where {R <: Real, T <: RealOrComplex{R}}
    z = g.a .* x .+ g.b
    v = prox!(y, g.f, z, (g.a .* g.a) .* gamma)
    y .-= g.b
    y ./= g.a
    return v
end

function prox_naive(g::PrecomposeDiagonal, x::AbstractArray{T}, gamma::Union{R, AbstractArray{R}}=R(1)) where {R <: Real, T <: RealOrComplex{R}}
    z = g.a .* x .+ g.b
    y, fy = prox_naive(g.f, z, (g.a .* g.a) .* gamma)
    return (y .- g.b)./g.a, fy
end

fun_name(f::PrecomposeDiagonal) = string("Precomposition by affine diagonal mapping of ", fun_name(f.f))
fun_dom(f::PrecomposeDiagonal) = fun_dom(f.f)
fun_expr(f::PrecomposeDiagonal) = "x ↦ f(diag(a)*x + b)"
fun_params(f::PrecomposeDiagonal) = string("f(x) = ", fun_expr(f.f), ", a = ", length(f.a) == 1 ? string(f.a[1]) : string(typeof(f.a)), ", b = ", length(f.b) == 1 ? string(f.b[1]) : string(typeof(f.b)))


text:
@testset "dual_sets.jl" begin
    greater_than_Float64 = MOI.GreaterThan(0.0)
    greater_than_Int64 = MOI.GreaterThan(0)
    less_than_Float64 = MOI.LessThan(0.0)
    less_than_Int64 = MOI.LessThan(0)
    nonpositives_3 = MOI.Nonpositives(3)
    nonnegatives_3 = MOI.Nonnegatives(3)
    zeros_3 = MOI.Zeros(3)
    nonpositives_4 = MOI.Nonpositives(4)
    nonnegatives_4 = MOI.Nonnegatives(4)
    zeros_4 = MOI.Zeros(4)

    # GreaterThan
    @test MOI.dual_set(MOI.GreaterThan(0.0)) == greater_than_Float64
    @test MOI.dual_set(MOI.GreaterThan(1.0)) == greater_than_Float64
    @test MOI.dual_set(MOI.GreaterThan(-1.0)) == greater_than_Float64
    @test MOI.dual_set(MOI.GreaterThan(0)) == greater_than_Int64
    @test MOI.dual_set(MOI.GreaterThan(1)) == greater_than_Int64
    @test MOI.dual_set(MOI.GreaterThan(-1)) == greater_than_Int64

    # LessThan
    @test MOI.dual_set(MOI.LessThan(0.0)) == less_than_Float64
    @test MOI.dual_set(MOI.LessThan(1.0)) == less_than_Float64
    @test MOI.dual_set(MOI.LessThan(-1.0)) == less_than_Float64
    @test MOI.dual_set(MOI.LessThan(0)) == less_than_Int64
    @test MOI.dual_set(MOI.LessThan(1)) == less_than_Int64
    @test MOI.dual_set(MOI.LessThan(-1)) == less_than_Int64

    # EqualTo
    @test MOI.dual_set(MOI.EqualTo(0.0)) == nothing
    @test MOI.dual_set(MOI.EqualTo(1.0)) == nothing
    @test MOI.dual_set(MOI.EqualTo(-1.0)) == nothing
    @test MOI.dual_set(MOI.EqualTo(0)) == nothing
    @test MOI.dual_set(MOI.EqualTo(1)) == nothing
    @test MOI.dual_set(MOI.EqualTo(-1)) == nothing

    # Nonpositives
    @test MOI.dual_set(nonpositives_3) == nonpositives_3
    @test MOI.dual_set(nonpositives_4) == nonpositives_4

    # Nonnegatives
    @test MOI.dual_set(nonnegatives_3) == nonnegatives_3
    @test MOI.dual_set(nonnegatives_4) == nonnegatives_4

    # Zeros
    @test MOI.dual_set(zeros_3) == MOI.Reals(3)
    @test MOI.dual_set(zeros_4) == MOI.Reals(4)

    #SOC
    soc = MOI.SecondOrderCone(2)
    soc3 = MOI.SecondOrderCone(3)
    @test MOI.dual_set(soc) == soc
    @test MOI.dual_set(soc) != soc3
    @test MOI.dual_set(soc3) == soc3

    #RSOC
    rsoc = MOI.RotatedSecondOrderCone(2)
    rsoc3 = MOI.RotatedSecondOrderCone(3)
    @test MOI.dual_set(rsoc) == rsoc
    @test MOI.dual_set(rsoc) != rsoc3
    @test MOI.dual_set(rsoc3) == rsoc3

    #PSD
    psd = MOI.PositiveSemidefiniteConeTriangle(2)
    psd3 = MOI.PositiveSemidefiniteConeTriangle(3)
    @test MOI.dual_set(psd) == psd
    @test MOI.dual_set(psd) != psd3
    @test MOI.dual_set(psd3) == psd3

    # Exponential
    exp = MOI.ExponentialCone()
    dual_exp = MOI.DualExponentialCone()
    @test MOI.dual_set(exp) == dual_exp
    @test MOI.dual_set(exp) != exp
    @test MOI.dual_set(dual_exp) == exp
    @test MOI.dual_set(dual_exp) != dual_exp

    # Power
    pow = MOI.PowerCone(0.3)
    pow4 = MOI.PowerCone(0.4)
    dual_pow = MOI.DualPowerCone(0.3)
    @test MOI.dual_set(pow) == dual_pow
    @test MOI.dual_set(pow) != pow
    @test MOI.dual_set(dual_pow) == pow
    @test MOI.dual_set(dual_pow) != pow4
    @test MOI.dual_set(dual_pow) != dual_pow
end


text:
function mm_dd3(t::Triangulation, simplices, convex_coeffs)

    # Some constants used throughout the funciton
    n_simplices = size(t.simplex_inds, 1)
    dim = size(t.points, 2)

    #=
    # update number of points in case a regular grid of points was employed
    # (if so, because the number of subsimplices generated by the
    # shape-preserving splitting depends only on the dimension of the space,
    # there will be more points than we asked for).
    =#
    n_coeffs = maximum(size(convex_coeffs))

    # Pre-allocated arrays (SizedArrays, for efficiency)
    pt          = Size(dim)(zeros(Float64, dim))
    s_arr       = Size((dim+1)^2)(zeros(Float64, (dim+1)^2))
    signs       = Size(dim + 1)(zeros(Float64, dim + 1))

    # Re-arrange simplices so that look-up is a bit more efficient
    simplex_arrs = Vector{Array{Float64, 2}}(n_simplices)
    imsimplex_arrs = Vector{Array{Float64, 2}}(n_simplices)
    for i in 1:n_simplices
        simplex_arrs[i] = t.points[t.simplex_inds[i, :], :]
        imsimplex_arrs[i] = t.impoints[t.simplex_inds[i, :], :]
    end

    # The Markov matrix
    M = zeros(Float64, n_simplices, n_simplices)

    for i in 1:n_simplices
        inds = potentially_intersecting_simplices2(t, i)
        is = imsimplex_arrs[i]

        @views  for k in 1:n_coeffs
            pt = convex_coeffs[:, k].' * is

            #@show i, k, pt
            for j in inds
                sx = simplices[:, j]
                if contains_point_lessalloc!(signs, s_arr, sx, pt, dim)
                    M[j, i] += 1.0
                end
            end

        end
    end
    return M.' ./ n_coeffs
end


text:
include("am.jl")

using Plots
gr()

function SolveEquilibrium(r0; verbose=false)
    am = Household(r = r0,
                          α = 0.33,
                          δ = 0.05,
                          σ = 3/2,
                          κ = 1.00,
                          ν = 2/3,
                          β = 0.97,
                          MC = MarkovChain([0.5 0.5; 0.5 0.5], [0.95; 1.05]),
                          b = 0.0,
                          a_max = 45,
                          a_size = 250)

    @unpack R, α, δ, r, r_ss, w, b, σ, β, Π, z_vals, a_size, a_vals, z_size, duc, duh, duc_inv, duh_inv, u = am
    V, g_a, g_c, g_h = initialize(am)

    g_c = compute_fixed_point(c -> EGMOperator_Labor(c, am),
                            g_c,
                            max_iter=1000,
                            err_tol=1e-4,
                            verbose=false,
                            print_skip = 100)
    Random.seed!(42)
    a_sim, K, L, r1 = SimulationAssets(am, g_c, labor=true)
    residual = (r1 - r0)/r0
    if verbose
        println("Solving Aiyagari Model for r = $r:")
        println("Endogeous Grid Method. Grid size = $a_size")
        println("Given r0 = $r0, then K = $K and L = $L, which implies r1 = $r1")
        println("Diff = $residual")
    end
    return residual, am, a_sim, g_c
end

r_star = find_zero(SolveEquilibrium, (0.01, 0.05))
_, am, a_sim, g_c = SolveEquilibrium(r_star)
histogram(a_sim, nbins = 100, leg = false, normed = true, xlabel = "assets")
plot(am.a_vals, g_c)


text:
using DataFrames, CSV, StatsPlots

ProjDir = @__DIR__

arm_log_0 = CSV.read(joinpath(ProjDir, "results", "arm_log_0_df.csv"), DataFrame)
arm_log_1 = CSV.read(joinpath(ProjDir, "results", "arm_log_1_df.csv"), DataFrame)
intel_log_0 = CSV.read(joinpath(ProjDir, "results", "intel_log_0_df.csv"), DataFrame)
intel_log_1 = CSV.read(joinpath(ProjDir, "results", "intel_log_1_df.csv"), DataFrame)
intel_tbb_log_0 = CSV.read(joinpath(ProjDir, "results", "intel_tbb_log_0_df.csv"), DataFrame)
intel_tbb_log_1 = CSV.read(joinpath(ProjDir, "results", "intel_tbb_log_1_df.csv"), DataFrame)

# Sort on (:num_threads,:num_chains0 instead on (:num_threads, :num_cpp_chains)
#sort(arm_log_1, [:num_threads, :num_chains])

function timings_plot(df::DataFrame; model = "1", ylim=(0, 250))
    colors = [:darkred, :darkblue, :darkgreen, :black]
    fig1 = plot(; xlim=(0, 9), ylim=ylim,
    xlab="c++ num_threads", ylab="elapsed time [s]")
    for (indx, i) in enumerate([1, 2, 4])
        dft = df[df.num_julia_chains .== i .&& df.num_cpp_chains .== 1, :]
        println(dft)
        marksym = :cross
        plot!(dft.num_threads, dft.median;
            marker=marksym, lab="$(i) julia]_chains")
        title!("$(model) Julia chains")
    end
    fig2 = plot(; xlim=(0, 9), ylim=ylim,
    xlab="c++ num_threads", ylab="elapsed time [s]")
    for (indx, i) in enumerate([1, 2, 4])
        dft = df[df.num_cpp_chains .== i .&& df.num_julia_chains .== 1, :]
        println(dft)
        marksym = :cross
        plot!(dft.num_threads, dft.median;
            marker=marksym, lab="$(i) cpp_chains")
        title!("$(model) C++chains")
    end
    plot(fig1, fig2, layout=(1, 2))
end

f0 = timings_plot(arm_log_0; model="ARM log_0")
savefig(joinpath(ProjDir, "graphs", "arm_log_0.png"))

f1 = timings_plot(arm_log_1; model="ARM log_1")
savefig(joinpath(ProjDir, "graphs", "arm_log_1.png"))

f3 = timings_plot(intel_tbb_log_0; model="Intel TBB log_0",
    ylim=(0, 250))
savefig(joinpath(ProjDir, "graphs", "intel_tbb_log_0.png"))

f4 = timings_plot(intel_tbb_log_1; model="Intel TBB log_1")
savefig(joinpath(ProjDir, "graphs", "intel_tbb_log_1.png"))

f5 = timings_plot(intel_log_0; model="Intel log_0",
    ylim=(0, 250))
savefig(joinpath(ProjDir, "graphs", "intel_log_0.png"))

f6 = timings_plot(intel_log_1; model="Intel log_1",
    ylim=(0, 250))
savefig(joinpath(ProjDir, "graphs", "intel_log_1.png"))


text:
module SDPA

using CxxWrap

# This 'using' is required to suppress a warning about SDPA not having Libdl in its
# dependencies (Libdl is used by BinaryProvider), e.g.: bicycle1885/CodecZlib.jl#26.
using Libdl

import SDPA_jll: libsdpawrap_path, libsdpawrap
@wrapmodule(libsdpawrap_path)

function __init__()
    @initcxx
end

include("blockdiag.jl")
include("blockmat.jl")
include("options.jl")
include("MOI_wrapper.jl")

end # module


text:
include("build_h.jl")

if isfile("version_git.jl")
    include("version_git.jl")
end

module common

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
# This file is not an active part of the package. This is the code
# that uses the Clang.jl package to wrap Sundials using the headers.

ENV["JULIAHOME"] = "/Users/dls/jl/julia/"

using Clang.cindex
using Clang.wrap_c

if (!has(ENV, "JULIAHOME"))
  error("Please set JULIAHOME variable to the root of your julia install")
end

clang_includes = map(x->joinpath(ENV["JULIAHOME"], x), [
  "deps/llvm-3.3/build_Release/Release/lib/clang/3.3/include",
  "deps/llvm-3.3/include",
  "deps/llvm-3.3/include",
  "deps/llvm-3.3/build_Release/include/",
  "deps/llvm-3.3/include",
])

clang_extraargs = ["-D", "__STDC_LIMIT_MACROS", "-D", "__STDC_CONSTANT_MACROS"]

header_path = "/Users/dls/jl/openssl-1.0.1e/include/openssl/"
headers_to_wrap = map(x -> joinpath(header_path, x), ["aes.h", "des.h", "md5.h", "md4.h", "mdc2.h", "sha.h", "rand.h", "ecdh.h", "ecdsa.h"])
@show headers_to_wrap

wc = wrap_c.init(CommonFile="libcrypto_common.jl",
                 ClangArgs=clang_extraargs,
                 ClangIncludes=clang_includes,
                 header_wrapped= (th, h) -> contains(headers_to_wrap, h),
                 header_library= h -> "libcrypto",
                 header_outputfile= h -> joinpath("src", "gen", last(split(h, "/")) * ".jl"))

wc.options.wrap_structs = true

wrap_c.wrap_c_headers(wc, map(ascii, headers_to_wrap))


text:
# ======================================================================
# Logarithmic norms of matrices (a.k.a matrix measures)
#
# Reference:
#
# C. Desoer and M. Vidyasagar, Feedback Systems: Input-Output Properties.
# Philadelphia, PA: Society for Industrial and Applied Mathematics,
# Jan. 2009, ser. Classics in Applied Mathematics.
# ======================================================================

# logarithmic norm (also known as matrix measure) for commonly used p-norms
function logarithmic_norm(A::AbstractMatrix, p::Real=Inf)
    if p == Inf
        return logarithmic_norm_inf(A)
    elseif p == 1
        return logarithmic_norm_1(A)
    elseif p == 2
        return logarithmic_norm_2(A)
    else
        throw(ArgumentError("logarithmic norm only implemented for p = 1, 2, or Inf, got p = $p"))
    end
end

# max_j a_jj + ∑_{i ≠ j} |a_ij|
function logarithmic_norm_1(A::AbstractMatrix{N}) where {N}
    out = -Inf
    @inbounds for j in 1:size(A, 2)
        α = A[j, j]
        for i in 1:size(A, 1)
            if i ≠ j
                α += abs(A[i, j])
            end
        end
        if α > out
            out = α
        end
    end
    return out
end

# max_i a_ii + ∑_{j ≠ i} |a_ij|
function logarithmic_norm_inf(A::AbstractMatrix{N}) where {N}
    out = -Inf
    @inbounds for i in 1:size(A, 1)
        α = A[i, i]
        for j in 1:size(A, 2)
            if i ≠ j
                α += abs(A[i, j])
            end
        end
        if α > out
            out = α
        end
    end
    return out
end

# max_j  1/2 * λⱼ(A + A^T)
function logarithmic_norm_2(A::AbstractMatrix)
    B = A + A'
    λ = eigvals(B)
    return maximum(λ) / 2
end

# ========================================================================
# Specialized functions on LazySets.Arrays.SingleEntryVector
# TODO refactor to LazySets.jl
# ========================================================================

# difference between single entry vectors (SEV); it is type unstable because
# it may return either a SEV or a regular Vector depending on the indices
# but it is much faster than the fallback if the indices match
function minus(e1::SingleEntryVector{N}, e2::SingleEntryVector{N}) where {N}
    e1.n == e2.n || throw(DimensionMismatch("dimensions must match, but they are $(length(e1)) and $(length(e2)) respectively"))

    if e1.i == e2.i
        return SingleEntryVector(e1.i, e1.n, e1.v - e2.v)
    else
        out = zeros(N, e1.n)
        @inbounds begin
            out[e1.i] = e1.v
            out[e2.i] = -e2.v
        end
        return out
    end
end

# addition between single entry vectors (SEV); it is type unstable because
# it may return either a SEV or a regular Vector depending on the indices
# but it is much faster than the fallback if the indices match
function plus(e1::SingleEntryVector{N}, e2::SingleEntryVector{N}) where {N}
    e1.n == e2.n || throw(DimensionMismatch("dimensions must match, but they are $(length(e1)) and $(length(e2)) respectively"))

    if e1.i == e2.i
        return SingleEntryVector(e1.i, e1.n, e1.v + e2.v)
    else
        out = zeros(N, e1.n)
        @inbounds begin
            out[e1.i] = e1.v
            out[e2.i] = e2.v
        end
        return out
    end
end

# norm of the difference of two SEV ||x - y|| in the (vector p norm)
function normdiff(e1::SingleEntryVector{N}, e2::SingleEntryVector{N}, p::Real=2) where {N}
    e1.n == e2.n || throw(DimensionMismatch("dimensions must match, but they are $(length(e1)) and $(length(e2)) respectively"))

    if e1.i == e2.i
        δ = e1.v - e2.v
        return abs(δ)

    else
        a = abs(e1.v)
        b = abs(e2.v)
        if isinf(p)
            return max(a, b)
        else
            s = a^p + b^p
            return s^(1/p)
        end
    end
end

# fallback
function normdiff(x::AbstractVector{N}, y::AbstractVector{N}, p::Real=2) where {N}
    norm(x - y, p)
end


text:
@testset "TODO: fix this code" begin
    Class                   = DiagnosisClassification.Class
    generate_ccs9_to_icd9   = DiagnosisClassification.generate_ccs9_to_icd9
    construct_icd9          = DiagnosisClassification.construct_icd9
    construct_ccs9          = DiagnosisClassification.construct_ccs9
    generate_ccs10_to_icd10 = DiagnosisClassification.generate_ccs10_to_icd10
    construct_icd10         = DiagnosisClassification.construct_icd10
    construct_ccs10         = DiagnosisClassification.construct_ccs10

    @testset "generate_ccs9_to_icd9" begin
        ccs9_to_icd9_struct     = generate_ccs9_to_icd9()
        ccs9_to_icd9            = ccs9_to_icd9_struct.ccs9_to_icd9
        ccs9_to_ccs_description = ccs9_to_icd9_struct.ccs9_to_ccs_description
        @test ccs9_to_icd9            isa Dict{Class, Set{Class}};
        @test ccs9_to_ccs_description isa Dict{Class, String};
        @test !isempty(ccs9_to_icd9);
        @test !isempty(ccs9_to_ccs_description);
        @test ccs9_to_ccs_description[construct_ccs9("1")] == "Tuberculosis"
        @test ccs9_to_ccs_description[construct_ccs9("1")] != "Septicemia"
        @test construct_icd9("011.00")  in ccs9_to_icd9[construct_ccs9("1")];
        @test !(construct_icd9("031.2") in ccs9_to_icd9[construct_ccs9("1")]);
    end

    @testset "generate_ccs10_to_icd10" begin
        ccs10_to_icd10_struct    = generate_ccs10_to_icd10()
        ccs10_to_icd10           = ccs10_to_icd10_struct.ccs10_to_icd10
        ccs10_to_ccs_description = ccs10_to_icd10_struct.ccs10_to_ccs_description
        @test ccs10_to_icd10           isa Dict{Class, Set{Class}};
        @test ccs10_to_ccs_description isa Dict{Class, String};
        @test !isempty(ccs10_to_icd10);
        @test !isempty(ccs10_to_ccs_description);
        @test ccs10_to_ccs_description[construct_ccs10("INF001")] == "Tuberculosis"
        @test ccs10_to_ccs_description[construct_ccs10("INF001")] != "Septicemia"
        @test construct_icd10("A15.7") in ccs10_to_icd10[construct_ccs10("INF001")];
        @test !(construct_icd10("A31.2") in ccs10_to_icd10[construct_ccs10("INF001")]);
    end
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ cf07621a-9a65-11ec-1fc0-51fbadbb086f
using Pkg; Pkg.activate()

# ╔═╡ 345a90f7-f8ac-43b2-b5bd-0282d9447556
begin
	using DifferentialEquations
	using ModelingToolkit
	using DynamicalSystems
	using ForwardDiff
	using LinearAlgebra
	using Graphs
	using Plots
	using PlutoUI
	using LaTeXStrings
end

# ╔═╡ eaf9dda2-d898-4771-9f86-ca4be174712c
md"""
# Dynamical Systems in Julia
In this tutorial, we'll looking at specifying, simulating and analyzing dynamical systems in Julia. This is one area where Julia especially shines, thanks in no small part to the `DifferentialEquations` package, which pretty much provides the largest and fastest suite of ODE solvers on the market today.
"""

# ╔═╡ eb6aed5a-0f95-4a76-8a99-b8719984bbe5
md"## Defining systems symbolically"

# ╔═╡ 0e28ede5-0af6-49a7-b95d-7094e320eee4
md"""
We'll start by using ModelingToolkit to define the famous Lorenz system. The process starts by a symbolic description of the problem.
"""

# ╔═╡ 6ccfc87f-40e7-4ee4-baae-ee41b81efe67
begin
	@parameters t σ ρ β
	@variables x(t) y(t) z(t)
	D = Differential(t)

	eqs = [
		D(x) ~ σ*(y-x),
       	D(y) ~ x*(ρ-z)-y,
       	D(z) ~ x*y - β*z
	]

	@named lorenz = ODESystem(eqs)
end

# ╔═╡ 3e4ec586-482f-4c1b-a8e7-4ae88e48fa38
md"""
So far, the `lorenz` is purely symbolic, which has the added benefit of being able to render the equations in LaTeX as above! The next step before simulating it is defining values for our initial condition and parameters.
"""

# ╔═╡ ce0b65b5-e57a-4def-acb6-892c49f182d9
begin
	u0 = [
		x => 0.0,
      	y => 0.5,
      	z => 0.5
	]

	p  = [
		ρ => 28.0,
      	σ => 10.0,
      	β => 8/3
	]
end

# ╔═╡ 1a7e8273-8069-42a8-ab7b-62f5fe1ade67
md"""
Simulating our system is a two step process. After defining our timespan of integration, we first create an `ODEProblem`, we describes an ODE + initial conditions and parameters, and then pass it to the `solve` function to compute the solution of our problem. The second argument to solve is the solver we want to use (`Tsit5` is a 4-5 order Runge-Kutta scheme which is generally a good default choice for non-stiff problems)
"""

# ╔═╡ 1d77a2be-773c-4edd-9d68-1ed84fb5dadf
begin
	tspan = (0.0,100.0)
	prob = ODEProblem(lorenz,u0,tspan,p,jac=true)
	sol = solve(prob,Tsit5())
end

# ╔═╡ 7c712f73-3931-4cf4-82d8-9d9e2263df14
plot(sol, vars=(x,y,z))

# ╔═╡ 432fa70e-9225-4e76-aec3-c27d31887fb9
md"""
NB. There are also packages for modeling specific classes of dynamical systems, such as `Catalyst` for chemical reaction networks, or `PowerDynamics` for power grids.
"""

# ╔═╡ 60063cc3-e71f-461b-be0a-491b1b566718
md"## Defining Systems by hand"

# ╔═╡ ea2eeac2-5514-4fe8-8921-6235cdb0d1ad
md"""
While defining systems symbolically is cool, this can be tedious for large systems, and sometimes you just want full control over how your ODE function is defined. To that end, we will now define another famous system, the Brusselator by hand and then build a network reaction diffusion system with it.

The Brusselator equations are as follows,
```math
\begin{align*}
\dot{x} &= a + x^2y -bx -x \\
\dot{y} &= bx - x^2y
\end{align*}
```
"""

# ╔═╡ 2bf13604-71c6-474f-90de-5b270c5c0708
md"""
We can construct ODEProblems by providing a function that computes the derivative, just like we do in Matlab. The functions we can use come in two flavours. The first type just takes as input the state `u`, parameters stored in `p` and time `t` (for non autonomous systems), and returns the new state.
"""

# ╔═╡ 743b31cd-6b2a-472c-8a70-a3955a2c26b6
function brusselator(u, p, t)
	a, b = p
	[
		a + u[1]^2 * u[2] - b*u[1] - u[1],
		b*u[1] - u[1]^2 * u[2]
	]
end

# ╔═╡ 03bd6ff6-107f-4b04-9cb8-fe5c28570e58
md"""
However, it is always recommended to use the second form, which also takes as input a preallocated array for the output and modifies it. Doing this is more efficient in terms of memory allocation, which is actually an important source of slowdowns.
"""

# ╔═╡ 08fc4616-d2af-471b-adc0-195bb566334f
function brusselator!(du, u, p, t)
	a, b = p
	du[1] = a + u[1]^2 * u[2] - b*u[1] - u[1]
	du[2] = b*u[1] - u[1]^2 * u[2]
end

# ╔═╡ d0c1e91e-d1df-4a12-a0e2-8513eeee7ba6
begin
	u0_b = [1.0, 1.0]
	p_b = [1.0, 3.0]
end

# ╔═╡ 8623c9b7-5724-4678-9b53-4312389a276d
begin
	prob_b = ODEProblem(brusselator!, u0_b, (0.0, 50.0), p_b)
	sol_b = solve(prob_b, Tsit5())
end

# ╔═╡ da0220ea-99aa-4a0d-8826-f0e461fd503e
plot(sol_b)

# ╔═╡ 60d1ea9b-b935-4e88-9e92-07497f06aeaa
md"""
Now that we have our brusselator function, let's use it on a network to make a reaction diffusion system. Such a system is of the form

```math
\dot{u_i} = f(u_i) + D \sum_{j=1}^N A_{ij}(u_j - u_i),
```
where ``u_i\in\mathbb{R}^m`` is the state of node ``i``, subject to local dynamics ``f:\mathbb{R}^m\rightarrow \mathbb{R}^m``, ``A\in\mathbb{R}^{N\times N}`` is the adjacency matrix of the graph (i.e. ``A_{ij}=1`` iff there is a link between nodes ``i`` and ``j``), and ``D`` is a diagonal matrix of diffusion coefficients.

We can rewrite this system in vector form (glossing over some details) as

```math
\dot{\mathbf{u}} = F(\mathbf{u}) + \mathbf{D} L \mathbf{u},
```
where ``L`` is the *Laplacian matrix* of the graph.

Reaction-Diffusion systems can take more general forms, but we'll stick to this one for today.
"""

# ╔═╡ efbb46be-c442-4ab3-85a0-b10cbc3a2a9e
md"""
Now, we need to create a function for our reaction-diffusion system.
"""

# ╔═╡ 8b2ec98a-b37f-4bd0-9fc8-9f3dcf9550c9
function brusselator_rd(u, p, t)
	a, b, L, D = p
	du = similar(u) # creates an array of same shape and element type as u
	for i in 1:size(du, 1)
		du[i,:] .= brusselator(u[i,:], [a, b], t)
	end
	du[:,1] .-= D[1] .* (L*u[:,1])
	du[:,2] .-= D[2] .* (L*u[:,2])
	return du
end

# ╔═╡ 25c22a1f-6ae4-4b72-8e4c-de844e2699d1
function brusselator_rd!(du, u, p, t)
	# NB. the commented out parts should be equivalent to what is used, but for some reason it is not
	a, b, L, D = p
	du .= 0.0
	for i in 1:size(du, 1)
		#brusselator!(du[i,:], u[i,:], [a, b], t)
		du[i,:] .+= brusselator(u[i,:], [a, b], t)
	end
	#mul!(du[:,1], L, u[:,1], -D[1], 0.0) # du[:,1] .+= D[1]*L*u[:,1]
	#mul!(du[:,2], L, u[:,2], -D[2], 0.0)
	du[:,1] .-= D[1] .* (L*u[:,1])
	du[:,2] .-= D[2] .* (L*u[:,2])
end

# ╔═╡ 582b74c7-aaf6-4a8e-901f-489199ccc0d6
md"""
Now, we need a network to run our system on. To do this, we can use the `Graphs` package.
"""

# ╔═╡ 3e575a45-50d6-4297-84e5-6699e36d00dd
g = erdos_renyi(10, 0.4)

# ╔═╡ ed1de34a-a42e-474d-ba36-411163a343f8
md"""
We can use the function `laplacian_matrix` to get the Laplacian matrix of the graph. The nice thing about this function is that it returns a sparse matrix, which uses optimized methods for matrix multiplications.
"""

# ╔═╡ 947c169c-a900-4728-ad29-5200b52f9599
laplacian_matrix(g)

# ╔═╡ 04bbc594-17c9-4b41-a599-d41e2e541c9f
begin
	u_rd_0 = 5*rand(nv(g), 2)
end

# ╔═╡ 29a0f8c8-5529-496c-99f1-b436d15f2517
md"""
Here's a nice feature of Pluto. It provides allows interactively modifying variables with html elements (sliders, buttons, ...) via the `PlutoUI` package.

Let's use it to make sliders for the parameters of our system.
"""

# ╔═╡ f346a0b3-4d1d-49ad-94ec-cc42ee75a05b
md"""
a: $(@bind a Slider(0.0:0.01:5.0, default=1.0, show_value=true))

b: $(@bind b Slider(0.0:0.01:5.0, default=3.0, show_value=true))

D₁: $(@bind D_1 Slider(0.0:0.01:10.0, default=1.0, show_value=true))

D₂ $(@bind D_2 Slider(0.0:0.01:10.0, default=1.0, show_value=true))
"""

# ╔═╡ 8a1c55e5-1008-44ee-b853-c53c2dbf3bb2
p_rd = [a, b, laplacian_matrix(g), [D_1, D_2]];

# ╔═╡ 2ee789c5-ae38-4e85-95e1-920a580df058
begin
	prob_rd = ODEProblem(brusselator_rd!, u_rd_0, (0.0, 50.0), p_rd)
	sol_rd = solve(prob_rd, Tsit5())
	nothing
end

# ╔═╡ 9d1f3408-e0f5-4093-bb6b-862485e99476
begin
	local p1 = plot(sol_rd, vars=(0,1:nv(g)), label="", ylabel=L"x_i")
	local p2 = plot(sol_rd, vars=(0,(nv(g)+1:2nv(g))), label="", ylabel=L"y_i")
	plot(p1, p2)
end

# ╔═╡ a4d4a26f-ff25-4bf0-b8ab-2db56131ee3e
md"""
### Exercise
Implement your favourite dynamical system and play around with it's parameters.
"""

# ╔═╡ 91f1ce37-0d68-4ed2-9359-dbcfd06f6cdd
md"""
## Intermezzo: using `DynamicalSystems` to study chaotic Systems

Before moving to the last part of this tutorial, let us make a detour to look at the `DynamicalSystems` package, which provides a library of famous dynamical systems, as well as generic tools for ananlyzing them, particularly chaotic systems.
"""

# ╔═╡ 4328e814-82a8-4e7d-8b8a-44a9d5ba17c5
ds = Systems.logistic()

# ╔═╡ 7b11ac9b-a470-4059-bbad-ba6bc443c68b
md"""
Here we load the well known logistic map from the library of dynamical systems, and generate its orbit diagram for varying values of the `r` parameter. The magic function is `orbitdiagram`, which should work for pretty much any discrete map.
"""

# ╔═╡ a5629c28-53bf-4b19-a49f-ab26df1d1b46
begin
	i = 1
	pvalues = 1:0.005:4
    ics = [rand() for m in 1:10]
    n = 2000
    Ttr = 2000
    p_index = 1

    output = orbitdiagram(ds, i, p_index, pvalues; n = n, Ttr = Ttr)
end

# ╔═╡ 74adb98a-1ad5-4181-a3be-ec5fe2565a35
begin
	local L = length(pvalues)
	local x = Vector{Float64}(undef, n*L)
	local y = copy(x)
	for j in 1:L
    	x[(1 + (j-1)*n):j*n] .= pvalues[j]
  		y[(1 + (j-1)*n):j*n] .= output[j]
	end
	scatter(x,y,
		xlabel="r",
		ylabel="x",
		label="",
		markersize=0.8,
		color=:black,
		markeralpha=0.05,
		markerstrokealpha=0.0
	)
end

# ╔═╡ c38e202f-a2d2-4449-8308-df22774b7af1
md"""
NB. While `DynamicalSystems` provides a bunch of cool tools, a lot of them seem to inexplicably break when used outside of the examples. For example, try changing the above code to use values of `r` above 4.
"""

# ╔═╡ d5ab9847-ace4-4a5b-8c81-d2eb8b56c896
md"""
## Analyzing the stability of our system

Reaction-Diffusion systems like the one we just saw are interesting because of their asymptotic behaviour.  Because the laplacian always has a zero eigenvalue associated to an eigenvector whose entries are all the same, a reaction-diffusion system always admits a homogeneous solution where the nodes are all synchronised.
By tuning the diffusion coefficients, it's possible to (de)synchronize the states of the nodes.

The standard tool for analysing this phenomenon is the Master Stability Function [Pecora-Caroll. 1998].

The main idea consists in expressing the system in the basis of eigenvectors of the laplacian matrix. By considering a small perturbation around the synchronized solution, and expressing its linearized equations in the eigenbasis of ``L``, one can obtain a set of decoupled equations involving the laplacian eigenvalues and the jacobian of the local dynamics.

Skipping over the details, we can check that the synchronized solution is linearly stable by checking that the following ODE is stable for every eigenvalue ``\lambda`` of the laplacian matrix

```math
\dot{\xi} = [\mathcal{J}f(u_*) - \lambda D] \xi
```
where ``\mathcal{J}f`` denotes the jacobian of ``f``.

This can be done, for example, by computing the Lyapunov exponent of this equation, thereby yielding a function ``MSF(\lambda)`` called the *Master Stability Function* of the system.
"""

# ╔═╡ 575e3106-704e-4bea-8155-7ee6e6c9da5d
md"First, let's compute the eigenvalues of the laplacian matrix. We don't actually need to compute the eigenvectors."

# ╔═╡ bff34b78-f7ad-452e-9b95-a498f3e377f5
λ = laplacian_spectrum(g)

# ╔═╡ bba5eaa1-16f2-401d-8cf3-f05e9761583d
md"""
Next, let's start by writing a function that takes the local dynamics as input, and generates a function for the MSF ODE.

Notice the use, of `ForwardDiff.jacobian` to automagically compute the Jacobian matrix using Automatic Differentiation (this is different from Finite Differences or Symbolic methods)
"""

# ╔═╡ ddbb84d1-945a-4184-850b-31ee3540e960
function msf_ode_fun(f)
	function msf_ode!(du, u, p, t)
		p_f, D, λ = p
		du[:,1] .= f(u[:,1], p_f, t)
		J = ForwardDiff.jacobian(x -> f(x,p_f,t), u[:,1])
		du[:,2] .= J*u[:,2] - λ*D*u[:,2]
	end
	return msf_ode!
end

# ╔═╡ f0215a78-146d-4009-ae73-8a2f35072db9
md"""
Now, we can move on to computing the MSF proper. To do this, we solve the MSF ODE, then estimating the Lyapunov exponent from the time series.

NB. Recall that the Lyapunov exponent for a dynamical sytem is defined by
```math
\lambda = \lim_{t\rightarrow \infty} \lim_{\xi_0 \rightarrow 0} \frac{1}{t} \ln \frac{\|\xi(t)\|}{\|\xi_0\|},
```
where ``\xi(0) = \xi_0`` is a small perturbation applied to the orbit of interest (here, the periodic solution).
"""

# ╔═╡ 84d86795-e2c4-4c2c-888b-ae999bcda2dd
function master_stability_function(f, u_0, p_f, D, t; δ=0.1, k=50, Δt=0.1)
	msf_ode! = msf_ode_fun(f)
	n = length(u_0)
	u_msf_0 = hcat(u_0, δ*randn(n))
	function msf(λ; u_msf_0=u_msf_0, p_f=p_f, D=D, t=t, Δt=Δt)
		p_msf = [p_f, D, λ]
		prob = ODEProblem(msf_ode!, u_msf_0, (0.0,t), p_msf)
		sol = solve(prob, Tsit5(), saveat=0.0:Δt:t, abstol=1e-10, reltol=1e-10)
		ξ_0 = norm(sol.u[1][:,2])
		l = length(sol.u)
		ξ_t = norm([sol.u[k][:,2] for k in l-k:l])
		return log(ξ_t/ξ_0)/sol.t[end]
	end
	return msf
end

# ╔═╡ 8c592351-4da7-4bed-a99c-595156c1ed8b
begin
	u0_msf = sol_b.u[end]
	p_f = [1.0, 3.0]
	D_msf = Diagonal([0.2,5.0])
end

# ╔═╡ e3411cdb-acf7-40da-9f5e-6852c1370797
msf = master_stability_function(brusselator, u0_msf, p_f, D_msf, 200.0)

# ╔═╡ 5fbfe4aa-f2e5-4845-92e7-9de00063059f
msf(0.0)

# ╔═╡ ffecadd1-95b5-412f-86e4-e4dee733cab1
md"""
Note that this implementation is not super fast, as it takes about a minute to run the next cell. There's lots of room for improvements though, and in the mean time, we still have a generic funcion for computing the MSF.
"""

# ╔═╡ bad93338-72e3-46fc-888f-e507c62cefc6
begin
	local λs = 0.0:0.01:7.0
	plot(λs, msf.(λs),
		label="",
		xlabel="-Re λ"
	)
	#plot!(λs, msf_ds.(λs), label="ds")
	plot!([0.0,7.0],[0.0,0.0], linestyle=:dash, label="", color=:black)
	scatter!(λ, msf.(λ), label="")
end

# ╔═╡ 63c14c9e-808e-464b-b5cc-25921aa34f56
md"""
Below, we examine the trajectory of the MSF equations for λ = 0. Note that this is quite ill-behaved (notice the log scale), so it's no wonder this is tricky to get right.
"""

# ╔═╡ 8f8a9dec-4de2-4fc7-9d87-a3638ec47ee3
begin
	f! = msf_ode_fun(brusselator)
	u0_ = hcat(u0_msf, 0.1*randn(2))
	p_msf = [p_f, D_msf, 0.0]
	prob_msf = ODEProblem(f!, u0_, (0.0,100.0), p_msf)
	sol_msf = solve(prob_msf, Rodas5(), abstol=1e-10, reltol=1e-10)
end

# ╔═╡ e85b9b6d-736d-446a-acb1-b6ed0ccdd907
begin
	local foo = (t,u,v) -> (t,norm([u;v]))
	plot(sol_msf, vars=(foo,0,3,4),
		yscale=:log10,
		label="",
		ylabel=L"\ln \|\xi\|"
	)
end

# ╔═╡ Cell order:
# ╟─eaf9dda2-d898-4771-9f86-ca4be174712c
# ╠═cf07621a-9a65-11ec-1fc0-51fbadbb086f
# ╠═345a90f7-f8ac-43b2-b5bd-0282d9447556
# ╟─eb6aed5a-0f95-4a76-8a99-b8719984bbe5
# ╟─0e28ede5-0af6-49a7-b95d-7094e320eee4
# ╠═6ccfc87f-40e7-4ee4-baae-ee41b81efe67
# ╟─3e4ec586-482f-4c1b-a8e7-4ae88e48fa38
# ╠═ce0b65b5-e57a-4def-acb6-892c49f182d9
# ╟─1a7e8273-8069-42a8-ab7b-62f5fe1ade67
# ╠═1d77a2be-773c-4edd-9d68-1ed84fb5dadf
# ╠═7c712f73-3931-4cf4-82d8-9d9e2263df14
# ╟─432fa70e-9225-4e76-aec3-c27d31887fb9
# ╟─60063cc3-e71f-461b-be0a-491b1b566718
# ╟─ea2eeac2-5514-4fe8-8921-6235cdb0d1ad
# ╟─2bf13604-71c6-474f-90de-5b270c5c0708
# ╠═743b31cd-6b2a-472c-8a70-a3955a2c26b6
# ╟─03bd6ff6-107f-4b04-9cb8-fe5c28570e58
# ╠═08fc4616-d2af-471b-adc0-195bb566334f
# ╠═d0c1e91e-d1df-4a12-a0e2-8513eeee7ba6
# ╠═8623c9b7-5724-4678-9b53-4312389a276d
# ╠═da0220ea-99aa-4a0d-8826-f0e461fd503e
# ╟─60d1ea9b-b935-4e88-9e92-07497f06aeaa
# ╟─efbb46be-c442-4ab3-85a0-b10cbc3a2a9e
# ╠═8b2ec98a-b37f-4bd0-9fc8-9f3dcf9550c9
# ╠═25c22a1f-6ae4-4b72-8e4c-de844e2699d1
# ╟─582b74c7-aaf6-4a8e-901f-489199ccc0d6
# ╠═3e575a45-50d6-4297-84e5-6699e36d00dd
# ╟─ed1de34a-a42e-474d-ba36-411163a343f8
# ╠═947c169c-a900-4728-ad29-5200b52f9599
# ╠═04bbc594-17c9-4b41-a599-d41e2e541c9f
# ╟─29a0f8c8-5529-496c-99f1-b436d15f2517
# ╟─f346a0b3-4d1d-49ad-94ec-cc42ee75a05b
# ╠═8a1c55e5-1008-44ee-b853-c53c2dbf3bb2
# ╠═2ee789c5-ae38-4e85-95e1-920a580df058
# ╠═9d1f3408-e0f5-4093-bb6b-862485e99476
# ╟─a4d4a26f-ff25-4bf0-b8ab-2db56131ee3e
# ╟─91f1ce37-0d68-4ed2-9359-dbcfd06f6cdd
# ╠═4328e814-82a8-4e7d-8b8a-44a9d5ba17c5
# ╟─7b11ac9b-a470-4059-bbad-ba6bc443c68b
# ╠═a5629c28-53bf-4b19-a49f-ab26df1d1b46
# ╠═74adb98a-1ad5-4181-a3be-ec5fe2565a35
# ╟─c38e202f-a2d2-4449-8308-df22774b7af1
# ╟─d5ab9847-ace4-4a5b-8c81-d2eb8b56c896
# ╟─575e3106-704e-4bea-8155-7ee6e6c9da5d
# ╠═bff34b78-f7ad-452e-9b95-a498f3e377f5
# ╟─bba5eaa1-16f2-401d-8cf3-f05e9761583d
# ╠═ddbb84d1-945a-4184-850b-31ee3540e960
# ╟─f0215a78-146d-4009-ae73-8a2f35072db9
# ╠═84d86795-e2c4-4c2c-888b-ae999bcda2dd
# ╠═8c592351-4da7-4bed-a99c-595156c1ed8b
# ╠═e3411cdb-acf7-40da-9f5e-6852c1370797
# ╠═5fbfe4aa-f2e5-4845-92e7-9de00063059f
# ╟─ffecadd1-95b5-412f-86e4-e4dee733cab1
# ╠═bad93338-72e3-46fc-888f-e507c62cefc6
# ╟─63c14c9e-808e-464b-b5cc-25921aa34f56
# ╠═8f8a9dec-4de2-4fc7-9d87-a3638ec47ee3
# ╠═e85b9b6d-736d-446a-acb1-b6ed0ccdd907

max_stars_repo_name: csimal/Julia-for-Dynamical-Systems

max_issues_repo_name: csimal/Julia-for-Dynamical-Systems

max_forks_repo_name: csimal/Julia-for-Dynamical-Systems


text:
using TropicalYao.Reversible
using Test
using LinearAlgebra
using Yao, NiLang

@testset "new instr" begin
    g4 = Diagonal(i_gvb_tensor(1.5))
    reg = ArrayReg(randn(1<<12) .|> Tropical)
    S = stack4reg(reg, 0)
    s1 = i_instruct!(copy(vec(reg.state)), g4, (3, 7), (), (), S)[1]
    nreg = copy(reg) |> put(12, (3, 7)=>matblock(g4))
    @test statevec(nreg) ≈ s1

    g4 = randn(4, 4) .|> Tropical
    reg = ArrayReg(randn(1<<12) .|> Tropical)
    S = stack4reg(reg, 1)
    s1 = i_instruct!(copy(vec(reg.state)), g4, (3, 7), (), (), S)[1]
    nreg = copy(reg) |> put(12, (3, 7)=>matblock(g4))
    @test statevec(nreg) ≈ s1
end

@testset "stack" begin
    a = one(Tropical{Float64})
    b = Tropical(randn())
    @instr TropicalYao.Reversible.unsafe_store!(a, b)
    @test a == b
    @instr (~TropicalYao.Reversible.unsafe_store!)(a, b)
    @test a == one(Tropical{Float64})
    @test check_inv(TropicalYao.Reversible.unsafe_store!, (a, b))
end


text:
mutable struct Tile
    id::Int64
    data::Array{Int64,2}
end
Base.show(io::IO, z::Tile) = print(io, z.id)

function rotate(tile::Tile)
    new_tile = Tile(tile.id, rotl90(tile.data))
    return new_tile
end

function rotate!(tile::Tile)
    tile.data = rotl90(tile.data)
    return tile
end

function flip_y(tile::Tile)
    new_tile = Tile(tile.id, reverse(tile.data, dims=2))
    return new_tile
end

function flip_y!(tile::Tile)
    tile.data = reverse(tile.data, dims=2)
    return tile
end

function flip_x(tile::Tile)
    new_tile = Tile(tile.id, reverse(tile.data, dims=1))
    return new_tile
end

function flip_x!(tile::Tile)
    tile.data = reverse(tile.data, dims=1)
    return tile
end

function flip_xy(tile::Tile)
    new_tile = Tile(tile.id, reverse(reverse(tile.data, dims=2), dims=1))
    return new_tile
end

function flip_xy!(tile::Tile)
    tile.data = reverse(reverse(tile.data, dims=2), dims=1)
    return tile
end

function sides(tile::Tile)::Array{Array{Int64,1}}
    arr = tile.data
    return [arr[:, 1], arr[end, :], arr[:, end], arr[1, :]] # t, r, b, l
end

function side(tile::Tile, n::Int64)::Array{Int64,1}
    return sides(tile)[n]
end

function parse_side(side::Union{Array{Int64,1},Nothing})::Int64
    if side === nothing
        return -1
    end
    return parse(Int, join(side), base = 2)
end

function find_sides(tiles::Array{Tile,1})
    all_sides = []
    for t in tiles
        append!(all_sides, sides(t))
        append!(all_sides, reverse.(sides(t)))
    end
    return all_sides
end

function find_corners(tiles::Array{Tile,1})::Array{Tile,1}
    all_sides = find_sides(tiles)

    counter = Dict()
    for s in all_sides
        if haskey(counter, s)
            counter[s] += 1
        else
            counter[s] = 1
        end
    end

    corners = filter(tile -> count(s -> counter[s] == 1, sides(tile))
    == 2, tiles)

    return corners
end

function draw(tile::Tile)
    for i in 1:size(tile.data, 2)
        println(join(map(x -> x == 1 ? "#" : ".", tile.data[i, :])))
    end
end

function draw(image_map::Array{Tile,2})
    for i in 1:size(image_map, 2)
        for r in 1:size(image_map[1,1].data, 2)
            for t in image_map[:, i]
                print(join(map(x -> x == 1 ? "#" : ".", t.data[r, :])), " ")
            end
            println()
        end
        println()
    end
end

tiles = open("image_tiles.txt") do f
    map(m -> Tile(parse(Int, m.captures[1]), hcat(map(r -> map(x -> x == '#' ? 1 : 0, collect(r)), split(strip(m.captures[2]), "\n"))...)), eachmatch(r"Tile (\d+):\n((?:[#\.]+\n)+)", read(f, String)))
end

println(prod(map(x -> x.id, find_corners(tiles))))

dim = Int(sqrt(length(tiles)))
image = Array{Tile}(undef, dim, dim)
# pick first corners
# start to assemble from top to bottom, left to right
# orient tiles to match

function find_tile(s::Array{Int64,1}, tiles::Array{Tile,1}, n::Int64)
    tile = nothing
    for t in tiles
        # println("looking at $(t)")
        if s in find_sides([t])
            # println("Tile $(t) matches side $(parse_side(s))")
            tile = t
            break
        end
    end
    if tile === nothing
        error("no matching tile for $(parse_side(s)) found")
    end

    # rotate tile to match side n
    while side(tile, n) != s
        # println("rotating")
        rotate!(tile)
        # println(parse_side.(sides(tile)))
        if side(tile, n) == s
            break
        end
        for op in [:flip_x, :flip_y, :flip_xy]
            # println("flipping $(op)")
            if side(eval(op)(tile), n) == s
                eval(Symbol(op, "!"))(tile)
                break
            end
        end
    end

    return tile
end

function rotate_corner!(tile::Tile, tiles::Array{Tile,1})::Tile
    bag = filter(t -> t != tile, tiles)
    all_sides = find_sides(bag)

    while !(side(tile, 2) in all_sides && side(tile, 3) in all_sides)
        rotate!(tile)
    end

    return tile
end

function find_monsters(tile::Tile)::Int64
    monsters = 0
    monster = [
        "^..................#.",
        "^#....##....##....###",
        "^.#..#..#..#..#..#..."]
    for r in 1:(size(tile.data, 1)-length(monster)+1)
        for c in 1:(size(tile.data, 2)-length(monster[1]))
            # println("($(r),$(c))")
            if occursin(Regex(monster[1]), join(map(x -> x == 1 ? "#" : ".",tile.data[r,c:end])))
                # println("head found")
                if occursin(Regex(monster[2]), join(map(x -> x == 1 ? "#" : ".",tile.data[r+1,c:end])))
                    # println("body found")
                    if occursin(Regex(monster[3]), join(map(x -> x == 1 ? "#" : ".",tile.data[r+2,c:end])))
                        # println("legs found")
                        monsters += 1
                    end
                end
            end
        end
    end

    return monsters
end

left = nothing
for x in 1:dim
    top = nothing
    for y in 1:dim
        global left
        # println("($(x),$(y)): top: $(parse_side(top)), left: $(parse_side(left))")
        if y == 1 && x == 1 # set corner tile
            tile = find_corners(tiles)[1]
            # println("($(x),$(y)): setting corner tile $(tile)")
            image[1,1] = tile
            rotate_corner!(tile, tiles)
            top = side(tile, 3) # bottom
            left = side(tile, 2) # right
        elseif y == 1 # new column, match left
            tile = find_tile(left, tiles, 4)
            # println("($(x),$(y)): found tile $(tile) matching left $(parse_side(left))")
            image[y,x] = tile
            top = side(tile, 3)
            left = side(tile, 2)
        else # match top
            tile = find_tile(top, tiles, 1)
            # println("($(x),$(y)): found tile $(tile) matching top $(parse_side(top)) == $(parse_side(side(tile,1)))")
            image[y,x] = tile
            # no rotaton necessary
            top = side(tile, 3)
        end
        filter!(x -> x.id != tile.id, tiles)
        # println("$(length(tiles)) tiles left")
    end
end
# println(image)
# draw(image)

# remove borders (strip sides)
for i in eachindex(image)
    image[i].data = image[i].data[2:end-1, 2:end-1]
end
# draw(image)
# remove gaps (make one huge tile?)
arr = []
for i in 1:size(image, 2)
    global arr
    local row
    for r in 1:size(image[1,1].data, 2)
        row = []
        for t in image[i, :]
            row = vcat(row, t.data[:, r])
        end
        push!(arr, row)
    end
end
image = Tile(0, hcat(arr...))

monsters = 0
for i in 1:4
    global monsters, image
    rotate!(image)
    monsters = [find_monsters(image), find_monsters(flip_x(image)), find_monsters(flip_y(image)), find_monsters(flip_xy(image))]
    if maximum(monsters) > 0
        break
    end
end

println(sum(image.data) - 15*maximum(monsters))


text:
"""
    ui_update(id, value, session)

Manually update a UI variable/element `id` to `value`. Use inside server-side functions for
_other_ outputs when you want calculations to have side effects in the UI.

`session` is the optional `session` argument to server-side functions that provides stateful
variables. See the guides on [server-side state](@ref g06a-server-side-state) and
[side effects](@ref g06b-side-effects) for more information.
"""
function update_output(id, value, session)
    response = try
        response = render_response(value)
        JSON.json(Dict("id" => id, "value" => response))
    catch e
        msg = sprint(showerror, e)
        @error string(msg, "\nError occurred trying to `update_output` for `$id` to `$value`")

        JSON.json(Dict("matte_error_msg" => msg))
    end
    Genie.WebChannels.message(session.genie_wsclient, response)
end


text:
function checkCorner(arr, corneridx)
    #1 = up left, 2 = up right, three = down left, four = bottom right
    islowpoint = false
    maxrow = size(arr, 1)
    maxcol = size(arr[1], 1)
    if corneridx == 1
        val = arr[1][1]
        if val < arr[1][2] && val < arr[2][1]
            islowpoint = true
        end
    elseif corneridx == 2
        val = arr[1][maxcol]
        if val < arr[1][maxcol-1] && val < arr[2][maxcol]
            islowpoint = true
        end
    elseif corneridx == 3
        val = arr[maxrow][1]
        if val < arr[maxrow-1][1] && val < arr[maxrow][2]
            islowpoint = true
        end
    elseif corneridx == 4
        val = arr[maxrow][maxcol]
        if val < arr[maxrow-1][maxcol] && val < arr[maxrow][maxcol-1]
            islowpoint = true
        end
    end
    return islowpoint
end

function checkOnBoundary(arr, i, j)
     #1 = top row, 2 = left col, three = right col, four = bottom row
     islowpoint = false
     maxrow = size(arr, 1)
     maxcol = size(arr[1], 1)
     if (i,j) == (1,1)
        islowpoint = checkCorner(arr, 1)
     elseif (i,j) == (1, maxcol)
        islowpoint = checkCorner(arr, 2)
     elseif (i,j) == (maxrow, 1)
        islowpoint = checkCorner(arr, 3)
     elseif (i,j) == (maxrow, maxcol)
        islowpoint = checkCorner(arr, 4)
     elseif i == 1
        val = arr[i][j]
        if val < arr[i][j-1] && val < arr[i][j+1] && val < arr[i+1][j]
            islowpoint = true
        end
    elseif j == 1
        val = arr[i][j]
        if val < arr[i-1][j] && val < arr[i+1][j] && val < arr[i][j+1]
            islowpoint = true
        end
    elseif j == maxcol
        val = arr[i][j]
        if val < arr[i-1][maxcol] && val < arr[i+1][maxcol] && val < arr[i][maxcol - 1]
            islowpoint = true
        end
    elseif i == maxrow
        val = arr[i][j]
        if val < arr[maxrow][j-1] && val < arr[maxrow][j+1] && val < arr[maxrow-1][j]
            islowpoint = true
        end
    end
    return islowpoint
end

function checkPoint(arr, i, j)
    islowpoint = false
    maxrow = size(arr, 1)
    maxcol = size(arr[1], 1)
    if i == 1 || i == maxrow || j == 1 || j == maxcol
        islowpoint = checkOnBoundary(arr, i, j)
    else
        val = arr[i][j]
        if val < arr[i][j-1] && val < arr[i][j+1] && val < arr[i-1][j] && val < arr[i+1][j]
            islowpoint = true
        end
    end
    return islowpoint
end

function solve(filename)
    values = readlines(filename*".txt")
    values = split.(values, "")
    intvalues = Vector{Int}[]
    for line in values
        convline = parse.(Int, line)
        intvalues = push!(intvalues, convline)
    end
    #println(intvalues)
    res = 0
    for (i, line) in enumerate(intvalues)
        for (j, el) in enumerate(line)
            if checkPoint(intvalues, i, j)
                res += 1 + el
            end
        end
    end
    return res
end

println("Sum of lowpoints: "*string(solve("smalltest")))
println("Sum of lowpoints: "*string(solve("input")))

function getNeighbors(arr, (i, j))
    maxrow = size(arr, 1)
    maxcol = size(arr[1], 1)
    neighboridxs = Tuple[]
    if i == 1 || i == maxrow || j == 1 || j == maxcol
        if (i,j) == (1,1)
            neighboridxs = push!(neighboridxs, (1,2))
            neighboridxs = push!(neighboridxs, (2,1))
        elseif (i,j) == (1, maxcol)
            neighboridxs = push!(neighboridxs, (1,maxcol-1))
            neighboridxs = push!(neighboridxs, (2,maxcol))
        elseif (i,j) == (maxrow, 1)
            neighboridxs = push!(neighboridxs, (maxrow-1,1))
            neighboridxs = push!(neighboridxs, (maxrow,2))
        elseif (i,j) == (maxrow, maxcol)
            neighboridxs = push!(neighboridxs, (maxrow-1,maxcol))
            neighboridxs = push!(neighboridxs, (maxrow,maxcol-1))
        elseif i == 1
            neighboridxs = push!(neighboridxs, (i, j-1))
            neighboridxs = push!(neighboridxs, (i, j+1))
            neighboridxs = push!(neighboridxs, (i+1, j))
        elseif j == 1
            neighboridxs = push!(neighboridxs, (i-1, j))
            neighboridxs = push!(neighboridxs, (i, j+1))
            neighboridxs = push!(neighboridxs, (i+1, j))
        elseif j == maxcol
            neighboridxs = push!(neighboridxs, (i-1, maxcol))
            neighboridxs = push!(neighboridxs, (i+1, maxcol))
            neighboridxs = push!(neighboridxs, (i, maxcol - 1))
        elseif i == maxrow
            neighboridxs = push!(neighboridxs, (maxrow, j-1))
            neighboridxs = push!(neighboridxs, (maxrow, j+1))
            neighboridxs = push!(neighboridxs, (maxrow-1, j))
        end
    else
        neighboridxs = push!(neighboridxs, (i, j-1))
        neighboridxs = push!(neighboridxs, (i, j+1))
        neighboridxs = push!(neighboridxs, (i+1, j))
        neighboridxs = push!(neighboridxs, (i-1, j))
    end
    return neighboridxs
end

function getBasin(arr, lowpointcoords)
    neighbors = getNeighbors(arr, lowpointcoords)
    basinidxs = []
    for idx in neighbors
        if arr[idx[1]][idx[2]] != 9
            basinidxs = push!(basinidxs, idx)
        end
    end
    for idx in basinidxs
        i = idx[1]
        j = idx[2]
        candidates = getNeighbors(arr, idx)
        for c in candidates
            ci = c[1]
            cj = c[2]
            if arr[i][j] < arr[ci][cj] && arr[ci][cj] != 9
                basinidxs = push!(basinidxs, c)
            end
        end
    end
    basinidxs = push!(basinidxs, lowpointcoords)
    return unique!(basinidxs)
end

function solve2(filename)
    values = readlines(filename*".txt")
    values = split.(values, "")
    intvalues = Vector{Int}[]
    for line in values
        convline = parse.(Int, line)
        intvalues = push!(intvalues, convline)
    end
    lowpointidx = Tuple[]
    for (i, line) in enumerate(intvalues)
        for (j, el) in enumerate(line)
            if checkPoint(intvalues, i, j)
                lowpointidx = push!(lowpointidx, (i, j))
            end
        end
    end
    basins = []
    basinlengths = []
    for point in lowpointidx
        basin = getBasin(intvalues, point)
        blength = length(basin)
        basins = push!(basins, basin)
        basinlengths = push!(basinlengths, blength)
    end
    sortedlengths = sort(basinlengths, rev=true)
    res = sortedlengths[1] * sortedlengths[2] * sortedlengths[3]
    return res
end

println("Product of basin lengths: "*string(solve2("smalltest")))
println("Product of basin lengths "*string(solve2("input")))