text:
# # Rate Constants, T-dependence
# Here, we will just briefly summarize how temperature affects enzymatic rate constants and how much uncertainty there actually is in literature (in terms of how to best define them and how plastic some of the variables are)

## Add usual tools we use:
#using Revise
using BenchmarkTools

# Use PyPlot to plot figures
using PyPlot

#----------------------------------------------------------------------------

## load Photosynthesis module:
using Land.Photosynthesis
#----------------------------------------------------------------------------

## Specify Field Type
const FT = Float32
#----------------------------------------------------------------------------

# ## T-dependence of V$_{c,max}$
# In literature, there are different implementations of the temperature dependence of enzymatic reaction rates. Common among all of them is a typical Arrhenius formulation with Activation energy $E_a$, which leads to a temperature dependence of a quantity $V$ given the standard value defined at 25°C as $T_{ref}$.
#
# $$V(T)= V(T_{ref})\, \underbrace{\exp\left(\frac{E_a}{RT_{ref}}-\frac{E_a}{RT}\right)}_{\text{Activation}}$$
#
# Other formulations add de-activation of proteins due to denaturalization at higher temperatures:
# $$V(T)= V(T_{ref})\, \exp\left(\frac{E_a}{RT_{ref}}-\frac{E_a}{RT}\right) \underbrace{\frac{1+\exp\left((\Delta ST_{ref}-H_d)/RT_{ref}\right)}{1+\exp\left((\Delta ST-H_d)/RT\right)}}_{\text{de-activation}} $$
# which includes an entropy term $\Delta S$ and the energy for deactivation $H_d$.
#
# To illustrate the differences, we show two different implementations, one only using the activation part based on Bernacchi et al 2001[^1] and one with the typical CLM5 implementation.
#
# [^1]: Bernacchi, C.J., Singsaas, E.L., Pimentel, C., Portis Jr, A.R. and Long, S.P., 2001. Improved temperature response functions for models of Rubisco‐limited photosynthesis. Plant, Cell & Environment, 24(2), pp.253-259.
#

## This looks a bit more tedious here than it needs to but in reality

vcmax_Bernacchi = Float32[]
vcmax_CLM = Float32[]

## Define T-range
Tleaf = collect(FT,260.0:1.0:315.0)

## Run through temperatures and save Vcmax values:
td_vc_bernacchi = Photosynthesis.VcmaxTDBernacchi(FT)
td_vc_clm       = Photosynthesis.VcmaxTDCLM(FT)
for T in Tleaf
    _Vcmax = Photosynthesis.photo_TD_from_val(td_vc_bernacchi, FT(100.0), T)
    push!(vcmax_Bernacchi, _Vcmax)
    _Vcmax = Photosynthesis.photo_TD_from_val(td_vc_clm, FT(100.0), T)
    push!(vcmax_CLM, _Vcmax)
end
#----------------------------------------------------------------------------

figure()
plot(Tleaf .- 273.15, vcmax_Bernacchi, label="Bernacchi", lw=2)
plot(Tleaf .- 273.15, vcmax_CLM,  label="CLM5", lw=2)
ylabel("Vcmax (µmol/m²/s)")
xlabel("Leaf Temperature (°C)")
legend()
gcf()
#----------------------------------------------------------------------------

# ### CO$_2$ compensation point $\Gamma_\star$

## Here, we only have one implementation:
Γ_CLM = Float32[]
Tleaf = collect(FT,260.0:1.0:315.0)
td_gamma_clm = Photosynthesis.ΓStarTDCLM(FT)

for T in Tleaf
    _ΓStar = Photosynthesis.photo_TD_from_set(td_gamma_clm, T)
    push!(Γ_CLM, _ΓStar)
end
#----------------------------------------------------------------------------
figure()
plot(Tleaf .- 273.15, Γ_CLM, label="Γ_CLM")
ylabel("Γ⋆ (Pa)")
xlabel("Leaf Temperature (°C)")
legend()
gcf()
#----------------------------------------------------------------------------

# ## T-dependence of J$_{max}$

Jmax_Bernacchi = Float32[]
Jmax_CLM = Float32[]
Tleaf = collect(FT,260.0:1.0:315.0)

td_j_bernacchi = Photosynthesis.JmaxTDBernacchi(FT)
td_j_clm = Photosynthesis.JmaxTDCLM(FT)

for T in Tleaf
    _Jmax = Photosynthesis.photo_TD_from_val(td_j_bernacchi, FT(100.0), T)
    push!(Jmax_Bernacchi, _Jmax)
    _Jmax = Photosynthesis.photo_TD_from_val(td_j_clm, FT(100.0), T)
    push!(Jmax_CLM, _Jmax)
end
#----------------------------------------------------------------------------
figure()
plot(Tleaf .- 273.15, Jmax_Bernacchi, label="Bernacchi")
plot(Tleaf .- 273.15, Jmax_CLM,  label="CLM5")
ylabel("Jmax (µmol/m2/s)")
xlabel("Leaf Temperature (°C)")
legend()
gcf()
#----------------------------------------------------------------------------

# ---
# ## Summary
# The way land surface models implement temperature variations in V$_{c,max}$ is still rather variable and there is no clear consensus in what model is best. In addition, gas exchange datasets performed at different temperatures might actually alias different confounding factors into the V$_{c,max}$ determination, so that the actually fitted V$_{c,max}$ strongly depends on the model formulation. E.g. mesophyll conductance is usually ignore in deriving V$_{c,max}$ but might be temperature dependent, which can cause T-dependent errors in V$_{c,max}$.
#
# That said, having the option to use different model forumations and actually using the ones that are implemented in land surface models directly when fitting leaf level measurements seems to be a prudent way to improve models in the future. More leaf level data might help in the future, especially looking at temperature dependencies, mesophyll conductance, and fluorescence yields. So many chemical and phsyical aspects change with temperature that it is hard to isolate impacting factors. Now, at least, we can run the models directly in the REPL, which is great for prototyping anc validating the model on the individual module level without any effort in terms of additional coding.

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
using Test
using TestSetExtensions
using Qaintessent
using StatsBase
using LinearAlgebra

@testset ExtendedTestSet "util pauli_vector" begin
    n1, n2, n3 = sample(1:15, 3, replace=true)
    @test all(Qaintessent.pauli_vector(n1, n2, n3) .≈ matrix(X) * n1 + matrix(Y) * n2 + matrix(Z) * n3)
end

@testset ExtendedTestSet "util intlog2" begin
    @test_throws ErrorException("Logarithm of 0 is undefined") Qaintessent.intlog2(0)
    @test_throws ErrorException("Logarithm of negative number is undefined") Qaintessent.intlog2(-30)
    N = rand(1:50)
    @test Qaintessent.intlog2(2^N) == N
    N = rand(1:200)
    @test Qaintessent.intlog2(N) <= log(2, N)
    @test Qaintessent.intlog2(N) >= (log(2, N) - 1)
end

@testset ExtendedTestSet "util binary_digits" begin
    n1 = abs(rand(Int, 1)[])
    # n1 = 5
    M = Qaintessent.intlog2(n1) + 1
    bd = Qaintessent.binary_digits(M, n1)
    x = 0
    count = 0

    for i in bd
        x+= i * 2^count
        count += 1
    end

    @test x ≈ n1
    @test Qaintessent.binary_to_int(bd) ≈ n1
end

@testset ExtendedTestSet "util binary_digits!" begin
    # n1 = abs(rand(Int, 1)[])
    n1 = 5
    M = Qaintessent.intlog2(n1) + 1
    m = BitArray{1}(undef, M)
    bd = Qaintessent.binary_digits!(m, n1)
    x = 0
    count = 0

    for i in bd
        x+= i * 2^count
        count += 1
    end

    @test x ≈ n1
end

@testset ExtendedTestSet "util quaternary_digits" begin
    n1 = abs(rand(Int, 1)[])
    M = Qaintessent.intlog2(n1)
    qd = Qaintessent.quaternary_digits(M, n1)
    x = 0
    count = 0
    for i in qd
        x+= i * 4^count
        count += 1
    end
    @test x ≈ n1
end

@testset ExtendedTestSet "util quaternary_digits!" begin
    n1 = abs(rand(Int, 1)[])
    M = Qaintessent.intlog2(n1)
    m = Vector{Int}(undef, M)
    qd = Qaintessent.quaternary_digits!(m, n1)
    x = 0
    count = 0
    for i in qd
        x+= i * 4^count
        count += 1
    end
    @test x ≈ n1
end

@testset ExtendedTestSet "util gramm_schmidt!" begin
    N = rand(4:10)

    @testset "util gramm_schmidt! complex" begin
        a = rand(ComplexF64, (N,N))

        Qaintessent.gramm_schmidt!(a)

        for i in 1:N
            for j in (i+1):N
                @test norm(dot(a[:, i], a[:, j])) < 1e-12
            end
            @test dot(a[:, i], a[:, i]) ≈ 1
        end
    end

    @testset "util gramm_schmidt! float" begin
        a = rand(Float64, (N,N))

        Qaintessent.gramm_schmidt!(a)

        for i in 1:N
            for j in (i+1):N
                @test norm(dot(a[:, i], a[:, j])) < 1e-12
            end
            @test dot(a[:, i], a[:, i]) ≈ 1
        end
    end
end

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

@testset ExtendedTestSet "util reduced density matrix" begin
    N = 4
    ψ = randn(ComplexF64, 2^N)
    χ = randn(ComplexF64, 2^N)
    # full density matrix |ψ⟩⟨χ|
    ρ = reshape(kron(conj(χ), ψ), 2^N, 2^N)

    id = Matrix{ComplexF64}(I, 2, 2)
    A = randn(ComplexF64, 2, 2)
    B = randn(ComplexF64, 2, 2)
    @testset ExtendedTestSet "reduced density matrix correctness" begin
        @test sum(kron(A, B) .* rdm(N, (4, 2), ψ, χ)) ≈ sum(kron(B, id, A, id) .* ρ)
    end

    @testset ExtendedTestSet "reduced density matrix exceptions" begin
        @test_throws ErrorException("Need at least one wire to act on.") rdm(N, (), ψ, χ)
        @test_throws ErrorException("Number of gate wires cannot be larger than total number of wires.") rdm(N, (1, 2, 3, 4, 5), ψ, χ)
        @test_throws ErrorException("Wire indices must be unique.") rdm(N, (2, 2), ψ, χ)
        @test_throws ErrorException("Wire index cannot be smaller than 1.") rdm(N, (-1, 2), ψ, χ)
        @test_throws ErrorException("Wire index cannot be larger than total number of wires.") rdm(N, (5, 1), ψ, χ)
    end
end

ext: jl


text:
using BAT
using Test
using Statistics
using StatsBase
using Distributions
using HypothesisTests

@testset "Gaussian Shell Distribution" begin
    #Tests different ways of instantiate GaussianShell Distribution
    @test @inferred(BAT.GaussianShell(r = 3., w = 0.1, n = 3)) isa BAT.GaussianShell
    @test @inferred(BAT.GaussianShell()) isa BAT.GaussianShell

    gs = BAT.GaussianShell(r = 3., w = 0.1, n = 3)

    @test @inferred(Base.length(gs)) == 3
    @test @inferred(Base.eltype(gs)) == Float64

    @test @inferred(StatsBase.params(gs)) == (3.0, 0.1, [0.0, 0.0, 0.0])

    @test isapprox(@inferred(Statistics.cov(gs)), [3.02718 0.00166016 0.0113455;
    0.00166016 3.01293 0.00453035; 0.0113455   0.00453035  3.01046], atol = 1e-1)

    @test isapprox(@inferred(Distributions._logpdf(gs, [0., 0., 0.])), -453.34571, atol = 1e-5)

    #In the limit when r is close to zero, the distribution should be a Gaussian
    gs = BAT.GaussianShell(r = 5e-10, w = 1., n = 1)
    ks_test = HypothesisTests.ExactOneSampleKSTest(rand(gs, 10^6)[:], Normal(0, 1.))
    @test pvalue(ks_test) > 0.05

end


text:
abstract type Variable end
abstract type Model end
const MaxOrder = 16

