text:
using Test
using TestSetExtensions
using LinearAlgebra
using Random
using RandomMatrices
using Qaintessent

##==----------------------------------------------------------------------------------------------------------------------
import LinearAlgebra.BLAS: @blasfunc

using LinearAlgebra: BlasInt
for (s, elty) in (("dlarfg_", Float64),
                  ("zlarfg_", ComplexF64))
    @eval begin
        function RandomMatrices.larfg!(n::Int, α::Ptr{$elty}, x::Ptr{$elty}, incx::Int, τ::Ptr{$elty})
	    ccall((@blasfunc($s), LAPACK.liblapack), Nothing,
		  (Ptr{Int}, Ptr{$elty}, Ptr{$elty}, Ptr{Int}, Ptr{$elty}),
                  Ref(n), α, x, Ref(incx), τ)
        end
    end
end

@testset ExtendedTestSet "compilediagonal unitaries helper functions" begin
    @testset "greyencode" begin
        g = Qaintessent.greyencode.(collect(0:15))
        @test all(g .== [0, 1, 3, 2, 6, 7, 5, 4, 12, 13, 15, 14, 10, 11, 9, 8])
    end

    @testset "svalue" begin
        nref = abs(rand(Int, 1)[1])%10
        b = collect(Qaintessent.svalue(nref))
        n = 0
        for i in b
            n += 1 << (i-1)
        end
        @test n == nref
    end

    @testset "fillψ!" begin
        N = 2
        d = exp.(im .* rand(Float64, 2^N))
        l = 2^(N-1)
        ψ = zeros(Float64, l-1)
        Qaintessent.fillψ!(d, ψ, l)
        ψref = zeros(Float64, l-1)
        for i in StepRange(1,1,l-1)
            ψref[i] = imag(log(d[2i-1]*d[2i+2]/(d[2i]*d[2i+1])))
        end
        @test all(ψref .≈ ψ)
    end

    @testset "ηcol" begin
        N = 4
        ref = [[-1, 0, 0, 0], [1, -1, 0, 0], [0, 1, -1, 0], [0, 0, 1, -1], [0, 0, 0, 1]]
        for i in 0:N
            @test all(ref[i+1] .== Qaintessent.ηcol(N, i))
        end
    end

    @testset "flip_state" begin
        N = 8
        @test all([1, 0, -1, 0, 1, 0, -1] .== Qaintessent.flip_state(3, N))
    end
end

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "compile2qubit unitaries helper functions" begin
    @testset "decomposeSO4" begin
        N = 2
        E = 1/sqrt(2) .* [1 im 0 0; 0 0 im 1; 0 0 im -1; 1 -im 0 0]

        M = Stewart(Float64, 2^N)
        Q, _ = qr(M)
        Q = Matrix{ComplexF64}(Q)
        if det(Q) == -1
            Q = Q*Q
        end
        A, B = Qaintessent.decomposeSO4(Q*Q)
        @test E'*kron(A, B)*E ≈ Q*Q
    end
end

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "general compile unitaries helper functions" begin
    @testset "QR decomposition" begin
        N = 2
        U, _ = qr(Matrix(Stewart(ComplexF64, 2^N)))
        U = Matrix(U)
        Uref = deepcopy(U)

        @testset "unblocked QR decomposition" begin
            QR, τ = Qaintessent.qr_unblocked(deepcopy(U))
            R = diag(QR)
            Q = Matrix{ComplexF64}(I, (2^N, 2^N))
            Id = Matrix{ComplexF64}(I, (2^N, 2^N))

            for i in 1:size(U)[1]-1
                u = pushfirst!(QR[i+1:2^N, i], 1)
                u = u ./ norm(u)
                H = deepcopy(Id)
                H[i:2^N, i:2^N] = H[i:2^N, i:2^N] - 2*u*u'
                Q = H*Q
                U = H*U
            end
            d = diag(U)

            @test diagm(d) ≈ U
            @test diag(U) ≈ R
            @test Q*Uref ≈ U
            @test inv(Q)*U ≈ Uref
        end

        U, _ = qr(Matrix(Stewart(ComplexF64, 2^N)))
        U = Matrix(U)
        Uref = deepcopy(U)

        @testset "unblocked! QR decomposition" begin
            QR = deepcopy(U)
            τ = Qaintessent.qr_unblocked!(QR)
            R = diag(QR)
            Q = Matrix{ComplexF64}(I, (2^N, 2^N))
            Id = Matrix{ComplexF64}(I, (2^N, 2^N))

            for i in 1:size(U)[1]-1
                u = pushfirst!(QR[i+1:2^N, i], 1)
                u = u ./ norm(u)
                H = deepcopy(Id)
                H[i:2^N, i:2^N] = H[i:2^N, i:2^N] - 2*u*u'
                Q = H*Q
                U = H*U
            end
            d = diag(U)

            @test diagm(d) ≈ U
            @test diag(U) ≈ R
            @test Q*Uref ≈ U
            @test inv(Q)*U ≈ Uref
        end
    end

    @testset "stateprep" begin
        @inline function allequal(x)
            length(x) < 2 && return true
            e1 = x[1]
            i = 2
            @inbounds for i=2:length(x)
                x[i] == e1 || return false
            end
            return true
        end

        N = 3
        ψ = rand(ComplexF64, 2^N)
        ψ = ψ ./ norm(ψ)
        angles = exp.(im.*angle.(ψ))
        ψ = ψ ./ angles
        ψref = deepcopy(ψ)
        ϕ = zeros(ComplexF64, (2^N))
        ϕ[1] = 1

        for j in 1:N-1
            cg = Qaintessent.stateprep(ψ[1:2^(j-1):end], N, j)
            ψ = apply(ψ, cg)
        end
        θ = real(atan(-ψ[2^(N-1)+1]./ψ[1]).*2)

        if !isnan(θ)
            cg = CircuitGate((N,), RyGate(θ))
            ψ = apply(ψ, cg)
        end

        @test ψ ≈ ϕ
    end

    @testset "inverseM" begin
        N = 4
        mref = [ 0.125   0.125   0.125   0.125   0.125   0.125   0.125   0.125;
                0.125  -0.125   0.125  -0.125   0.125  -0.125   0.125  -0.125;
                0.125  -0.125  -0.125   0.125   0.125  -0.125  -0.125   0.125;
                0.125   0.125  -0.125  -0.125   0.125   0.125  -0.125  -0.125;
                0.125   0.125  -0.125  -0.125  -0.125  -0.125   0.125   0.125;
                0.125  -0.125  -0.125   0.125  -0.125   0.125   0.125  -0.125;
                0.125  -0.125   0.125  -0.125  -0.125   0.125  -0.125   0.125;
                0.125   0.125   0.125   0.125  -0.125  -0.125  -0.125  -0.125
                ]
        @test Qaintessent.inverseM(N) ≈ mref
    end
end

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "unitary compilation" begin
    @testset "general compile" begin
        N = 6
        U, _ = qr(Matrix(Stewart(ComplexF64, 2^N)))
        U = Matrix(U)
        M = Stewart(ComplexF64, 2^N)

        cgs = unitary2circuit(deepcopy(U), N)

        ψ = rand(ComplexF64, 2^N)
        ψ_ref = U*ψ
        ψ_compiled = apply(ψ, cgs)

        @test ψ_ref'*M*ψ_ref ≈ ψ_compiled'*M*ψ_compiled
    end

    @testset "general compile exceptions" begin
        N = 6
        U = rand(ComplexF64, (2^N, 2^N))
        @test_throws ErrorException("Only unitary matrices can be compiled into a quantum circuit") unitary2circuit(U)

        U = rand(ComplexF64, (2^N, 2^(N-1)))
        @test_throws ErrorException("Only unitary matrices can be compiled into a quantum circuit") unitary2circuit(U)
    end

    @testset "compile 1 qubit standard unitary" begin
        N = 1
        M = Stewart(ComplexF64, 2)
        random_θ = rand(Float64, 3)
        for gate in [X, Y, Z, HadamardGate(), TGate(), SGate(), RxGate(random_θ[1]), RyGate(random_θ[2]), RzGate(random_θ[3])]
            U = matrix(gate)

            cgs = unitary2circuit(U)

            ψ = rand(ComplexF64, 2^N)

            ψ_ref = U*ψ
            ψ_compiled = apply(ψ, cgs)
            @test ψ_ref'*M*ψ_ref ≈ ψ_compiled'*M*ψ_compiled
        end
    end

    @testset "compile 1 qubit random unitary" begin
        N = 1
        U, _ = qr(Stewart(ComplexF64, 2^N))
        U = Matrix(U)
        M = Stewart(ComplexF64, 2)

        cgs = unitary2circuit(U, N)

        ψ = rand(ComplexF64, 2^N)

        ψ_ref = U*ψ
        ψ_compiled = apply(ψ, cgs)

        @test ψ_ref'*M*ψ_ref ≈ ψ_compiled'*M*ψ_compiled
    end

    @testset "compile diagonal unitaries" begin
        N = 5
        U = diagm(exp.(im .* rand(Float64, 2^N)))
        M = Stewart(ComplexF64, 2^N)

        cgs = unitary2circuit(U, N)

        ψ = rand(ComplexF64, 2^N)

        ψ_ref = U*ψ
        ψ_compiled = apply(ψ, cgs)

        @test ψ_ref'*M*ψ_ref ≈ ψ_compiled'*M*ψ_compiled
    end

    @testset "compile 2 qubit unitaries" begin
        N = 2
        U, _ = qr(Matrix(Stewart(ComplexF64, 2^N)))
        U = Matrix(U)
        M = Stewart(ComplexF64, 2^N)

        cgs = unitary2circuit(U, N)
        ψ = rand(ComplexF64, 2^N)

        ψ_ref = U*ψ
        ψ_compiled = apply(ψ, cgs)

        @test ψ_ref'*M*ψ_ref ≈ ψ_compiled'*M*ψ_compiled
    end

    @testset "compile standard 2 qubit unitaries" begin
        N = 2
        M = Stewart(ComplexF64, 2^N)
        l = deepcopy(M)
        for _ in 1:400
            random_θ = rand(Float64, 3)
            for gate in [X, Y, Z, TGate(), SGate(), RxGate(random_θ[1]), RyGate(random_θ[2]), RzGate(random_θ[3])]

                nU = Matrix(sparse_matrix(circuit_gate(1, gate, 2)))
                cgs = unitary2circuit(nU, N)

                ψ = rand(ComplexF64, 2^N)
                ψ_ref = nU*ψ
                ψ_compiled = apply(ψ, cgs)

                @test isapprox(ψ_ref'*(M*ψ_ref), ψ_compiled'*(M*ψ_compiled), rtol=1e-5, atol=1e-5)
            end
        end
    end
end


text:
using Dash
using DashTextareaAutocomplete

const WORDLIST = ["apple", "application", "apartment",
                  "boat", "banana", "boss",
                  "coast", "code", "cat"]

app = dash(;external_stylesheets=["https://codepen.io/chriddyp/pen/bWLwgP.css"])

app.layout = html_div(className = "container") do
    dashtextareaautocomplete(id = "input",
                             # value = "initial value",
                             placeholder = "Type something!",
                             wordList = WORDLIST,
                             # common options with their default values
                             # triggerChar = ":",
                             # minChar = 1,
                             style = (minHeight = "100px",),
                             dropdownStyle = (maxHeight = "120px", overflow = "auto")),
    html_div(id = "output"),
    html_button("CLICK", id = "btn"),
    html_div(id = "output2")
end

callback!(app, Output("output", "children"),
               Input("input", "value")) do value
    return "You have entered $value"
end

callback!(app, Output("output2", "children"),
               Input("btn", "n_clicks"),
               State("input", "value")) do clicks, value
    return "After click, you have: $value"
end

run_server(app, "0.0.0.0", 8050, debug=true)

ext: jl


text:
ia> [RpcApi.SimpleOrganism(0x0000000000014444,:(function (o::Creature.Organism,)
     local var_1::Int8 = 61
     Creature.stepRight(o)
     Creature.eatDown(o,Int(var_1))
     Creature.eatUp(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,234471,Dict{Int16,Int16}(),[19,99],4),RpcApi.SimpleOrganism(0x0000000000014252,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.eatUp(o,Int(var_1))
     Creature.stepDown(o)
     Creature.eatDown(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,189356,Dict{Int16,Int16}(),[3,58],4),RpcApi.SimpleOrganism(0x0000000000012f53,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.eatRight(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,163929,Dict{Int16,Int16}(),[0,35],2),RpcApi.SimpleOrganism(0x000000000000a639,:(function (o::Creature.Organism,)
     local var_1::Int8 = 41
     Creature.stepLeft(o)
     return true
 end),[100,300,95,0,1,1,1],1,400,1,160499,Dict{Int16,Int16}(),[23,14],2),RpcApi.SimpleOrganism(0x00000000000146e2,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.stepUp(o)
     Creature.stepDown(o)
     return true
 end),[100,300,95,0,1,1,1],1,400,1,133996,Dict{Int16,Int16}(),[0,35],3),RpcApi.SimpleOrganism(0x0000000000014963,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     if var_1 >= var_1
     end
     return true
 end),[100,300,95,0,1,1,1],1,400,1,121367,Dict{Int16,Int16}(),[21,41],2),RpcApi.SimpleOrganism(0x000000000001262e,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     var_1 = var_1 - var_1
     Creature.eatUp(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,120371,Dict{Int16,Int16}(),[0,35],3),RpcApi.SimpleOrganism(0x000000000000c00d,:(function (o::Creature.Organism,)
     Creature.stepDown(o)
     return true
 end),[100,300,95,0,1,1,1],1,400,1,116664,Dict{Int16,Int16}(),[4,58],1),RpcApi.SimpleOrganism(0x000000000001429e,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.stepLeft(o)
     Creature.eatLeft(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,114691,Dict{Int16,Int16}(),[0,36],3),RpcApi.SimpleOrganism(0x000000000001439c,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.stepUp(o)
     Creature.eatLeft(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,111356,Dict{Int16,Int16}(),[0,36],3),RpcApi.SimpleOrganism(0x0000000000014266,:(function (o::Creature.Organism,)
     local var_1::Int8 = 61
     Creature.eatDown(o,Int(var_1))
     Creature.stepDown(o)
     var_1 = var_1 / var_1
     return true
 end),[100,300,95,0,1,1,1],1,400,1,109727,Dict{Int16,Int16}(),[20,98],4),RpcApi.SimpleOrganism(0x0000000000014415,:(function (o::Creature.Organism,)
     function func_4(var_2::Int64=6870157196147773660,var_3::Int16=-29673)
         return var_2
     end
     local var_1::Int8 = 99
     Creature.stepUp(o)
     Creature.eatLeft(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,108763,Dict{Int16,Int16}(),[0,36],4),RpcApi.SimpleOrganism(0x000000000001521a,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.eatLeft(o,Int(var_1))
     Creature.stepLeft(o)
     Creature.eatUp(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,107001,Dict{Int16,Int16}(),[0,39],4),RpcApi.SimpleOrganism(0x0000000000014491,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.eatLeft(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,100090,Dict{Int16,Int16}(),[0,36],2),RpcApi.SimpleOrganism(0x00000000000144da,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     function func_4(var_2::Int16=5846,var_3::Int8=59)
         Creature.stepDown(o)
         return var_2
     end
     for i = 1:var_1
     end
     Creature.eatLeft(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,99262,Dict{Int16,Int16}(),[0,36],5),RpcApi.SimpleOrganism(0x00000000000129e7,:(function (o::Creature.Organism,)
     function func_3(var_2::String="LK2ZxJSv")
         return var_2
     end
     local var_1::Int8 = 99
     Creature.stepDown(o)
     Creature.eatLeft(o,Int(var_1))
     Creature.eatUp(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,97758,Dict{Int16,Int16}(),[0,36],5),RpcApi.SimpleOrganism(0x00000000000145bd,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.eatDown(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,96823,Dict{Int16,Int16}(),[0,36],2),RpcApi.SimpleOrganism(0x000000000001452d,:(function (o::Creature.Organism,)
     Creature.eatLeft(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,96760,Dict{Int16,Int16}(),[0,36],1),RpcApi.SimpleOrganism(0x0000000000014521,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.eatLeft(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,95022,Dict{Int16,Int16}(),[0,36],2),RpcApi.SimpleOrganism(0x000000000001432f,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.stepRight(o)
     return true
 end),[100,300,95,0,1,1,1],1,400,1,94093,Dict{Int16,Int16}(),[0,37],2),RpcApi.SimpleOrganism(0x000000000001463e,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.stepLeft(o)
     Creature.eatUp(o,Int(var_1))
     Creature.stepUp(o)
     return true
 end),[100,300,95,0,1,1,1],1,400,1,89478,Dict{Int16,Int16}(),[1,29],4),RpcApi.SimpleOrganism(0x000000000001504c,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.stepRight(o)
     Creature.eatRight(o,Int(var_1))
     Creature.stepDown(o)
     return true
 end),[100,300,95,0,1,1,1],1,400,1,88248,Dict{Int16,Int16}(),[67,80],4),RpcApi.SimpleOrganism(0x00000000000144cb,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     local var_2::Bool = false
     Creature.eatLeft(o,Int(var_1))
     Creature.stepLeft(o)
     Creature.stepRight(o)
     return true
 end),[100,300,95,0,1,1,1],1,400,1,87933,Dict{Int16,Int16}(),[0,36],5),RpcApi.SimpleOrganism(0x000000000001468b,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.stepUp(o)
     Creature.stepUp(o)
     Creature.eatLeft(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,86576,Dict{Int16,Int16}(),[0,36],4),RpcApi.SimpleOrganism(0x000000000001434d,:(function (o::Creature.Organism,)
     local var_1::Int8 = 41
     Creature.eatRight(o,Int(var_1))
     Creature.eatDown(o,Int(var_1))
     var_1 = var_1 + var_1
     Creature.eatDown(o,Int(var_1))
     Creature.stepDown(o)
     return true
 end),[100,300,95,0,1,1,1],1,400,1,83536,Dict{Int16,Int16}(),[55,94],6),RpcApi.SimpleOrganism(0x0000000000014411,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     var_1 = var_1 / var_1
     var_1 = var_1 + var_1
     return true
 end),[100,300,95,0,1,1,1],1,400,1,82340,Dict{Int16,Int16}(),[0,37],3),RpcApi.SimpleOrganism(0x0000000000014bdb,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.stepUp(o)
     Creature.stepLeft(o)
     Creature.eatUp(o,Int(var_1))
     Creature.eatRight(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,81699,Dict{Int16,Int16}(),[3,28],5),RpcApi.SimpleOrganism(0x00000000000147f3,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.stepUp(o)
     Creature.eatLeft(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,78967,Dict{Int16,Int16}(),[0,36],3),RpcApi.SimpleOrganism(0x00000000000144bf,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.stepLeft(o)
     return true
 end),[100,300,95,0,1,1,1],1,400,1,77628,Dict{Int16,Int16}(),[0,36],2),RpcApi.SimpleOrganism(0x00000000000143c8,:(function (o::Creature.Organism,)
     local var_1::Int8 = 61
     Creature.stepLeft(o)
     var_1 = var_1 / var_1
     return true
 end),[100,300,95,0,1,1,1],1,400,1,75354,Dict{Int16,Int16}(),[1,95],3),RpcApi.SimpleOrganism(0x0000000000014fea,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.eatUp(o,Int(var_1))
     Creature.eatRight(o,Int(var_1))
     Creature.stepLeft(o)
     Creature.stepRight(o)
     return true
 end),[100,300,95,0,1,1,1],1,400,1,74932,Dict{Int16,Int16}(),[13,28],5),RpcApi.SimpleOrganism(0x0000000000014216,:(function (o::Creature.Organism,)
     local var_1::Int8 = 41
     local var_2::Bool = true
     Creature.eatLeft(o,Int(var_1))
     Creature.eatUp(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,74553,Dict{Int16,Int16}(),[21,13],4),RpcApi.SimpleOrganism(0x000000000001494b,:(function (o::Creature.Organism,)
     Creature.eatLeft(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,74524,Dict{Int16,Int16}(),[0,36],1),RpcApi.SimpleOrganism(0x0000000000014d26,:(function (o::Creature.Organism,)
     local var_1::Int8 = 74
     Creature.stepUp(o)
     Creature.eatUp(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,71398,Dict{Int16,Int16}(),[33,9],3),RpcApi.SimpleOrganism(0x0000000000014934,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.eatLeft(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,71008,Dict{Int16,Int16}(),[0,36],2),RpcApi.SimpleOrganism(0x0000000000014d09,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.stepRight(o)
     Creature.eatRight(o,Int(var_1))
     Creature.stepRight(o)
     if var_1 === var_1
     end
     return true
 end),[100,300,95,0,1,1,1],1,400,1,69482,Dict{Int16,Int16}(),[51,67],5),RpcApi.SimpleOrganism(0x00000000000148ec,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.stepUp(o)
     Creature.eatLeft(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,69300,Dict{Int16,Int16}(),[0,36],3),RpcApi.SimpleOrganism(0x00000000000149f2,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.stepUp(o)
     return true
 end),[100,300,95,0,1,1,1],1,400,1,67718,Dict{Int16,Int16}(),[0,36],2),RpcApi.SimpleOrganism(0x0000000000014956,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.stepLeft(o)
     Creature.eatLeft(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,67197,Dict{Int16,Int16}(),[0,36],3),RpcApi.SimpleOrganism(0x000000000000eca6,:(function (o::Creature.Organism,)
     local var_1::Int8 = 17
     Creature.stepDown(o)
     Creature.eatDown(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,67084,Dict{Int16,Int16}(),[77,59],3),RpcApi.SimpleOrganism(0x0000000000014606,:(function (o::Creature.Organism,)
     Creature.stepRight(o)
     Creature.eatRight(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,65455,Dict{Int16,Int16}(),[0,37],2),RpcApi.SimpleOrganism(0x0000000000008cc4,:(function (o::Creature.Organism,)
     local var_1::Int8 = 41
     return true
 end),[100,300,95,0,1,1,1],1,400,1,64154,Dict{Int16,Int16}(),[22,14],1),RpcApi.SimpleOrganism(0x000000000001423c,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.stepUp(o)
     return true
 end),[100,300,95,0,1,1,1],1,1374,1,63775,Dict{Int16,Int16}(),[0,36],2),RpcApi.SimpleOrganism(0x0000000000014978,:(function (o::Creature.Organism,)
     local var_1::Int8 = 41
     Creature.stepLeft(o)
     Creature.eatLeft(o,Int(var_1))
     Creature.eatLeft(o,Int(var_1))
     Creature.stepDown(o)
     return true
 end),[100,300,95,0,1,1,1],1,400,1,62588,Dict{Int16,Int16}(),[23,100],5),RpcApi.SimpleOrganism(0x0000000000014b13,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.eatDown(o,Int(var_1))
     Creature.eatLeft(o,Int(var_1))
     Creature.eatUp(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,61452,Dict{Int16,Int16}(),[0,36],4),RpcApi.SimpleOrganism(0x0000000000014631,:(function (o::Creature.Organism,)
     local var_2::Int16 = -11035
     local var_1::Int8 = 99
     Creature.stepUp(o)
     Creature.stepRight(o)
     Creature.eatDown(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,61337,Dict{Int16,Int16}(),[0,37],5),RpcApi.SimpleOrganism(0x0000000000014b3e,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.stepRight(o)
     Creature.eatRight(o,Int(var_1))
     var_1 = var_1 / var_1
     Creature.stepLeft(o)
     return true
 end),[100,300,95,0,1,1,1],1,400,1,61226,Dict{Int16,Int16}(),[15,68],5),RpcApi.SimpleOrganism(0x000000000001471a,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     local var_2::Int64 = -7709863351867858983
     Creature.eatUp(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,61188,Dict{Int16,Int16}(),[0,37],3),RpcApi.SimpleOrganism(0x000000000001446c,:(function (o::Creature.Organism,)
     local var_1::Int8 = 99
     Creature.stepUp(o)
     for i = 1:var_1
     end
     return true
 end),[100,300,95,0,1,1,1],1,400,1,58846,Dict{Int16,Int16}(),[0,37],3),RpcApi.SimpleOrganism(0x0000000000014d5e,:(function (o::Creature.Organism,)
     local var_1::Int8 = 41
     Creature.eatLeft(o,Int(var_1))
     Creature.stepDown(o)
     Creature.eatLeft(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],1,400,1,58123,Dict{Int16,Int16}(),[23,99],4)]

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
"Algorithms for discovering sparse bases for the symmetry-allowed space of exchange couplings."

# ##### Math behind `basis_for_symmetry_allowed_couplings()` #####
#
# A crystal spacegroup consists of a set of symmetry operations. Each symop is
# defined by an orthogonal (rotation or reflection) matrix R and a translation
# vector T. The symop transforms every atom position `x` to a new one `x′ = Rx +
# T`, but leaves the crystal as a whole invariant. Consider some bond b=(i,j),
# defined as an ordered pair of atoms, not necessarily in the same crystal unit
# cell. Suppose this bond carries an exchange interaction of the form `S_iᵀ J
# S_j`. After applying a symop defined by (R, T), the bond b=(i,j) maps to a
# transformed bond b′=(i′,j′). This new bond must carry an exchange interaction
# of the form `S_{i′}ᵀ J′ S_{j′}`, where J′ is related to J, as we will now
# show. Besides mapping the atom j to a new position j′, the symop also
# transforms its spin vector from `S_j` to `S′_{j′} = R S_j`. Similarly, `S_iᵀ`
# maps to `S′_{i′}ᵀ = S_iᵀ Rᵀ`. The energy along a bond, however, is invariant
# under the symmetry transformation, so we require `S_iᵀ J S_j = S′_{i′}ᵀ J′
# S′_{j′}`, or
#
#   `S_iᵀ J S_j = S_iᵀ Rᵀ J′ R S_j`.
#
# This equation must be valid for arbitrary `S_i` and `S_j`, which implies `J =
# Rᵀ J′ R`, or equivalently,
#
#   J′ = R J Rᵀ
#
# because the matrix R is orthogonal. The conclusion is that specifying the
# exchange matrix J for one bond implicitly constrains the exchange matrices J′
# for all symmetry equivalent bonds.
#
# Often a symmetry operation will map a bond b into itself, or into the same
# bond but with atom positions reversed. The existence of such symops constrains
# the space of allowed exchange matrices J for the bond b. Specifically, we
# require
#
#   (1)  J = R J Rᵀ   or   (2)  Jᵀ = R J Rᵀ
#
# for every symop `s = (R, T)` that maps `b` into `b` (constraint 1), or into
# `reverse(b)` (constraint 2). The intersection of all such constraints defines
# the symmetry-allowed space of exchange matrices J. The allowed matrices can be
# expressed as a linear combination of basis matrices because all constraints
# are linear in J.
#
# It is convenient to express the constraints (1) or (2) in the form `F J = 0`,
# where F denotes the linear operator such that `F J = R J Rᵀ - J` or `F J = R J
# Rᵀ - Jᵀ`. Here, we are viewing the 3×3 matrix J as a flattened 9-dimensional
# vector, and F as a 9x9 matrix. In the first case, `F = R⊗R - I`. In the second
# case, we should replace `I` by the suitable transpose operation,
# `transpose_op_3x3`, defined explicitly below.
#
# Any linear combination of "vectors" (3x3 matrices) in the null space of F,
# i.e. `F J = 0`, satisfies the symmetry constraint (1) or (2) for the given
# symop. The singular value decomposition
#
#     F = U Σ Vᵀ
#
# can be used to produce an orthogonal basis for this null space. It is spanned
# by columns v of V corresponding to the zero singular values. The space spanned
# by v equivalently represented as a projection matrix,
#
#     P = v vᵀ
#
# When there are multiple constraints (F1, F2, ... Fn) we should take the
# intersection of the spanned spaces of P1, ... Pn. To calculate this
# intersection, form the product of the projectors:
#
#     P = P1 P2 ... Pn
#
# The allowable J values correspond to the eigenvectors of P with eigenvalue 1.
# An orthogonal basis for this space can again be calculated with an SVD.

# A 9x9 matrix that, when applied to a flattened 3x3 matrix (viewed as a
# 9-dimensional vector), generates the transposed 3x3 matrix in flattened form.
const transpose_op_3x3 = SMatrix{9, 9, Float64}([
    1 0 0  0 0 0  0 0 0
    0 0 0  1 0 0  0 0 0
    0 0 0  0 0 0  1 0 0

    0 1 0  0 0 0  0 0 0
    0 0 0  0 1 0  0 0 0
    0 0 0  0 0 0  0 1 0

    0 0 1  0 0 0  0 0 0
    0 0 0  0 0 1  0 0 0
    0 0 0  0 0 0  0 0 1
])

# Returns a projection operator P that maps to zero any symmetry-unallowed
# coupling matrix J. The space spanned by the eigenvectors of P with eigenvalue
# 1 represents the allowed coupling matrices J.
function projector_for_symop(cryst::Crystal, s::SymOp, parity::Bool)
    # Cartesian-space rotation operator corresponding to `s`
    R = cryst.lat_vecs * s.R * inv(cryst.lat_vecs)

    # Constraint is modeled as `F J = 0`
    F = kron(R, R) - (parity ? SMatrix{9, 9, Float64}(I) : transpose_op_3x3)

    # Orthogonal column vectors that span the null space of F
    v = nullspace(F; atol=1e-12)

    # Projector onto the null space of F
    P = SMatrix{9, 9, Float64}(v * v')
    return P
end

# Return an operator P that implicitly gives the space of symmetry allowed
# coupling matrices for bond b. Specifically, x is an allowed coupling if and
# only if it is an eigenvector of P with eigenvalue 1, i.e., `P x = x`.
function symmetry_allowed_couplings_operator(cryst::Crystal, b::BondRaw)
    P = SMatrix{9, 9, Float64}(I)
    for (s, parity) in symmetries_between_bonds(cryst, b, b)
        P = P * projector_for_symop(cryst, s, parity)
    end
    # Allowed coupling matrices J are simultaneously eigenvectors for all
    # projectors above, with eigenvalue 1.
    return P
end

# Check whether a coupling matrix J is consistent with symmetries of a bond
function is_coupling_valid(cryst::Crystal, b::BondRaw, J::Mat3)
    for (s, parity) in symmetries_between_bonds(cryst, b, b)
        R = cryst.lat_vecs * s.R * inv(cryst.lat_vecs)
        # TODO use symprec to handle case where matrix R is not precise
        if norm(R*J*R' - (parity ? J : J')) > 1e-12
            return false
        end
    end
    return true
end

function is_coupling_valid(cryst::Crystal, b::Bond, J::Mat3)
    return is_coupling_valid(cryst, BondRaw(cryst, b), J)
end

# Orthonormal basis of 3x3 symmetric matrices
const sym_basis = begin
    b = [diagm([1, 0, 0]),
         diagm([0, 1, 0]),
         diagm([0, 0, 1]),
         [0 1 0
          1 0 0
          0 0 0]/√2,
         [0 0 1
          0 0 0
          1 0 0]/√2,
         [0 0 0
          0 0 1
          0 1 0]/√2,]
    SMatrix{9, 6, Float64}(hcat(reshape.(b, 9)...))
end

# Orthonormal basis of 3x3 antisymmetric matrices
const asym_basis = begin
    b = [[ 0  1  0
          -1  0  0
           0  0  0]/√2,
         [ 0  0 -1
           0  0  0
           1  0  0]/√2,
         [ 0  0  0
           0  0  1
           0 -1  0]/√2]
    SMatrix{9, 3, Float64}(hcat(reshape.(b, 9)...))
end

@assert sym_basis * sym_basis' + asym_basis * asym_basis' ≈ I

# Given an m×n matrix A with empty nullspace, linearly combine the n columns to
# make them sparser. Solution was proposed here:
# https://math.stackexchange.com/q/4227648/660903 . Closely related to finding
# reduced row echolon form.
function sparsify_columns(A; atol)
    if size(A, 2) <= 1
        return A
    else
        # By assumption, the n columns of A are linearly independent
        @assert isempty(nullspace(A; atol))
        # Since row rank equals column rank, it should be possible to find n
        # linearly independent rows in A. Store these independent rows of A as
        # column vectors in indep_rows.
        indep_rows = zeros(Float64, size(A, 2), 0)
        for row in eachrow(A)
            # Add `row` to list if linearly independent with existing ones
            if isempty(nullspace(hcat(indep_rows, row); atol))
                indep_rows = hcat(indep_rows, row)
            end
        end
        return A * inv(indep_rows)'
    end
end

const _basis_elements_by_priority = [1, 5, 9, 8, 3, 4]

function _score_basis_matrix(J)
    return findfirst(i -> abs(J[i]) > 1e-12, _basis_elements_by_priority)
end

"""    basis_for_symmetry_allowed_couplings(cryst::Crystal, b::Bond)

Returns a list of ``3×3`` matrices that form a linear basis for the
symmetry-allowed coupling matrices associated with bond `b`.
"""
function basis_for_symmetry_allowed_couplings(cryst::Crystal, b::BondRaw)
    # Expected floating point precision for 9x9 matrix operations
    atol = 1e-12

    P = symmetry_allowed_couplings_operator(cryst, b)
    # Any solution to the original symmetry constraints `R J Rᵀ = J` or `R J Rᵀ
    # = Jᵀ` decomposes into purely symmetric/antisymmetric solutions. Therefore
    # we can pick a basis that separates into symmetric and antisymmetric parts.
    # We will do so by decomposing P. By construction, P = P_sym+P_asym.
    P_sym  = P *  sym_basis *  sym_basis'
    P_asym = P * asym_basis * asym_basis'

    acc_sym = SVector{9, Float64}[]
    acc_asym = SVector{9, Float64}[]

    # If any "reference" basis vectors are eigenvalues of P_sym with eigenvalue
    # 1, use them as outputs, and remove them from P_sym
    for x in eachcol(sym_basis)
        if isapprox(P_sym*x, x; atol)
            push!(acc_sym, x)
            P_sym = P_sym * (I - x*x')
        end
    end
    # Same for P_asym
    for x in eachcol(asym_basis)
        if isapprox(P_asym*x, x; atol)
            push!(acc_asym, x)
            P_asym = P_asym * (I - x*x')
        end
    end

    # Search for eigenvectors of P_sym with eigenvalue 1. These provide an
    # orthonormal basis for symmetric couplings.
    v = nullspace(P_sym-I; atol)
    v = sparsify_columns(v; atol)
    append!(acc_sym, eachcol(v))
    # Same for P_asym
    v = nullspace(P_asym-I; atol)
    v = sparsify_columns(v; atol)
    append!(acc_asym, eachcol(v))

    # Sort basis elements according to the indices where the nonzero elements
    # first appear
    sort!(acc_sym;  by=_score_basis_matrix)
    sort!(acc_asym; by=_score_basis_matrix)

    acc = [acc_sym; acc_asym]
    return map(acc) do x
        # Normalize each basis vector so that its maximum component is 1. The
        # shift by atol avoids unnecessary sign change in the case where the
        # maximum magnitude values of x appear symmetrically as ±c for some c.
        _, i = findmax(abs.(x.+atol))
        x = x / x[i]

        # Reinterpret as 3x3 matrix
        x = Mat3(reshape(x, 3, 3))

        # Double check that x indeed satifies the necessary symmetries
        @assert is_coupling_valid(cryst, b, x)

        return x
    end
end

function basis_for_symmetry_allowed_couplings(cryst::Crystal, b::Bond)
    return basis_for_symmetry_allowed_couplings(cryst, BondRaw(cryst, b))
end

"""
    all_symmetry_related_couplings_for_atom(cryst::Crystal, i::Int, b::Bond, J)

Given a reference bond `b` and coupling matrix `J` on that bond, return a list
of symmetry-equivalent bonds (constrained to start from atom `i`), and a
corresponding list of symmetry-transformed coupling matrices.
"""
function all_symmetry_related_couplings_for_atom(cryst::Crystal, i::Int, b_ref::Bond, J_ref)
    J_ref = Mat3(J_ref)
    @assert is_coupling_valid(cryst, b_ref, J_ref)

    bs = Bond[]
    Js = Mat3[]

    for b in all_symmetry_related_bonds_for_atom(cryst, i, b_ref)
        push!(bs, b)
        (s, parity) = first(symmetries_between_bonds(cryst, BondRaw(cryst, b), BondRaw(cryst, b_ref)))
        R = cryst.lat_vecs * s.R * inv(cryst.lat_vecs)
        push!(Js, R * (parity ? J_ref : J_ref') * R')
    end

    return (bs, Js)
end

"""
    all_symmetry_related_couplings(cryst::Crystal, b::Bond, J)

Given a reference bond `b` and coupling matrix `J` on that bond, return a list
of symmetry-equivalent bonds and a corresponding list of symmetry-transformed
coupling matrices.
"""
function all_symmetry_related_couplings(cryst::Crystal, b_ref::Bond, J_ref)
    J_ref = Mat3(J_ref)

    bs = Bond[]
    Js = Mat3[]

    for i in eachindex(cryst.positions)
        (bs_i, Js_i) = all_symmetry_related_couplings_for_atom(cryst, i, b_ref, J_ref)
        append!(bs, bs_i)
        append!(Js, Js_i)
    end

    return (bs, Js)
end


text:
module Normalize

using LinearAlgebra
using Optim, NLSolversBase
using Random, Statistics, StatsBase
using SpecialFunctions: loggamma

const ∞ = Inf

"""
    logmean(x;ϵ=1)

Compute the geometric mean: ``\exp\\left(\\langle \\log\\left(x + \\epsilon \\right) \\rangle\\right) - \\epsilon``
"""
const logmean(x;ϵ=1) = exp.(mean(log.(x.+ϵ)))-ϵ

"""
    logmean(x;ϵ=1)

Compute the geometric variance: ``exp\\left(\\langle \\log\\left(x + \\epsilon \\right)^2 \\rangle_c\\right) - \\epsilon``
"""
const logvar(x;ϵ=1)  = exp.(var(log.(x.+ϵ)))-ϵ

"""
    negativebinomial(count, depth)

Compute the log likelihood of a negative binomial generalized linear model (GLM) with log link function for `count` of a single gene.
The sequencing depth for each sequenced cell is assumed to be the only confounding variables.
"""
function negativebinomial(count, depth)
    loglikelihood = function(Θ)
        Θ₁,Θ₂,Θ₃ = Θ
        return -sum(
            loggamma(n+Θ₃)
          - loggamma(n+1)
          - loggamma(Θ₃)
          + n*(Θ₁+Θ₂*d)
          + Θ₃*log(Θ₃)
          - (n+Θ₃)*log(exp(Θ₁+Θ₂*d)+Θ₃)
          for (n,d) ∈ zip(count,depth)
        )
    end

    μ  = logmean(count)
    Θ₀ = [log(μ), 1, logvar(count)/μ - 1]
    if Θ₀[end] < 0 || isinf(Θ₀[end]) || isnan(Θ₀[end])
        Θ₀[end] = 1
    end

    return (
        Θ₀ = Θ₀,
        likelihood = TwiceDifferentiable(loglikelihood, Θ₀; autodiff=:forward),
        constraint = TwiceDifferentiableConstraints([-∞,-∞,0],[+∞,+∞,+∞]),
        residual   = function(Θ)
            Θ₁,Θ₂,Θ₃ = Θ
            μ = @. exp(Θ₁ + Θ₂*depth)
            σ = @. √(μ + μ^2/Θ₃)
            z = @. (count - μ) / σ

            z[ z .< -5 ] .= -5
            z[ z .> +5 ] .= +5

            return z
        end
    )
end

"""
    gamma(count, depth)

Compute the log likelihood of a gamma distributed generalized linear model (GLM) with log link function for the estimated mean `count` of a single gene.
The sequencing depth for each sequenced cell is assumed to be the only confounding variables.
"""
function gamma(count, depth)
    loglikelihood = function(Θ)
        Θ₁,Θ₂,Θ₃ = Θ
        return -sum(
            let
                α = Θ₃*exp(Θ₁+Θ₂*d)
                α*log(Θ₃) + (α-1)*log(x) - (Θ₃*x) - loggamma(α)
            end for (x,d) ∈ zip(count,depth)
        )
    end

    μ  = logmean(count; ϵ = 1e-10)
    Θ₀ = [log(μ), 1.25, μ/logvar(count; ϵ=1e-10)]
    if Θ₀[end] < 0 || isinf(Θ₀[end]) || isnan(Θ₀[end])
        Θ₀[end] = 1
    end

    return (
        Θ₀ = Θ₀,
        likelihood = TwiceDifferentiable(loglikelihood, Θ₀; autodiff=:forward),
        constraint = TwiceDifferentiableConstraints([-∞,-∞,0],[+∞,+∞,+∞]),
        residual   = function(Θ)
            Θ₁, Θ₂, Θ₃ = Θ
            μ = @. exp(Θ₁ + Θ₂*depth)
            σ = @. sqrt(μ / Θ₃)
            z = @. (count - μ) / σ

            z[ z .< -5 ] .= -5
            z[ z .> +5 ] .= +5

            return z
        end
    )
end

"""
    prior(params)

Estimate a generalized normal distribution to `params` by maximum likelihood estimation.
In practice, used to compute the empirical prior for overdispersion factor `Θ₃` in the negative binomial.
"""
function prior(params)
    objective = function(Θ)
        Θ₁, Θ₂, Θ₃ = Θ
        return sum(
            loggamma(1/Θ₃) + log(Θ₂) .+ (abs.(params .- Θ₁)./ Θ₂).^Θ₃ .- log(Θ₃)
        )
    end

    Θ₁ = mean(params)
    Θ₂ = std(params)
    Θ₀ = [Θ₁,Θ₂,2]

    likelihood = TwiceDifferentiable(f,Θ₀;autodiff=:forward)
    constraint = TwiceDifferentiableConstraints([0,0,0],[+∞,+∞,+∞])
    hyperparam = optimize(likelihood, constraint, Θ₀, IPNewton())

    return Optim.minimizer(hyperparam)
end

"""
    FitType = NamedTuple{
    (
     :likelihood,
     :parameters,
     :uncertainty,
     :residual
    ),
    Tuple{Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}}

Stores the result of MLE fit of _one_ gene.
"""
const FitType = NamedTuple{
    (
     :likelihood,
     :parameters,
     :uncertainty,
     :residual
    ),
    Tuple{Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}}
}

"""
    fit(stochastic, count, depth)

Fit a generative model `stochastic` to gene expression `count` data, assuming confounding sequencing `depth`.
`stochastic` can be either [`negativebinomial`](@ref) or [`gamma`](@ref).
"""
function fit(stochastic, count, depth)
    model = stochastic(count, depth)
    param = optimize(model.likelihood, model.constraint, model.Θ₀, IPNewton())

    return (
        likelihood  = Optim.minimum(param),
        parameters  = Optim.minimizer(param),
        uncertainty = diag(inv(hessian!(model.likelihood, Optim.minimizer(param)))),
        residual    = model.residual(Optim.minimizer(param))
    )
end

"""
    bootstrap(count, depth; stochastic=negativebinomial, samples=50)

Empirically verify the MLE fit of `count`, using a GLM model generated by `stochastic` with confounding `depth` variables by bootstrap.
One third of cells are removed and the parameters are re-estimated with the remaining cells.
This process is repeated `samples` times.
The resultant distribution of estimation is returned.
"""
function bootstrap(count, depth; stochastic=negativebinomial, samples=50)
    N = length(depth)

    Θ₁ = Array{Float64}(undef,samples)
    Θ₂ = Array{Float64}(undef,samples)
    Θ₃ = Array{Float64}(undef,samples)
    δL = Array{Float64}(undef,samples)

    for n in 1:samples
        ι = randperm(N)[1:2*N÷3]
        f = fit(stochastic, count[ι],depth[ι])

        δL[n] = f.likelihood
        Θ₁[n], Θ₂[n], Θ₃[n] = f.parameters
    end

    return Θ₁, Θ₂, Θ₃, δL
end

"""
    glm(data; stochastic=negativebinomial, ϵ=1)

Fit a generalized linear model (GLM) to the matrix `data`.
Genes are assumed to be on rows, cells over columns.
The underlying generative model is passed by `stochastic`.
"""
function glm(data; stochastic=negativebinomial, ϵ=1)
    average(x) = logmean(x; ϵ=ϵ)
    depth = map(eachcol(data)) do col
        col |> vec |> average
    end

    fits = Array{NamedTuple}(undef,size(data,1))
    Threads.@threads for (i,gene) in collect(enumerate(eachrow(data)))
        fits[i] = fit(stochastic,vec(gene),depth)
    end

    return (
        likelihood = map((f)->f.likelihood,  fits),
        residual   = Matrix(reduce(hcat, map((f)->f.residual, fits))'),

        Θ₁  = map((f)->f.parameters[1],  fits),
        Θ₂  = map((f)->f.parameters[2],  fits),
        Θ₃  = map((f)->f.parameters[3],  fits),

        δΘ₁ = map((f)->f.uncertainty[1], fits),
        δΘ₂ = map((f)->f.uncertainty[2], fits),
        δΘ₃ = map((f)->f.uncertainty[3], fits),
    )
end

end

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
export AbstractThresholdLink
export cdf, invcdf
export NormalThresholdLink, StudentThresholdLink

"Abstract type for link functions between probabilities and values on the whole real line"
abstract type AbstractThresholdLink{TUnit, TReal} end

"""
    cdf(atl::AbstractThresholdLink)

Get the link function from real line to probability for the link specification `atl`
"""
function cdf(atl::AbstractThresholdLink) end

"""
    invcdf(atl::AbstractThresholdLink)

Get the link function from probability to real line for the link specification `atl`
"""
function invcdf(atl::AbstractThresholdLink) end

"Normal distribution link specification using StatsFuns implementation of normal cdf and invcdf"
struct NormalThresholdLink{T <: AbstractFloat} <: AbstractThresholdLink{T, T}
end

"""
    cdf(ntl::NormalThresholdLink)

Get the link function from real line to probability for `ntl` which is always the normal cdf
"""
function cdf(ntl::NormalThresholdLink)
    return normcdf
end

"""
    invcdf(ntl::NormalThresholdLink)

Get the link function from probability to real line for `ntl` which is always the normal inverse cdf
"""
function invcdf(ntl::NormalThresholdLink)
    return norminvcdf
end

"Student's distribution link specification using StatsFuns implementation of Student's cdf and invcdf"
struct StudentThresholdLink <: AbstractThresholdLink{Float64, Float64}
    df::Union{Float64, Int64}

    """
        StudentThresholdLink(df::Union{Float64, Int64})

    Construct a StudentThresholdLink with `df` degrees of freedom
    """
    function StudentThresholdLink(df::Union{Float64, Int64})
        if df <= 0
            throw(DomainError(df, "Degrees of freedom was not positive"))
        end
        return new(df)
    end
end

"""
    cdf(stl::StudentThresholdLink)

Get the link function from real line to probability for `stl` which is the Student's cdf with degrees
    of freedom specified by `stl`
"""
function cdf(stl::StudentThresholdLink)
    return q -> tdistcdf(stl.df, q)
end

"""
    invcdf(stl::StudentThresholdLink)

Get the link function from probability to real line for `stl` which is the Student's inverse cdf with degrees
    of freedom specified by `stl`
"""
function invcdf(stl::StudentThresholdLink)
    return p -> tdistinvcdf(stl.df, p)
end


text:
using Pollen
using DataAugmentation
using FilePathsBase

project = Pollen.documentationproject(DataAugmentation)
Pollen.fullbuild(project, Pollen.FileBuilder(Pollen.HTML(), p"dev/"))


text:
using Plots, SimpleDrawing

export life_run, life_movie

"""
`life_run(A)` runs the game of life starting with the matrix `A`. It is
visualized on the screen. It will run forever unless:
* there are no live cells remaining,
* the current state is the same as the previous state, or
* the current state is the same as the state two steps earlier.
If `life_run` halts, it returns the number of iterations.

Optional named arguments:
* `pause=0.0`: number of extra seconds between frames
* `wrap=false`: determine if the board wraps (is toroidal)
* `counter=false`: show the step number below the image
* `max_steps`: maximum number of steps
"""
function life_run(A::Matrix{Int};
        pause=0.0,
        wrap::Bool=false,
        max_steps::Int=typemax(Int),
        counter::Bool=false
    )
    @assert _life_check(A) "Matrix may only contain 0s and 1s"
    step = 0
    back_1 = 0*A    # previous A
    back_2 = 0*A    # previous previous A

    while true
        p = my_spy(A)
        if counter
            plot!(xlabel=string(step))
        end
        display(p)
        sleep(pause)

        back_2 = back_1
        back_1 = A
        A = life_step(A,wrap)

        if A==back_1 || A==back_2 || sum(A) == 0
            return step
        end
        if step >= max_steps
            return step
        end
        step += 1

    end
end

"""
`life_movie(A)` creates a GIF animation of the game of life.

Optional named arguments:
* `file_name="life.gif"`: name of file to hold the animation
* `rate=5`: frame rate (frames per second)
* `wrap=false`: determine if the board wraps (is toroidal)
* `max_steps`: maximum number of steps; if omitted, runs forever
"""
function life_movie(A::Matrix{Int};
        file_name="life.gif",
        wrap::Bool=false,
        max_steps::Int=typemax(Int),
        rate = 5
    )
    @assert _life_check(A) "Matrix may only contain 0s and 1s"
    step = 0
    back_1 = 0*A    # previous A
    back_2 = 0*A    # previous previous A

    movie = @animate for k=1:max_steps
        print(".")
        p = my_spy(A)
        back_2 = back_1
        back_1 = A
        A = life_step(A,wrap)

        if A==back_1 || A==back_2 || sum(A) == 0
            break
        end
        if step >= max_steps
            break
        end
        step += 1

    end
    println()
    G = gif(movie, file_name, fps = rate)
    return G
end


text:
include("simulate.jl")

using Ipopt
using LazyArrays
using GenericLinearAlgebra
using BenchmarkTools
using SparseArrays

function assembleFdiff(x::AbstractArray{T}, u::AbstractArray{S}, j::Vector{Joint}, GravityInInertial:: Vector{Float64}) where {T<:Real, S<:Real}
    maxBodyID = j[end].RB2.bodyID
    F = Vector(undef,7*(maxBodyID-1))

    for i=1:length(j)
        k = j[i].RB2.bodyID - 1
        F[7*(k-1)+1:7*k] =
        genExtF(x,j[i].RB2,u[:,k+1],GravityInInertial)

        if j[i].type == "Spring"
            F_b1, F_b2 = ForceConSpr(x,j[i]); # Forces acting on the bodies connected by the spring. Note that a spring does not constrain 2 bodies any more than they previously were, i.e., degrees of freedom remain the same.

            b1_id = j[i].RB1.bodyID -1; b2_id = j[i].RB2.bodyID - 1;
            if j[i].RB1.m == 0.0 # First body is the inertial frame
                F[7*(b2_id-1)+1:7*b2_id] += F_b2;
            else # First body is not the inertial frame
                F[7*(b1_id-1)+1:7*b1_id] += F_b1;
                F[7*(b2_id-1)+1:7*b2_id] += F_b2;
            end
        end
    end
    return F
end
## Ipopt directly

# function trim_kronLazy(j::Vector{Joint},GravityInInertial::Vector{Float64};
#   ix::Vector{Integer}=zeros(Integer,20*(nB-1)+1),
#   iy::Vector{Integer}=zeros(Integer,7*(nB-1)))
#
#   nB = length(j) + 1; # number of bodies in the joint tree
#   # ix: 0 indicates that the corresponding variable is free during optimization
#
#   # iy: 0 indicates that the corresponding second derivative term needs to be zero
#   # (specify indices for iy starting from second body)
#   # iy = [1;4;5] means that the ẍ. β̈1 and β̈2 are free for j[1].RB2
#
#   # Find which variables are free for optimization
#   freeVarInd = findall(x->x==0,ix); # indices of free variables
#
#   # find which acceleration-level terms are desired to be 0/constrained
#   freeAccInd = findall(x->x==1,iy) # indices of free accelerations
#
#   # Initial guesses
#   # Initial guess is provided initial location
#   x0 = []; # Initial guess
#   for k = 1:length(j)
#     append!(x0,j[k].RB2.x)
#   end
#   extFListCopy = externalFTotal(0.0,j);
#   u0 = genU(extFListCopy); # Initial guess for control input
#
#   n_xu = length(freeVarInd)-1; # Number of optimization variables
#   lb = [-Inf*ones(n_xu);0.0]; ub = [Inf*ones(n_xu);Inf] # Bounds on optimization variables
#   # lb[4] = 0.0;
#   # lb[14+4] = 0.0;
#
#   m = 1; # Number of nonlinear constraints (xdot = f(x,u)) for second order variables
#   lb_g = [-Inf]; ub_g = [0.0];

j = qRotor_j;
x0Orig, u0Orig = getXU_0(j);
nB = length(j) + 1;
ix= zeros(Integer,20*(nB-1)+1);
# WITH U=0
ix[14*(nB-1)+1:20*(nB-1)] .= 1;
freeVarInd = findall(x->x==0,ix); # indices of free variables
# extract states and control from optimization variables
function getXU(z::AbstractArray{T}) where T<:Real
  z0 = [x0;u0[:,2:end][:];0.0];
  z0[freeVarInd] .= 0;

  Mz = zeros(length(z0), length(freeVarInd))
  for (i,ind) in enumerate(freeVarInd)
  Mz[ind,i] = 1;
  end

  zFull = z0 + Mz*z;

  X = zFull[1:14*(nB-1)];
  U = zFull[14*(nB-1)+1:end-1]
  return X,U
end

## Constraint functions
function qddFn(xuVar)
  x = xuVar[1:14*nB]; u = reshape(xuVar[14*nB+1:end],(6,nB));
  xdot = fxdot(x,u,j,g);
  y = Vector{Float64}(undef,7*nB); # qdd
  for i=1:nB
    y[7*(i-1)+1:7*(i)] = xdot[7*(2*i-1)+1:14*i]
  end
  deleteat!(y,1:7); # inertial frame has no acceleration
  # deleteat!(y,freeAccInd); # free acceleration variables
  return sparse(y)
end

function quatNormConstr(xuVar)
  x = xuVar[1:14*nB]; u = reshape(xuVar[14*nB+1:end],(6,nB));
  q_constr = Vector(undef, nB)
  for i=1:nB
    q_constr[i] = norm(x[14*(i-1)+4:14*(i-1)+7])^2 - 1;
  end
  return q_constr
end

function jointLocConstr(xuVar)
  x = xuVar[1:14*nB]; u = reshape(xuVar[14*nB+1:end],(6,nB));
  jLoc_constr = Matrix(undef,3,length(j))
  for i=1:length(j)
    b1_id = j[i].RB1.bodyID; b2_id = j[i].RB2.bodyID;
    x1 = x[(b1_id-1)*14+1:b1_id*14]
    x2 = x[(b2_id-1)*14+1:b2_id*14]
    j_b1 = x1[1:3] + transpose(quat2dcm(x1[4:7]))*j[i].pos1;
    j_b2 = x2[1:3] + transpose(quat2dcm(x2[4:7]))*j[i].pos2;
    jLoc_constr[:,i] = j_b1 - j_b2;
  end
  return jLoc_constr
end

function _nonlinear_constraints(z::AbstractArray{T}) where T<:Real
  xVar, uVar = getXU(z)
  x = [Vector{T}(j[1].RB1.x); xVar]
  u = [zeros(T,6) reshape(uVar,(6, nB-1))]
  xuVar = [x;u[:]]; # Vector of x,u

  # Dynamics constraint
  y = qddFn(xuVar);
  # Quaternion norm constraint
  q_constr = quatNormConstr(xuVar);
  # JointLocation constraint  (translation constraint)
  jLoc_constr = jointLocConstr(xuVar);

  constr = [y;q_constr;jLoc_constr[:]];

  J = sum(constr.*constr) - z[end];
  return J
end

## Helper functions for derivatives
function kronSum(A,B)
  # kronecker sum of 2 arrays
  out = kronProd(A, Matrix(I,(size(B)))) + kronProd(Matrix(I,(size(A))),B)
  return sparse(out)
end

function kronProd(A,B)
  # K = ApplyArray(kron,A[:,:],B[:,:]);
  # # out = Matrix{Float64}(undef,size(K)...); copyto!(out,K);
  # out = Matrix{Float64}(undef,size(K)...); copyto!(out,K)
  # return sparse(out)

  K = kron(sparse(A),sparse(B));
  return sparse(K)
end

function createCommMat(M)
  # commutation matrix
  r,m = size(M);
  return createCommMat(r,m);
  # K = spzeros(m*r,m*r);
  # Ir = sparse(I,r,r); Im = sparse(I,m,m);
  # for i=1:r
  #   for j=1:m
  #       ei = Ir[:,i]; # ith-canonical unit vector of dimension r
  #       ej = Im[:,j];# jth-canonical unit vector of dimension m
  #       K += kronProd(ei*permutedims(ej), ej*permutedims(ei));
  #   end
  # end
  # return sparse(K)
end

function createCommMat(r,m)
  # commutation matrix for a matrix of size r,ms
  K = spzeros(m*r,m*r);
  Ir = sparse(1.0I,r,r); Im = sparse(1.0I,m,m);
  for i=1:r
    for j=1:m
        ei = @view Ir[:,i]; # ith-canonical unit vector of dimension r
        ej = @view Im[:,j];# jth-canonical unit vector of dimension m
        K += kronProd(ei*permutedims(ej), ej*permutedims(ei));
    end
  end
  return sparse(K)
end

function fdJ_kronAB(A,B, dA, dB)
  n,q = size(A[:,:]); p,r = size(B[:,:]);
  # Iq = 1.0I(q); Ip = 1.0I(p); In = 1.0I(n); Ir = 1.0I(r);
  Iq = sparse(1.0I,q,q); Ip = sparse(1.0I,p,p); In = sparse(1.0I,n,n); Ir = sparse(1.0I,r,r);
  Inq = sparse(1.0I,n,q); Ipr = sparse(1.0I,p,r)
  Krn = sparse(createCommMat(r,n));
  Ay1 = sparse(kronProd(Krn,Ip)); Ay2 = sparse(kronProd(In,vec(B)));
  # println("=========")
  # println("typeof(Ay1) = ", typeof(Ay1))
  # println("typeof(Ay2) = ", typeof(Ay2))
  # println("Ay1 = ", Ay1);
  # println("=========")
  # println("Ay2 = ", Ay2);

  Ay = sparse(kronProd(Iq, Ay1*Ay2));
  Bx1 = sparse(kronProd(Iq,Krn)); Bx2 = sparse(kronProd(vec(A),Ir));
  Bx = sparse(kronProd(Bx1*Bx2, Ip));
  out = sparse(Ay*dA + Bx*dB)
  return out
end

xv = x0Orig;
uv = u0Orig;
xu_var = [xv;uv[:]];

# function gh_nc(xv, uv)
xu_var = [xv;uv[:]];
function Mdiff1(x)
  # First derivatives of mass matrix and powers (1/2, -1/2, -1)
  M = Array(assembleM(x,j));
  Ms = (real(sqrt(M)));
  M_inv = (inv(M))
  Ms_inv = (real(sqrt(M_inv)));

  dM = sparse(ForwardDiff.jacobian(z -> assembleM(z,j),x))
  dM_inv = sparse(-(kronProd(permutedims(M_inv),M_inv))*dM)
  dMs = sparse(inv(Array(kronSum(permutedims(Ms), Ms)))*dM)
  dMs_inv = sparse(-kronProd(permutedims(Ms_inv),Ms_inv)*(dMs));
  return sparse(M), sparse(Ms), sparse(Ms_inv), dM, dMs, dMs_inv, sparse(M_inv), dM_inv
end
M, Ms, Ms_inv, dM, dMs, dMs_inv, M_inv, dM_inv = Mdiff1(xv);
# @btime Mdiff1(xv);

function Abdiff1(x)
  A, b = Ab_VecOfMat(x,j);
  dA = ForwardDiff.jacobian(z->Ab_VecOfMat(z,j)[1],x);
  db = ForwardDiff.jacobian(z->Ab_VecOfMat(z,j)[2],x);
  # return A, dA, b, db
  return sparse(A), sparse(dA), sparse(b), sparse(db)
end
A, dA, b, db = Abdiff1(xv);
# @btime Abdiff1(xv);

function Gpdiff1(x)
  function t1diff1(x)
    # t1 = A'*(AM)^{-1}*A';
    sz_A1 = size(A,1)
    t1 = sparse(permutedims(A)*inv(Array(A*M_inv*permutedims(A))))

    dA_tr = sparse(createCommMat(A)*dA);

    Y = A*M_inv*permutedims(A); invY = inv(Array(Y));
    dY1 = sparse(kronProd(A*M_inv,1.0I(sz_A1))*dA);
    dY2 = sparse(kronProd(A,A)*dM_inv);
    dY3 = sparse(kronProd(1.0I(sz_A1), A*M_inv)*dA_tr);
    dY = dY1 + dY2 + dY3;
    dY_inv = sparse(-kronProd(permutedims(invY),invY)*dY);

    dt1_1 = sparse(kronProd(permutedims(invY),1.0I(size(A,2)))*dA_tr);
    dt1_2 = sparse(kronProd(1.0I(sz_A1),permutedims(A))*dY_inv;)
    dt1 = sparse(dt1_1 + dt1_2);
    return t1, dt1
  end
  G = Array(A*Ms_inv);
  Gp = sparse(permutedims(G)*((G*permutedims(G))\I(size(G,1))))

  t1, dt1 = t1diff1(x)
  dGp1 = sparse(kronProd(permutedims(t1),1.0I(size(M,1)))*dMs_inv)

  sz_t2 =size(t1,2)
  t2 = sparse(kronProd(1.0I(sz_t2), Ms_inv))
  dGp2 = sparse(t2*dt1);

  dGp = sparse(dGp1 + dGp2)
  return Gp, dGp
end
Gp, dGp = Gpdiff1(xv);
# @btime Gpdiff1(xv);

function Fudiff1(x,u)
  Fu = assembleFdiff(x,u,j,g);
  dFu = ForwardDiff.jacobian(y->assembleFdiff(y,u,j,g),x);
  dFu_u = ForwardDiff.jacobian(y->assembleFdiff(x,y,j,g),u);

  # Fu = convert(SparseVector{Float64}, sparse(Fu));
  # dFu = convert(SparseMatrixCSC{Float64}, sparse(dFu));
  # dFu_u = convert(SparseMatrixCSC{Float64}, sparse(dFu_u));
  return (Fu), (dFu), (dFu_u)
end
Fu, dFu, dFu_u = Fudiff1(xv,uv);
# @btime Fudiff1(xv,uv);

Fu = convert(SparseVector{Float64}, sparse(Fu));
dFu = convert(SparseMatrixCSC{Float64}, sparse(dFu));
dFu_u = convert(SparseMatrixCSC{Float64}, sparse(dFu_u));

function hFn(x,u)
  h = b - A*M_inv*Fu;
  return sparse(h)
end

function hdiff1(x,u)
  # has derivatives wrt both x,u
  h = hFn(x,u);

  # dh
  term1_1 = sparse(kronProd(permutedims(M_inv*Fu), 1.0I(size(A,1))));
  # println("typeof(term1_1) = ", typeof(term1_1))
  # println("typeof(dA) = ", typeof(dA))
  term1 = sparse(term1_1*dA);
  term2 = sparse(kronProd(permutedims(Fu), A)*dM_inv);
  term3 = sparse(A*M_inv*dFu);
  dh_x = sparse(db - (term1 + term2 + term3));

  # dh_u
  dh_u = sparse(-A*M_inv*dFu_u);

  return sparse(h), sparse(dh_x), sparse(dh_u)
end
h, dh_x, dh_u = hdiff1(xv,uv);
# @btime hdiff1(xv,uv);

function Fcdiff1(x,u)
  Fc = Ms*Gp*h

  dFc1 = kronProd(permutedims(Gp*h),I(size(M,1)))*dMs
  dFc2 = kronProd(permutedims(h),Ms)*dGp
  dFc3 = Ms*Gp*dh_x;

  dFc = sparse(dFc1 + dFc2 + dFc3);

  dFc_u = Ms*Gp*dh_u;
  return sparse(Fc), sparse(dFc), sparse(dFc_u)
end
Fc, dFc, dFc_u = Fcdiff1(xv,uv);
# @btime Fcdiff1(xv,uv);

function accdiff1(x,u)
  # derivative of qdd w.r.t x
  dacc1 = sparse(kronProd(transpose(Fc+Fu),1.0I(size(M_inv,1)))*dM_inv)
  dacc2 = sparse(M_inv*(dFc + dFu))
  dqdd_x = sparse(dacc1 + dacc2)

  return (dqdd_x)
end
dqdd_x = accdiff1(xv,uv);
# @btime accdiff1(xv,uv);

# Constraint derivative functions for gradient

function dqddFn(xuVar)
  dqdd_u = M_inv*(dFc_u + dFu_u);
  out = [dqdd_x dqdd_u];
  return sparse(out)
end
qdd = qddFn(xu_var);
qdd = convert(SparseVector{Float64}, sparse(qdd));

# @btime qddFn(xu_var);
dqdd = dqddFn(xu_var);
# @btime dqddFn(xu_var);

function dynConstr_jac(xuVar)
  out = (2*permutedims(qdd)*dqdd)[:];
  return out
end

function quatNormConstr_diff(xuVar)
  f(y) = sum(quatNormConstr(y).*quatNormConstr(y));
  df(x) = ForwardDiff.gradient(y->f(y), x);
  grad_f = df(xuVar);
  d2f = ForwardDiff.jacobian(y->df(y), xuVar);
  return grad_f, d2f
end

function jointLocConstr_diff(xuVar)
  f(y) = sum((jointLocConstr(y)[:]).*(jointLocConstr(y)[:]));
  df(x) = ForwardDiff.gradient(y->f(y), x);
  grad_f = df(xuVar);
  d2f = ForwardDiff.jacobian(y->df(y), xuVar);
  return grad_f, d2f
end

function grad_nonlinear_constraints()
  dqdd_constr = dynConstr_jac(xu_var); # Dyn
  dquatNorm_constr,_ = quatNormConstr_diff(xu_var); # Quat
  djointLoc_constr, _ = jointLocConstr_diff(xu_var); # JointLoc

  dconstr_z = dqdd_constr + dquatNorm_constr + djointLoc_constr;

  # take care of free optimization variables
  freeVarInd_grad = 14 .+ freeVarInd[1:end-1];
  grad_constr_freeInd = dconstr_z[freeVarInd_grad];
  out = [grad_constr_freeInd;-1.0];
  return out
end
grad_nc = grad_nonlinear_constraints();
# @btime grad_nonlinear_constraints();

##  Derivatives for hessian
## Mdiff2
# dM_fn(y) = ForwardDiff.jacobian(z->assembleM(z,j),y);
# dM_xx = sparse(ForwardDiff.jacobian(z->dM_fn(z),xv));
#
# # function z1Fn(x)
# z1= (kronSum(permutedims(Ms), Ms));
# z1_inv = sparse(inv(Array(z1)));
# dMs_t = sparse(createCommMat(Ms)*dMs);
# dz1_1 = fdJ_kronAB(permutedims(Ms),1.0I(size(Ms,1)),dMs_t,zeros(size(dMs_t)));
# dz1_2 = fdJ_kronAB(1.0I(size(Ms,2)),Ms,zeros(size(dMs)),dMs);
# dz1 = dz1_1 + dz1_2;
#
# dz1_inv_1 = kronProd(permutedims(z1_inv),z1_inv);
# dz1_inv = -(kronProd(permutedims(z1_inv),z1_inv))*dz1;
#   # return z1_inv, dz1_inv
# # end
# # z1_inv, dz1_inv = z1Fn(x);
# println("=====z1Done=======");
# out1 = kronProd(permutedims(dM), 1.0I(size(z1_inv,1)))*dz1_inv;
# println("=====out1 Done=======");
# out2 = kronProd(1.0I(size(dM,2)), z1_inv)*dM_xx;
# dMs_xx = out1 + out2;

function Mdiff2(x)
  # dM_xx
  dM_fn(y) = ForwardDiff.jacobian(z->assembleM(z,j),y);
  dM_xx = sparse(ForwardDiff.jacobian(z->dM_fn(z),x));

  # dMs_xx
  function z1Fn(x)
    z1= kronSum(permutedims(Ms), Ms);
    z1_inv = sparse(inv(Array(z1)));
    dMs_t = sparse(createCommMat(Ms)*dMs);
    dz1_1 = fdJ_kronAB(permutedims(Ms), sparse(1.0I, size(Ms)...), dMs_t, spzeros(size(dMs_t)...));
    dz1_2 = fdJ_kronAB(sparse(1.0I, size(Ms)...), Ms, spzeros(size(dMs)...), dMs);
    dz1 = sparse(dz1_1 + dz1_2);
    dz1_inv = sparse(-(kronProd(permutedims(z1_inv),z1_inv))*dz1);
    return z1_inv, dz1_inv
  end
  z1_inv, dz1_inv = z1Fn(x);
  out1 = sparse(kronProd(permutedims(dM), 1.0I(size(z1_inv,1)))*dz1_inv);
  out2 = sparse(kronProd(1.0I(size(dM,2)), z1_inv)*dM_xx);
  dMs_xx = sparse(out1 + out2);

  # dMsInv_xx
  R = kronProd(permutedims(Ms_inv), Ms_inv);
  dMs_inv_t = sparse(createCommMat(Ms_inv)*dMs_inv);
  dR = fdJ_kronAB(permutedims(Ms_inv), Ms_inv, dMs_inv_t, dMs_inv);

  out1 = sparse(kronProd(permutedims(dMs),sparse(1.0I(size(R,1))))*dR);
  out2 = kronProd(sparse(1.0I(size(dMs,2))),R)*dMs_xx;
  dMsInv_xx = sparse(-(out1 + out2));

  # dMinv_xx
  Q = kronProd(permutedims(M_inv), M_inv);
  dM_inv_t = sparse(createCommMat(M_inv)*dM_inv);
  dQ = fdJ_kronAB(permutedims(M_inv), M_inv, dM_inv_t, dM_inv);
  out1 = sparse(kronProd(permutedims(dM), sparse(1.0I(size(Q,1))))*dQ);
  out2 = sparse(kronProd(sparse(1.0I(size(dM,2))), Q)*dM_xx);
  dMinv_xx = sparse(-(out1 + out2));

  return dM_xx, dMs_xx, dMsInv_xx, dMinv_xx
end
dM_xx, dMs_xx, dMsInv_xx, dMinv_xx = Mdiff2(xv);
# @btime Mdiff2(xv);

#
function Abdiff2(x)

  dAFn(y) = ForwardDiff.jacobian(z->Ab_VecOfMat(z,j)[1],y);
  cfg1 = ForwardDiff.JacobianConfig(dAFn, x, ForwardDiff.Chunk{1}());
  dA_xx = ForwardDiff.jacobian(dAFn,x, cfg1);
  dbFn(y) = ForwardDiff.jacobian(z->Ab_VecOfMat(z,j)[2],y);
  cfg1 = ForwardDiff.JacobianConfig(dbFn, x, ForwardDiff.Chunk{1}());
  db_xx = ForwardDiff.jacobian(dbFn,x, cfg1);
  return sparse(dA_xx), sparse(db_xx);
end
dA_xx, db_xx = Abdiff2(xv);
# @btime Abdiff2(xv);
# A, dA = Adiff1(xv);

function Gdiff2(x,u)
  G = sparse(A*Ms_inv);
  #dG
  dG1 = sparse(kronProd(permutedims(Ms_inv),1.0I(size(A,1)))*dA);
  dG2 = sparse(kronProd(1.0I(size(Ms_inv,2)),A)*dMs_inv);
  dG = sparse(dG1 + dG2);

  # dG_xx
  Q1 = kronProd(permutedims(Ms_inv), 1.0I(size(A,1)));
  dMs_inv_t = sparse(createCommMat(Ms_inv)*dMs_inv);
  dQ1 = fdJ_kronAB(permutedims(Ms_inv), 1.0I(size(A,1)),dMs_inv_t,zeros(size(A,1)^2,length(x)));

  Q2 = kronProd(1.0I(size(Ms_inv,2)), A);
  dQ2 = fdJ_kronAB(1.0I(size(Ms_inv,2)), A, zeros(size(dMs_inv)), dA);

  out1 = sparse(kronProd(permutedims(dA),1.0I(size(Q1,1)))*dQ1);
  out2 = sparse(kronProd(1.0I(size(dA,2)), Q1)*dA_xx);
  out3 = sparse(kronProd(permutedims(dMs_inv), 1.0I(size(Q2,1)))*dQ2);
  out4 = sparse(kronProd(1.0I(size(dMs_inv,2)), Q2)*dMsInv_xx);

  dG_xx = sparse(out1 + out2 + out3 + out4);
  return G, dG, dG_xx
end
G, dG, dG_xx = Gdiff2(xv,uv);
# @btime Gdiff2(xv,uv);

function Hdiff2(x,u)
  # H = G*G^T
  dG_t = sparse(createCommMat(G)*dG);

  H = sparse(G*permutedims(G));

  dH = sparse(kronProd(G,1.0I(size(G,1)))*dG + kronProd(1.0I(size(G,1)),G)*dG_t);

  P1 = kronProd(G,1.0I(size(G,1)));
  dP1 = fdJ_kronAB(G, 1.0I(size(G,1)), dG, zeros(size(G,1)^2,length(x)));

  P2 = kronProd(1.0I(size(G,1)), G);
  dP2 = fdJ_kronAB(1.0I(size(G,1)), G, zeros(size(G,1)^2,length(x)), dG);

  dG2_t = sparse(kronProd(1.0I(size(dG,2)), createCommMat(G))*dG_xx);

  out1 = sparse(kronProd(permutedims(dG), 1.0I(size(P1,1)))*dP1);
  out2 = sparse(kronProd(1.0I(size(dG,2)), P1)*dG_xx);
  out3 = sparse(kronProd(permutedims(dG_t), 1.0I(size(P2,1)))*dP2);
  out4 = sparse(kronProd(1.0I(size(dG_t,2)), P2)*dG2_t);

  dH_xx = sparse(out1 + out2 + out3 + out4);
  return H, dH, dH_xx
end
H, dH, dH_xx = Hdiff2(xv,uv);
# @btime Hdiff2(xv,uv);

function Gpdiff2(x,u)
  dG_t = sparse(createCommMat(G)*dG);
  dG2_t = sparse(kronProd(1.0I(size(dG,2)), createCommMat(G))*dG_xx);

  H_inv = sparse(inv(Array(H)));
  dH_inv = sparse(-kronProd(permutedims(H_inv), H_inv)*dH);
  dH_inv_t = sparse(createCommMat(H_inv)*dH_inv);

  W = kronProd(permutedims(H_inv), H_inv);
  dW = fdJ_kronAB(permutedims(H_inv), H_inv, dH_inv_t, dH_inv);
  h1 = sparse(kronProd(permutedims(dH),1.0I(size(W,1)))*dW);
  h2 = sparse(kronProd(1.0I(size(dH,2)),W)*dH_xx);
  dH_inv2 = sparse(-(h1 + h2));

  V1 = kronProd(permutedims(H_inv), 1.0I(size(G,2)));
  dV1 = fdJ_kronAB(permutedims(H_inv), 1.0I(size(G,2)), dH_inv_t, zeros(size(G,2)^2, length(x)));

  V2 = kronProd(1.0I(size(H_inv,2)), permutedims(G));
  dV2 = fdJ_kronAB(1.0I(size(H_inv,2)),permutedims(G), zeros(size(dH_inv_t)), dG_t);

  out1 = sparse(kronProd(permutedims(dG_t), 1.0I(size(V1,1)))*dV1);
  out2 = sparse(kronProd(1.0I(size(dG_t,2)), V1)*dG2_t);
  out3 = sparse(kronProd(permutedims(dH_inv), 1.0I(size(V2,1)))*dV2);
  out4 = sparse(kronProd(1.0I(size(dH_inv,2)), V2)*dH_inv2);

  dGp_xx = sparse(out1 + out2 + out3 + out4);

  return dGp_xx;
end
dGp_xx = Gpdiff2(xv,uv);
# @btime Gpdiff2(xv,uv);

function Fudiff2(x,u)
  dFu_xx = ForwardDiff.jacobian(y->Fudiff1(y,u)[2],x);
  dFu_ux = ForwardDiff.jacobian(y->Fudiff1(x,y)[2],u);
  dFu_xu = ForwardDiff.jacobian(y->Fudiff1(y,u)[3],x);
  dFu_uu = ForwardDiff.jacobian(y->Fudiff1(x,u)[3],u);
  return sparse(dFu_xx), sparse(dFu_ux), sparse(dFu_xu), sparse(dFu_uu)
end
dFu_xx, dFu_ux, dFu_xu, dFu_uu = Fudiff2(xv,uv);
# @btime Fudiff2(xv,uv);

Fu = convert(SparseVector{Float64}, sparse(Fu));
dFu = convert(SparseMatrixCSC{Float64}, sparse(dFu));
dFu_u = convert(SparseMatrixCSC{Float64}, sparse(dFu_u));
# #
dFu_xx = convert(SparseMatrixCSC{Float64}, dFu_xx);
dFu_ux = convert(SparseMatrixCSC{Float64}, dFu_ux);
dFu_xu = convert(SparseMatrixCSC{Float64}, dFu_xu);
dFu_uu = convert(SparseMatrixCSC{Float64}, dFu_uu);

function hdiff2(x,u)
  # dh_xx
  # terms 1,2
  T1 = kronProd(permutedims(M_inv*Fu), 1.0I(size(A,1)))
  dMinvFu = sparse(kronProd(permutedims(Fu), 1.0I(size(M_inv,1)))*dM_inv + M_inv*dFu);
  dMinvFu_t = dMinvFu;

  dT1 = fdJ_kronAB(permutedims(M_inv*Fu), Matrix(1.0I(size(A,1))), dMinvFu_t, Matrix(zeros(size(A,1)^2, length(x))));
  term1 = sparse(kronProd(permutedims(dA), 1.0I(size(T1,1)))*dT1)
  term2 = sparse(kronProd(1.0I(size(dA,2)), T1)*dA_xx);

  # terms 3,4
  T2 = kronProd(permutedims(Fu), A);
  dT2 = fdJ_kronAB(permutedims(Fu), A, dFu, dA);
  term3 = sparse(kronProd(permutedims(dM_inv), 1.0I(size(T2,1)))*dT2);
  term4 = sparse(kronProd(1.0I(size(dM_inv,2)), T2)*dMinv_xx);

  # terms 5,6
  T3 = sparse(A*M_inv);
  dT3 = sparse(kronProd(permutedims(M_inv), 1.0I(size(A,1)))*dA + kronProd(1.0I(size(M_inv,2)), A)*dM_inv);
  term5 = sparse(kronProd(permutedims(dFu), 1.0I(size(A,1)))*dT3);
  term6 = sparse(kronProd(1.0I(length(x)), T3)*dFu_xx);

  dh_xx = sparse(db_xx - (term1 + term2 + term3 + term4 + term5 + term6));

  # dh_ux
  E1 = kronProd(permutedims(M_inv*Fu), 1.0I(size(A,1)));
  dE1_u = fdJ_kronAB(permutedims(M_inv*Fu), 1.0I(size(A,1)), M_inv*dFu_u, zeros(size(A,1)^2, length(u)));

  E2 = kronProd(permutedims(Fu), A);
  dE2_u = fdJ_kronAB(permutedims(Fu), A, dFu_u, zeros(length(A), length(u)));

  term1 = sparse(kronProd(permutedims(dA),1.0I(size(E1,1)))*dE1_u);
  term2 = sparse(kronProd(permutedims(dM_inv), 1.0I(size(E2,1)))*dE2_u);
  term3 = sparse(kronProd(1.0I(size(dFu,2)), A*M_inv)*dFu_ux);

  dh_ux = sparse(-(term1 + term2 + term3));

  # dh_xu
  out1 = sparse(kronProd(permutedims(M_inv*dFu_u), 1.0I(size(A,1)))*dA);
  out2 = sparse(kronProd(permutedims(dFu_u), A)*dM_inv);
  out3 = sparse(kronProd(1.0I(length(u)), A*M_inv)*dFu_xu);
  dh_xu = sparse(-(out1 + out2 + out3));

  # dh_uu
  dh_uu = sparse(-kronProd(1.0I(length(u)), A*M_inv)*dFu_uu);

  return dh_xx, dh_ux, dh_xu, dh_uu
end
dh_xx, dh_ux, dh_xu, dh_uu = hdiff2(xv,uv);
h, dh_x, dh_u = hdiff1(xv,uv);

# @btime hdiff2(xv,uv);
# println("typeof(dh_xx) = ", typeof(dh_xx))
# println("typeof(dh_ux) = ", typeof(dh_ux))
# println("typeof(dh_xu) = ", typeof(dh_xu))
# println("typeof(dh_uu) = ", typeof(dh_uu))
# h, dh_x, dh_u = hdiff1(xv,uv);

function Fcdiff2(x,u)
  # dFc_xx
  function jacT1(x)
    D1 = kronProd(permutedims(Gp*h), 1.0I(size(Ms,1)));
    dGph = sparse(kronProd(permutedims(h), 1.0I(size(Gp,1)))*dGp + Gp*dh_x);
    dGph_t = dGph;
    dD1 = fdJ_kronAB(permutedims(Gp*h), 1.0I(size(Ms,1)), dGph_t, zeros(size(dMs)));

    dMs_t = sparse(createCommMat(Ms)*dMs);

    dT1_1 = sparse(kronProd(permutedims(dMs), 1.0I(size(D1,1)))*dD1);
    dT1_2 = sparse(kronProd(1.0I(size(dMs,2)), D1)*dMs_xx);

    dT1 = sparse(dT1_1 + dT1_2);
    return dT1
  end

  function jacT2(x)
    D2 = kronProd(permutedims(h), Ms);
    dD2 = fdJ_kronAB(permutedims(h), Ms, dh_x, dMs);
    # dGp_xx = Gpdiff2(x,u);

    dT2_1 = sparse(kronProd(permutedims(dGp),1.0I(size(D2,1)))*dD2);
    dT2_2 = sparse(kronProd(1.0I(size(dGp,2)), D2)*dGp_xx);
    dT2 = sparse(dT2_1 + dT2_2);
    return dT2
  end

  function jacT3(x)
    dT3_1 = sparse(kronProd(permutedims(Gp*dh_x), 1.0I(size(Ms,1)))*dMs);
    dT3_2 = sparse(kronProd(permutedims(dh_x), Ms)*dGp);
    dT3_3 = sparse(kronProd(1.0I(size(dh_x,2)), Ms*Gp)*dh_xx);
    dT3 = sparse(dT3_1 + dT3_2 + dT3_3);

    return dT3
  end
  dFc_xx = sparse(jacT1(x) + jacT2(x) + jacT3(x));

  # dFc_ux
  # # first term
  D1 = kronProd(permutedims(Gp*h), 1.0I(size(Ms,1)));
  dGph_u = sparse(Gp*dh_u);
  dGph_u_t = dGph_u;
  dD1_u = fdJ_kronAB(permutedims(Gp*h), 1.0I(size(Ms,1)), dGph_u_t, zeros(size(dMs,1),length(u)));
  out1 = sparse(kronProd(permutedims(dMs), 1.0I(size(D1,1)))*dD1_u);
  # # second term
  D2 = sparse(kronProd(permutedims(h), Ms)*dGp);
  dD2_u = fdJ_kronAB(permutedims(h), Ms, dh_u, zeros(size(dMs,1), length(u)));
  out2 = sparse(kronProd(permutedims(dGp), 1.0I(size(D2,1)))*dD2_u);
  # # third term
  out3 = sparse(kronProd(1.0I(size(dh_x,2)), Ms*Gp)*dh_ux);
  dFc_ux = sparse(out1 + out2 + out3);

  # dFc_xu
  out1 = sparse(kronProd(permutedims(Gp*dh_u), 1.0I(size(Ms,1)))*dMs);
  out2 = sparse(kronProd(permutedims(dh_u), Ms)*dGp);
  out3 = sparse(kronProd(1.0I(length(u)), Ms*Gp)*dh_xu);
  dFc_xu = sparse(out1 + out2 + out3);

  # dFc_uu
  dFc_uu = sparse(kronProd(1.0I(size(dh_u,2)), Ms*Gp)*dh_uu);

  return dFc_xx, dFc_ux, dFc_xu, dFc_uu
end
dFc_xx, dFc_ux, dFc_xu, dFc_uu = Fcdiff2(xv,uv);
# @btime Fcdiff2(xv,uv);
# println("typeof(dFc_xx) = ", typeof(dFc_xx))
# println("typeof(dFc_ux) = ", typeof(dFc_ux))
# println("typeof(dFc_xu) = ", typeof(dFc_xu))
# println("typeof(dFc_uu) = ", typeof(dFc_uu))

function hessqdd(x,u)
  T1 = kronProd(permutedims(Fc + Fu), 1.0I(size(M_inv,1)));
  dT1 = fdJ_kronAB(permutedims(Fc + Fu), 1.0I(size(M_inv,1)), (dFc + dFu), zeros(size(dM_inv)));
  T2 = sparse(dFc + dFu);
  dT2 = sparse(dFc_xx + dFu_xx);

  out1 = sparse(kronProd(permutedims(dM_inv), 1.0I(size(T1,1)))*dT1);
  out2 = sparse(kronProd(1.0I(size(dM_inv,2)), T1)*dMinv_xx);
  out3 = sparse(kronProd(permutedims(T2), 1.0I(size(M_inv,1)))*dM_inv);
  out4 = sparse(kronProd(1.0I(size(T2,2)), M_inv)*dT2);

  dqdd_xx = sparse(out1 + out2 + out3 + out4);
  return dqdd_xx
end
dqdd_xx = hessqdd(xv,uv);
# @btime hessqdd(xv,uv);
# println("typeof(dqdd_xx) = ", typeof(dqdd_xx))

function cross_qdd_ux(x,u)
  T1 = kronProd(permutedims(Fc + Fu), 1.0I(size(M_inv,1)));
  dT1_u = fdJ_kronAB(permutedims(Fc + Fu), 1.0I(size(M_inv,1)), dFc_u + dFu_u, zeros(size(M_inv,1)^2, length(u)));

  T2 = sparse(dFc + dFu);
  dT2_u = sparse(dFc_ux + dFu_ux);

  out1 = sparse(kronProd(permutedims(dM_inv), 1.0I(size(T1,1)))*dT1_u);
  out2 = sparse(kronProd(1.0I(length(x)), M_inv)*dT2_u);
  dqdd_ux = sparse(out1+out2);
  return dqdd_ux
end
dqdd_ux = cross_qdd_ux(xv,uv);
# @btime cross_qdd_ux(xv,uv);
# println("typeof(dqdd_ux) = ", typeof(dqdd_ux))

function cross_qdd_xu(x,u)
  out1 = sparse(kronProd(permutedims(dFc_u + dFu_u), 1.0I(size(M_inv,1)))*dM_inv);
  out2 = sparse(kronProd(1.0I(length(u)), M_inv)*(dFc_xu + dFu_xu));
  dqdd_xu = sparse(out1 + out2);
  return dqdd_xu
end
dqdd_xu = cross_qdd_xu(xv,uv);
# @btime cross_qdd_xu(xv,uv);
# println("typeof(dqdd_xu) = ", typeof(dqdd_xu))

function hessqdd_u(x,u)
  dqdd_uu = sparse(kronProd(1.0I(length(u)), M_inv)*(dFc_uu + dFu_uu));
  return dqdd_uu
end
dqdd_uu = hessqdd_u(xv,uv);
# @btime hessqdd_u(xv,uv);
# println("typeof(dqdd_uu) = ", typeof(dqdd_uu))

function hessDyn(xuVar)
  x = xuVar[1:14*nB]; u = reshape(xuVar[14*nB+1:end], (6,nB));

  y = qdd;
  dy = dqdd;

  d2y_xx = dqdd_xx;#  hessqdd(x,u);
  d2y_xu = dqdd_xu;# cross_qdd_xu(x,u);
  d2y_ux = dqdd_ux;# cross_qdd_ux(x,u);
  d2y_uu = dqdd_uu;# hessqdd_u(x,u);

  d2y = [d2y_xx d2y_ux; d2y_xu d2y_uu];

  out1 = sparse(kronProd(1.0I(length(xuVar)), permutedims(y))*d2y);
  out2 = sparse(permutedims(dy)*dy);
  out = sparse(2*(out1 + out2)); # hessian of dynamics cost = sum(y.*y)
  return out;
end
dqdd_xuVar = hessDyn(xu_var);
# @btime hessDyn(xu_var);
# println("typeof(dqdd_xuVar) = ", typeof(dqdd_xuVar))

function hess_nonlinear_constraints()
  zLen = 71;

  hess_dyn = hessDyn(xu_var); # hessian of dynamics cost
  _, hess_qNorm = quatNormConstr_diff(xu_var); # hessian of quatNormConstr cost
  _, hess_jLoc = jointLocConstr_diff(xu_var); # hessian of jointLocConstr cost
  # println()
  # println("typeof(hess_dyn) = ", typeof(hess_dyn))
  # println("typeof(hess_qNorm) = ", typeof(hess_qNorm))
  # println("typeof(hess_jLoc) = ", typeof(hess_jLoc))

  hess_constr = hess_dyn +  hess_qNorm + hess_jLoc;

  # taking only free variables into consideration
  freeVarInd_hess = 14 .+ freeVarInd[1:end-1];
  hess_constr_free = hess_constr[freeVarInd_hess, freeVarInd_hess];

  out = zeros(typeof(hess_constr_free[1]),(zLen, zLen));
  out[1:zLen-1, 1:zLen - 1] = hess_constr_free;
  return out
end
hess_nc = hess_nonlinear_constraints();

#   return grad_nc, hess_nc;
# end

# grd, hss = gh_nc(xv, uv);
# @btime gh_nc(xv,uv);

# @btime hess_nonlinear_constraints();

## testing the trim function
# Ay1 = [1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0];

# Ay2 = [1.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 1.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 1.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 1.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 1.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 1.0];

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
import RecursiveArrayTools, DiffEqBase

# for compatibility with OrdinaryDiffEq
# Based on ApproxFun definitions
#  https://github.com/SciML/RecursiveArrayTools.jl/blob/6e779acb321560c75e27739a89ae553cd0f332f1/src/init.jl#L8-L12
# and
#  https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/16?u=simonbyrne

# This one is annoying
#  https://github.com/SciML/OrdinaryDiffEq.jl/blob/181dcf265351ed3c02437c89a8d2af3f6967fa85/src/initdt.jl#L80
Base.any(f, field::Field) = any(f, parent(field))

Base.muladd(x, field::Field) = muladd.(x, field)

Base.similar(field::F, ::Type{F}) where {F <: Field} = similar(field)

Base.vec(field::Field) = vec(parent(field))

RecursiveArrayTools.recursive_unitless_eltype(field::Field) = typeof(field)

RecursiveArrayTools.recursive_unitless_bottom_eltype(field::Field) =
    RecursiveArrayTools.recursive_unitless_bottom_eltype(parent(field))

RecursiveArrayTools.recursive_bottom_eltype(field::Field) =
    RecursiveArrayTools.recursive_bottom_eltype(parent(field))

function RecursiveArrayTools.recursivecopy!(dest::F, src::F) where {F <: Field}
    copy!(parent(dest), parent(src))
    return dest
end

RecursiveArrayTools.recursivecopy(field::Field) = copy(field)
# avoid call to deepcopy
RecursiveArrayTools.recursivecopy(a::AbstractArray{<:Field}) =
    map(RecursiveArrayTools.recursivecopy, a)

@inline RecursiveArrayTools.unpack(
    bc::Base.Broadcast.Broadcasted{
        RecursiveArrayTools.ArrayPartitionStyle{Style},
    },
    i,
) where {Style <: Fields.FieldStyle} = Base.Broadcast.broadcasted(
    bc.f,
    RecursiveArrayTools.unpack_args(i, bc.args)...,
)

function DiffEqBase.calculate_residuals!(
    out::Fields.Field,
    ũ::Fields.Field,
    u₀::Fields.Field,
    u₁::Fields.Field,
    α,
    ρ,
    internalnorm,
    t,
)
    DiffEqBase.calculate_residuals!(
        parent(out),
        parent(ũ),
        parent(u₀),
        parent(u₁),
        α,
        ρ,
        internalnorm,
        t,
    )
end

@inline function calculate_residuals!(out, ũ, u₀, u₁, α, ρ, internalnorm, t)
    @. out = calculate_residuals(ũ, u₀, u₁, α, ρ, internalnorm, t)
    nothing
end

# Play nice with DiffEq ArrayPartition
DiffEqBase.UNITLESS_ABS2(field::Field) =
    sum(DiffEqBase.UNITLESS_ABS2, parent(field))
DiffEqBase.recursive_length(field::Field) = 1
DiffEqBase.NAN_CHECK(field::Field) = DiffEqBase.NAN_CHECK(parent(field))


text:
module XCB

include("../deps/deps.jl")
include("include_h.jl")
include("include_manual.jl")

#this should work at least

# #include <unistd.h>      /* pause() */

# #include <xcb/xcb.h>

# int
# main ()
# {
#   xcb_connection_t *c;
#   xcb_screen_t     *screen;
#   xcb_window_t      win;

#   /* Open the connection to the X server */
#   c = xcb_connect (NULL, NULL);

#   /* Get the first screen */
#   screen = xcb_setup_roots_iterator (xcb_get_setup (c)).data;

#   /* Ask for our window's Id */
#   win = xcb_generate_id(c);

#   /* Create the window */
#   xcb_create_window (c,                             /* Connection          */
#                      XCB_COPY_FROM_PARENT,          /* depth (same as root)*/
#                      win,                           /* window Id           */
#                      screen->root,                  /* parent window       */
#                      0, 0,                          /* x, y                */
#                      150, 150,                      /* width, height       */
#                      10,                            /* border_width        */
#                      XCB_WINDOW_CLASS_INPUT_OUTPUT, /* class               */
#                      screen->root_visual,           /* visual              */
#                      0, NULL);                      /* masks, not used yet */

#   /* Map the window on the screen */
#   xcb_map_window (c, win);

#   /* Make sure commands are sent before we pause, so window is shown */
#   xcb_flush (c);

#   pause ();    /* hold client until Ctrl-C */

#   return 0;
# }

# ...so

function connect(displayname::AbstractString,screenp::Int32)
	c = ccall((:xcb_connect, _jl_libxcb), Ptr{Void},
                (Ptr{UInt8},Ptr{Int32}), bytestring(displayname), &screenp)
    xcb_connection_t(c)
end

function connect()
  c = ccall((:xcb_connect, _jl_libxcb), Ptr{Void},
                (Ptr{UInt8},Ptr{Int32}), C_NULL, C_NULL)
    xcb_connection_t(c)
end

function get_setup(c::xcb_connection_t)
	ccall((:xcb_get_setup, _jl_libxcb), Ptr{xcb_setup_t},
                (xcb_connection_t,), c)
end

function setup_roots_iterator(s::Ptr{xcb_setup_t})
	s1 = ccall((:xcb_setup_roots_iterator, _jl_libxcb), xcb_screen_iterator_t,
                (Ptr{xcb_setup_t},), s)
  xcb_screen_iterator_t(s1)
end

function generate_id(c::xcb_connection_t)
	ccall((:xcb_generate_id, _jl_libxcb), xcb_window_t,
                (xcb_connection_t,), c)
end

function map_window(c::xcb_connection_t,w::xcb_window_t)
  ccall((:xcb_map_window, _jl_libxcb), Void,
                (xcb_connection_t, xcb_window_t), c, w)
  nothing
end

function flush(c::xcb_connection_t)
  ccall((:xcb_flush, _jl_libxcb), Void,
                (xcb_connection_t,), c)
  nothing
end

# xcb_void_cookie_t
# xcb_create_window (xcb_connection_t *c  /**< */,
#                    uint8_t           depth  /**< */,
#                    xcb_window_t      wid  /**< */,
#                    xcb_window_t      parent  /**< */,
#                    int16_t           x  /**< */,
#                    int16_t           y  /**< */,
#                    uint16_t          width  /**< */,
#                    uint16_t          height  /**< */,
#                    uint16_t          border_width  /**< */,
#                    uint16_t          _class  /**< */,
#                    xcb_visualid_t    visual  /**< */,
#                    uint32_t          value_mask  /**< */,
#                    const uint32_t   *value_list  /**< */);

function create_window(c::xcb_connection_t,
                   depth::uint8_t,
                   wid::xcb_window_t,
                   parent::xcb_window_t,
                   x::int16_t,
                   y::int16_t,
                   width::uint16_t,
                   height::uint16_t,
                   border_width::uint16_t,
                   _class::uint32_t,
                   visual::xcb_visualid_t,
                   value_mask::uint32_t,
                   value_list::Ptr{Void})
	ccall((:xcb_create_window, _jl_libxcb), Void,
		(xcb_connection_t, UInt8, xcb_window_t, xcb_window_t, Int16, Int16,
      UInt16, UInt16, UInt16, UInt32,
			xcb_visualid_t, UInt32, Ptr{Void}), c,depth,wid,parent,x,y,width,height,border_width,
		_class,visual,value_mask,value_list)
	nothing
	end

function test0()

	c = connect()

    s0 = get_setup(c)

	s = setup_roots_iterator(s0)

	screen = unsafe_load(s.data)

	win = generate_id(c)

    a = UInt32(0)
    create_window(c,
        XCB_COPY_FROM_PARENT,win,screen.root,
        Int16(0),Int16(0),UInt16(150),UInt16(150),UInt16(10),
        XCB_WINDOW_CLASS_INPUT_OUTPUT,screen.root_visual,UInt32(0),pointer_from_objref(a))

    map_window(c,win)

    flush(c)
    c,screen,win
end

function screen_allowed_depths_iterator(screen::xcb_screen_t)
    ccall((:xcb_screen_allowed_depths_iterator, _jl_libxcb), xcb_depth_iterator_t,
                (xcb_screen_t,), screen)
end

function depth_next(d::xcb_depth_iterator_t)
    ccall((:xcb_depth_next, _jl_libxcb), Void,
                (Ptr{xcb_depth_iterator_t},), &d)
end

function screen_next(s::xcb_screen_iterator_t)
    ccall((:xcb_screen_next, _jl_libxcb), Void,
                (Ptr{xcb_screen_iterator_t},), &s)
end

function get_visual(c::xcb_connection_t, visual::xcb_visualid_t)

    screen_iter = setup_roots_iterator(get_setup(c))
    display(screen_iter)

    while screen_iter.rem > 0
        depth_iter = screen_allowed_depths_iterator(unsafe_load(screen_iter.data))

        while depth_iter.rem > 0
            display(depth_iter)

            depth_next(depth_iter)
        end
        display(screen_iter)
        screen_next(screen_iter)
    end

    0

end
# static xcb_visualtype_t *find_visual(xcb_connection_t *c, xcb_visualid_t visual)
# {
#   xcb_screen_iterator_t screen_iter = xcb_setup_roots_iterator(xcb_get_setup(c));

#   for (; screen_iter.rem; xcb_screen_next(&screen_iter)) {
#     xcb_depth_iterator_t depth_iter = xcb_screen_allowed_depths_iterator(screen_iter.data);
#     for (; depth_iter.rem; xcb_depth_next(&depth_iter)) {
#       xcb_visualtype_iterator_t visual_iter = xcb_depth_visuals_iterator(depth_iter.data);
#       for (; visual_iter.rem; xcb_visualtype_next(&visual_iter))
#         if (visual == visual_iter.data->visual_id)
#           return visual_iter.data;
#     }
#   }

#   return NULL;
# }

#  To get the xcb_visualtype_t structure, it's a bit less easy. You have to get the xcb_screen_t structure that you want, get its root_visual member, then iterate over the xcb_depth_ts and the xcb_visualtype_ts, and compare the xcb_visualid_t of these xcb_visualtype_ts: with root_visual:

# xcb_connection_t *c;
# xcb_screen_t     *screen;
# int               screen_nbr;
# xcb_visualid_t    root_visual = { 0 };
# xcb_visualtype_t  *visual_type = NULL;    /* the returned visual type */

# /* you init the connection and screen_nbr */

# screen = screen_of_display (c, screen_nbr);
# if (screen) {
#   xcb_depth_iterator_t depth_iter;

#   depth_iter = xcb_screen_allowed_depths_iterator (screen);
#   for (; depth_iter.rem; xcb_depth_next (&depth_iter)) {
#     xcb_visualtype_iterator_t visual_iter;

#     visual_iter = xcb_depth_visuals_iterator (depth_iter.data);
#     for (; visual_iter.rem; xcb_visualtype_next (&visual_iter)) {
#       if (screen->root_visual == visual_iter.data->visual_id) {
#         visual_type = visual_iter.data;
#         break;
#       }
#     }
#   }
# }

# /* visual_type contains now the visual structure, or a NULL visual structure if no screen is found */

end # module


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Libxc"
version = v"5.1.7"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://gitlab.com/libxc/libxc/-/archive/$(version)/libxc-$(version).tar.gz",
                  "2d82b7bcfd8749490f6bb0906acf99fbf03050696dd2213da4b7a7600fc14328"),
]

# Bash recipe for building across all platforms
# Notes:
#   - Autotools fully supported upstream, but Windows builds only work with CMake
#   - 3rd and 4th derivatives (KXC, LXC) not built since gives a binary size of ~200MB
script = raw"""
cd $WORKSPACE/srcdir/libxc-*/

if [[ "${target}" = *-mingw* ]]; then
    mkdir libxc_build
    cd libxc_build
    cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
        -DCMAKE_BUILD_TYPE=Release -DENABLE_FORTRAN=ON -DENABLE_XHOST=OFF -DBUILD_SHARED_LIBS=ON \
        -DDISABLE_VXC=OFF -DDISABLE_FXC=OFF -DDISABLE_KXC=ON -DDISABLE_LXC=ON ..
else
    autoreconf -vi
    export CFLAGS="$CFLAGS -std=c99"
    export FCFLAGS="-pipe -O3"
    ./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
        --disable-static --enable-shared \
        --enable-vxc=yes --enable-fxc=yes --enable-kxc=no --enable-lxc=no
fi

make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_gfortran_versions(supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct("libxc", :libxc)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"5", julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# # Statistical plots

using RDatasets: dataset
using SplitApplyPlot, CairoMakie
mpg = dataset("ggplot2", "mpg");
mpg.IsAudi = mpg.Manufacturer .== "audi"

data(mpg) *
    mapping(:Displ, :Hwy, col=:IsAudi => nonnumeric) *
    visual(QQPlot, qqline=:fit) |> draw
AbstractPlotting.save("qqplot.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](qqplot.svg)


text:
@testset "Forecast data" begin
    ps_dict = PowerSystems.parsestandardfiles(joinpath(MATPOWER_DIR, "case5_re.m"))
    da_time_series = PowerSystems.read_data_files(joinpath(FORECASTS_DIR, "5bus_ts");
                                                  REGEX_FILE=r"da_(.*?)\.csv")
    rt_time_series = PowerSystems.read_data_files(joinpath(FORECASTS_DIR, "5bus_ts");
                                                  REGEX_FILE=r"rt_(.*?)\.csv")
    ps_dict = PowerSystems.assign_ts_data(ps_dict,rt_time_series)

    buses, generators, storage, branches, loads, loadZones, shunts, services =
        PowerSystems.ps_dict2ps_struct(ps_dict)
    sys_5 = PowerSystems.System(buses, generators, loads, branches, storage,
                                     ps_dict["baseMVA"])

    forecast_gen = PowerSystems.make_forecast_dict(da_time_series["gen"],
                                                   Dates.Day(1), 24, generators)
    forecast_load = PowerSystems.make_forecast_dict(da_time_series, Dates.Day(1), 24, loads)
    forecast_gen = PowerSystems.make_forecast_dict(da_time_series["gen"], Dates.Day(1), 24,
                                                   generators)
    forecast_load = PowerSystems.make_forecast_dict(da_time_series, Dates.Day(1), 24, loads)
    forecast_gen = PowerSystems.make_forecast_dict(da_time_series["gen"], Dates.Day(1), 24,
                                                   generators);
    forecast_load = PowerSystems.make_forecast_dict(da_time_series, Dates.Day(1), 24,
                                                    loads);

    ps_dict = PowerSystems.parsestandardfiles(joinpath(MATPOWER_DIR, "RTS_GMLC.m"))
    da_time_series = PowerSystems.read_data_files(joinpath(FORECASTS_DIR,
                                                           "RTS_GMLC_forecasts");
                                                  REGEX_FILE=r"DAY_AHEAD(.*?)\.csv")
    rt_time_series = PowerSystems.read_data_files(joinpath(FORECASTS_DIR,
                                                           "RTS_GMLC_forecasts");
                                                  REGEX_FILE=r"REAL_TIME(.*?)\.csv")
    ps_dict = PowerSystems.assign_ts_data(ps_dict,rt_time_series)

    buses, generators, storage, branches, loads, loadZones, shunts, services =
        PowerSystems.ps_dict2ps_struct(ps_dict)
    sys_5 = PowerSystems.System(buses, generators, loads, branches, storage,
                                     ps_dict["baseMVA"])

    forecast_gen = PowerSystems.make_forecast_dict(da_time_series["gen"],
                                                   Dates.Day(1), 24, generators);
    forecast_load = PowerSystems.make_forecast_dict(rt_time_series, Dates.Day(1),
                                                    288, loads, loadZones)

    forecast_struct = PowerSystems.make_forecast_array(forecast_load)
end


text:
using ProgressMeter, Plots
include("util.jl")
const N = 50+2

import Base.zero
function zero(::Type{Char})::Char
    return '.'
end

function setup()::Array{Char,2}
    lines = getlines("inputs/18.txt")
    retval = zeros(Char,N,N)
    for lnum=1:length(lines)
        for cnum=1:length(lines[lnum])
            retval[lnum+1,cnum+1] = lines[lnum][cnum]
        end
    end
    return retval
end

function localstep(smallgrid::Array{Char,2})::Char
    if smallgrid[2,2] == '.'
        numtrees = 0
        for cell in smallgrid
            if cell == '|'
                numtrees += 1
            end
        end
        if numtrees >= 3
            return '|'
        else
            return '.'
        end
    elseif smallgrid[2,2] == '|'
        numlumberyards = 0
        for cell in smallgrid
            if cell == '#'
                numlumberyards += 1
            end
        end
        if numlumberyards >= 3
            return '#'
        else
            return '|'
        end
    elseif smallgrid[2,2] == '#'
        numlumberyards = 0
        numtrees = 0
        for cell in smallgrid
            if cell == '#'
                numlumberyards += 1
            elseif cell == '|'
                numtrees += 1
            end
        end
        if numlumberyards >= 2 && numtrees >= 1
            return '#'
        else
            return '.'
        end
    end
end

function timestep(grid::Array{Char,2})::Array{Char,2}
    retval = zeros(Char,N,N)
    @sync for col in 2:N-1
        @async for row in 2:N-1
            smallgrid = grid[row-1:row+1,col-1:col+1]
            retval[row,col] = localstep(smallgrid)
        end
    end
    return retval
end

function print_grid(grid)
    s = ""
    for row in 2:N-1
        s *= join(grid[row,2:N-1])
        s *= "\n"
    end
    println(s)
end

function run_10_mins(grid)
    for i in 1:10
        grid = timestep(grid)
        print_grid(grid)
    end
    return grid
end

function num_char(grid, char)::Int64
    c = 0
    for x in grid
        if x == char
            c += 1
        end
    end
    return c
end

function part_1(grid)::Int64
    numwood = num_char(grid, '|')
    numlumber = num_char(grid, '#')
    #println(numwood)
    #println(numlumber)
    return numwood*numlumber
end

function get_scores(grid)
    scores = []
    @showprogress for i in 1:10000
        grid = timestep(grid)
        push!(scores, part_1(grid))
    end
    return scores
end

grid = setup()

scores = get_scores(grid)
part_1(grid)

cycle = scores[5013:5040]
function get_cycle_data()
    data = []
    for d in 1:300
        data = vcat(data, cycle)
    end
    return data
end
data = get_cycle_data()
plot(scores[5000:5100])
plot!(data[5000:5100])

function get_value_at(x::Int64, cycle)
    _x = x % length(cycle)
    return cycle[_x]
end

function part_2(cycle)
    return get_value_at(1000000000, cycle)
end

println(part_2(cycle))