"""
mutable struct Configuration

    Struct that contains everything needed for MC.

    There are three different componenets:

 # Members

 ## Static parameters

 - `seed`: seed to initialize random numebr generator, also serves as the unique pid of the configuration

 - `rng`: a MersenneTwister random number generator, seeded by `seed`

 - `para`: user-defined parameter, set to nothing if not needed

 - `totalStep`: the total number of updates for this configuration

 - `var`: TUPLE of variables, each variable should be derived from the abstract type Variable, see variable.jl for details). Use a tuple rather than a vector improves the performance.

 ## integrand properties

- `neighbor::Vector{Vector{Int}}` : vectors that indicates the neighbors of each integrand. e.g., ([2, ], [1, ]) means the neighbor of the first integrand is the second one, while the neighbor of the second integrand is the first.
    There is a MC update proposes to jump from one integrand to another. If these two integrands' degrees of freedom are very different, then the update is unlikely to be accepted. To avoid this problem, one can specify neighbor to guide the update.

    By default, we assume the N integrands are in the increase order, meaning the neighbor will be set to ([N+1, 2], [1, 3], [2, 4], ..., [N-1,], [1, ]), where the first N entries are for diagram 1, 2, ..., N and the last entry is for the normalization diagram. Only the first diagram is connected to the normalization diagram.

 - `dof::Vector{Vector{Int}}`: degrees of freedom of each integrand, e.g., [[0, 1], [2, 3]] means the first integrand has zero var#1 and one var#2; while the second integrand has two var#1 and 3 var#2.

 - `observable`: observables that is required to calculate the integrands, will be used in the `measure` function call.
    It is either an array of any type with the common operations like +-*/^ defined.

 - `reweight`: reweight factors for each integrands. The reweight factor of the normalization diagram is assumed to be 1. Note that you don't need to explicitly add the normalization diagram.

 - `visited`: how many times this integrand is visited by the Markov chain.

 ## current MC state

 - `step`: the number of MC updates performed up to now

 - `curr`: the current integrand, initialize with 1

 - `norm`: the index of the normalization diagram. `norm` is larger than the index of any user-defined integrands

 - `normalization`: the accumulated normalization factor. Physical observable = Configuration.observable/Configuration.normalization.

 - `absWeight`: the abolute weight of the current integrand. User is responsible to initialize it after the contructor is called.

 - `propose/accept`: array to store the proposed and accepted updates for each integrands and variables.
    Their shapes are (number of updates X integrand number X max(integrand number, variable number).
    The last index will waste some memory, but the dimension is small anyway.
"""
mutable struct Configuration{V,P,O}
    ########### static parameters ###################
    seed::Int # seed to initialize random numebr generator, also serves as the unique pid of the configuration
    rng::MersenneTwister # random number generator seeded by seed
    para::P
    totalStep::Int64
    var::V

    ########### integrand properties ##############
    neighbor::Vector{Vector{Int}}
    dof::Vector{Vector{Int}} # degrees of freedom
    observable::O  # observables for each integrand
    reweight::Vector{Float64}
    visited::Vector{Float64}

    ############# current state ######################
    step::Int64
    curr::Int # index of current integrand
    norm::Int # index of the normalization diagram
    normalization::Float64 # normalization factor for observables
    absWeight::Float64 # the absweight of the current diagrams. Store it for fast updates

    propose::Array{Float64,3} # updates index, integrand index, integrand index
    accept::Array{Float64,3} # updates index, integrand index, integrand index

    """
    Configuration(totalStep, var::V, dof, obs::O; para::P=nothing, state=nothing, reweight=nothing, seed=nothing, neighbor=Vector{Vector{Int}}([])) where {V,P,O}

    Create a Configuration struct

 # Arguments

 ## Static parameters

 - `totalStep`: the total number MC steps of each block (one block, one configuration)

 - `var`: TUPLE of variables, each variable should be derived from the abstract type Variable, see variable.jl for details). Use a tuple rather than a vector improves the performance.

 - `dof::Vector{Vector{Int}}`: degrees of freedom of each integrand, e.g., [[0, 1], [2, 3]] means the first integrand has zero var#1 and one var#2; while the second integrand has two var#1 and 3 var#2.

 - `obs`: observables that is required to calculate the integrands, will be used in the `measure` function call
    It is either an array of any type with the common operations like +-*/^ defined.

 - `para`: user-defined parameter, set to nothing if not needed

 - `reweight`: reweight factors for each integrands. If not set, then all factors will be initialized with one.

 - `seed`: seed to initialize random numebr generator, also serves as the unique pid of the configuration. If it is nothing, then use RandomDevice() to generate a random seed in [1, 1000_1000]

- `neighbor::Vector{Vector{Int}}` : vectors that indicates the neighbors of each integrand. e.g., ([2, ], [1, ]) means the neighbor of the first integrand is the second one, while the neighbor of the second integrand is the first.
    There is a MC update proposes to jump from one integrand to another. If these two integrands' degrees of freedom are very different, then the update is unlikely to be accepted. To avoid this problem, one can specify neighbor to guide the update.

    By default, we assume the N integrands are in the increase order, meaning the neighbor will be set to ([N+1, 2], [1, 3], [2, 4], ..., [N-1,], [1, ]), where the first N entries are for diagram 1, 2, ..., N and the last entry is for the normalization diagram. Only the first diagram is connected to the normalization diagram.
"""
    function Configuration(totalStep, var::V, dof, obs::O; para::P=nothing, reweight=nothing, seed=nothing, neighbor=Vector{Vector{Int}}([])) where {V,P,O}
        @assert totalStep > 0 "Total step should be positive!"
        # @assert O <: AbstractArray "observable is expected to be an array. Noe get $(typeof(obs))."
        @assert V <: Tuple{Vararg{Variable}} || V <: Tuple{Variable} "Configuration.var must be a tuple of Variable to maximize efficiency. Now get $(typeof(V))"
        Nv = length(var) # number of variables

        ################# integrand initialization #########################
        @assert typeof(dof) == Vector{Vector{Int}} "Configuration.dof should be with a type of Vector{Vector{Int}} to avoid mistakes. Now get $(typeof(dof))"
        # add normalization diagram to dof
        dof = deepcopy(dof) # don't modify the input dof
        push!(dof, zeros(Int, length(var))) # add the degrees of freedom for the normalization diagram

        Nd = length(dof) # number of integrands + renormalization diagram
        @assert Nd > 1 "At least one integrand is required."
        # make sure dof has the correct size that matches var and neighbor
        for nv in dof
            @assert length(nv) == Nv "Each element of `dof` should have the same dimension as `var`"
        end

        if neighbor == []
            # By default, only the order-1 and order+1 diagrams are considered to be the neighbors
            # Nd is the normalization diagram, by default, it only connects to the first diagram
            neighbor = Vector{Vector{Int}}([[d - 1, d + 1] for d = 1:Nd])
            neighbor[1] = (Nd == 2 ? [2,] : [Nd, 2]) # if Nd=2, then 2 must be the normalization diagram
            neighbor[end] = [1,] # norm to the first diag
            (Nd >= 3) && (neighbor[end-1] = [Nd - 2,]) # last diag to the second last, possible only for Nd>=3
        end
        @assert typeof(neighbor) == Vector{Vector{Int}} "Configuration.neighbor should be with a type of Vector{Vector{Int}} to avoid mistakes. Now get $(typeof(neighbor))"
        @assert Nd == length(neighbor) "$Nd elements are expected for neighbor=$neighbor"

        ############# initialize reweight factors ########################
        if isnothing(reweight)
            reweight = [1.0 for d = 1:Nd] # the last element is for the normalization diagram
        else
            push!(reweight, 1.0)
        end
        @assert Nd == length(reweight) "reweight vector size is wrong! Note that the last element in reweight vector is for the normalization diagram."

        if isnothing(seed)
            seed = rand(Random.RandomDevice(), 1:1000000)
        end
        rng = MersenneTwister(seed)

        curr = 1 # set the current diagram to be the first one
        norm = Nd
        # a small initial absweight makes the initial configuaration quickly updated,
        # so that no error is caused even if the intial absweight is wrong,
        absweight = 1.0e-10
        normalization = 1.0e-10

        # visited[end] is for the normalization diagram
        visited = zeros(Float64, Nd) .+ 1.0e-8  # add a small initial value to avoid Inf when inverted

        # propose and accept shape: number of updates X integrand number X max(integrand number, variable number)
        # the last index will waste some memory, but the dimension is small anyway
        propose = zeros(Float64, (2, Nd, max(Nd, Nv))) .+ 1.0e-8 # add a small initial value to avoid Inf when inverted
        accept = zeros(Float64, (2, Nd, max(Nd, Nv)))

        return new{V,P,O}(seed, rng, para, totalStep, var,  # static parameters
            collect(neighbor), collect(dof), obs, collect(reweight), visited, # integrand properties
            0, curr, norm, normalization, absweight, propose, accept  # current MC state
        )
    end
end

function reset!(config, reweight=nothing)
    if typeof(config.observable) <: AbstractArray
        fill!(config.observable, zero(eltype(config.observable))) # reinialize observable
    else
        config.observable = zero(config.observable)
    end
    if isnothing(reweight) == false
        fill!(reweight, 1.0)
    end
    config.curr = 1
    config.normalization = 1.0e-10
    fill!(config.visited, 1.0e-8)
    fill!(config.propose, 1.0e-8)
    fill!(config.accept, 1.0e-10)
end

mutable struct FermiK{D} <: Variable
    # data::Vector{MVector{D,Float64}}
    data::Matrix{Float64}
    # data::Vector{Vector{Float64}}
    kF::Float64
    δk::Float64
    maxK::Float64
    offset::Int
    function FermiK(dim, kF, δk, maxK, size=MaxOrder; offset=0)
        @assert offset + 1 < size
        k = zeros(dim, size) .+ kF / sqrt(dim)
        # k0 = MVector{dim,Float64}([kF for i = 1:dim])
        # k0 = @SVector [kF for i = 1:dim]
        # k = [k0 for i = 1:size]
        return new{dim}(k, kF, δk, maxK, offset)
    end
end

Base.getindex(Var::FermiK{D}, i::Int) where {D} = Var.data[:, i]
function Base.setindex!(Var::FermiK{D}, v, i::Int) where {D}
    view(Var.data, :, i) .= v
end
Base.lastindex(Var::FermiK{D}) where {D} = size(Var.data)[2] # return index, not the value

mutable struct RadialFermiK <: Variable
    data::Vector{Float64}
    kF::Float64
    δk::Float64
    offset::Int
    function RadialFermiK(kF=1.0, δk=0.01, size=MaxOrder; offset=0)
        @assert offset + 1 < size
        k = [kF * (i - 0.5) / size for i = 1:size] #avoid duplication
        return new(k, kF, δk, offset)
    end
end

mutable struct BoseK{D} <: Variable
    data::Vector{SVector{D,Float64}}
    maxK::Float64
end

mutable struct Tau <: Variable
    data::Vector{Float64}
    λ::Float64
    β::Float64
    offset::Int
    function Tau(β=1.0, λ=0.5, size=MaxOrder; offset=0)
        @assert offset + 1 < size
        t = [β * (i - 0.5) / size for i = 1:size] #avoid duplication
        return new(t, λ, β, offset)
    end
end

mutable struct Continuous <: Variable
    data::Vector{Float64}
    λ::Float64
    lower::Float64
    range::Float64
    offset::Int
    function Continuous(bound, λ=nothing, size=MaxOrder; offset=0)
        lower, upper = bound
        @assert offset + 1 < size
        @assert upper > lower
        @assert isnothing(λ) || (0 < λ < (upper - lower))
        t = [lower + (upper - lower) * (i - 0.5) / size for i = 1:size] #avoid duplication

        if isnothing(λ)
            λ = (upper - lower) / 2.0
        end

        return new(t, λ, lower, upper - lower, offset)
    end
end

mutable struct Angle <: Variable
    data::Vector{Float64}
    λ::Float64
    offset::Int
    function Angle(λ=0.5, size=MaxOrder; offset=0)
        @assert offset + 1 < size
        theta = [π * (i - 0.5) / size for i = 1:size] #avoid dulication
        return new(theta, λ, offset)
    end
end

mutable struct TauPair <: Variable
    data::Vector{MVector{2,Float64}}
    λ::Float64
    β::Float64
    offset::Int
    function TauPair(β=1.0, λ=0.5, size=MaxOrder; offset=0)
        @assert offset + 1 < size
        t = [@MVector [β * (i - 0.4) / size, β * (i - 0.6) / size] for i = 1:size] #avoid duplication
        return new(t, λ, β, offset)
    end
end

mutable struct Discrete <: Variable
    data::Vector{Int}
    lower::Int
    upper::Int
    size::Int
    offset::Int
    function Discrete(lower, upper, size=MaxOrder; offset=0)
        d = [i for i = 1:size] #avoid dulication
        @assert offset + 1 < size
        @assert upper > lower
        return new(d, lower, upper, upper - lower + 1, offset)
    end
end

# mutable struct ContinuousND{D} <: Variable
#     data::Vector{Float64}
#     lower::Vector{Float64}
#     range::Vector{Float64}
#     offset::Int
#     function ContinuousND{dim}(lower, upper, size=MaxOrder; offset=0) where {dim}
#         if lower isa Number
#             lower = ones(Float64, dim) * lower
#         else
#             @assert length(lower) == dim && eltype(lower) isa Number
#         end
#         if upper isa Number
#             upper = ones(Float64, dim) * upper
#         else
#             @assert length(upper) == dim && eltype(upper) isa Number
#         end
#         @assert offset + 1 < size
#         @assert all(x -> x > 0, upper .- lower)
#         println(lower, ", ", upper)

#         ######## deterministic initialization #####################
#         t = []
#         for i in 1:size
#             for d in 1:dim
#                 # the same value should not appear twice!
#                 init = lower[d] + (upper[d] - lower[d]) * ((i - 1) * dim + d - 0.5) / (size * dim)
#                 @assert lower[d] <= init <= upper[d]
#                 append!(t, init)
#             end
#         end

#         return new{dim}(t, lower, upper .- lower, offset)
#     end
# end

Base.getindex(Var::Variable, i::Int) = Var.data[i]
function Base.setindex!(Var::Variable, v, i::Int)
    Var.data[i] = v
end
Base.firstindex(Var::Variable) = 1 # return index, not the value
Base.lastindex(Var::Variable) = length(Var.data) # return index, not the value

# struct Uniform{T,D} <: Model
#     lower::T
#     upper::T
#     function Uniform{T,D}(lower, upper) where {T<:Number,D}
#         return new{T,D}(lower, upper)
#     end
# end

# mutable struct Var{T,D,M} <: Variable
#     data::T
#     model::M
#     offset::Int
#     function Var(model{type,D}::Model, size=MaxOrder; offset=0) where {type<:Number,D}
#         # lower, upper = model.lower, model.upper
#         # k = zeros(type, dim, size + offset)
#         # for i in 1:size
#         #     for d in 1:dim
#         #         init = lower[d] + (upper[d] - lower[d]) * ((i - 1) * dim + d - 0.5) / (size * dim)
#         #         @assert lower[d] <= init <= upper[d]
#         #         k[d, i+offset] = init
#         #     end
#         # end
#         if D == 1
#             data = zeros(type, size)
#         else
#             data = zeros(type, (D, size))
#         end
#         return new{typeof(data),D,typeof(model)}(data, model, offset)
#     end
# end

# Base.getindex(var::Var{T,1,M}, i::Int) where {T,M} = var.data[i]
# function Base.setindex!(var::Var{T,1,M}, v, i::Int) where {T,M}
#     var.data[i] = v
# end
# Base.lastindex(var::Var{T,1,M}) where {T,M} = length(var.data) # return index, not the value

# Base.getindex(var::Var{T,D,M}, i::Int) where {T,D,M} = var.data[:, i]
# function Base.setindex!(var::Var{T,D,M}, v, i::Int) where {T,M}
#     var.data[:, i] = v
# end
# Base.lastindex(var::Var{T,D,M}) where {T,D,M} = size(var.data)[2] # return index, not the value


text:
"""
    GeneralizedGamma(a,d,p)

The *GeneralizedGamma distribution* with parameter a, d and p has probability density
function

```math
f(x; a, d, p) = \\frac{p}{a^p\\Gamma(d/p)} x^{d-1} e^{-\\left(\\frac{x}{a}\\right)^p},
\\quad x > 0
```

```julia
GeneralizedGamma()          # GeneralizedGamma distribution GeneralizedGamma(1, 1, 1)
GeneralizedGamma(a)         # GeneralizedGamma distribution GeneralizedGamma(a, 1, 1)

params(d)        # Get the parameters, i.e. (a, d, p)

S = rand(GeneralizedGamma(2, 3, 4), 5000); ##Sample from GeneralizedGamma(2, 3, 4)
fit_mle(GeneralizedGamma, S)  #MLE
```

External links

* [Generalized gamma distribution on Wikipedia](https://en.wikipedia.org/wiki/Generalized_gamma_distribution)
"""

struct GeneralizedGamma{T<:Real} <: ContinuousUnivariateDistribution
    a::T
    d::T
    p::T
    GeneralizedGamma{T}(a, d, p) where {T} = new{T}(a, d, p)
end

function GeneralizedGamma(a::T, d::T, p::T; check_args = true) where {T<:Real}
    check_args &&
        @check_args(GeneralizedGamma, a > zero(a) && d > zero(d) && p > zero(p))
    return GeneralizedGamma{T}(a, d, p)
end

GeneralizedGamma(a::Real, d::Real, p::Real) =
    GeneralizedGamma(promote(a, d, p)...)
GeneralizedGamma(a::Integer, d::Integer, p::Integer) =
    GeneralizedGamma(float(a), float(d), float(p))
GeneralizedGamma(a::T) where {T<:Real} = GeneralizedGamma(a, one(T), one(T))
GeneralizedGamma() = GeneralizedGamma(1.0, 1.0, 1.0, check_args = false)

@distr_support GeneralizedGamma 0.0 Inf

#### Conversions
convert(::Type{GeneralizedGamma{T}}, a::S, d::S, p::S) where {T<:Real,S<:Real} =
    GeneralizedGamma(T(a), T(d), T(p))
convert(
    ::Type{GeneralizedGamma{T}},
    d::GeneralizedGamma{S},
) where {T<:Real,S<:Real} =
    GeneralizedGamma(T(d.a), T(d.d), T(d.p), check_args = false)

#### Parameters

#shape(d::GeneralizedGamma) = d.p ##or d.d??
#scale(d::GeneralizedGamma) = d.a
#rate(d::GeneralizedGamma) = 1 / d.a

params(d::GeneralizedGamma) = (d.a, d.d, d.p)
partype(::GeneralizedGamma{T}) where {T} = T

#### Statistics

function mean(d::GeneralizedGamma)
    ip = 1 / d.p
    d.a * gamma((d.d + 1) * ip) / gamma(d.d * ip)
end
function var(d::GeneralizedGamma)
    ip = 1 / d.p
    iΓdp = 1 / gamma(d.d * ip)
    Γd1p = gamma((d.d + 1) * ip)
    Γd2p = gamma((d.d + 2) * ip)
    d.a^2 * (Γd2p * iΓdp - (Γd1p * iΓdp)^2)
end

function skewness(d::GeneralizedGamma)
    (a, d, p) = params(d)
    ip = 1 / p
    iΓdp = 1 / gamma(d * ip)
    Γd1p = gamma((d + 1) * ip)
    Γd2p = gamma((d + 2) * ip)
    Γd3p = gamma((d + 3) * ip)

    μ = a * Γd1p * iΓdp
    σ2 = a^2 * (Γd2p * iΓdp - (Γd1p * iΓdp)^2)
    EX3 = a^3 * Γd3p * iΓdp

    (EX3 - 3 * μ * σ2 - μ^3) / sqrt(σ2)^3
end

function kurtosis(d::GeneralizedGamma)
    (a, d, p) = params(d)
    ip = 1 / p
    iΓdp = 1 / gamma(d * ip)
    Γd1p = gamma((d + 1) * ip)
    Γd2p = gamma((d + 2) * ip)
    Γd3p = gamma((d + 3) * ip)
    Γd4p = gamma((d + 4) * ip)

    μ = a * Γd1p * iΓdp
    EX2 = a^2 * Γd2p * iΓdp
    EX3 = a^3 * Γd3p * iΓdp
    EX4 = a^4 * Γd4p * iΓdp
    (EX4 - 4 * EX3 * μ + 6 * EX2 * μ^2 - 3 * μ^4) / (EX2 - μ^2)^2 - 3
end

function mode(d::GeneralizedGamma)
    (a, d, p) = params(d)
    ip = 1 / p
    d >= 1 ? a * ((d - 1) * ip)^ip :
    error("GeneralizedGamma has no mode when d < 1")
end

function entropy(d::GeneralizedGamma)
    (a, d, p) = params(d)
    ip = 1 / p
    log(a * ip) + loggamma(d * ip) + d * ip - (d - 1) * ip * digamma(d * ip)
end

function logpdf(d::GeneralizedGamma{T}, x::Real) where {T<:Real}
    if x >= 0
        (a, d, p) = params(d)
        return log(p) - d * log(a) + (d - 1) * log(x) - (x / a)^p -
               loggamma(d / p)
    else
        return -T(Inf)
    end
end

pdf(d::GeneralizedGamma{T}, x::Real) where {T<:Real} = exp(logpdf(d, x))

# function cdf(d::GeneralizedGamma{T}, x::Real) where {T<:Real}
#     if x >= 0
#         (a, d, p) = params(d)
#         Γinc = gamma_inc(d / p, (x / a)^p, 0)
#         return Γinc[1] / (Γinc[1] + Γinc[2])
#     else
#         return zero(T)
#     end
# end
cdf(d::GeneralizedGamma{T}, x::Real) where {T<:Real} =
    cdf(Gamma(d.d / d.p, 1), (x / d.a)^d.p)
ccdf(d::GeneralizedGamma{T}, x::Real) where {T<:Real} =
    ccdf(Gamma(d.d / d.p, 1), (x / d.a)^d.p)

function quantile(d::GeneralizedGamma{T}, q::Real) where {T<:Real}
    (a, d, p) = params(d)
    ip = 1 / p
    a * quantile(Gamma(d * ip, 1), q)^ip
end

#function mgf(d::GeneralizedGamma{T}, t::Real) where {T<:Real}
#        g(x) = exp(logpdf(d, x) + t*x)
#        return quadgk(g, 0.0, Inf)
#end

# cf(d::GeneralizedGamma, t::Real)

function gradlogpdf(d::GeneralizedGamma{T}, x::Real) where {T<:Real}
    (a, d, p) = params(d)
    insupport(GeneralizedGamma, x) ? (-x^p * a^(-p) * p + d - 1) / x : zero(T)
end

#### Sampling
function rand(rng::AbstractRNG, d::GeneralizedGamma)
    (a, d, p) = params(d)
    if p == 1
        return rand(rng, Gamma(d, a))
    end
    ip = 1 / p
    return a * rand(rng, Gamma(d * ip, 1))^ip
end

#### Fit model
function fit_mle(
    ::Type{<:GeneralizedGamma},
    x::AbstractArray{T};
    maxiter::Int = 5000,
    tol::Float64 = 1e-6,
) where {T<:Real}
    n = length(x)
    logx = zeros(T, n)
    Σlogx = 0.0
    for i = 1:n
        logx[i] = log(x[i])
        Σlogx += logx[i]
    end
    ##Initial values
    pu, pl = 0.0, 0.0
    for i = 1.0:2.0:5000.0
        MM = gengamma_MLE_update_function(i, n, x, logx, Σlogx)
        if MM[4] < 0
            pu = i
            pl = i - 2

            break
        end
    end
    if gengamma_MLE_update_function(pu * 1.05, n, x, logx, Σlogx)[5]
        return error("K is zero.") # For some extreme conditions.
    end
    #print((pl,pu))
    if pu == 0
        return error("p is very large.")
    end
    pl = max(pl, 1.0e-7)
    # if isnan(gengamma_MLE_update_function(pl, n, x, logx, Σlogx)[4])
    #     return error("p is very small.")
    # end
    #if isnan(gengamma_MLE_update_function(pu, n, x, logx, Σlogx)[4])
    #    return error("p is very large.")
    #end
    a = 0.0
    d = 0.0
    p = 0.0
    p0 = (pl + pu) * 0.5
    for i = 1:maxiter
        #println(p0)
        if p0 < pl || p0 > pu || isnan(p0)
            p0 = pl + (pu - pl) * rand()
        end
        a, d, p = gengamma_MLE_update_function(p0, n, x, logx, Σlogx)
        if abs(p0 - p) <= tol
            return GeneralizedGamma(a, d, p)
        end
        p0 = p
    end
    @warn("Maximum iteration reached.")
    return GeneralizedGamma(a, d, p)
end

function gengamma_MLE_update_function(p, n, x, logx, Σlogx)
    Σxp = 0.0
    ∂Σxp = 0.0
    ∂2Σxp = 0.0
    ∂3Σxp = 0.0
    ip = 1 / p
    for i = 1:n
        xp = x[i]^p
        log_x = logx[i]
        Σxp += xp
        ∂Σxp += xp * log_x
        ∂2Σxp += xp * log_x^2
        ∂3Σxp += xp * log_x^3
    end
    K1 = n * ∂Σxp - Σlogx * Σxp
    K = 1 / K1
    ap = n^2 * ip * K
    a = ap^(-ip) ##exp(-ip*log(ap))
    log_ap = log(ap)
    log_a = -ip * log_ap
    ia = 1 / a
    d = n * Σxp * K
    dip = d * ip
    pΣxp = p * Σxp
    loggamma_dip = loggamma(d * ip)
    #ld = -n * log(a) + Σlogx - n * digamma_dip * ip
    ##Derivatives
    ∂ip = -ip^2
    ∂K1 = n * ∂2Σxp - Σlogx * ∂Σxp
    ∂K = -∂K1 * K^2
    ∂ap = n^2 * (∂ip * K + ip * ∂K)
    ∂log_ap = ∂ap / ap

    a_temp = -(∂ip * log_ap + ip * ∂log_ap)
    ∂a = a_temp * a
    ∂log_a = ∂a * ia
    #∂ia = (∂ip * log(ap) + ip * ∂ap / ap) * ia
    ∂d = n * (∂Σxp * K + Σxp * ∂K)
    #∂pΣxp = Σxp + p * ∂Σxp
    ∂dip = (∂d * ip + ∂ip * d)
    digamma_dip = digamma(dip)

    ∂loggamma_dip = ∂dip * digamma_dip
    #∂digamma_dip = trigamma(d * ip) * (∂d * ip + ∂ip * d)
    #∂ld = -n * ∂a * ia - n * (∂digamma_dip * ip + digamma_dip * ∂ip)
    ##Second derivatives
    ∂2ip = -2 * ∂ip * ip
    ∂2K1 = n * ∂3Σxp - Σlogx * ∂2Σxp
    ∂2K = -∂2K1 * K^2 - 2 * K * ∂K * ∂K1
    ∂2ap = n^2 * ((∂2ip * K + ∂ip * ∂K) + (∂ip * ∂K + ip * ∂2K))
    ∂2log_ap = (∂2ap * ap - ∂ap^2) / ap^2
    ∂a_temp =
        -((∂2ip * log_ap + ∂ip * ∂log_ap) + (∂ip * ∂log_ap + ip * ∂2log_ap))
    ∂2a = ∂a_temp * a + a_temp * ∂a
    ∂2log_a = (∂2a * a - ∂a^2) * ia^2

    ∂2d = n * ((∂2Σxp * K + ∂Σxp * ∂K) + (∂Σxp * ∂K + Σxp * ∂2K))
    ∂digamma_dip = ∂dip * trigamma(dip)
    ∂2dip = (∂2d * ip + ∂d * ∂ip + ∂2ip * d + ∂ip * ∂d)
    ∂2loggamma_dip = ∂2dip * digamma_dip + ∂dip * ∂digamma_dip
    ### log likelihood as a function of p. (a and d have closed form depend on p)
    llike =
        n * log(p) - n * d * log_a + (d - 1) * Σlogx - Σxp * ap -
        n * loggamma_dip

    ∂llike =
        n * ip - n * (∂d * log_a + d * ∂log_a) + ∂d * Σlogx -
        (∂Σxp * ap + ∂ap * Σxp) - n * ∂loggamma_dip

    ∂2llike =
        n * ∂ip - n * (∂2d * log_a + 2 * ∂d * ∂log_a + d * ∂2log_a) +
        ∂2d * Σlogx - (∂2Σxp * ap + 2 * ∂Σxp * ∂ap + Σxp * ∂2ap) -
        n * ∂2loggamma_dip

    return a, d, p - ∂llike / ∂2llike, ∂llike, ∂K == 0, llike
end

# Functions h(α) and h'(α). We will solve h(α) = 0 using Newton's method.
# α = p
# function gengamma_h_dh(
#     α::Float64,
#     n::Int64,
#     X::AbstractArray{T},
#     logX::AbstractArray{T},
#     slogX::Float64,
# ) where {T<:Real}
#     iα1 = 1 / (α + 1)
#     Sx = 0.0
#     Sxlog = 0.0
#     dSxlog = 0.0
#     for i = 1:n
#         xα = X[i]^α
#         Sx += xα
#         Sxlog += xα * logX[i]
#         dSxlog += xα * logX[i]^2
#     end
#     k = (n - 1.5 - α * iα1) / α - log(α) * α * iα1^2
#     dk = (α - 1) * log(α) * iα1^3 - (n - 1.5) / α^2
#     dSx = Sxlog
#     A = k * Sx
#     B = n * Sxlog - Sx * slogX
#     ϕ = A / B
#     dA = dk * Sx + k * dSx
#     dB = n * dSxlog - dSx * slogX
#     dϕ = (dA * B - A * dB) / B^2
#
#     H = n * digamma(n * ϕ) + α * slogX - 1 / ϕ - n * digamma(ϕ) - n * log(Sx)
#
#     dH =
#         n^2 * dϕ * trigamma(n * ϕ) + slogX + dϕ / ϕ^2 - n * dϕ * trigamma(ϕ) -
#         n * dSx / Sx
#     return H, dH, ϕ, Sx
# end

# # function sumlogpdf(x, n, slogx, a, d, p)
# #     a1 = 1 / a
# #     n * log(p) + (d - 1) * slogx - sum((x * a1) .^ p) + n * d * log(a1) -
# #     n * log(gamma(d / p))
# # end


text:
"""
    schottky_genus_4(τ)

Compute the value of the genus 4 Schottky-Igusa polynomial, given as input a 4x4 matrix in the Siegel upper half space.
"""
function schottky_genus_4(τ::Array{<:Number})
    R = RiemannMatrix(τ, siegel=false, ϵ=1.0e-10, nderivs=0);
    m1 = [[1,0,1,0], [1,0,1,0]];
    m2 = [[0,0,0,1], [1,0,0,0]];
    m3 = [[0,0,1,1], [1,0,1,1]];
    n1 = [[0,0,0,1], [1,1,1,0]];
    n2 = [[0,0,1,1], [0,0,0,1]];
    n3 = [[0,0,1,0], [1,0,1,1]];
    p1 = 1;
    p2 = 1;
    p3 = 1;
    z = [0; 0; 0; 0];
    for i1=0:1, i2=0:1, i3=0:1
        n = i1*n1 + i2*n2 + i3*n3;
        c1 = theta(z, R, char=remainder_char(m1+n));
        c2 = theta(z, R, char=remainder_char(m2+n));
        c3 = theta(z, R, char=remainder_char(m3+n));
        p1 *= c1;
        p2 *= c2;
        p3 *= c3;
    end
    schottky_poly = p1^2 + p2^2 + p3^2 - 2*p1*p2 - 2*p1*p3 - 2*p2*p3;
    return schottky_poly
end

"""
    random_nonschottky_genus_4(tol=0.1, trials=100)

Find a random 4x4 matrix in the Siegel upper half space which is not in the Schottky locus, up to the input tolerance and number of trials.
"""
function random_nonschottky_genus_4(tol::Real=0.1, trials::Integer=100)
    t = 0; # largest value of schottky polynomial
    i = 0; # counter for number of trials
    max_matrix = rand(4,4); # stores best non-schottky candidate
    while t < tol && i < trials
        τ = random_siegel(4);
        s = schottky_genus_4(τ);
        if abs(s) > t
            t = abs(s)
            max_matrix = τ;
        end
        i += 1;
    end
    return [max_matrix, t];
end


text:
using Profile, Plots

include("simulate.jl")
include("serverbuffer.jl")
include("projectscenarios.jl")

function do_experiments(;max_time = 10.0^5, λ = [0.2, 0.4, 0.6, 0.8, 1, 1.5, 2, 3], params_all = [scenario1, scenario2, scenario3, scenario4, scenario5])

    # 10^5 because it'll take too long to do 10^7. Each sim with the latter time runs for ~30 secs, which I'd say is acceptable, but not acceptable enough for me to test it half an hour before deadline
    mean_sum = 0
    orbit_sum = 0
    last_time = .0

    # Adds sum of jobs in system to running total
    function record_total_sum(time::Float64, state::DiscreteState)
        mean_sum += state.jobs_in_system*(time-last_time)
    end

    # Adds sum of orbiting jobs to running total
    function record_total_orbit(time::Float64, state::DiscreteState)
        orbit_sum += state.moving_jobs_num*(time-last_time)
    end

    # Wraps all our recordings
    function record_all(time::Float64, state::DiscreteState)
        record_total_sum(time, state)
        record_total_orbit(time, state)
        last_time = time
    end

    current_scenario = 1

    #Iterate through all our records, give back a chart
    for params in params_all
        mean_vec = []
        orbit_vec = []
        for i in λ
            mean_sum = 0
            orbit_sum = 0
            last_time = .0
            new_params = NetworkParameters(params.L, params.gamma_scv, i, params.η, params.μ_vector, params.P, params.Q, params.p_e, params.K)
            state = DiscreteState(new_params)
            @time simulate(state, generate_timed_event(.0, state, ArrivalEvent(state.last_job)), max_time=max_time, callback=record_all)
            push!(mean_vec, mean_sum/max_time)
            push!(orbit_vec, orbit_sum/mean_sum)
        end
        print(mean_vec)
        print(orbit_vec)
        mean_plot = plot(λ, mean_vec, title="Scenario $current_scenario: Mean Jobs in System", xlabel="λ", ylabel="Total Items",legend=false)
        orbit_plot = plot(λ, mean_vec, title="Scenario $current_scenario: Mean Orbiting Jobs in System", xlabel="λ", ylabel="Orbiting Items",legend=false)
        display(plot(mean_plot, orbit_plot, layout = (2, 1),legend=false, size = (600, 900)))
        current_scenario += 1
    end
end

do_experiments()


text:
module NormalCopula
    using Distributions
    using LinearAlgebra

    export normcopula_density, normcopula_logdensity

    const stdnormal = Normal()
    function normcopula_density(u, R)
        pvec = similar(u)
        for i = 1:length(u)
            pvec[i] = quantile(stdnormal, u[i])
        end
        RmI = inv(R)
        for i = 1:length(u)
            RmI[i,i] -= 1.
        end
        return 1/sqrt(det(R)) * exp(-0.5 * dot(pvec, RmI*pvec))
    end

    function normcopula_logdensity(u, R)
        pvec = similar(u)
        for i = 1:length(u)
            pvec[i] = quantile(stdnormal, u[i])
        end
        RmI = inv(R)
        for i = 1:length(u)
            RmI[i,i] -= 1.
        end
        return -log(sqrt(det(R))) + -0.5 * dot(pvec, RmI*pvec)
    end

end # module

max_stars_repo_name: chrished/NormalCopula


text:
using FocusedBlindDecon
using Test
using LinearAlgebra
using Misfits

include("adj_tests.jl")
include("fbd_tests.jl")


text:
using SparseDiffTools, CuArrays, Test, LinearAlgebra
using ArrayInterface: allowed_getindex, allowed_setindex!
function f(dx,x)
    dx[2:end-1] = x[1:end-2] - 2x[2:end-1] + x[3:end]
    allowed_setindex!(dx,-2allowed_getindex(x,1) + allowed_getindex(x,2),1)
    allowed_setindex!(dx,-2allowed_getindex(x,30) + allowed_getindex(x,29),30)
    nothing
end

_J1 = similar(rand(30,30))
_denseJ1 = cu(collect(_J1))
x = cu(rand(30))
CuArrays.allowscalar(false)
forwarddiff_color_jacobian!(_denseJ1, f, x)
forwarddiff_color_jacobian!(_denseJ1, f, x, sparsity = _J1)
forwarddiff_color_jacobian!(_denseJ1, f, x, colorvec = repeat(1:3,10), sparsity = _J1)
_Jt = similar(Tridiagonal(_J1))
@test_broken forwarddiff_color_jacobian!(_denseJ1, f, x, colorvec = repeat(1:3,10), sparsity = _Jt)


text:
import DataFrames
import SQLite

Base.eval(:(have_color = true))

function gameFileSetup(gameFile, typeOfGame, cheating, time_limit, time_add)

    db= SQLite.DB(gameFile)

    SQLite.execute!(db,"create table meta(key,value);")

    if typeOfGame == "S"
        SQLite.execute!(db,"insert into meta values ('type','standard');")
    elseif typeOfGame == "M"
        SQLite.execute!(db,"insert into meta values ('type','mini');")
    elseif typeOfGame == "C"
        SQLite.execute!(db,"insert into meta values ('type','chu');")
    elseif typeOfGame == "T"
        SQLite.execute!(db,"insert into meta values ('type','ten');")
    end

    if cheating == "T"
        SQLite.execute!(db,"insert into meta values ('legality','cheating');")
    else
        SQLite.execute!(db,"insert into meta values ('legality','legal');")
    end

    if time_limit == 0
        SQLite.execute!(db,"insert into meta values ('timed','no');")
        SQLite.execute!(db,"insert into meta (key) values ('time_add');")
        SQLite.execute!(db,"insert into meta (key) values ('sente_time');")
        SQLite.execute!(db,"insert into meta (key) values ('gote_time');")
    else
        SQLite.execute!(db,"insert into meta values ('timed','yes');")
        SQLite.execute!(db,"insert into meta values ('time_add','$time_add');")
        SQLite.execute!(db,"insert into meta values ('sente_time','$time_limit');")
        SQLite.execute!(db,"insert into meta values ('gote_time','$time_limit');")
    end

    SQLite.execute!(db,"insert into meta values ('seed','$(round(Int64,time()*1000))');")

    SQLite.execute!(db,"create table moves (move_number integer primary key ,move_type ,sourcex,sourcey,targetx,targety,option,i_am_cheating,targetx2,targety2,targetx3,targety3);")

    return db

end

function gameSetup(gameFile)
    # metaTable is a DataFrame holding SQLite meta table data
    metaTable = SQLite.query(gameFile, "select * from meta")

    gameType = get(metaTable[2][1]) # standard, mini, chu, or ten shogi
    if gameType == "ten"
        gameType = "tenjiku"
    end
    is_legal = get(metaTable[2][2]) == "yes"? true:false
    is_timed = get(metaTable[2][3]) == "yes"? true:false

    seed = 0
    try
        seed = parse(Int64, get(metaTable[2][7]))
    catch e
        seed = parse(Float64, get(metaTable[2][7])) # seed - unix time since epoch
        seed_str = string(seed)
        seed = parse(Int64, string(seed_str[1]) * seed_str[3:length(seed_str)-2])
    end

    # movesTable is a DataFrame holding SQLite moves table data
    movesTable = SQLite.query(gameFile, "select * from moves")

    # the number of moves that have been played so far in the game
    movesPlayed = length(movesTable[1])
    return metaTable, movesTable, gameType, seed, is_legal, is_timed
end

function getTimeData(metaTable)
    time_add = parse(Float64, get(metaTable[2][4]))
    sente_time = parse(Float64, get(metaTable[2][5]))
    gote_time = parse(Float64, get(metaTable[2][6]))
    return time_add, sente_time, gote_time
end

function boardSetup(gameType)
    # board size: 16x16 for tenjiku, 12x12 for chu, 9x9 for standard, 5x5 for mini
    if gameType == "tenjiku"
        boardLength = 16
    elseif gameType == "chu"
        boardLength = 12
    elseif gameType == "standard"
        boardLength = 9
    else
        boardLength = 5
    end

    # promotion ranks: 5 for tenjiku, 4 for chu, 3 for standard, 1 for mini
    if gameType == "tenjiku"
        numOfPromotionRanks = 5
    elseif gameType == "chu"
        numOfPromotionRanks = 4
    elseif gameType == "standard"
        numOfPromotionRanks = 3
    else
        numOfPromotionRanks = 1
    end

    # populate a board of empty pieces
    # syntax is board[y][x] : 1st index is y coordinate, 2nd index is x coordinate
    # board[1][1] is the top left corner
    board = [ [ empty for i = 1:boardLength ] for j = 1:boardLength ]

    if gameType == "tenjiku"

        # sente/black/odd pieces
        board[1][16] = piece("lance", "lance", 16, 1, 1, false)
        board[16][16] = piece("lance", "lance", 16, 16, 1, false)
        board[2][16] = piece("knight", "knight", 16, 2, 1, false)
        board[15][16] = piece("knight", "knight", 16, 15, 1, false)
        board[3][16] = piece("leopard", "leopard", 16, 3, 1, false) # ferocious leopard
        board[14][16] = piece("leopard", "leopard", 16, 14, 1, false) # ferocious leopard
        board[4][16] = piece("iron", "iron", 16, 4, 1, false) # iron general
        board[13][16] = piece("iron", "iron", 16, 13, 1, false) # iron general
        board[5][16] = piece("copper", "copper", 16, 5, 1, false) # copper general
        board[12][16] = piece("copper", "copper", 16, 12, 1, false) # copper general
        board[6][16] = piece("silver", "silver", 16, 6, 1, false) # silver general
        board[11][16] = piece("silver", "silver", 16, 11, 1, false) # silver general
        board[7][16] = piece("gold", "gold", 16, 7, 1, false) # gold general
        board[10][16] = piece("gold", "gold", 16, 10, 1, false) # gold general
        board[8][16] = piece("elephant", "elephant", 16, 8, 1, false) # drunk elephant
        board[9][16] = piece("king", "king", 16, 9, 1, false)
        board[1][15] = piece("chariot", "chariot", 15, 1, 1, false) # reverse chariot
        board[16][15] = piece("chariot", "chariot", 15, 16, 1, false) # reverse chariot
        board[3][15] = piece("csoldier", "csoldier", 15, 3, 1, false) # chariot soldier
        board[4][15] = piece("csoldier", "csoldier", 15, 4, 1, false) # chariot soldier
        board[13][15] = piece("csoldier", "csoldier", 15, 13, 1, false) # chariot soldier
        board[14][15] = piece("csoldier", "csoldier", 15, 14, 1, false) # chariot soldier
        board[6][15] = piece("tiger", "tiger", 15, 6, 1, false) # blind tiger
        board[11][15] = piece("tiger", "tiger", 15, 11, 1, false) # blind tiger
        board[7][15] = piece("phoenix", "phoenix", 15, 7, 1, false)
        board[8][15] = piece("queen", "queen", 15, 8, 1, false)
        board[9][15] = piece("lion", "lion", 15, 9, 1, false)
        board[10][15] = piece("kirin", "kirin", 15, 10, 1, false)
        board[1][14] = piece("ssoldier", "ssoldier", 14, 1, 1, false) # side soldier
        board[16][14] = piece("ssoldier", "ssoldier", 14, 16, 1, false) # side soldier
        board[2][14] = piece("vsoldier", "vsoldier", 14, 2, 1, false) # vertical soldier
        board[15][14] = piece("vsoldier", "vsoldier", 14, 15, 1, false) # vertical soldier
        board[3][14] = piece("bishop", "bishop", 14, 3, 1, false)
        board[14][14] = piece("bishop", "bishop", 14, 14, 1, false)
        board[4][14] = piece("horse", "horse", 14, 4, 1, false) # dragon horse
        board[13][14] = piece("horse", "horse", 14, 13, 1, false) # dragon horse
        board[5][14] = piece("dragon", "dragon", 14, 5, 1, false) # dragon king
        board[12][14] = piece("dragon", "dragon", 14, 12, 1, false) # dragon king
        board[6][14] = piece("buffalo", "buffalo", 14, 6, 1, false) # water buffalo
        board[11][14] = piece("buffalo", "buffalo", 14, 11, 1, false) # water buffalo
        board[7][14] = piece("demon", "demon", 14, 7, 1, false) # fire demon
        board[10][14] = piece("demon", "demon", 14, 10, 1, false) # fire demon
        board[8][14] = piece("eagle", "eagle", 14, 8, 1, false) # free eagle
        board[9][14] = piece("hawk", "hawk", 14, 9, 1, false) # lion hawk
        board[1][13] = piece("smover", "smover", 13, 1, 1, false) # side mover
        board[16][13] = piece("smover", "smover", 13, 16, 1, false) # side mover
        board[2][13] = piece("vmover", "vmover", 13, 2, 1, false) # vertical mover
        board[15][13] = piece("vmover", "vmover", 13, 15, 1, false) # vertical mover
        board[3][13] = piece("rook", "rook", 13, 3, 1, false)
        board[14][13] = piece("rook", "rook", 13, 14, 1, false)
        board[4][13] = piece("falcon", "falcon", 13, 4, 1, false) # horned falcon
        board[13][13] = piece("falcon", "falcon", 13, 13, 1, false) # horned falcon
        board[5][13] = piece("soaring", "soaring", 13, 5, 1, false) # soaring eagle
        board[12][13] = piece("soaring", "soaring", 13, 12, 1, false) # soaring eagle
        board[6][13] = piece("bgeneral", "bgeneral", 13, 6, 1, false) # bishop general
        board[11][13] = piece("bgeneral", "bgeneral", 13, 11, 1, false) # bishop general
        board[7][13] = piece("rgeneral", "rgeneral", 13, 7, 1, false) # rook general
        board[10][13] = piece("rgeneral", "rgeneral", 13, 10, 1, false) # rook general
        board[8][13] = piece("vice", "vice", 13, 8, 1, false) # vice general
        board[9][13] = piece("great", "great", 13, 9, 1, false) # great general
        board[1][12] = piece("pawn", "pawn", 12, 1, 1, false)
        board[2][12] = piece("pawn", "pawn", 12, 2, 1, false)
        board[3][12] = piece("pawn", "pawn", 12, 3, 1, false)
        board[4][12] = piece("pawn", "pawn", 12, 4, 1, false)
        board[5][12] = piece("pawn", "pawn", 12, 5, 1, false)
        board[6][12] = piece("pawn", "pawn", 12, 6, 1, false)
        board[7][12] = piece("pawn", "pawn", 12, 7, 1, false)
        board[8][12] = piece("pawn", "pawn", 12, 8, 1, false)
        board[9][12] = piece("pawn", "pawn", 12, 9, 1, false)
        board[10][12] = piece("pawn", "pawn", 12, 10, 1, false)
        board[11][12] = piece("pawn", "pawn", 12, 11, 1, false)
        board[12][12] = piece("pawn", "pawn", 12, 12, 1, false)
        board[13][12] = piece("pawn", "pawn", 12, 13, 1, false)
        board[14][12] = piece("pawn", "pawn", 12, 14, 1, false)
        board[15][12] = piece("pawn", "pawn", 12, 15, 1, false)
        board[16][12] = piece("pawn", "pawn", 12, 16, 1, false)
        board[5][11] = piece("dog", "dog", 11, 5, 1, false)
        board[12][11] = piece("dog", "dog", 11, 12, 1, false)

        # gote/white/even pieces
        board[1][1] = piece("lance", "lance", 1, 1, 0, false)
        board[16][1] = piece("lance", "lance", 1, 16, 0, false)
        board[2][1] = piece("knight", "knight", 1, 2, 0, false)
        board[15][1] = piece("knight", "knight", 1, 15, 0, false)
        board[3][1] = piece("leopard", "leopard", 1, 3, 0, false) # ferocious leopard
        board[14][1] = piece("leopard", "leopard", 1, 14, 0, false) # ferocious leopard
        board[4][1] = piece("iron", "iron", 1, 4, 0, false) # iron general
        board[13][1] = piece("iron", "iron", 1, 13, 0, false) # iron general
        board[5][1] = piece("copper", "copper", 1, 5, 0, false) # copper general
        board[12][1] = piece("copper", "copper", 1, 12, 0, false) # copper general
        board[6][1] = piece("silver", "silver", 1, 6, 0, false) # silver general
        board[11][1] = piece("silver", "silver", 1, 11, 0, false) # silver general
        board[7][1] = piece("gold", "gold", 1, 7, 0, false) # gold general
        board[10][1] = piece("gold", "gold", 1, 10, 0, false) # gold general
        board[8][1] = piece("king", "king", 1, 8, 0, false) # king
        board[9][1] = piece("elephant", "elephant", 1, 9, 0, false) # drunk elephant
        board[1][2] = piece("chariot", "chariot", 2, 1, 0, false) # reverse chariot
        board[16][2] = piece("chariot", "chariot", 2, 16, 0, false) # reverse chariot
        board[3][2] = piece("csoldier", "csoldier", 2, 3, 0, false) # chariot soldier
        board[4][2] = piece("csoldier", "csoldier", 2, 4, 0, false) # chariot soldier
        board[13][2] = piece("csoldier", "csoldier", 2, 13, 0, false) # chariot soldier
        board[14][2] = piece("csoldier", "csoldier", 2, 14, 0, false) # chariot soldier
        board[6][2] = piece("tiger", "tiger", 2, 6, 0, false) # blind tiger
        board[11][2] = piece("tiger", "tiger", 2, 11, 0, false) # blind tiger
        board[7][2] = piece("kirin", "kirin", 2, 7, 0, false)
        board[8][2] = piece("lion", "lion", 2, 8, 0, false)
        board[9][2] = piece("queen", "queen", 2, 9, 0, false)
        board[10][2] = piece("phoenix", "phoenix", 2, 10, 0, false)
        board[1][3] = piece("ssoldier", "ssoldier", 3, 1, 0, false) # side soldier
        board[16][3] = piece("ssoldier", "ssoldier", 3, 16, 0, false) # side soldier
        board[2][3] = piece("vsoldier", "vsoldier", 3, 2, 0, false) # vertical soldier
        board[15][3] = piece("vsoldier", "vsoldier", 3, 15, 0, false) # vertical soldier
        board[3][3] = piece("bishop", "bishop", 3, 3, 0, false)
        board[14][3] = piece("bishop", "bishop", 3, 14, 0, false)
        board[4][3] = piece("horse", "horse", 3, 4, 0, false) # dragon horse
        board[13][3] = piece("horse", "horse", 3, 13, 0, false) # dragon horse
        board[5][3] = piece("dragon", "dragon", 3, 5, 0, false) # dragon king
        board[12][3] = piece("dragon", "dragon", 3, 12, 0, false) # dragon king
        board[6][3] = piece("buffalo", "buffalo", 3, 6, 0, false) # water buffalo
        board[11][3] = piece("buffalo", "buffalo", 3, 11, 0, false) # water buffalo
        board[7][3] = piece("demon", "demon", 3, 7, 0, false) # fire demon
        board[10][3] = piece("demon", "demon", 3, 10, 0, false) # fire demon
        board[8][3] = piece("hawk", "hawk", 3, 8, 0, false) # lion hawk
        board[9][3] = piece("eagle", "eagle", 3, 9, 0, false) # free eagle
        board[1][4] = piece("smover", "smover", 4, 1, 0, false) # side mover
        board[16][4] = piece("smover", "smover", 4, 16, 0, false) # side mover
        board[2][4] = piece("vmover", "vmover", 4, 2, 0, false) # vertical mover
        board[15][4] = piece("vmover", "vmover", 4, 15, 0, false) # vertical mover
        board[3][4] = piece("rook", "rook", 4, 3, 0, false)
        board[14][4] = piece("rook", "rook", 4, 14, 0, false)
        board[4][4] = piece("falcon", "falcon", 4, 4, 0, false) # horned falcon
        board[13][4] = piece("falcon", "falcon", 4, 13, 0, false) # horned falcon
        board[5][4] = piece("soaring", "soaring", 4, 5, 0, false) # soaring eagle
        board[12][4] = piece("soaring", "soaring", 4, 12, 0, false) # soaring eagle
        board[6][4] = piece("bgeneral", "bgeneral", 4, 6, 0, false) # bishop general
        board[11][4] = piece("bgeneral", "bgeneral", 4, 11, 0, false) # bishop general
        board[7][4] = piece("rgeneral", "rgeneral", 4, 7, 0, false) # rook general
        board[10][4] = piece("rgeneral", "rgeneral", 4, 10, 0, false) # rook general
        board[8][4] = piece("great", "great", 4, 8, 0, false) # great general
        board[9][4] = piece("vice", "vice", 4, 9, 0, false) # vice general
        board[1][5] = piece("pawn", "pawn", 5, 1, 0, false)
        board[2][5] = piece("pawn", "pawn", 5, 2, 0, false)
        board[3][5] = piece("pawn", "pawn", 5, 3, 0, false)
        board[4][5] = piece("pawn", "pawn", 5, 4, 0, false)
        board[5][5] = piece("pawn", "pawn", 5, 5, 0, false)
        board[6][5] = piece("pawn", "pawn", 5, 6, 0, false)
        board[7][5] = piece("pawn", "pawn", 5, 7, 0, false)
        board[8][5] = piece("pawn", "pawn", 5, 8, 0, false)
        board[9][5] = piece("pawn", "pawn", 5, 9, 0, false)
        board[10][5] = piece("pawn", "pawn", 5, 10, 0, false)
        board[11][5] = piece("pawn", "pawn", 5, 11, 0, false)
        board[12][5] = piece("pawn", "pawn", 5, 12, 0, false)
        board[13][5] = piece("pawn", "pawn", 5, 13, 0, false)
        board[14][5] = piece("pawn", "pawn", 5, 14, 0, false)
        board[15][5] = piece("pawn", "pawn", 5, 15, 0, false)
        board[16][5] = piece("pawn", "pawn", 5, 16, 0, false)
        board[5][6] = piece("dog", "dog", 6, 5, 0, false)
        board[12][6] = piece("dog", "dog", 6, 12, 0, false)

    elseif gameType == "chu" # game setup for chu shogi

        # sente/black/odd pieces
        board[4][8] = piece("cobra", "cobra", 8, 4, 1, false) # AKA go-between
        board[9][8] = piece("cobra", "cobra", 8, 9, 1, false) # AKA go-between
        board[1][9] = piece("pawn", "pawn", 9, 1, 1, false)
        board[2][9] = piece("pawn", "pawn", 9, 2, 1, false)
        board[3][9] = piece("pawn", "pawn", 9, 3, 1, false)
        board[4][9] = piece("pawn", "pawn", 9, 4, 1, false)
        board[5][9] = piece("pawn", "pawn", 9, 5, 1, false)
        board[6][9] = piece("pawn", "pawn", 9, 6, 1, false)
        board[7][9] = piece("pawn", "pawn", 9, 7, 1, false)
        board[8][9] = piece("pawn", "pawn", 9, 8, 1, false)
        board[9][9] = piece("pawn", "pawn", 9, 9, 1, false)
        board[10][9] = piece("pawn", "pawn", 9, 10, 1, false)
        board[11][9] = piece("pawn", "pawn", 9, 11, 1, false)
        board[12][9] = piece("pawn", "pawn", 9, 12, 1, false)
        board[1][10] = piece("smover", "smover", 10, 1, 1, false)
        board[2][10] = piece("vmover", "vmover", 10, 2, 1, false)
        board[3][10] = piece("rook", "rook", 10, 3, 1, false)
        board[4][10] = piece("horse", "horse", 10, 4, 1, false) # dragon horse
        board[5][10] = piece("dragon", "dragon", 10, 5, 1, false) # dragon king
        board[6][10] = piece("queen", "queen", 10, 6, 1, false)
        board[7][10] = piece("lion", "lion", 10, 7, 1, false)
        board[8][10] = piece("dragon", "dragon", 10, 8, 1, false) # dragon king
        board[9][10] = piece("horse", "horse", 10, 9, 1, false) # dragon horse
        board[10][10] = piece("rook", "rook", 10, 10, 1, false)
        board[11][10] = piece("vmover", "vmover", 10, 11, 1, false) # vertical mover
        board[12][10] = piece("smover", "smover", 10, 12, 1, false) # side mover
        board[1][11] = piece("chariot", "chariot", 11, 1, 1, false) # reverse chariot
        board[3][11] = piece("bishop", "bishop", 11, 3, 1, false)
        board[5][11] = piece("tiger", "tiger", 11, 5, 1, false) # blind tiger
        board[6][11] = piece("phoenix", "phoenix", 11, 6, 1, false)
        board[7][11] = piece("kirin", "kirin", 11, 7, 1, false)
        board[8][11] = piece("tiger", "tiger", 11, 8, 1, false) # blind tiger
        board[10][11] = piece("bishop", "bishop", 11, 10, 1, false)
        board[12][11] = piece("chariot", "chariot", 11, 12, 1, false) # reverse chariot
        board[1][12] = piece("lance", "lance", 12, 1, 1, false)
        board[2][12] = piece("leopard", "leopard", 12, 2, 1, false) # ferocious leopard
        board[3][12] = piece("copper", "copper", 12, 3, 1, false) # copper general
        board[4][12] = piece("silver", "silver", 12, 4, 1, false) # silver general
        board[5][12] = piece("gold", "gold", 12, 5, 1, false) # gold general
        board[6][12] = piece("elephant", "elephant", 12, 6, 1, false) # drunk elephant
        board[7][12] = piece("king", "king", 12, 7, 1, false)
        board[8][12] = piece("gold", "gold", 12, 8, 1, false) # gold general
        board[9][12] = piece("silver", "silver", 12, 9, 1, false) # silver general
        board[10][12] = piece("copper", "copper", 12, 10, 1, false) # copper general
        board[11][12] = piece("leopard", "leopard", 12, 11, 1, false) # ferocious leopard
        board[12][12] = piece("lance", "lance", 12, 12, 1, false)

        # gote/white/even pieces
        board[4][5] = piece("cobra", "cobra", 5, 4, 0, false) # AKA go-between
        board[9][5] = piece("cobra", "cobra", 5, 9, 0, false) # AKA go-between
        board[1][4] = piece("pawn", "pawn", 4, 1, 0, false)
        board[2][4] = piece("pawn", "pawn", 4, 2, 0, false)
        board[3][4] = piece("pawn", "pawn", 4, 3, 0, false)
        board[4][4] = piece("pawn", "pawn", 4, 4, 0, false)
        board[5][4] = piece("pawn", "pawn", 4, 5, 0, false)
        board[6][4] = piece("pawn", "pawn", 4, 6, 0, false)
        board[7][4] = piece("pawn", "pawn", 4, 7, 0, false)
        board[8][4] = piece("pawn", "pawn", 4, 8, 0, false)
        board[9][4] = piece("pawn", "pawn", 4, 9, 0, false)
        board[10][4] = piece("pawn", "pawn", 4, 10, 0, false)
        board[11][4] = piece("pawn", "pawn", 4, 11, 0, false)
        board[12][4] = piece("pawn", "pawn", 4, 12, 0, false)
        board[1][3] = piece("smover", "smover", 3, 1, 0, false) # side mover
        board[2][3] = piece("vmover", "vmover", 3, 2, 0, false) # vertical mover
        board[3][3] = piece("rook", "rook", 3, 3, 0, false)
        board[4][3] = piece("horse", "horse", 3, 4, 0, false) # dragon horse
        board[5][3] = piece("dragon", "dragon", 3, 5, 0, false) # dragon king
        board[6][3] = piece("lion", "lion", 3, 6, 0, false)
        board[7][3] = piece("queen", "queen", 3, 7, 0, false)
        board[8][3] = piece("dragon", "dragon", 3, 8, 0, false) # dragon king
        board[9][3] = piece("horse", "horse", 3, 9, 0, false) # dragon horse
        board[10][3] = piece("rook", "rook", 3, 10, 0, false)
        board[11][3] = piece("vmover", "vmover", 3, 11, 0, false) # vertical mover
        board[12][3] = piece("smover", "smover", 3, 12, 0, false) # side mover
        board[1][2] = piece("chariot", "chariot", 2, 1, 0, false) # reverse chariot
        board[3][2] = piece("bishop", "bishop", 2, 3, 0, false)
        board[5][2] = piece("tiger", "tiger", 2, 5, 0, false) # blind tiger
        board[6][2] = piece("kirin", "kirin", 2, 6, 0, false)
        board[7][2] = piece("phoenix", "phoenix", 2, 7, 0, false)
        board[8][2] = piece("tiger", "tiger", 2, 8, 0, false) # blind tiger
        board[10][2] = piece("bishop", "bishop", 2, 10, 0, false)
        board[12][2] = piece("chariot", "chariot", 2, 12, 0, false) # reverse chariot
        board[1][1] = piece("lance", "lance", 1, 1, 0, false)
        board[2][1] = piece("leopard", "leopard", 1, 2, 0, false) # ferocious leopard
        board[3][1] = piece("copper", "copper", 1, 3, 0, false) # copper general
        board[4][1] = piece("silver", "silver", 1, 4, 0, false) # silver general
        board[5][1] = piece("gold", "gold", 1, 5, 0, false) # gold general
        board[6][1] = piece("king", "king", 1, 6, 0, false)
        board[7][1] = piece("elephant", "elephant", 1, 7, 0, false) # drunk elephant
        board[8][1] = piece("gold", "gold", 1, 8, 0, false) # gold general
        board[9][1] = piece("silver", "silver", 1, 9, 0, false) # silver general
        board[10][1] = piece("copper", "copper", 1, 10, 0, false) # copper general
        board[11][1] = piece("leopard", "leopard", 1, 11, 0, false) #ferocious leopard
        board[12][1] = piece("lance", "lance", 1, 12, 0, false)

    elseif gameType == "standard" # game setup for standard shogi

        # sente/black/odd player
        board[1][7] = piece("pawn", "pawn", 7, 1, 1, false)
        board[2][7] = piece("pawn", "pawn", 7, 2, 1, false)
        board[3][7] = piece("pawn", "pawn", 7, 3, 1, false)
        board[4][7] = piece("pawn", "pawn", 7, 4, 1, false)
        board[5][7] = piece("pawn", "pawn", 7, 5, 1, false)
        board[6][7] = piece("pawn", "pawn", 7, 6, 1, false)
        board[7][7] = piece("pawn", "pawn", 7, 7, 1, false)
        board[8][7] = piece("pawn", "pawn", 7, 8, 1, false)
        board[9][7] = piece("pawn", "pawn", 7, 9, 1, false)
        board[2][8] = piece("rook", "rook", 8, 2, 1, false)
        board[8][8] = piece("bishop", "bishop", 8, 8, 1, false)
        board[1][9] = piece("lance", "lance", 9, 1, 1, false)
        board[2][9] = piece("knight", "knight", 9, 2, 1, false)
        board[3][9] = piece("silver", "silver", 9, 3, 1, false) # silver general
        board[4][9] = piece("gold", "gold", 9, 4, 1, false) # gold general
        board[5][9] = piece("king", "king", 9, 5, 1, false)
        board[6][9] = piece("gold", "gold", 9, 6, 1, false) # gold general
        board[7][9] = piece("silver", "silver", 9, 7, 1, false) # silver general
        board[8][9] = piece("knight", "knight", 9, 8, 1, false)
        board[9][9] = piece("lance", "lance", 9, 9, 1, false)

        # gote/white/even player
        board[1][3] = piece("pawn", "pawn", 3, 1, 0, false)
        board[2][3] = piece("pawn", "pawn", 3, 2, 0, false)
        board[3][3] = piece("pawn", "pawn", 3, 3, 0, false)
        board[4][3] = piece("pawn", "pawn", 3, 4, 0, false)
        board[5][3] = piece("pawn", "pawn", 3, 5, 0, false)
        board[6][3] = piece("pawn", "pawn", 3, 6, 0, false)
        board[7][3] = piece("pawn", "pawn", 3, 7, 0, false)
        board[8][3] = piece("pawn", "pawn", 3, 8, 0, false)
        board[9][3] = piece("pawn", "pawn", 3, 9, 0, false)
        board[2][2] = piece("bishop", "bishop", 2, 2, 0, false)
        board[8][2] = piece("rook", "rook", 2, 8, 0, false)
        board[1][1] = piece("lance", "lance", 1, 1, 0, false)
        board[2][1] = piece("knight", "knight", 1, 2, 0, false)
        board[3][1] = piece("silver", "silver", 1, 3, 0, false) # silver general
        board[4][1] = piece("gold", "gold", 1, 4, 0, false) # gold general
        board[5][1] = piece("king", "king", 1, 5, 0, false)
        board[6][1] = piece("gold", "gold", 1, 6, 0, false) # gold general
        board[7][1] = piece("silver", "silver", 1, 7, 0, false) # silver general
        board[8][1] = piece("knight", "knight", 1, 8, 0, false)
        board[9][1] = piece("lance", "lance", 1, 9, 0, false)

    else # game setup for minishogi

        # sente/black/odd player
        board[5][4] = piece("pawn", "pawn", 4, 5, 1, false)
        board[1][5] = piece("rook", "rook", 5, 1, 1, false)
        board[2][5] = piece("bishop", "bishop", 5, 2, 1, false)
        board[3][5] = piece("silver", "silver", 5, 3, 1, false) # silver general
        board[4][5] = piece("gold", "gold", 5, 4, 1, false) # gold general
        board[5][5] = piece("king", "king", 5, 5, 1, false)

        # gote/white/even player
        board[1][2] = piece("pawn", "pawn", 2, 1, 0, false)
        board[1][1] = piece("king", "king", 1, 1, 0, false)
        board[2][1] = piece("gold", "gold", 1, 2, 0, false) # gold general
        board[3][1] = piece("silver", "silver", 1, 3, 0, false) # silver general
        board[4][1] = piece("bishop", "bishop", 1, 4, 0, false)
        board[5][1] = piece("rook", "rook", 1, 5, 0, false)

    end
    return board, boardLength, numOfPromotionRanks
end

max_stars_repo_name: DarienFlamont/Shogi

max_issues_repo_name: DarienFlamont/Shogi

max_forks_repo_name: DarienFlamont/Shogi


text:
const FlatlandGraph = DataDiGraph{Int,NTuple{4,Int},Nothing,Float64,Matrix{UInt16}}

function flatland_graph(pyenv::Py, agents::Vector{Agent})
    # Find stations
    initial_positions = unique(agent.initial_position for agent in agents)
    target_positions = unique(agent.target_position for agent in agents)

    # Retrieve grid
    width = pyconvert(Int, pyenv.width)
    height = pyconvert(Int, pyenv.height)
    grid = pyconvert(Matrix{UInt16}, pyenv.rail.grid)

    # Initialize g
    g = DataDiGraph{Int}(;
        VL=NTuple{4,Int}, VD=Nothing, ED=Float64, graph_data=grid
    )

    # Create vertices from non empty grid cells
    for i in 1:height, j in 1:width
        cell = grid[i, j]
        if cell > 0
            transition_map = bitstring(cell)
            # Real vertices
            for direction in CARDINAL_POINTS
                if direction_exists(transition_map, direction)
                    label = (i, j, direction, REAL)
                    add_vertex!(g, label)
                end
            end
            # Departure vertices
            if (i, j) in initial_positions
                for direction in CARDINAL_POINTS
                    if direction_exists(transition_map, direction)
                        label = (i, j, direction, DEPARTURE)
                        add_vertex!(g, label)
                    end
                end
            end
            # Arrival vertices
            if (i, j) in target_positions
                label = (i, j, NO_DIRECTION, ARRIVAL)
                add_vertex!(g, label)
            end
        end
    end

    # Create out edges for every vertex
    for v in vertices(g)
        label_s = get_label(g, v)
        (i, j, direction, kind) = label_s
        if kind == REAL  # from real vertices
            cell = grid[i, j]
            transition_map = bitstring(cell)
            # to themselves
            add_edge!(g, label_s, label_s, 1.0)
            # to other real vertices
            for out_direction in CARDINAL_POINTS
                if transition_exists(transition_map, direction, out_direction)
                    i2, j2 = neighbor_cell(i, j, out_direction)
                    label_d = (i2, j2, out_direction, REAL)
                    add_edge!(g, label_s, label_d, 1.0)
                end
            end
            # to arrival vertices
            if (i, j) in target_positions
                label_d = (i, j, NO_DIRECTION, ARRIVAL)
                add_edge!(g, label_s, label_d, 0.0)
            end
        elseif kind == DEPARTURE  # from departure vertices
            # to themselves
            add_edge!(g, label_s, label_s, 1.0)
            # to real vertices
            label_d = (i, j, direction, REAL)
            add_edge!(g, label_s, label_d, 0.0)
        elseif kind == ARRIVAL  # from arrival vertices
            # to themselves
            add_edge!(g, label_s, label_s, 0.0)
        end
    end
    return g
end

get_grid(g::FlatlandGraph) = get_data(g)
get_height(g::FlatlandGraph) = size(get_grid(g), 1)
get_width(g::FlatlandGraph) = size(get_grid(g), 2)


text:
using SimGBS

## parameters
### Step One: Generate GBS Fragments
genofile = "ref.fa.gz"
re = [SimGBS.ApeKI]; # specify the restriction enzyme to be used for virtual digestion
useChr = [1]; # specify either the number of chromosome(s) or a set of chromosome(s) to be used for simulating GBS data
useChrLen = Array{Float64}(undef, 0); # specify the length of chromosome(s) in cM to be simulated, otherwise using the entire chromosome (Array{Float64}(undef, 0) means empty set)
lower = 65; # lower bound of fragment size selection
upper = 195; # upper bound of fragment size selection
winSize = 1000000 # window size to be used to sample SNP density

### Step Two: Define Population Structure
numFounders = 100; # number of founders in the base population
endSize = 1000; # number of individuals to end up in the changingPopSize step
numGenCha = 20; # number of generations for changingPopSize function
numGenCon = 50; # number of generations for constantPopSize function
numGenFinal = 4; # number of final generations to be used to select individual
numInd = 96; # number of individuals to be simulated
useWeights = Array{Float64}(undef, 0) # weights of each contributing genetarion in the final population composition
usePedigree = false;  # false if you don't use pedigree, otherwise specify the pedigree file to be used
pedFile = "newPed.ped"; # file stores the pedigree (default = "sim.ped")
pedOutput = false; # true if print out pedigree (default = false)

### Step Three: Simulate GBS Process
totalQTL = 100 # total (across all chromosomes) number of QTLs to be simulated
totalSNP = 0 # total (across all chromosomes) number of QTLs to be simulated, set [totalSNP = 0] if sampling based on density
muDensity = -2.0 # location parameter of log-Laplace distribution (for sampling SNP density)
sigmasqDensity = 0.001 # scale parameter of log-Laplace distribution (for sampling SNP density)
winSize = 1000000 # window size to be used to sample SNP density
muAlleleFreq = 0.5 # mean of sampled allele frequency
sigmasqAlleleFreq = 0.0061 # variance of sampled allele frequency
meanDepth = 20.0 # expected sequencing depth
barcodeFile = "GBS_Barcodes.txt" # file contains GBS barcodes
useChr = [1]; # specify either the number of chromosome(s) or a set of chromosome(s) to be used for simulating GBS data

### miscellaneous
plotOutput = false; # true if plots are required
writeOutput = true; # true if outputs in text file are required
onlyOutputGBS = true; # true if only GBS data is required

## Run SimGBS
## Note: Please maker sure that your working directory is same as where you've downloaed input files!
## You can use pwd() to check your current working directory
## or cd() to change your working directory

### Step One: Generate GBS Fragments
@time digestGenome(
    genofile,
    re,
    useChr,
    useChrLen,
    lower,
    upper,
    winSize,
    plotOutput,
    writeOutput,
);

### Step Two: Define Population Structure
@time definePopulation(
    numFounders,
    endSize,
    numGenCha,
    numGenCon,
    numGenFinal,
    numInd,
    useWeights,
    usePedigree,
    pedFile,
    pedOutput,
);

### Step Three: Simulate GBS Process
@time GBS(
    totalQTL,
    totalSNP,
    muDensity,
    sigmasqDensity,
    winSize,
    muAlleleFreq,
    sigmasqAlleleFreq,
    re,
    meanDepth,
    barcodeFile,
    useChr,
    plotOutput,
    writeOutput,
    onlyOutputGBS,
);


text:
using OptimizationProblems

## All pbs dimension > 1 and not in CUTEst
Problems = [brownden, cliff, clplatea, clplateb, clplatec, nasty,
palmer1c, palmer1d, palmer2c, palmer3c, palmer4c, palmer5c, palmer5d,
palmer6c, palmer7c, palmer8c,
beale, brownbs,
chnrosnb_mod,
errinros_mod, fletcbv3_mod, genrose_nash,
meyer3, scosine]

algo = CGlin
# algo = CRlin

D = ["model     algo nvar   f(x)    f(x0)   ‖g(x)‖  ‖g(x0)‖   #f  #g  #Hv  #it #s_it #vs_it #rej_it optimal"]

for problem in Problems
    model = MathProgNLPModel(problem(), name=string(problem))
    L = Linesearch(model, algo)
    S = @sprintf("%5s %5s %4d %5s %5s %5s %5s %4d %4d %4d %4d %5s", L[1], L[2], L[3], L[4], L[5], L[6], L[7], L[8], L[9], L[10], L[11], L[12])
    D = vcat(D, S)
    reset!(model)
end
writedlm("PbsCG.txt", D)


text:
using StatsPlots
using IndexedTables

# This demonstrates the PROCESS by plotting a correlation matrix I THINK.
include(joinpath(WATERFALL_DIR,"src","includes.jl"))

begin value=:Value; label=:Label; sample=:Sample end
distribution=:normal
fuzziness=(0.01,0.3)
numsample=1000
ncor=2
kwargs = (label=label, distribution=distribution, fuzziness=fuzziness,
    numsample=numsample,
    ncor=ncor,
)

df = CSV.read(joinpath(DATA_DIR,"pvrd2-investment-metrics.csv"), DataFrames.DataFrame)
df[!,:Label] .= string.(titlecase.(df[:,:Process])," (",df[:,:Step],")")
df[!,:Units] .= "Efficiency (%)"

label=:Process
df = DataFrames.transform(fuzzify(df; kwargs...))
dfw = DataFrames.unstack(DataFrames.select(df, [sample,label,value]), label, value)
steps = size(dfw,2)-1

cols = propertynames(dfw)[2:end]
vcorr = random_correlation(steps, ncor)
# @df iris corrplot(cols(1:4), grid = false)

StatsPlots.gr(size=(800,800))

iis = sum(hcat([vcat(sum(LinearAlgebra.UnitLowerTriangular(vcorr); dims=ii)...) for ii in 1:2]...); dims=2).!==2.0
iis = unique([[ii for ii in collect(2:steps+1) .* iis if ii!==0]; steps+1])

@df dfw StatsPlots.corrplot(cols(iis),
    grid=false,
    compact=true,
    label =["$i" for i=iis],
)

# using StatsPlots
# using IndexedTables

# import RDatasets

# df = DataFrames.DataFrame(a = 1:10, b = 10 .* rand(10), c = 10 .* rand(10))
# @df df plot(:a, [:b :c], colour = [:red :blue])
# @df df scatter(:a, :b, markersize = 4 .* log.(:c .+ 0.1))
# t = table(1:10, rand(10), names = [:a, :b]) # IndexedTable
# @df t scatter(2 .* :b)

# @df df plot(:a, cols(2:3), colour = [:red :blue])

# # SCHOOL
# # school = RDatasets.dataset("mlmRev","Hsb82")
# # @df school density(:MAch, group = (:Sx, :Sector), legend = :topleft)


text:
export Addition, +, -

"""
Description:

    An addition constraint factor node.

    f(out,in1,in2) = δ(in1 + in2 - out)

Interfaces:

    1. out
    2. in1
    3. in2

Construction:

    Addition(out, in1, in2, id=:some_id)
"""
mutable struct Addition <: DeltaFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function Addition(out, in1, in2; id=generateId(Addition))
        @ensureVariables(out, in1, in2)
        self = new(id, Array{Interface}(undef, 3), Dict{Int,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:in1] = self.interfaces[2] = associate!(Interface(self), in1)
        self.i[:in2] = self.interfaces[3] = associate!(Interface(self), in2)

        return self
    end
end

slug(::Type{Addition}) = "+"

function +(in1::Variable, in2::Variable)
    out = Variable()
    Addition(out, in1, in2)
    return out
end

function +(in1::Variable, in2)
    @ensureVariables(in2)
    in1 + in2
end

+(in1, in2::Variable) = +(in2, in1)

"""
    -(in1::Variable, in2::Variable)

A subtraction constraint based on the addition factor node.

"""
function -(in1::Variable, in2::Variable)
    out = Variable()
    Addition(in1, out, in2)
    return out
end

function -(in1::Variable, in2)
    @ensureVariables(in2)
    in1 - in2
end

function -(in1, in2::Variable)
    @ensureVariables(in1)
    in1 - in2
end


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 = "MMseqs2"
version = v"13"

# MMseqs2 seem to use as versioning scheme of "major version + first 5
# characters of the tagged commit"
# https://github.com/soedinglab/MMseqs2/releases
version_commitprefix = "45111"

# Possible build variants
# - OpenMP (default)
# - MPI (-DHAVE_MPI=1)
# - single-threaded (-DREQUIRE_OPENMP=0)

# TODO
# - built-in zstd, replace with Zstd_jll? (-DUSE_SYSTEM_ZSTD=1 cmake option)
# - aarch64 macos: single-threaded, openmp not used
# - os-specific build script examples under util/build_{osx,windows}

# Build failures
# - windows:
#   - cmake can't find openmp, this check can be avoided by passing -DREQUIRE_OPENMP=0 to cmake
#   - compile error afterwards
#     error: ‘posix_memalign’ was not declared in this scope
#     (and more following errors)
# - macos aarch64
#   - g++-11 build fails (tries to compile for x86 simd)
#   - compiles with clang, but openmp doesn't work
# - i686: compile error due to bitwidth issues, haven't investigated more
# - powerpc build fails with g++-7.x (tries to compile for x86 simd),
#   works with g++-8.x and above

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/soedinglab/MMseqs2/archive/refs/tags/$(version.major)-$(version_commitprefix).tar.gz",
                  "6444bb682ebf5ced54b2eda7a301fa3e933c2a28b7661f96ef5bdab1d53695a2"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd MMseqs2-*/

# patch CMakeLists.txt so it doesn't set -march unnecessarily on ARM
atomic_patch -p1 ../patches/arm-simd-march-cmakefile.patch

# macos x86_64: use gcc/g++ so we can use openmp
if [[ "${target}" == x86_64-apple-darwin* ]]; then
    CMAKE_TARGET_TOOLCHAIN="${CMAKE_TARGET_TOOLCHAIN/%.cmake/_gcc.cmake}"
    echo "[INFO] setting CMAKE_TARGET_TOOLCHAIN = ${CMAKE_TARGET_TOOLCHAIN}"
fi

# macos aarch64: disable openmp as it doesn't work with clang
EXTRA_CMAKE_FLAGS=
if [[ "${target}" == aarch64-apple-darwin* ]]; then
    EXTRA_CMAKE_FLAGS="-DREQUIRE_OPENMP=0"
fi

# architecture extensions
ARCH_FLAGS=
if [[ "${target}" == x86_64-* || "${target}" == i686-* ]]; then
    ARCH_FLAGS="-DHAVE_SSE2=1 -DHAVE_SSE4_1=1 -DHAVE_AVX2=1"
elif [[ "${target}" == powerpc64le-* ]]; then
    ARCH_FLAGS="-DHAVE_POWER8=1 -DHAVE_POWER9=1"
elif [[ "${target}" == aarch64-* ]]; then
    ARCH_FLAGS="-DHAVE_ARM8=1"
fi

mkdir build
cd build
cmake .. \
    -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=RELEASE \
    -DNATIVE_ARCH=0 ${EXTRA_CMAKE_FLAGS} ${ARCH_FLAGS}
make -j${nproc}
make install

install_license ../LICENSE.md
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true, exclude = p -> Sys.iswindows(p) || arch(p) == "i686")
# expand cxxstring abis on platforms where we use g++
platforms = expand_cxxstring_abis(platforms; skip = p -> Sys.isfreebsd(p) || (Sys.isapple(p) && arch(p) == "aarch64"))

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

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

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
fname = joinpath(dirname(@__FILE__), "..", "data", "test-4d.dat")
sname = joinpath(dirname(@__FILE__), "..", "data", "tmp", "same.dat")

x, y, z, s, t = read_dat(fname)

@fact size(x) --> (30,)
@fact size(y) --> (36,)
@fact size(z) --> (28,)
@fact size(s) --> (30,36,28,2)
@fact size(t) --> (2,)

@fact maximum(x) --> 72.5
@fact maximum(y) --> 71.220001
@fact maximum(z) --> 76.809998
@fact maximum(s) --> 0.067409396
@fact maximum(t) --> 0.24

@fact minimum(x) --> -72.5
@fact minimum(y) --> -103.779999
@fact minimum(z) --> -58.189999
@fact minimum(s) --> 0.0
@fact minimum(t) --> 0.12

write_dat(sname, x, y, z, s[:,:,:,:], t)

x2, y2, z2, s2, t2 = read_dat(sname)

@fact x --> x2
@fact y --> y2
@fact z --> z2
@fact s --> s2
@fact t --> t2

fname = joinpath(dirname(@__FILE__), "..", "data", "test-3d.dat")
x, y, z, s, t = read_dat(fname)

@fact size(x) --> (30,)
@fact size(y) --> (36,)
@fact size(z) --> (28,)
@fact size(s) --> (30,36,28,1)
@fact size(t) --> (1,)

@fact maximum(x) --> 72.5
@fact maximum(y) --> 71.220001
@fact maximum(z) --> 76.809998
@fact maximum(s) --> 33.2692985535
@fact maximum(t) --> 0

@fact minimum(x) --> -72.5
@fact minimum(y) --> -103.779999
@fact minimum(z) --> -58.189999
@fact minimum(s) --> -7.5189352036
@fact minimum(t) --> 0


text:
# We use nHidden as a vector, containing the number of hidden units in each layer
using DataFrames
# Function that returns total number of parameters
function NeuralNet_nParams_class(d,nHidden,k)

	# Connections from inputs to first hidden layer
	nParams = d*nHidden[1]

	# Connections between hidden layers
	for h in 2:length(nHidden)
		nParams += nHidden[h-1]*nHidden[h]
	end

	# Connections from last hidden layer to output
	nParams += nHidden[end]*k

end

# Compute squared error and gradient
# for a single training example (x,y)
# (x is assumed to be a column-vector)
function NeuralNet_backprop_class(bigW,x,y,nHidden,k)
	d = length(x)
	nLayers = length(nHidden)

	#### Reshape 'bigW' into vectors/matrices
	W1 = reshape(bigW[1:nHidden[1]*d],nHidden[1],d)
	ind = nHidden[1]*d
	Wm = fill([],(nLayers-1))
	for layer in 2:nLayers
		Wm[layer-1] = reshape(bigW[ind+1:ind+nHidden[layer]*nHidden[layer-1]],nHidden[layer],nHidden[layer-1])
		ind += nHidden[layer]*nHidden[layer-1]
	end
	w = reshape(bigW[ind+1:end],k,nHidden[nLayers])

	### Bin y variables

	#### Define activation function and its derivative
	h(z) = tanh.(z)
	dh(z) = (sech.(z)).^2

	#### Forward propagation
	z = fill([],(nLayers))
	z[1] = W1*x
	for layer in 2:nLayers
		z[layer] = Wm[layer-1]*h(z[layer-1])
	end
	yhat,_ = findmax(w*h(z[end]))
	yhat=yhat-1

	f,r = softmax(w,x,y,k)

	#### Backpropagation
	dr = r
	err = dr

	# Output weights
	Gout = err*h(z[end])

	Gm = fill([],(nLayers-1))
	if nLayers > 1
		# Last Layer of Hidden Weights
		backprop = err*(dh(z[end]).*w)
		Gm[end] = backprop*h(z[end-1])'

		# Other Hidden Layers
		for layer in nLayers-2:-1:1
			backprop = (Wm[layer+1]'*backprop).*dh(z[layer+1])
			Gm[layer] = backprop*h(z[layer])'
		end

		# Input Weights
		backprop = (Wm[1]'*backprop).*dh(z[1])
		G1 = backprop*x'
	else
		# Input weights
		print(size(err),size(x'))
		G1 = err*(dh(z[1]).*w)*x'
	end

	#### Put gradients into vector
	g = zeros(size(bigW))
	g[1:nHidden[1]*d] = G1
	ind = nHidden[1]*d
	for layer in 2:nLayers
		g[ind+1:ind+nHidden[layer]*nHidden[layer-1]] = Gm[layer-1]
		ind += nHidden[layer]*nHidden[layer-1]
	end
	g[ind+1:end] = Gout

	return (f,g)
end

# Computes predictions for a set of examples X

function NeuralNet_predict_class(bigW,Xhat,nHidden,k)
	(t,d) = size(Xhat)
	nLayers = length(nHidden)

	#### Reshape 'bigW' into vectors/matrices
	W1 = reshape(bigW[1:nHidden[1]*d],nHidden[1],d)
	ind = nHidden[1]*d
	Wm = fill([],(nLayers-1))
	for layer in 2:nLayers
		Wm[layer-1] = reshape(bigW[ind+1:ind+nHidden[layer]*nHidden[layer-1]],nHidden[layer],nHidden[layer-1])
		ind += nHidden[layer]*nHidden[layer-1]
	end
	w = reshape(bigW[ind+1:end],k,nHidden[nLayers])

	#### Define activation function and its derivative
	h(z) = tanh.(z)
	dh(z) = (sech.(z)).^2

	#### Forward propagation on each example to make predictions
	yhat = zeros(t,1)
	for i in 1:t
		# Forward propagation
		z = fill([],(nLayers))
		z[1] = W1*Xhat[i,:]
		for layer in 2:nLayers
			z[layer] = Wm[layer-1]*h(z[layer-1])
		end
		yhat[i] = w'*h(z[end])
	end
	return yhat
end

function softmax(w,X,y,k)
	d=length(X)

	yunique=unique(y)
	f=0.0
	g=zeros(k,d)
	for i in 1:n
		print(yunique.==y[i])
		for j in length(yunique)
			if yunique[j]==y[i]
				s=j
			end
		end
		s=firstindex((yunique.==y[i]),1)
		print(s)
        wyi=transpose(w[s,:])
        f+=log(sum(exp.((w*X[i,:]))))-wyi*X[i,:]

        for r in 1:k
            if r==s
                indic=1.0
            else
                indic=0.0
			end

            for j in 1:d
                g[r,j]+=X[i,j]*((exp.((transpose(w[r,:])*X[i,:]))/sum(exp.((w*X[i,:]))))-indic)
	      	end
	  	end
	end

    return f,g
end


text:
using CoolDicts
using Test

cooldict = (
    empty = CoolDict(),
    one_arg = CoolDict(:a=>5),
    two_arg = CoolDict(:a=>5, :b=>2),
    vector = CoolDict([:a=>5, :b=>2]),
    hetero = CoolDict(:a=>5, :b=>[5, 6]),
    dict = CoolDict(Dict([:a=>5, :b=>2])),
    cooldict = CoolDict(CoolDict(:a=>5)),
    kwargs = CoolDict(a=5, b=2),
    nt = CoolDict((a=5, b=2)),
    empty_explicit = CoolDict{Float64}(),
    one_arg_explicit = CoolDict{Float64}(:a=>5),
    dict_explicit = CoolDict{Float64}(Dict([:a=>5, :b=>2])),
    convert = CoolDict{Float64}(CoolDict(:a=>5, :b=>2)),
    kwargs_explicit = CoolDict{Float64}(a=5, b=2),
    nt_explicit = CoolDict{Float64}((a=5, b=2)),
)

@testset "Constructors" begin
    @test cooldict.empty isa CoolDict{Any}
    @test cooldict.one_arg isa CoolDict{Int}
    @test cooldict.two_arg isa CoolDict{Int}
    @test cooldict.hetero isa CoolDict{Any}
    @test cooldict.kwargs isa CoolDict{Int}
    @test cooldict.nt isa CoolDict{Int}
    @test cooldict.empty_explicit isa CoolDict{Float64}
    @test cooldict.one_arg_explicit isa CoolDict{Float64}
    @test cooldict.dict_explicit isa CoolDict{Float64}
    @test cooldict.convert isa CoolDict{Float64}
    @test cooldict.kwargs_explicit isa CoolDict{Float64}
    @test cooldict.nt_explicit isa CoolDict{Float64}
end

@testset "Interface" begin
    @testset "length" begin
        @test length(cooldict.empty) == 0
        @test length(cooldict.two_arg) == 2
    end

    @testset "iterate" begin
        out_pairs = Any[]
        out_key_val = Any[]

        for pair in cooldict.vector
            push!(out_pairs, pair)
        end
        @test out_pairs[1] == (:a => 5)

        for (key, val) in cooldict.vector
            push!(out_key_val, (key, val))
        end
        @test out_key_val[1] == (:a, 5)
    end

    @testset "{get|set}{index|property}" begin
        cooldict_vector = deepcopy(cooldict.vector)
        cooldict_vector.a = 100
        cooldict_vector[:b] = 200

        @test cooldict.vector[:a] == cooldict.vector.a == 5
        @test cooldict_vector.a == 100
        @test cooldict_vector.b == 200
    end

    @testset "keys, propertynames, values" begin
        @test collect(keys(cooldict.hetero)) == [:a, :b]
        @test propertynames(cooldict.hetero) == (:a, :b)
        @test collect(values(cooldict.hetero)) == [5, [5, 6]]
    end

    @testset "convert" begin
        @testset "Dict" begin
            @test convert(Dict, cooldict.vector) == Dict(:a=>5, :b=>2)
            @test Dict(cooldict.vector) == Dict(:a=>5, :b=>2)
        end

        @testset "NamedTuple" begin
            @test convert(NamedTuple, cooldict.vector) == (a=5, b=2)
            @test NamedTuple(cooldict.vector) == (a=5, b=2)
        end
    end
end


text:
using OrdinaryDiffEq
using DiffEqSensitivity
using Zygote
using Plots

function sir_ode!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    N = S+I+R
    @inbounds begin
        du[1] = -β*c*I/N*S
        du[2] = β*c*I/N*S - γ*I
        du[3] = γ*I
    end
    nothing
end;

δt = 1.0
tmax = 40.0
tspan = (0.0,tmax)
t = 0.0:δt:tmax
num_timepoints = length(t);

u0 = [990.0,10.0,0.0] # S,I,R
num_states = length(u0);

p = [0.05,10.0,0.25]; # β,c,γ
num_params = length(p);

prob_ode = ODEProblem(sir_ode!,u0,tspan,p);

sim_ode = (u0,p)-> solve(prob_ode,Tsit5(),u0=u0,p=p,saveat=t,sensealg=QuadratureAdjoint());

sol_ode = sim_ode(u0,p);

du0,dp = Zygote.jacobian(sim_ode,u0,p);

dβ = reshape(dp[:,1],(num_states,:))' # as β is the first parameter
dc = reshape(dp[:,2],(num_states,:))' # c is 2nd parameter
dγ = reshape(dp[:,3],(num_states,:))' # γ is 3rd parameter
dI₀ = reshape(du0[:,2],(num_states,:))'; # I₀ is the 2nd initial condition

plot(sol_ode.t,
     Array(sol_ode(t))',
     labels = ["S" "I" "R"],
     xlabel = "Time",
     ylabel = "Number")

l = @layout [a b; c d]
pl1 = plot(t,dβ,xlabel="Time",ylabel="dp",label=["S" "I" "R"],title="Sensitivity to β")
pl2 = plot(t,dc,xlabel="Time",ylabel="dp",label=["S" "I" "R"],title="Sensitivity to c")
pl3 = plot(t,dγ,xlabel="Time",ylabel="dp",label=["S" "I" "R"],title="Sensitivity to γ")
pl4 = plot(t,dI₀,xlabel="Time",ylabel="dp",label=["S" "I"  "R"],title="Sensitivity to I₀")
plot(pl1,pl2,pl3,pl4,layout=l)


text:
using Documenter, ZXCalculus
using Multigraphs

makedocs(;
    modules = [ZXCalculus],
    format = Documenter.HTML(prettyurls = !("local" in ARGS)),
    pages = ["Home" => "index.md",
        "Tutorials" => "tutorials.md",
        "Examples" => "examples.md",
        "APIs" => "api.md"
    ],
    repo = "https:/github.com/QuantumBFS/ZXCalculus.jl",
    sitename = "ZXCalculus.jl",
)

deploydocs(repo = "github.com/QuantumBFS/ZXCalculus.jl.git")


text:
using FranklinUtils

# ----------------------------------- #
# Academic blocks // General elements #
# ----------------------------------- #

@env function section(md; name="", class="wg-$name",rowclass="")
    id = Franklin.refstring(name)
    return html("""

""") * md * html("""

""")
end

@lx function sectionheading(title; class="")
    return html("""

$title

        """)
end

# Insert an image with given class, alt, src
@lx img(src; class="", alt="") = html("""""")

# ---------------------------------------- #
# Academic blocks // Landing page elements #
# ---------------------------------------- #

# Portrait block with a few optional fields: name, job title, social buttons
@lx function portrait(; name="", job="", link="", linkname="",
                     twitter="", gscholar="", github="", linkedin="")
    io = IOBuffer()
    write(io, html("
"))
    isempty(name) || write(io, html("
$name
"))
    isempty(job) || write(io, html("
$job
"))
    if !isempty(link)
        if isempty(linkname)
            write(io, html("""
$link
"""))
        else
            write(io, html("""
$linkname
"""))
        end
    end
    if !all(isempty, (twitter, gscholar, github, linkedin))
        write(io, html("
"))
        isempty(twitter) || write(io, html("""

"""))
        isempty(gscholar) || write(io, html("""

"""))
        isempty(github) || write(io, html("""

"""))
        isempty(linkedin) || write(io, html("""

"""))
        write(io, html("
"))
    end
    write(io, html("
"))
    return String(take!(io))
end

# Biography block with optional resume link
@env function biography(md; resume="")
    io = IOBuffer()
    write(io, html("""
Biography
""") * md)
    isempty(resume) || write(io, html("""

Download my resumé.
"""))
    return String(take!(io))
end

# Short CV block with a column for interests and one for education
@lx function shortcv(; interests=nothing, education=nothing)
    io = IOBuffer()
    write(io, html("""
"""))
    if !isnothing(interests)
        write(io, html("""

Interests

"""))
        for i in interests
            write(io, html("""
$i
"""))
        end
        write(io, html("""

"""))
    end
    if !isnothing(education)
        write(io, html("""

Education

"""))
        for (course, school) in education
            write(io, html("""

$course

$school

                """))
        end
        write(io, html("""

"""))
    end
    write(io, html("""
""")) # end row
    return String(take!(io))
end

# skill featurette
@lx function skill(name, sub=""; img="", fa="",
                   imgstyle="display:inline-block; width:56px;",
                   fastyle="")
    illustration = ""
    if !isempty(img)
        illustration = """"""
    elseif !isempty(fa)
        illustration = """"""
    end

    return """

          """ * illustration * """

$name

          $(ifelse(isempty(sub), "", "
$sub
"))

""" |> html
end

# experience cards
@lx function experience(; title="", company="", descr="",
                          from="", to="", location="", active=false,
                          first=active, last=false)
    fill = ifelse(active, "exp-fill", "")
    # elements for the vertical bar with filled/unfilled pill
    # they are assembled depending on 'first' so that they can connect.
    pill = """

"""
    vbar = """

"""
    vspace = """

"""

    return html("""

            $(ifelse(first, vspace * pill * vbar, ifelse(last, vbar * pill * vspace, vbar*pill*vbar)))

$title

$company

$from – $to
                  $location

""") * descr * html("""

""")
end

@lx function certificate(; title="", meta="", metalink="", date="", descr="",
                           cert="See certificate", certlink="")
    origin = ifelse(isempty(metalink), meta, """
        $meta
        """)
    certificate = ifelse(isempty(certlink), "", """
        $cert
        """)
    description = ifelse(isempty(descr), "", """

$(Franklin.fd2html(descr, internal=true))

        """)
    return html("""

$title

              $origin

              $date

            $description
            $certificate

        """)
end

# -------------------- #
# List of recent posts #
# -------------------- #

function all_posts()
    # expected file structure here is posts/YYYY/MM/name-of-post.md
    posts = Pair{String,Date}[]
    for (root, _, files) in walkdir(joinpath(Franklin.FOLDER_PATH[], "posts"))
        for file in files
            endswith(file, ".md") || continue
            ppath = joinpath(root, file)
            endswith(ppath, joinpath("posts", "index.md")) && continue
            tmp, fn = splitdir(ppath)
            tmp, mm = splitdir(tmp)
            tmp, yy = splitdir(tmp)
            rpath = joinpath("posts", yy, mm, splitext(fn)[1])
            date = pagevar(rpath, "pubdate")
            isnothing(date) && (date = Date("$yy-$mm-01"))
            push!(posts, rpath => date)
        end
    end
    # sort by chron order, most recent first
    return sort(posts, by=(e->e.second), rev=true)
end

function show_posts(posts; byyear=false)
    isempty(posts) && return ""
    curyear = year(posts[1].second)
    io = IOBuffer()
    byyear && write(io, """

$curyear

        """)
    for post in posts
        if byyear && year(post.second) < curyear
            curyear = year(post.second)
            write(io, """

$curyear

                """)
        end
        rpath = post.first
        title = pagevar(rpath, "title")
        isnothing(title) && (title = "Untitled")
        summary = pagevar(rpath, "summary")
        isnothing(summary) && (summary = "")
        date = Dates.format(post.second, dateformat"u d, Y")
        imgpath = pagevar(rpath, "img")
        if isnothing(imgpath)
            imgpath = ""
        else
            if imgpath[1] != '/'
                imgpath = "$imgpath"
            end
        end
        write(io, """

$title

$summary

Published on $date.

""")
    end
    return String(take!(io))
end

function hfun_recentposts(params)
    n = parse(Int, params[1])
    allposts = all_posts()
    posts = allposts[1:min(n, length(allposts))]
    return show_posts(posts)
end

function hfun_allposts()
    return show_posts(all_posts(), byyear=true)
end

# --------------------------------- #
# Working with Javascript Libraries #
# --------------------------------- #

function env_mermaid(e, _)
    content = Franklin.content(e)
    return "@@mermaid ~~~$content~~~@@"
end

max_stars_repo_name: Pramodh-G/ft-academic-resume

max_issues_repo_name: Pramodh-G/ft-academic-resume

max_forks_repo_name: Pramodh-G/ft-academic-resume