text:
module GR_Spherical

using DifferentialEquations
using BoundaryValueDiffEq
using OrdinaryDiffEq
using Fun1d
using DataFrames
using CSV
using Plots
using Roots

using BenchmarkTools
using InteractiveUtils
using RecursiveArrayTools
using StaticArrays
using LinearAlgebra

using Profile

struct Param{T}
    rtmin::T
    rtmax::T
    drt::T
    Mtot::T
    grid::Grid{T}
    r::Function
    drdrt::Function
    d2rdrt::Function
    rsamp::Vector{T}
    drdrtsamp::Vector{T}
    d2rdrtsamp::Vector{T}
    state::ArrayPartition{T, NTuple{12, Vector{T}}}
    drstate::ArrayPartition{T, NTuple{12, Vector{T}}}
    dr2state::ArrayPartition{T, NTuple{12, Vector{T}}}
    dtstate::ArrayPartition{T, NTuple{12, Vector{T}}}
    dissipation::ArrayPartition{T, NTuple{12, Vector{T}}}
    temp::ArrayPartition{T, NTuple{12, Vector{T}}}
end

numvar = 12

@inline function Base.similar(::Type{ArrayPartition},::Type{T},size::Int) where T
    return ArrayPartition(
        similar(Vector{T}(undef,size)),similar(Vector{T}(undef,size)),similar(Vector{T}(undef,size)),
        similar(Vector{T}(undef,size)),similar(Vector{T}(undef,size)),similar(Vector{T}(undef,size)),
        similar(Vector{T}(undef,size)),similar(Vector{T}(undef,size)),similar(Vector{T}(undef,size)),
        similar(Vector{T}(undef,size)),similar(Vector{T}(undef,size)),similar(Vector{T}(undef,size))
    )
end

function printlogo()
    # Just a fancy ASCII logo for the program
    println(
"\n",
"           /\\\\\\\\\\\\\\\\\\       /\\\\\\\\\\\\\\\\\\\\\\\\\n",
" ________/\\\\\\\\\\\\\\\\\\\\\\\\\\\\___/ /\\\\\\\\/ / /\\\\\\\\_____________________________\n",
"  ______//\\\\\\\\\\/ / / / /___\\/ /\\\\\\\\/_/ / /\\\\\\\\___________________________\n",
"   _____//\\\\\\\\/_/_/_/_/_____\\/ /\\\\\\\\ \\/ / /\\\\\\\\___________________________\n",
"    ____/ /\\\\\\\\____/\\\\\\\\\\\\___\\/ /\\\\\\\\\\\\\\\\\\\\\\\\\\\\_______in___________________\n",
"     ___\\/ /\\\\\\\\__/ /  /\\\\\\___\\/ /\\\\\\\\\\\\\\\\\\\\\\\\/________Spherical____________\n",
"      ___\\/ /\\\\\\\\_\\/__/ /\\\\\\___\\/ /\\\\\\\\ / /\\\\\\\\__________Symmetry____________\n",
"       ___\\/ /\\\\\\\\____\\/ /\\\\\\___\\/ /\\\\\\\\_/ //\\\\\\\\_____________________________\n",
"        ___\\/ /\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\___\\/ /\\\\\\\\\\// //\\\\\\\\____________________________\n",
"         ___\\/ / /\\\\\\\\\\\\\\\\\\  /____\\/ /\\\\\\\\_\\// //\\\\\\\\_____by Conner Dailey_____\n",
"          ___\\/_/ / / / / /_/______\\/ /  /___\\// /  /______________________________\n",
"           _____\\/_/_/_/_/__________\\/__/______\\/__/______________________________\n"
)

end

function sample!(f::Vector{T}, grid::Grid{S}, fun) where {S,T}

    f .= T[fun(location(grid, n)) for n in 1:(grid.ncells + 4)]

end

function init!(state::ArrayPartition, param) where T

    ############################################
    # Specifies the Initial Conditions
    ############################################

    αreg,A,βr,Br,χ,γtrrreg,γtθθreg,Arrreg,Kreg,Γr,𝜙,K𝜙 = state.x

    grid = param.grid
    drt = spacing(grid)
    r = param.r
    drdrt = param.drdrt
    d2rdrt = param.d2rdrt
    rtmin = param.rtmin
    rtmax = param.rtmax

    n = grid.ncells + 4
    m = 1.
    rtspan = (rtmin,rtmax)

    # Initial conditions for Schwarzschild metric (Kerr-Schild Coordinates)

    # Mass (no real reason not to use 1 here)
    #M = 1

    fα(M,rt) = real((1+2*M/(r(rt))+0im)^(-1/2))
    fA(rt) = 0.
    fβr(M,rt) = (2*M/r(rt))*fα(M,rt)^2
    fBr(rt) = 0.
    fχ(rt) = 1.
    fγtrr(M,rt) = 1+2*M/r(rt)
    fγtθθ(rt) = r(rt)^2
    fArr(M,∂M,rt) = (4/3)*(r(rt)*(M+r(rt))*∂M-M*(3*M+2*r(rt)))/real(((r(rt)^5)*(r(rt)+2*M)+0im)^(1/2))
    fK(M,∂M,rt) = (2*M*(3*M+r(rt))+2*r(rt)*∂M*(M+r(rt)))/real((r(rt)*(r(rt)+2*M)+0im)^(3/2))
    fΓr(M,∂M,rt) = (r(rt)*∂M-2*r(rt)-5*M)/(r(rt)+2*M)^2

    r0 = 10.
    σr = 0.5
    #Amp = 1.
    Amp = 0.

    f𝜙(rt) = Amp*(1/r(rt))*exp(-(1/2)*((r(rt)-r0)/σr)^2)
    f∂𝜙(rt) = Amp*exp(-(1/2)*((r(rt)-r0)/σr)^2)*(r(rt)*r0-r(rt)^2-σr^2)/(r(rt)^2*σr^2)
    fK𝜙(rt) = 0.

    fρ(M,rt) = (2*fK𝜙(rt)^2 + (1/2)*(fχ(rt)/fγtrr(M,rt))*f∂𝜙(rt)^2
        + (1/2)*m^2*f𝜙(rt)^2)

    function f∂rtM(M,rt)
         if r(rt) < 2*M
             return 0.
         else
             4*pi*(r(rt)^2)*fρ(M,rt)*drdrt(rt)
         end
    end

    function f∂M(M,rt)
         if r(rt) < 2*M
             return 0.
         else
             4*pi*(r(rt)^2)*fρ(M,rt)
         end
    end

    fαreg(M,rt) = 2*M
    fγtrrreg(M,rt) = 2*M
    fArrreg(M,∂M,rt) = sqrt(r(rt)^5)*fArr(M,∂M,rt)
    fKreg(M,∂M,rt) = sqrt(r(rt)^3)*fK(M,∂M,rt)

    # Constraint Equations

    function constraintSystem(M, param, rt)
        f∂rtM(M,rt)
    end

    atol = 1e-15

    BVP = ODEProblem(constraintSystem, 1., rtspan, param)
    M = solve(BVP, Tsit5(), abstol=atol, dt=drt, adaptive=false)

    ∂M(rt) = f∂M(M(rt),rt)

    sample!(αreg, grid, rt -> fαreg(M(rt),rt) )
    sample!(A, grid, fA)
    sample!(βr, grid, rt -> fβr(M(rt),rt) )
    sample!(Br, grid, fBr)
    sample!(χ, grid, fχ)
    sample!(γtrrreg, grid, rt -> fγtrrreg(M(rt),rt) )
    sample!(γtθθreg, grid, rt -> 0)
    sample!(Arrreg, grid, rt -> fArrreg(M(rt),∂M(rt),rt) )
    sample!(Kreg, grid, rt -> fKreg(M(rt),∂M(rt),rt) )
    sample!(Γr, grid, rt -> fΓr(M(rt),∂M(rt),rt))
    sample!(𝜙, grid, f𝜙)
    sample!(K𝜙, grid, fK𝜙)

end

@inline function deriv!(df::Vector{T}, f::Vector{T}, n::Int64, dx::T) where T

    # @inbounds @fastmath @simd

    df[1] = T((-25. *f[1] + 48. *f[2] - 36. *f[3] + 16. *f[4] - 3. *f[5])/(12. *dx))

    df[2] = T((-3. *f[1] - 10. *f[2] + 18. *f[3] - 6. *f[4] + f[5])/(12. *dx))

    for i in 3:(n - 2)
        df[i] = T(((f[i-2] - 8. *f[i-1] + 8. *f[i+1] - f[i+2])/(12. *dx)))
    end

    df[n-1] = T(-(-3. *f[n] - 10. *f[n-1] + 18. *f[n-2] - 6. *f[n-3] + f[n-4])/(12. *dx))

    df[n] = T(-(-25. *f[n] + 48. *f[n-1] - 36. *f[n-2] + 16. *f[n-3] - 3. *f[n-4])/(12. *dx))

end

@inline function dissipation!(df::Vector{T}, f::Vector{T},drdrt::Vector{T}, n::Int64) where T

    ############################################
    # Calculates the numerical dissipation terms
    #
    # Finite differencing fails to model high
    # frequency modes in the system, and these
    # modes can lead to instabilities. Adding
    # numerical dissipation terms damps these
    # high frequency modes by subtracting off
    # a high order derivative from each dynamical
    # variable in the system.
    ############################################

    # @simd @inbounds @fastmath

    df[1:2] .= 0.

    for i in 3:(n - 2)
        df[i] = T((f[i-2] - 4. *f[i-1] + 6. *f[i] - 4. *f[i+1] + f[i+2])/(drdrt[i]))
    end

    df[(n-1):n] .= 0.

end

function rhs_test(dtstate::ArrayPartition{T, NTuple{12, Vector{T}}},regstate::ArrayPartition{T, NTuple{12, Vector{T}}}, param::Param{T}, t) where T

    ############################################
    # Caculates the right hand ride of the
    # evolved variables
    #
    # This is the main meat of the program.
    # This function contains all of the boundary
    # conditions, coordinate conversions,
    # spatial derivative calculations,
    # evolution equations, and numerical
    # dissipation. Each time the Julia DiffEq
    # Solver moves one time step, it calls
    # this function to calculate the new
    # values of the evolved variables.
    ############################################

    # Unpack the parameters

    m = 1.
    Mtot = 1.

    grid = param.grid
    drt = param.drt
    r = param.rsamp
    drdrt = param.drdrtsamp
    d2rdrt = param.d2rdrtsamp
    rtmin = param.rtmin
    rtmax = param.rtmax

    state = param.state
    drstate = param.drstate
    dr2state = param.dr2state
    dissipation = param.dissipation
    dtstate2 = param.dtstate
    temp = param.temp

    n = grid.ncells + 4

    # Copy the state into the parameters
    # so that it can be changed

    #######################
    # Attention!
    #
    # Do not do the following:
    # state .= regstate
    #
    # This results in an intense slowdown
    # Do instead:
    for i in 1:12
        state.x[i] .= regstate.x[i]
    end

    # Give names to individual variables

    α,A,βr,Br,χ,γtrr,γtθθ,Arr,K,Γr,𝜙,K𝜙 = state.x
    ∂α,∂A,∂βr,∂Br,∂χ,∂γtrr,∂γtθθ,∂Arr,∂K,∂Γr,∂𝜙,∂K𝜙 = drstate.x
    ∂2α,∂2A,∂2βr,∂2Br,∂2χ,∂2γtrr,∂2γtθθ,∂2Arr,∂2K,∂2Γr,∂2𝜙,∂2K𝜙 = dr2state.x
    ∂tα,∂tA,∂tβr,∂tBr,∂tχ,∂tγtrr,∂tγtθθ,∂tArr,∂tK,∂tΓr,∂t𝜙,∂tK𝜙 = dtstate.x
    ∂4α,∂4A,∂4βr,∂4Br,∂4χ,∂4γtrr,∂4γtθθ,∂4Arr,∂4K,∂4Γr,∂4𝜙,∂4K𝜙 = dissipation.x

    # Calculate first derivatives

    deriv!(∂α,α,n,drt)
    deriv!(∂A,A,n,drt)
    deriv!(∂βr,βr,n,drt)
    deriv!(∂Br,Br,n,drt)
    deriv!(∂χ,χ,n,drt)
    deriv!(∂γtrr,γtrr,n,drt)
    deriv!(∂γtθθ,γtθθ,n,drt)
    deriv!(∂Arr,Arr,n,drt)
    deriv!(∂K,K,n,drt)
    deriv!(∂Γr,Γr,n,drt)
    deriv!(∂𝜙,𝜙,n,drt)
    deriv!(∂K𝜙,K𝜙,n,drt)

    # Calculate second derivatives

    deriv!(∂2α,∂α,n,drt)
    deriv!(∂2βr,∂βr,n,drt)
    deriv!(∂2χ,∂χ,n,drt)
    deriv!(∂2γtrr,∂γtrr,n,drt)
    deriv!(∂2γtθθ,∂γtθθ,n,drt)
    deriv!(∂2𝜙,∂𝜙,n,drt)

    # Convert between computational rt coordinnate
    # and trasditional r coordinate

    ∂α ./= drdrt
    ∂A ./= drdrt
    ∂βr ./= drdrt
    ∂Br ./= drdrt
    ∂χ ./= drdrt
    ∂γtrr ./= drdrt
    ∂γtθθ ./= drdrt
    ∂Arr ./= drdrt
    ∂K ./= drdrt
    ∂Γr ./= drdrt
    ∂𝜙 ./= drdrt
    ∂K𝜙 ./= drdrt

    @. ∂2α = (∂2α - d2rdrt*∂α)/(drdrt^2)
    @. ∂2βr = (∂2βr - d2rdrt*∂βr)/(drdrt^2)
    @. ∂2χ = (∂2χ - d2rdrt*∂χ)/(drdrt^2)
    @. ∂2γtrr = (∂2γtrr - d2rdrt*∂γtrr)/(drdrt^2)
    @. ∂2γtθθ = (∂2γtθθ - d2rdrt*∂γtθθ)/(drdrt^2)
    @. ∂2𝜙 = (∂2𝜙 - d2rdrt*∂𝜙)/(drdrt^2)

    # Convert between regularized variables
    # and cannonical variables

    @. α = real((1 + α/r + 0im)^(-1/2))
    @. ∂α = α*(1 - (1 + ∂α)*α^2)/(2*r)
    @. ∂2α = (6*r*∂α^2 - 4*α*∂α - ∂2α*α^4)/(2*r*α)

    @. γtrr = γtrr/r + 1
    @. ∂γtrr = (1 - γtrr + ∂γtrr)/r
    @. ∂2γtrr = (∂2γtrr - 2*∂γtrr)/r

    @. γtθθ = (r^2)*(γtθθ + 1)
    @. ∂γtθθ = (2*γtθθ + ∂γtθθ*r^3)/r
    @. ∂2γtθθ = (4*∂γtθθ*r - 6*γtθθ + ∂2γtθθ*r^4)/(r^2)

    @. K = sqrt(r^(-3))*K
    @. ∂K = sqrt(r^(-3))*∂K - (3/2)*K/r

    @. Arr = sqrt(r^(-5))*Arr
    @. ∂Arr = sqrt(r^(-5))*∂Arr - (5/2)*Arr/r

    #########################################################
    # Evolution Equations
    #
    # This is the full suite of evolution equations
    # for GR in spherical symmetry in the BSSN framework.
    # I have tried to keep them looking as close to their
    # mathematically written counterpart as possible.
    #
    # They are written in the order they appear in the
    # reference (arXiv:0705.3845v2) except for the ∂tBr
    # equation since it contains a ∂tΓr term.
    #
    #########################################################

    # Lagrangian Gauge condition
    v = 1.

    @. ∂tχ = ((2/3)*K*α*χ - (1/3)*v*βr*χ*∂γtrr/γtrr - (2/3)*v*βr*χ*∂γtθθ/γtθθ
     - (2/3)*v*χ*∂βr + βr*∂χ)

    @. ∂tγtrr = (-2*Arr*α - (1/3)*v*βr*∂γtrr + βr*∂γtrr
     - (2/3)*v*γtrr*βr*∂γtθθ/γtθθ + 2*γtrr*∂βr - (2/3)*v*γtrr*∂βr)

    @. ∂tγtθθ = (Arr*γtθθ*α/γtrr - (1/3)*v*γtθθ*βr*∂γtrr/γtrr - (2/3)*v*βr*∂γtθθ
     + βr*∂γtθθ - (2/3)*v*γtθθ*∂βr)

    @. ∂tArr = (-2*α*(Arr^2)/γtrr + K*α*Arr - (1/3)*v*βr*Arr*∂γtrr/γtrr
     - (2/3)*v*βr*Arr*∂γtθθ/γtθθ - (2/3)*v*Arr*∂βr + 2*Arr*∂βr
     + (2/3)*α*χ*(∂γtrr/γtrr)^2 - (1/3)*α*χ*(∂γtθθ/γtθθ)^2
     - (1/6)*α*(∂χ^2)/χ - (2/3)*α*χ*γtrr/γtθθ + βr*∂Arr
     + (2/3)*α*χ*γtrr*∂Γr - (1/2)*α*χ*(∂γtrr/γtrr)*(∂γtθθ/γtθθ)
     + (1/3)*χ*∂γtrr*∂α/γtrr + (1/3)*χ*∂α*∂γtθθ/γtθθ - (1/6)*α*∂γtrr*∂χ/γtrr
     - (1/6)*α*∂γtθθ*∂χ/γtθθ - (2/3)*∂α*∂χ - (1/3)*α*χ*∂2γtrr/γtrr
     + (1/3)*α*χ*∂2γtθθ/γtθθ - (2/3)*χ*∂2α + (1/3)*α*∂2χ)

    @. ∂tK = ((3/2)*α*(Arr/γtrr)^2 + (1/3)*α*K^2 + βr*∂K
     + (1/2)*χ*∂γtrr*∂α/(γtrr^2) - χ*∂α*(∂γtθθ/γtθθ)/γtrr
     + (1/2)*∂α*∂χ/γtrr - χ*∂2α/γtrr)

    @. ∂tΓr = (-v*βr*((∂γtθθ/γtθθ)^2)/γtrr + α*Arr*(∂γtθθ/γtθθ)/(γtrr^2)
     - (1/3)*v*∂βr*(∂γtθθ/γtθθ)/γtrr + ∂βr*(∂γtθθ/γtθθ)/γtrr
     + βr*∂Γr + α*Arr*∂γtrr/(γtrr^3) - (4/3)*α*∂K/γtrr
     - 2*Arr*∂α/(γtrr^2) + (1/2)*v*∂βr*∂γtrr/(γtrr^2)
     - (1/2)*∂βr*∂γtrr/(γtrr^2) - 3*α*Arr*(∂χ/χ)/(γtrr^2)
     + (1/6)*v*βr*∂2γtrr/(γtrr^2) + (1/3)*v*βr*(∂2γtθθ/γtθθ)/γtrr
     + (1/3)*v*∂2βr/γtrr + ∂2βr/γtrr)

    #######################################################
    # Gauge Evolution

    @. ∂tα = -2*α*A
    @. ∂tA = ∂tK

    @. ∂tβr = (3/4)*Br
    @. ∂tBr = ∂tΓr

    # Gauge choices for the evolution of the
    # determinant of the conformal metric
    # (must have v = 1 to use this)

    # ∂tlnγt = temp.x[5]
    # ∂rt∂tlnγt = temp.x[6]
    #
    # ∂tlnγt .= 0
    #
    # #∂tlnγt = -8*pi*Sr.*real((γtθθ./γtrr .+ 0im).^(1/2))
    #
    # deriv!(∂rt∂tlnγt,∂tlnγt,n,drt)
    #
    # @. ∂r∂tlnγt = ∂rt∂tlnγt/drdrt
    #
    # # ∂tα = -(1/2)*α.*∂tlnγt
    # # ∂tβr = (χ./γtrr).*∂tlnγt
    #
    # @. ∂tχ += (1/3)*χ*∂tlnγt
    # @. ∂tγtrr += (1/3)*γtrr*∂tlnγt
    # @. ∂tγtθθ += (1/3)*γtθθ*∂tlnγt
    # @. ∂tArr += (1/3)*Arr*∂tlnγt
    # @. ∂tΓr += -(1/3)*Γr*∂tlnγt - (1/6)*(χ/γtrr)*∂r∂tlnγt

    #########################################################
    # Source Terms and Source Evolution
    #
    # This currently includes the addition of source terms
    # to GR that come from a Klein-Gordon scalar field
    #
    #########################################################

    # Klein-Gordon System

    @. ∂t𝜙 = βr*∂𝜙 - 2*α*K𝜙

    @. ∂tK𝜙 = (βr*∂K𝜙 + α*K*K𝜙 + (1/2)*m^2*α*𝜙 - (1/2)*χ.*∂α.*∂𝜙./γtrr
        - (1/2)*(α*χ/γtrr)*(∂2𝜙 + ∂𝜙*(∂γtθθ/γtθθ - (1/2)*∂γtrr/γtrr
        - (1/2)*(∂χ/χ))))

    ρ = temp.x[1]
    Sr = temp.x[2]
    S = temp.x[3]
    Srr = temp.x[4]

    @. ρ = 2*K𝜙^2 + (1/2)*(χ/γtrr)*∂𝜙^2 + (1/2)*(m^2)*𝜙^2
    #Lower Index
    @. Sr = 2*K𝜙*∂𝜙
    @. S = 6*K𝜙^2 - (1/2)*(χ/γtrr)*∂𝜙^2 - (3/2)*(m^2)*𝜙^2
    @. Srr = (γtrr/χ)*(2*K𝜙^2 + (1/2)*(χ/γtrr)*∂𝜙^2 - (1/2)*(m^2)*𝜙^2)

    @. ∂tArr += -8*pi*α*(χ*Srr - (1/3)*S*γtrr)
    @. ∂tK += 4*pi*α*(ρ + S)
    @. ∂tΓr += -16*pi*α*Sr/γtrr

    # fr = param.r
    #
    # fα(M,rt) = real((1+2*M/(fr(rt))+0im)^(-1/2))
    # fβr(M,rt) = (2*M/fr(rt))*fα(M,rt)^2
    # fγtrr(M,rt) = 1+2*M/fr(rt)
    # fγtθθ(rt) = fr(rt)^2
    # fArr(M,∂M,rt) = (4/3)*(fr(rt)*(M+fr(rt))*∂M-M*(3*M+2*fr(rt)))/real(((fr(rt)^5)*(fr(rt)+2*M)+0im)^(1/2))
    # fK(M,∂M,rt) = (2*M*(3*M+fr(rt))+2*fr(rt)*∂M*(M+fr(rt)))/real((fr(rt)*(fr(rt)+2*M)+0im)^(3/2))
    # fΓr(M,∂M,rt) = (fr(rt)*∂M-2*fr(rt)-5*M)/(fr(rt)+2*M)^2
    #
    # f∂α(M,rt) = M*real((fr(rt)*(fr(rt)+2*M+0im)^3)^(-1/2))
    # f∂βr(M,rt) = -2*M/(fr(rt)+2*M)^2
    # f∂γtrr(M,rt) = -2*M/(fr(rt)^2)
    # f∂γtθθ(rt) = 2*fr(rt)
    # f∂Arr(M,rt) = (4*M/3)*(15*M^2+15*M*fr(rt)+4*fr(rt)^2)/real(((fr(rt)^7)*((fr(rt)+2*M)^3)+0im)^(1/2))
    # f∂K(M,rt) = -2*M*(9*M^2+10*M*fr(rt)+2*fr(rt)^2)/real((fr(rt)*(fr(rt)+2*M)+0im)^(5/2))
    # f∂Γr(M,rt) = 2*(fr(rt)+3*M)/(fr(rt)+2*M)^3
    #
    # rt = sample(Float64, A.grid, rt->rt)
    #
    # # for i=1:2
    #
    # ∂tα[1:2] .= (α[1:2] .- fα.(1.,rt[1:2]))./r[1:2] + ∂α[1:2] - f∂α.(1.,rt[1:2])
    # ∂tA[1:2] .= (A[1:2] .- 0.)./r[1:2] + ∂A[1:2]
    # ∂tβr[1:2] .= (βr[1:2] .- fβr.(1.,rt[1:2]))./r[1:2] + ∂βr[1:2] - f∂βr.(1.,rt[1:2])
    # ∂tBr[1:2] .= (Br[1:2] .- 0.)./r[1:2] + ∂Br[1:2]
    # ∂tχ[1:2] .= (χ[1:2] .- 1.)./r[1:2] + ∂χ[1:2]
    # ∂tγtrr[1:2] .= (γtrr[1:2] .- fγtrr.(1.,rt[1:2]))./r[1:2] + ∂γtrr[1:2] - f∂γtrr.(1.,rt[1:2])
    # ∂tγtθθ[1:2] .= (γtθθ[1:2] .- fγtθθ.(rt[1:2]))./r[1:2] + ∂γtθθ[1:2] - f∂γtθθ.(rt[1:2])
    # ∂tArr[1:2] .= (Arr[1:2] .- fArr.(1.,0.,rt[1:2]))./r[1:2] + ∂Arr[1:2] - f∂Arr.(1.,rt[1:2])
    # ∂tK[1:2] .= (K[1:2] .- fK.(1.,0.,rt[1:2]))./r[1:2] + ∂K[1:2] - f∂K.(1.,rt[1:2])
    # ∂tΓr[1:2] .= (Γr[1:2] .- fΓr.(1.,0.,rt[1:2]))./r[1:2] + ∂Γr[1:2] - f∂Γr.(1.,rt[1:2])
    # ∂t𝜙[1:2] .= (𝜙[1:2] .- 0.)./r[1:2] + ∂𝜙[1:2]
    # ∂tK𝜙[1:2] .= (K𝜙[1:2] .- 0.)./r[1:2] + ∂K𝜙[1:2]

    ######################################################
    # Calculate the numerical dissipation

    dissipation!(∂4α,α,drdrt,n)
    dissipation!(∂4A,A,drdrt,n)
    dissipation!(∂4βr,βr,drdrt,n)
    dissipation!(∂4Br,Br,drdrt,n)
    dissipation!(∂4χ,χ,drdrt,n)
    dissipation!(∂4γtrr,γtrr,drdrt,n)
    dissipation!(∂4γtθθ,γtθθ,drdrt,n)
    dissipation!(∂4Arr,Arr,drdrt,n)
    dissipation!(∂4K,K,drdrt,n)
    dissipation!(∂4Γr,Γr,drdrt,n)
    dissipation!(∂4𝜙,𝜙,drdrt,n)
    dissipation!(∂4K𝜙,K𝜙,drdrt,n)

    # Magnitude of dissipation
    σ = 0.3

    @. ∂tα -= (1/16)*σ*∂4α/drt
    @. ∂tA -= (1/16)*σ*∂4A/drt
    @. ∂tβr -= (1/16)*σ*∂4βr/drt
    @. ∂tBr -= (1/16)*σ*∂4Br/drt
    @. ∂tχ -= (1/16)*σ*∂4χ/drt
    @. ∂tγtrr -= (1/16)*σ*∂4γtrr/drt
    @. ∂tγtθθ -= (1/16)*σ*∂4γtθθ/drt
    @. ∂tArr -= (1/16)*σ*∂4Arr/drt
    @. ∂tK -= (1/16)*σ*∂4K/drt
    @. ∂tΓr -= (1/16)*σ*∂4Γr/drt
    @. ∂t𝜙 -= (1/16)*σ*∂4𝜙/drt
    @. ∂tK𝜙 -= (1/16)*σ*∂4K𝜙/drt

    # Convert back to regularized variables
    # for the time derivatives

    @. ∂tα = -2*r*∂tα/(α^3)
    @. ∂tγtrr = r*∂tγtrr
    @. ∂tγtθθ = ∂tγtθθ/(r^2)
    @. ∂tArr = sqrt(r^5)*∂tArr
    @. ∂tK = sqrt(r^3)*∂tK

    # Specify the inner temporal boundary conditions

    for i in 1:12
        dtstate.x[i][1:2] .= 0.
    end

    # Specify the outer temporal boundary conditions

    for i in 1:12
        dtstate.x[i][(n-1):n] .= 0.
    end

    # Store the calculated state into the param
    # so that we can print it to the screen

    for i in 1:12
        dtstate2.x[i] .= dtstate.x[i]
    end

end

function rhs_all(regstate::ArrayPartition{Float64, NTuple{12, Vector{Float64}}}, param::Param{Float64}, t)

    n = param.grid.ncells + 4

    dtstate = similar(ArrayPartition,Float64,n)

    rhs_test(dtstate,regstate,param,t)

    return dtstate

end

function constraints(state::ArrayPartition,drstate::ArrayPartition,dr2state::ArrayPartition,param)

    ############################################
    # Caculates the constraints of the system
    #
    # These are outputed at every saved time step.
    # They should limit to zero as the spatial
    # resolution is increased. They are important to
    # monitor to make sure the physics of the system
    # is being properly modeled.
    ############################################

    # Unpack Variables

    α,A,βr,Br,χ,γtrr,γtθθ,Arr,K,Γr,𝜙,K𝜙 = state.x

    ∂α,∂A,∂βr,∂Br,∂χ,∂γtrr,∂γtθθ,∂Arr,∂K,∂Γr,∂𝜙,∂K𝜙 = drstate.x

    ∂2α,∂2A,∂2βr,∂2Br,∂2χ,∂2γtrr,∂2γtθθ,∂2Arr,∂2K,∂2Γr,∂2𝜙,∂2K𝜙 = dr2state.x

    m = 1.
    n = param.grid.ncells + 4
    drt = param.drt
    r = param.rsamp
    drdrt = param.drdrtsamp
    d2rdrt = param.d2rdrtsamp

    deriv!(∂χ,χ,n,drt)
    deriv!(∂γtrr,γtrr,n,drt)
    deriv!(∂γtθθ,γtθθ,n,drt)
    deriv!(∂Arr,Arr,n,drt)
    deriv!(∂K,K,n,drt)
    deriv!(∂Γr,Γr,n,drt)
    deriv!(∂𝜙,𝜙,n,drt)
    deriv!(∂K𝜙,K𝜙,n,drt)

    deriv!(∂2χ,∂χ,n,drt)
    deriv!(∂2γtrr,∂γtrr,n,drt)
    deriv!(∂2γtθθ,∂γtθθ,n,drt)
    deriv!(∂2𝜙,∂𝜙,n,drt)

    ∂χ ./= drdrt
    ∂γtrr ./= drdrt
    ∂γtθθ ./= drdrt
    ∂Arr ./= drdrt
    ∂K ./= drdrt
    ∂Γr ./= drdrt
    ∂𝜙 ./= drdrt
    ∂K𝜙 ./= drdrt

    @. ∂2χ = (∂2χ - d2rdrt*∂χ)/(drdrt^2)
    @. ∂2γtrr = (∂2γtrr - d2rdrt*∂γtrr)/(drdrt^2)
    @. ∂2γtθθ = (∂2γtθθ - d2rdrt*∂γtθθ)/(drdrt^2)
    @. ∂2𝜙 = (∂2𝜙 - d2rdrt*∂𝜙)/(drdrt^2)

    @. γtrr = γtrr/r + 1
    @. ∂γtrr = (1 - γtrr + ∂γtrr)/r
    @. ∂2γtrr = (∂2γtrr - 2*∂γtrr)/r

    @. γtθθ = (r^2)*(γtθθ + 1)
    @. ∂γtθθ = (2*γtθθ + ∂γtθθ*r^3)/r
    @. ∂2γtθθ = (4*∂γtθθ*r - 6*γtθθ + ∂2γtθθ*r^4)/(r^2)

    @. K = sqrt(r^(-3))*K
    @. ∂K = sqrt(r^(-3))*∂K - (3/2)*K/r

    @. Arr = sqrt(r^(-5))*Arr
    @. ∂Arr = sqrt(r^(-5))*∂Arr - (5/2)*Arr/r

    ρ = 2*K𝜙.^2 + (1/2)*(χ./γtrr).*∂𝜙.^2 + (1/2)*m^2*𝜙.^2
    Sr = 2*K𝜙.*∂𝜙

    # Constraint Equations

    𝓗 = (-(3/2)*(Arr./γtrr).^2 + (2/3)*K.^2 - (5/2)*((∂χ.^2)./χ)./γtrr
     + 2*∂2χ./γtrr + 2*χ./γtθθ - 2*χ.*(∂2γtθθ./γtθθ)./γtrr + 2*∂χ.*(∂γtθθ./γtθθ)./γtrr
     + χ.*(∂γtrr./(γtrr.^2)).*(∂γtθθ./γtθθ) - ∂χ.*∂γtrr./(γtrr.^2)
     + (1/2)*χ.*((∂γtθθ./γtθθ).^2)./γtrr - 16*pi*ρ)

    𝓜r = (∂Arr./γtrr - (2/3)*∂K - (3/2)*Arr.*(∂χ./χ)./γtrr
     + (3/2)*Arr.*(∂γtθθ./γtθθ)./γtrr - Arr.*∂γtrr./(γtrr.^2)
     - 8*pi*Sr)

    𝓖r = -(1/2)*∂γtrr./(γtrr.^2) + Γr + (∂γtθθ./γtθθ)./γtrr

    return (𝓗, 𝓜r, 𝓖r, ρ)

end

function horizon(T,state::ArrayPartition,param)

    ############################################
    # Caculates the apparent horizon
    #
    # Where the function crosses zero is the
    # apparent horizon of the black hole.
    ############################################

    #v = param[3]

    # Unpack Variables

    χ = state.χ
    γtrr = state.γtrr
    γtθθ = state.γtθθ
    Arr = state.Arr
    K = state.K

    # Gauge condition

    # if v == 1 # Lagrangian Condition
    #     γtθθreg[2] = ((-315*γtθθreg[3] + 210*γtθθreg[4] - 126*γtθθreg[5]
    #     + 45*γtθθreg[6] - 7*γtθθreg[7])/63)
    # end

    #r = sample(T, χ.grid, param[4])
    drdrt = param.drdrtsamp

    # Conversions from regularized variables to canonical variables

    # γtrr = γtrrreg./r .+ 1
    #
    # γtθθ = (r.^2).*(γtθθreg .+ 1)
    #
    # K = real((r .+ 0im).^(-3/2)).*Kreg
    #
    # Arr = real((r .+ 0im).^(-5/2)).*Arrreg

    # Intermediate calculations

    Kθθ = ((1/3)*γtθθ.*K - (1/2)*Arr.*γtθθ./γtrr)./χ

    grr =  γtrr./χ

    gθθ =  γtθθ./χ

    # Spatial Derivatives

    ∂rtgθθ = deriv(gθθ,4,1)

    # Coordinate transformations from computational rt coordinate
    # to physical r coordinate

    ∂gθθ = ∂rtgθθ./drdrt

    # Apparent horizon function

    Θ = (∂gθθ./gθθ)./real((grr .+ 0im).^(1/2)) - 2*Kθθ./gθθ

    # cross = GridFun(χ.grid, sign.(Θ))

    return Θ

end

# function crossings(Θ::GridFun)
#
#     GridFun(χ.grid, sign.(Θ))
#
#     for
#
# end

function custom_progress_message(dt,state::ArrayPartition{T, NTuple{numvar, Vector{T}}},param,t) where T

    ###############################################
    # Outputs status numbers while the program runs
    ###############################################

    dtstate = param.dtstate::ArrayPartition{T, NTuple{numvar, Vector{T}}}

    ∂tα,∂tA,∂tβr,∂tBr,∂tχ,∂tγtrr,∂tγtθθ,∂tArr,∂tK,∂tΓr,∂t𝜙,∂tK𝜙 = dtstate.x

    println("  ",
    #rpad(string(param[1]),6," "),
    #rpad(string(round(dt,digits=3)),10," "),
    rpad(string(round(t,digits=3)),10," "),
    rpad(string(round(maximum(abs.(∂tα)),digits=3)),12," "),
    rpad(string(round(maximum(abs.(∂tχ)),digits=3)),12," "),
    rpad(string(round(maximum(abs.(∂tγtrr)),digits=3)),14," "),
    rpad(string(round(maximum(abs.(∂tγtθθ)),digits=3)),14," "),
    rpad(string(round(maximum(abs.(∂tArr)),digits=3)),12," "),
    # rpad(string(round(maximum(abs.(derivstate.𝜙)),digits=3)),12," "),
    # rpad(string(round(maximum(abs.(derivstate.K𝜙)),digits=3)),14," ")
    rpad(string(round(maximum(abs.(∂tK)),digits=3)),12," "),
    rpad(string(round(maximum(abs.(∂tΓr)),digits=3)),14," ")
    )

    #PrettyTables.jl

    # param[1] += param[2]

    return

end

function solution_saver(T,grid,sol,param,folder)

    ###############################################
    # Saves all of the variables in nice CSV files
    # in the choosen data folder directory
    ###############################################

    vars = (["α","A","βr","Br","χ","γtrr","γtθθ","Arr","K","Γr","𝜙","K𝜙",
    "∂tα","∂tA","∂tβr","∂tBr","∂tχ","∂tγtrr","∂tγtθθ","∂tArr","∂tK","∂tΓr","∂t𝜙","∂tK𝜙",
    "H","Mr","Gr","ρ","appHorizon"])
    varlen = length(vars)
    #mkdir(string("data\\",folder))
    tlen = size(sol)[2]
    rlen = grid.ncells + 4
    r = param.rsamp
    #cons = Array{GridFun,2}(undef,tlen,4)
    #state = Array{ArrayPartition,1}(undef,tlen)
    #dstate = Array{ArrayPartition,1}(undef,tlen)
    #d2state = Array{ArrayPartition,1}(undef,tlen)
    #derivs = Array{ArrayPartition,1}(undef,tlen)
    #apphorizon = Array{GridFun,1}(undef,tlen)

    # for i in 1:tlen
    #
    #     αreg,A,βr,Br,χ,γtrrreg,γtθθreg,Arrreg,Kreg,Γr,𝜙,K𝜙 = sol[i].x
    #
    #     # Conversions from regularized variables to canonical variables
    #     α = real((1 .+ αreg./r .+ 0im).^(-1/2))
    #     γtrr = γtrrreg./r .+ 1
    #     γtθθ = (r.^2).*(γtθθreg .+ 1)
    #     K = sqrt.(r.^(-3)).*Kreg
    #     Arr = sqrt.(r.^(-5)).*Arrreg
    #
    #     state[i] = ArrayPartition(α,A,βr,Br,χ,γtrr,γtθθ,Arr,K,Γr,𝜙,K𝜙)
    #     derivs[i] = param.dtstate
    #     rhs_test(derivs[i],sol[i],param,0)
    #     #cons[i,1:4] .= constraints!(T,state[i],dstate[i],d2state[i],param)
    #     #apphorizon[i] = horizon(T,state[i],param)
    # end

    drstate = param.drstate
    dr2state = param.dr2state

    dtstate = [rhs_all(sol[i],param,0.) for i = 1:tlen]

    cons = [constraints(sol[i],drstate,dr2state,param) for i = 1:tlen]

    array = Array{T,2}(undef,tlen+1,rlen+1)

    array[1,1] = 0
    array[1,2:end] .= r

    for j = 1:varlen
        if j < 13
            for i = 2:tlen+1
                array[i,1] = sol.t[i-1]
                array[i,2:end] .= sol[i-1].x[j]
                # if j==1
                #     println(array[i,2:10])
                # end
            end
        elseif j < 25

            for i = 2:tlen+1

                array[i,1] = sol.t[i-1]
                array[i,2:end] .= dtstate[i-1].x[j-12]
            end
        elseif j < 29
            for i = 2:tlen+1
                #println(size(cons[i-1][j-24]))
                array[i,1] = sol.t[i-1]
                array[i,2:end] .= cons[i-1][j-24]
            end
        # elseif j == varlen
        #     for i = 2:tlen+1
        #         array[i,1] = sol.t[i-1]
        #         array[i,2:end] .= apphorizon[i-1]
        #     end
        end

        CSV.write(
        string("data/",folder,"/",vars[j],".csv"),
        DataFrame(array, :auto),
        header=false
        )

    end

    # for i = 2:tlen+1
    #     array[i,1] = sol.t[i-1]
    #     array[i,2:end] .= sol[i-1].x[1][:]
    #     println(sol[i-1].x[1][1:10])
    # end
    #
    # CSV.write(
    # string("data/",folder,"/",vars[1],".csv"),
    # DataFrame(array, :auto),
    # header=false
    # )

end

function main(points)

    ###############################################
    # Main Program
    #
    # Calls each of the above functions to run a
    # simulation. Sets up the numerical grid,
    # sets the gauge conditions, sets up a new
    # computational rt coordinate that makes the
    # physical r coordinate step larger near the
    # outer boundary, sets the initial conditions,
    # and finally runs the numerical DiffEq
    # package to run the time integration.
    #
    # All data is saved in the folder specified to
    # the solution_saver, each in their own CSV
    # file.
    ###############################################

    T = Float64
    rspan = T[1.,210.]
    rtspan = T[1.,21.]

    rtmin, rtmax = rtspan

    domain = Domain{S}(rtmin, rtmax)
    grid = Grid(domain, points)

    n = grid.ncells + 4

    drt = spacing(grid)
    dt = drt/4.

    tspan = T[0.,20. ]
    tmin, tmax = tspan

    printtimes = 1.

    v = 1.

    m = 0.

    Mtot = 1.

    f(b) = b*tan(rtspan[2]/b)-rspan[2]

    scale = find_zero(f, 0.64*rtspan[2])

    r(rt) = scale*tan(rt/scale)
    drdrt(rt) = sec(rt/scale)^2
    d2rdrt(rt) = (2/scale)*(sec(rt/scale)^2)*tan(rt/scale)

    # r(rt) = rt
    # drdrt(rt) = 1
    # d2rdrt(rt) = 0

    atol = eps(T)^(T(3) / 4)

    #alg = KuttaPRK2p5()
    alg = RK4()
    #alg = ImplicitMidpoint()

    #printlogo()

    custom_progress_step = round(Int, printtimes/dt)
    step_iterator = custom_progress_step

    regstate = similar(ArrayPartition,T,n)
    state = similar(ArrayPartition,T,n)
    drstate = similar(ArrayPartition,T,n)
    dr2state = similar(ArrayPartition,T,n)
    dtstate = similar(ArrayPartition,T,n)
    dissipation = similar(ArrayPartition,T,n)
    temp = similar(ArrayPartition,T,n)

    # regstate = similar(Array{T}(undef,numvar,n))
    # state = similar(Array{T}(undef,numvar,n))
    # drstate = similar(Array{T}(undef,numvar,n))
    # dr2state = similar(Array{T}(undef,numvar,n))
    # dtstate = similar(Array{T}(undef,numvar,n))
    # dissipation = similar(Array{T}(undef,numvar,n))
    # temp = similar(Array{T}(undef,numvar,n))

    #println("Defining Problem...")
    rsamp = similar(Vector{T}(undef,n))
    drdrtsamp = similar(Vector{T}(undef,n))
    d2rdrtsamp = similar(Vector{T}(undef,n))

    sample!(rsamp, grid, r)
    sample!(drdrtsamp, grid, drdrt)
    sample!(d2rdrtsamp, grid, d2rdrt)

    param = Param(rtmin,rtmax,drt,Mtot,grid,r,drdrt,d2rdrt,rsamp,drdrtsamp,d2rdrtsamp,state,drstate,dr2state,dtstate,dissipation,temp)
    #
    init!(regstate, param)

    # @btime rhs!($regstate,$regstate,$param, 0.0)
    # # @btime radial_derivatives!($regstate,$regstate,$regstate, $param, 0.0)
    # # @btime regular_to_canonical!($regstate,$regstate,$regstate, $param, 0.0)
    # # # #
    # return

    prob = ODEProblem(rhs_test, regstate, tspan, param)

    println("Starting Solution...")

    #return

    println("")
    println("| Time | max α'(t) | max χ'(t) | max γtrr'(t) | max γtθθ'(t) | max Arr'(t) | max K'(t) | max Γr'(t) |")
    println("|______|___________|___________|______________|______________|_____________|___________|____________|")
    println("")

    sol = solve(
        prob, alg,
        abstol = atol,
        dt = drt/4,
        adaptive = false,
        saveat = printtimes,
        alias_u0 = true,
        progress = true,
        progress_steps = custom_progress_step,
        progress_message = custom_progress_message
    )

    #@code_native deriv!([1.0,1.0,1.0,1.0,1.0,1.0], [1.0,1.0,1.0,1.0,1.0,1.0], 6, 1.)

    #@btime α,A,βr,Br,χ,γtrr,γtθθ,Arr,K,Γr,𝜙,K𝜙 = $state

    #@btime deriv!($dtstate[1],$regstate[1], $n, 1.)

    # rhs_test(dtstate,regstate, param, 0)
    #
    # Profile.clear_malloc_data()
    #
    # rhs_test(dtstate,regstate, param, 0)

    #@btime rhs_test($dtstate,$regstate, $param, 1.)

    #@code_warntype rhs_test(dtstate,regstate, param, 1.)

    # @profile rhs_test(dtstate,regstate, param, 1.)
    #
    # Profile.print()

    #@btime α,A,βr,Br,χ,γtrr,γtθθ,Arr,K,Γr,𝜙,K𝜙 = state

    #@btime radial_derivatives!($drstate, $dr2state, $regstate, $param, 0)

    # df = similar(Vector{T}(undef,n))
    # f = similar(Vector{T}(undef,n))
    #
    # a = GridFun(grid,Vector{T}(zeros(n)))
    # b = GridFun(grid,Vector{T}(zeros(n)))
    #
    # test = similar(VectorOfArray,T,n)
    # )

    #return

    solution_saver(T,grid,sol,param,"test")

    return

end

end

# k1 = zero(GBSSN_Variables{T,T}, grid)::GBSSN_Variables
# k2 = zero(GBSSN_Variables{T,T}, grid)::GBSSN_Variables
# k3 = zero(GBSSN_Variables{T,T}, grid)::GBSSN_Variables
# k4 = zero(GBSSN_Variables{T,T}, grid)::GBSSN_Variables
#
# temp = zero(GBSSN_Variables{T,T}, grid)::GBSSN_Variables
#
# saveat = dt
#
# trange = tmin:dt:tmax
#
# tlen = length(trange)
#
# save_range = tmin:saveat:tmax
#
# save_len = length(save_range)
#
# save = Array{GBSSN_Variables,1}(undef,save_len)
#
# save[1] = regstate
#
# for t = trange
#
#     rhs!(dtstate, regstate, param, t)
#
#     k1 .= dtstate
#
#     temp .= regstate.+dt.*k1./2
#
#     rhs!(dtstate, temp, param , t + dt/2)
#
#     k2 .= dtstate
#
#     temp .= regstate.+dt.*k2./2
#
#     rhs!(dtstate, temp, param , t + dt/2)
#
#     k3 .= dtstate
#
#     temp .= regstate.+dt.*k3
#
#     rhs!(dtstate, temp, param , t + dt)
#
#     k4 .= dtstate
#
#     @. regstate = regstate + dt*(k1+2*k2+2 *k3+k4)/6
#
#     #save
#
#     if t in save_range
#
#         i = Int(round((tmax-tmin)/saveat) + 1)
#
#         save[i] = regstate
#
#     end
#
# end
#
# function M_init(::Type{T}, grid::Grid, param) where {T}
#
#     ############################################
#     # Specifies the Initial Conditions
#     ############################################
#
#     n = grid.ncells + 4
#     domain = grid.domain
#     initgrid = grid
#     drt = spacing(grid)
#     r = param[4]
#     drdrt = param[5]
#     d2rdrt = param[6]
#     m = param[7]
#     rtspan = param[8]
#
#     num = 0
#
#     fρ(M,rt) = (2*fK𝜙(rt)^2 + (1/2)*(fχ(rt)/fγtrr(M,rt))*f∂𝜙(rt)^2
#         + (1/2)*m^2*f𝜙(rt)^2)
#
#     fSr(rt) = 2*fK𝜙(rt)*f∂𝜙(rt)
#
#     #f∂M(M,rt) = 4*pi*(r(rt)^2)*fρ(M,rt)
#
#     function f∂M(M,rt)
#          if rt < 2
#              return 0.
#          else
#              4*pi*(r(rt)^2)*fρ(M,rt)
#          end
#     end
#
#     function f𝓗(M,∂M,rt)
#          (-(3/2)*(fArr(M,∂M,rt)/fγtrr(M,rt))^2 + (2/3)*fK(M,∂M,rt)^2
#          - (5/2)*((f∂χ(rt)^2)/fχ(rt))/fγtrr(M,rt) + 2*f∂2χ(rt)/fγtrr(M,rt)
#          + 2*fχ(rt)/fγtθθ(rt) - 2*fχ(rt)*(f∂2γtθθ(rt)/fγtθθ(rt))/fγtrr(M,rt)
#          + 2*f∂χ(rt)*(f∂γtθθ(rt)/fγtθθ(rt))/fγtrr(M,rt)
#          + fχ(rt)*(f∂γtrr(M,∂M,rt)/(fγtrr(M,rt)^2))*(f∂γtθθ(rt)/fγtθθ(rt))
#          - f∂χ(rt)*f∂γtrr(M,∂M,rt)/(fγtrr(M,rt)^2)
#          + (1/2)*fχ(rt)*((f∂γtθθ(rt)/fγtθθ(rt))^2)/fγtrr(M,rt) - 16*pi*fρ(M,rt))
#     end
#
#     fαreg(M,rt) = 2*M
#     fγtrrreg(M,rt) = 2*M
#     fArrreg(M,∂M,rt) = real((r(rt)+ 0im)^(5/2))*fArr(M,∂M,rt)
#     fKreg(M,∂M,rt) = real((r(rt)+ 0im)^(3/2))*fK(M,∂M,rt)
#
#     # Constraint Equations
#
#     rtspan = (rtspan[1], rtspan[2])
#     #rtspan = (rtspan[2], 0.5)
#
#     function constraintSystem(M, param, rt)
#         f∂M(M,rt)
#     end
#
#     # function boundaryCondition!(residual, M, param, rt)
#     #     residual = M[end] - 1. #inner boundary condition
#     # end
#
#     atol = 1e-15
#
#     BVP = ODEProblem(constraintSystem, 1., rtspan, param)
#     M = solve(BVP, Tsit5(), abstol=atol, dt=drt, adaptive=false)
#
#     ∂M(rt) = f∂M(M(rt),rt)
#
#     # M(rt) = 1.
#     # ∂M(rt) = 0
#
#
#     state = GridFun(grid, M)
#
#     return
#
# end
#
# function M_rhs(M::GridFun, param, t)
#
#     #########################################################
#     # Source Terms and Source Evolution
#     #
#     # This currently includes the addition of source terms
#     # to GR that come from a Klein-Gordon scalar field
#     #
#     #########################################################
#
#     # Klein-Gordon System
#
#     ∂t𝜙 = βr.*∂𝜙 - 2*α.*K𝜙
#     ∂tK𝜙 = (βr.*∂K𝜙 + α.*K.*K𝜙 - (1/2)*α.*χ.*∂2𝜙./γtrr
#         + (1/4)*α.*χ.*∂γtrr.*∂𝜙./γtrr.^2 - (1/4)*α.*∂χ.*∂𝜙./γtrr
#         - (1/2)*χ.*∂α.*∂𝜙./γtrr - (1/2)*χ.*∂γtθθ.*∂𝜙./(γtrr.*γtθθ)
#         + (1/2)*∂χ.*∂𝜙./(γtrr) + (1/2)*m^2*𝜙)
#
#     ρ = 2*K𝜙.^2 + (1/2)*(χ./γtrr).*∂𝜙.^2 + (1/2)*m^2*𝜙.^2
#     #Lower Index
#     Sr = 2*γtrr.*K𝜙.*∂𝜙./χ
#     # S = 6*K𝜙.^2 - (1/2)*(χ./γtrr).*∂𝜙.^2 - (3/2)*m^2*𝜙.^2
#     # Srr = (γtrr./χ).*(2*K𝜙.^2 + (1/2)*(χ./γtrr).*∂𝜙.^2 - (1/2)*m^2*𝜙.^2)
#
#     # ∂tArr .+= -8*pi*α.*(χ.*Srr - (1/3)*S.*γtrr)
#     # ∂tK .+= 4*pi*α.*(ρ + S)
#     # ∂tΓr .+= -16*pi*α.*Sr./γtrr
#
#     # Inner temporal boundary Conditions
#
#     # ∂tα[1:2] .= 0.
#     # ∂tA[1:2] .= 0.
#     # ∂tβr[1:2] .= 0.
#     # ∂tBr[1:2] .= 0.
#     # ∂tχ[1:2] .= 0.
#     # ∂tγtrr[1:2] .= 0.
#     # ∂tγtθθ[1:2] .= 0.
#     # ∂tArr[1:2] .= 0.
#     # ∂tK[1:2] .= 0.
#     # ∂tΓr[1:2] .= 0.
#     # ∂t𝜙[1:2] .= 0.
#     # ∂tK𝜙[1:2] .= 0.
#
#     return GBSSN_Variables(∂tαreg,∂tA,∂tβr,∂tBr,∂tχ,∂tγtrrreg,∂tγtθθreg,∂tArrreg,∂tKreg,∂tΓr,∂t𝜙,∂tK𝜙)
#
# end


text:
using MAT

export read_datasets

abstract type Dataset end

abstract type LabelledDataset <: Dataset end

"""
$(TYPEDEF)

Dataset of images stored as a 4-dimensional array of size `(num_samples, image_height,
image_width, num_channels)`, with accompanying labels (sorted in the same order) of size
`num_samples`.
"""
struct LabelledImageDataset{T<:Real,U<:Integer} <: LabelledDataset
    images::Array{T,4}
    labels::Array{U,1}

    function LabelledImageDataset{T,U}(
        images::Array{T,4},
        labels::Array{U,1},
    )::LabelledImageDataset where {T<:Real,U<:Integer}
        (num_image_samples, image_height, image_width, num_channels) = size(images)
        (num_label_samples,) = size(labels)
        @assert num_image_samples == num_label_samples
        return new(images, labels)
    end
end

function LabelledImageDataset(
    images::Array{T,4},
    labels::Array{U,1},
)::LabelledImageDataset where {T<:Real,U<:Integer}
    LabelledImageDataset{T,U}(images, labels)
end

function num_samples(dataset::LabelledDataset)
    return length(dataset.labels)
end

function Base.show(io::IO, dataset::LabelledImageDataset)
    image_size = size(dataset.images[1, :, :, :])
    num_samples = MIPVerify.num_samples(dataset)
    min_pixel = minimum(dataset.images)
    max_pixel = maximum(dataset.images)
    min_label = minimum(dataset.labels)
    max_label = maximum(dataset.labels)
    num_unique_labels = length(unique(dataset.labels))
    print(
        io,
        "{LabelledImageDataset}",
        "\n    `images`: $num_samples images of size $image_size, with pixels in [$min_pixel, $max_pixel].",
        "\n    `labels`: $num_samples corresponding labels, with $num_unique_labels unique labels in [$min_label, $max_label].",
    )
end

"""
$(TYPEDEF)

Named dataset containing a training set and a test set which are expected to contain the
same kind of data.

$(FIELDS)
"""
struct NamedTrainTestDataset{T<:Dataset,U<:Dataset} <: Dataset
    """
    Name of dataset.
    """
    name::String
    """
    Training set.
    """
    train::T
    """
    Test set.
    """
    test::U
    # TODO (vtjeng): train and test should be the same type of struct (but might potentially have different parameters).
end

function Base.show(io::IO, dataset::NamedTrainTestDataset)
    #! format: off
    print(
        io,
        "$(dataset.name):",
        "\n  `train`: $(dataset.train)",
        "\n  `test`: $(dataset.test)",
    )
    #! format: on
end

"""
$(SIGNATURES)

Makes popular machine learning datasets available as a `NamedTrainTestDataset`.

# Arguments
* `name::String`: name of machine learning dataset. Options:
    * `MNIST`: [The MNIST Database of handwritten digits](http://yann.lecun.com/exdb/mnist/). Pixel values in original dataset are provided as uint8 (0 to 255), but are scaled to range from 0 to 1 here.
    * `CIFAR10`: [Labelled subset in 10 classes of 80 million tiny images dataset](https://www.cs.toronto.edu/~kriz/cifar.html). Pixel values in original dataset are provided as uint8 (0 to 255), but are scaled to range from 0 to 1 here.
"""
function read_datasets(name::String)::NamedTrainTestDataset
    name = lowercase(name)

    if name in ["mnist", "cifar10"]
        dir = joinpath("datasets", name)

        m_train = prep_data_file(dir, "$(name)_int_train.mat") |> matread
        train = LabelledImageDataset(m_train["images"] / 255, m_train["labels"][:])

        m_test = prep_data_file(dir, "$(name)_int_test.mat") |> matread
        test = LabelledImageDataset(m_test["images"] / 255, m_test["labels"][:])
        return NamedTrainTestDataset(name, train, test)
    else
        throw(DomainError("Dataset $name not supported."))
    end
end

ext: jl


text:
using Gadfly

"""
Coins Model Selection Demo
Attempt to determine if a coin is biased or not as we accumulate more and more data, (coin flips).
"""
function coinsModelSelDemo()
    θ = .7
    N = 5
    α = 1
    ε = 0.001

    αh = α
    αt = α
    ii = 0:(2^N-1)

    flips = hcat([digits(i, 2, N) for i=ii]...)
    Nh    = reshape(sum(flips, 1), length(ii))
    Nt    = N-Nh
    marglik = exp(lbeta(αh+Nh, αt+Nt) - lbeta(αh, αt))
    mle     = Nh./N
    loglik  = Nh.*log(mle + ε) + Nt.*log(1-mle + ε)
    logBF   = lbeta(αh+Nh, αt+Nt) - lbeta(αh, αt) - N*log(0.5)

    # Sort in order of number of heads
    idx = sortperm(Nh)
    Nh      = Nh[idx]
    marglik = marglik[idx]
    logBF   = logBF[idx]
    loglik  = loglik[idx]

    # Plot marginal likelihood
    p1 = plot(x=ii, y=marglik, Geom.line, Geom.point, Theme(default_color=colorant"blue", line_width=3px),
              yintercept=[(1/2)^N], Geom.hline(color=colorant"black", size=3px))

    # Plot Bayes factors
    p2 = plot(x=ii, y=exp(logBF), Geom.line, Geom.point, Theme(default_color=colorant"blue", line_width=3px))

    # Plot BIC
    BIC1 = loglik - 1
    p3 = plot(x=ii, y=BIC1, Geom.line, Geom.point, Theme(default_color=colorant"blue", line_width=3px))

    # Plot log of marginal likelihood
    p4 = plot(x=ii, y=log(marglik), Geom.line, Geom.point, Theme(default_color=colorant"blue", line_width=3px))

    [p1, p2, p3, p4]
end


text:
using Gurobi
using QuadraticToBinary

const GRB_OPTIMIZER = Gurobi.Optimizer()
MOI.set(GRB_OPTIMIZER, MOI.Silent(), true)
const GRB_CACHE = MOIU.UniversalFallback(MOIU.Model{Float64}())
const GRB_CACHED = MOIU.CachingOptimizer(GRB_CACHE, GRB_OPTIMIZER)
const GRB_BRIDGED = MOI.Bridges.full_bridge_optimizer(GRB_CACHED, Float64)

push!(solvers, (opt = GRB_BRIDGED, mode = BilevelJuMP.SOS1Mode()))
push!(solvers_sos, (opt = GRB_BRIDGED, mode = BilevelJuMP.SOS1Mode()))
push!(solvers_quad, (opt = GRB_BRIDGED, mode = BilevelJuMP.SOS1Mode()))
push!(solvers_sos_quad, (opt = GRB_BRIDGED, mode = BilevelJuMP.SOS1Mode()))

push!(solvers_sos_quad_bin, (
    opt = QuadraticToBinary.Optimizer{Float64}(GRB_BRIDGED),
    mode = BilevelJuMP.SOS1Mode()))

push!(solvers, (opt = GRB_BRIDGED, mode = BilevelJuMP.PositiveSOS1Mode()))
push!(solvers_sos, (opt = GRB_BRIDGED, mode = BilevelJuMP.PositiveSOS1Mode()))
push!(solvers_quad, (opt = GRB_BRIDGED, mode = BilevelJuMP.PositiveSOS1Mode()))
push!(solvers_sos_quad, (opt = GRB_BRIDGED, mode = BilevelJuMP.PositiveSOS1Mode()))


text:
#Benchmarking of FVGQ20 variants
#loads the files from the unit test folder, where they should be manually updated.
using DifferentiableStateSpaceModels, BenchmarkTools
using DifferentiableStateSpaceModels.Examples
isdefined(Main, :FVGQ20) || include(joinpath(pkgdir(DifferentiableStateSpaceModels),
                                             "test/generated_models/FVGQ20.jl"))
const FVGQ = BenchmarkGroup()
const FVGQ["fvgq"] = BenchmarkGroup()

# Prepare the model for the FVGQ20 observables
const m_fvgq = PerturbationModel(Main.FVGQ20)
# Basic Steady State
p_d_fvgq = (β = 0.998, h = 0.97, ϑ = 1.17, κ = 9.51, α = 0.21, θp = 0.82, χ = 0.63,
            γR = 0.77, γy = 0.19, γΠ = 1.29, Πbar = 1.01, ρd = 0.12, ρφ = 0.93, ρg = 0.95,
            g_bar = 0.3, σ_A = exp(-3.97), σ_d = exp(-1.51), σ_φ = exp(-2.36),
            σ_μ = exp(-5.43), σ_m = exp(-5.85), σ_g = exp(-3.0), Λμ = 3.4e-3, ΛA = 2.8e-3)
const p_f_fvgq = (δ = 0.025, ε = 10, ϕ = 0, γ2 = 0.001, Ω_ii = sqrt(1e-5))
const c_fvgq_1 = SolverCache(m_fvgq, Val(1), p_d_fvgq)
const c_fvgq_2 = SolverCache(m_fvgq, Val(2), p_d_fvgq)

# The derivatives caches need to be already filled in case they are done out of order
const c_fvgq_1_p = SolverCache(m_fvgq, Val(1), p_d_fvgq)
const c_fvgq_2_p = SolverCache(m_fvgq, Val(2), p_d_fvgq)
generate_perturbation(m_fvgq, p_d_fvgq, p_f_fvgq, Val(1); cache = c_fvgq_1_p)
generate_perturbation(m_fvgq, p_d_fvgq, p_f_fvgq, Val(2); cache = c_fvgq_2_p)

# Can't hurt to  call the derivatives to precompile as well, eventhough they will be written over
generate_perturbation_derivatives!(m_fvgq, p_d_fvgq, p_f_fvgq, c_fvgq_1_p)
generate_perturbation_derivatives!(m_fvgq, p_d_fvgq, p_f_fvgq, c_fvgq_2_p)

# fvgq
const FVGQ["fvgq"]["first_order"] = @benchmarkable generate_perturbation($m_fvgq, $p_d_fvgq,
                                                                         $p_f_fvgq, Val(1);
                                                                         cache = $c_fvgq_1)
const FVGQ["fvgq"]["first_order_p"] = @benchmarkable generate_perturbation_derivatives!($m_fvgq,
                                                                                $p_d_fvgq,
                                                                                $p_f_fvgq,
                                                                                $c_fvgq_1_p)
const FVGQ["fvgq"]["second_order"] = @benchmarkable generate_perturbation($m_fvgq,
                                                                          $p_d_fvgq,
                                                                          $p_f_fvgq, Val(2);
                                                                          cache = $c_fvgq_2)
const FVGQ["fvgq"]["second_order_p"] = @benchmarkable generate_perturbation_derivatives!($m_fvgq,
                                                                                $p_d_fvgq,
                                                                                $p_f_fvgq,
                                                                                $c_fvgq_2_p)
const FVGQ["fvgq"]["second_order_no_cache"] = @benchmarkable generate_perturbation($m_fvgq,
                                                                                $p_d_fvgq,
                                                                                $p_f_fvgq,
                                                                                Val(2)) # just to see how long the cache takes

# return for the suite
FVGQ

# Uncomment to run on its own, otherwise nests into full benchmarking
# results = run(fvgq; verbose = true)
# results["fvgq"] # etc.


text:
"""
    `get_align_strand(PAIRED_END,FLAG1,FLAG2)`

    Function that returns the strand of methylation call based on Bismark's logic. In single-end mode,
    OT and CTOT reads will both receive a FLAG of 0 while OB and CTOB both receive a FLAG of 16. In
    paired end mode:

                    Read 1       Read 2
        OT:         99          147
        OB:         83          163
        CTOT:      147           99
        CTOB:      163           83

    This table was extracted from
        https://github.com/FelixKrueger/Bismark/issues/151

    # Examples
    ```julia-repl
    julia> CpelNano.get_align_strand(true,UInt16(99),UInt16(147))
    "OT"
    ```
"""
function get_align_strand(pe::Bool, flag1::UInt16, flag2::UInt16)::String

    # Treat SE and PE separately
    if !pe
        if flag1 == 0
            s = "OT"
        elseif flag1 == 16
            s = "OB"
        else
            print_log("SE: was expecting FLAG 0 or and encountered $(flag1) instead.")
            print_log("Exiting julia ...")
            exit(1)
        end
    else
        if flag1 == 99 && flag2 == 147
            s = "OT"
        elseif flag1 == 83 && flag2 == 163
            s = "OB"
        elseif flag1 == 147 && flag2 == 99
            s = "CTOT"
        elseif flag1 == 163 && flag2 == 83
            s = "CTOB"
        else
            print_log("PE: unexpected flag combination. Expected 99/147, 147/99, 83/163, 163/83.")
            print_log("Exiting julia ...")
            exit(1)
        end
    end

    # Return strand
    return s

end
"""
    `order_bams(PAIRED_END,RECORDS)`

    Function that returns an AlignTemp object with R1 as the first record in the forward strand and
    with the methylation call strand taken from `get_align_strand`.

    # Examples
    ```julia-repl
    julia> CpelNano.order_bams(true,RECORDS)
    ```
"""
function order_bams(pe::Bool, records::Vector{BAM.Record})::AlignTemp

    # Check which record appears first in forward strand
    if BAM.position(records[1]) <= BAM.position(records[2])
        s = get_align_strand(pe, BAM.flag(records[1]), BAM.flag(records[2]))
        return AlignTemp(s, records[1], records[2])
    else
        s = get_align_strand(pe, BAM.flag(records[2]), BAM.flag(records[1]))
        return AlignTemp(s, records[2], records[1])
    end

end
"""
    `clean_records(PAIRED_END,RECORDS)`

    Function that takes a set of records and returns an AllAlignTemps object that contains all the
    properly aligned templates as an array of AlignTemp, which contains information about the
    methylation call strand as well as the relevant BAM records. In the PE case, R1 corresponds to
    the BAM record that appears before in the forward strand.

    # Examples
    ```julia-repl
    julia> CpelNano.clean_records(true,RECORDS)
    ```
"""
function clean_records(pe::Bool, records::Vector{BAM.Record})::AllAlignTemps

    # Initialize struct
    out = AllAlignTemps(pe, [])

    # Consider PE vs. SE
    if pe
        # Handle PE case, first retrieve unique template names
        temp_names = unique([BAM.tempname(x) for x in records])
        for temp_name in temp_names
            # Get records with same template name
            temp_recs = filter(x -> BAM.tempname(x) == temp_name, records)

            # There should only be two records with the same template name
            length(temp_recs) == 2 && push!(out.templates, order_bams(pe, temp_recs))
        end
    else
        # Handle SE case
        out.templates = [AlignTemp(get_align_strand(pe, BAM.flag(x), UInt16(0)), x, BAM.Record()) for x in records]
    end

    # Return struct
    return out

end
"""
    `try_olaps(READER,CHR,WINDOW)`

    Function that tries to find BAM records overlaping with `CHR` at positions `WINDOW`.

    # Examples
    ```julia-repl
    julia> CpelNano.try_olaps(reader,chr,win)
    ```
"""
function try_olaps(reader::BAM.Reader, chr::String, win::UnitRange{Int64})::Vector{BAM.Record}

    # NOTE: known bug that has to do w/ BioAlignments when close to end?
    records_olap =
    try
        collect(eachoverlap(reader, chr, minimum(win):maximum(win)))
    catch
        Array{BAM.Record,1}()
    end

    return records_olap

end
"""
    `split_bam(BAM,FASTA,WIN_SIZE)`

    Function that splits BAM file in equally size windows of size WIN_SIZE.

    # Examples
    ```julia-repl
    julia> CpelNano.split_bam(bam,fasta,win_size)
    ```
"""
function split_bam(bam::String, fasta::String, win_size::Int64)::Nothing

    # Output BAM files prefix
    bam_prefix = splitext(bam)[1]

    # Get FASTA info
    chr_names, chr_sizes = get_genome_info(fasta)

    # Open reader
    j = 1
    reader = open(BAM.Reader, bam, index=bam * ".bai")
    header = BAM.header(reader)

    for (i, chr) in enumerate(chr_names)

        # Print info
        print_log("Splitting chromosome $(chr)...")

        # Get partitions
        chr_part = 1:win_size:chr_sizes[i]

        # Get reads overlapping with each region
        for chr_st in chr_part

            # Get end
            chr_end = min(chr_sizes[i] - 200, chr_st + win_size)

            # Get overlap reads
            recs = try_olaps(reader, chr, chr_st:chr_end)

            # Open BAM writer with header
            bamw = BAM.Writer(BGZFStream(open(bam_prefix * "_$(j).bam", "w"), "w"), header)

            # Write records
            map(rec -> write(bamw, rec), recs)

            # Close writer
            close(bamw)

            # Increase counter
            j += 1

        end

    end

    # Close reader
    close(reader)

    # Return nothing
    return nothing

end
"""
    `get_calls_bam(BAM_PATH,CHR,FEAT_ST,FEAT_END,CPG_POS,PE,TRIM)`

    Function that reads in BAM file in `BAM_PATH` and returns methylation vectors for those records
    that overlap with (1-based) genomic coordinates `chr:chrSt-chrEnd` at `cpg_pos`. A trimming given
    by TRIM=(Rf_5p,Rf_3p,Rr_5p,Rr_3p) is applied to the reads. The information was taken from:

        XM: meth calls (https://github.com/FelixKrueger/Bismark/tree/master/Docs)
        XS: uniqueness (http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml)

    For info on OT, CTOT, OB, CTOB nomenclature see

        http://software.broadinstitute.org/software/igv/book/export/html/37.

    # Examples
    ```julia-repl
    julia> CpelNano.get_calls_bam(BAM_PATH,"chr1",30,80,[40,60],false,(0,0,0,0))
    ```
"""
function get_calls_bam(bam::String,chr::String,roi_st::Int64,roi_end::Int64,cpg_pos::Vector{Int64},
    pe::Bool,trim::NTuple{4,Int64})::Tuple{Int64,Vector{Vector{MethCallCpgGrp}}}

    # Init return object
    calls = Vector{Vector{MethCallCpgGrp}}()

    # Get records overlapping window if any
    reader =
    try
        open(BAM.Reader, bam, index=bam * ".bai")
    catch
        nothing
    end
    isnothing(reader) && return (0, Vector{MethCallCpgGrp}[])

    # If no error proceed
    chr_size = reader.refseqlens[findfirst(x -> x == chr, reader.refseqnames)]
    roi_end = min(roi_end, chr_size - 151)
    records_olap = try_olaps(reader, chr, roi_st:roi_end)
    length(records_olap) > 0 || return 0, calls
    close(reader)

    # Relevant flags in BAM file (both Bismark and Arioc)
    filter!(x -> (BAM.ismapped(x)) && (haskey(x, "XM")) && (!haskey(x, "XS")) && (BAM.flag(x) in
           FLAGS_ALLOWED) && (BAM.mappingquality(x) > THRESH_MAPQ),records_olap)

    # Clean records & organize them
    records_org = clean_records(pe, records_olap)

    # Loop over organized records
    for record in records_org.templates

        # Get methylation call and offset (depends on strand where call is made)
        if !(records_org.paired)

            # Obtain methylation call from single end
            meth_call = record.R1[:"XM"]
            meth_call = SubString(meth_call, (1 + trim[1]):(length(meth_call) - trim[2]))
            OFFSET = record.strand in ["OT","CTOT"] ? 1 : 2
            OFFSET -= BAM.position(record.R1) + trim[1]

        else

            # Obtain methylation call
            R1_call = record.R1[:"XM"]
            R2_call = record.R2[:"XM"]
            R1_call = SubString(R1_call, (1 + trim[1]):(length(R1_call) - trim[2]))
            R2_call = SubString(R2_call, (1 + trim[4]):(length(R2_call) - trim[3]))
            dist_st = abs(BAM.position(record.R2) + trim[4] - BAM.position(record.R1) - trim[1])
            meth_call = R1_call * "."^max(0, dist_st - length(R1_call))
            meth_call *= R2_call[max(1, (length(R1_call) + 1 - dist_st)):end]
            OFFSET = record.strand in ["OT","CTOT"] ? 1 : 2
            OFFSET -= BAM.position(record.R1) + trim[1]

        end

        # Cross positions of CpG sites if template contains CpG sites
        obs_cpgs = [x.offset for x in eachmatch(r"[zZ]", meth_call)] .- OFFSET
        length(obs_cpgs) > 0 || continue

        # If overlapping CpG sites then store and add to xobs
        olap_cpgs = findall(x -> x in obs_cpgs, cpg_pos)
        length(olap_cpgs) > 0 || continue
        x = zeros(Int64, length(cpg_pos))
        obs_cpgs = obs_cpgs[findall(x -> x in cpg_pos, obs_cpgs)] .+ OFFSET
        x[olap_cpgs] = reduce(replace, ["Z" => 1,"z" => -1], init=split(meth_call[obs_cpgs], ""))

        # Define call for read pair
        call = Vector{MethCallCpgGrp}()
        @inbounds for n = 1:length(x)

            # Init call for CpG site
            cpg_call = MethCallCpgGrp()

            if x[n] == -1
                # If call is unmethylated assign p(y=-1|x=-1)≈1
                cpg_call.obs = true
                cpg_call.log_pyx_u = log_pyx_right_x
                cpg_call.log_pyx_m = log_pyx_wrong_x
            elseif x[n] == 1
                # If call is methylated assign p(y=1|x=1)≈1
                cpg_call.obs = true
                cpg_call.log_pyx_u = log_pyx_wrong_x
                cpg_call.log_pyx_m = log_pyx_right_x
            else
                # nothing
            end

            # Push call
            push!(call, cpg_call)

        end

        # Add to set of observations
        push!(calls, call)

    end # end loop over templates sequenced

    # Return
    return length(calls), calls

end
"""
    `pmap_analyze_reg_bam(REG_INT,CHR,NANO,FA_REC,CONFIG)`

    Function that estimates ϕ for a given analysis region starting at REG_ST and stores MML, NME, as well as the
    parameter vector ϕ.

    # Examples
    ```julia-repl
    julia> CpelNano.pmap_analyze_reg_bam(reg_int,chr,nano,fa_rec,config)
    ```
"""
function pmap_analyze_reg_bam(reg_int::UnitRange{Int64}, chr::String, bam::String, fa_rec::FASTA.Record, config::CpelNanoConfig)::RegStruct

    ## Get info about analysis region

    # Initizalize empty region structure
    rs = RegStruct()

    # Store chromosomal info (1-based)
    rs.chr = chr
    rs.chrst = minimum(reg_int)
    rs.chrend = maximum(reg_int)

    # Print out region being processed
    print_log("Processing $(rs.chr):$(rs.chrst)-$(rs.chrend)")

    # Get CpG site positions (1-based)
    get_grp_info!(rs, fa_rec, 1)
    length(rs.cpg_pos) > 9 || return rs
    config.verbose && print_log("CpG sites: $(length(rs.cpg_pos))")
    config.verbose && print_log("CpG groups: $(length(rs.cpg_grps))")

    # Retrieve region calls
    rs.m, rs.calls = get_calls_bam(bam, chr, rs.chrst, rs.chrend, rs.cpg_pos, config.pe, config.trim)
    config.verbose && print_log("Number of reads: $(rs.m)")
    config.verbose && print_log("Average depth: $(get_ave_depth(rs))")
    config.verbose && print_log("Percentage groups observed: $(perc_gprs_obs(rs) * 100)%")
    # config.verbose && length(rs.calls)>0 && print_log("Calls: $(rs.calls[1])")
    # config.verbose && readline()

    # Skip region if not enough data
    get_ave_depth(rs) >= config.min_cov || return rs
    perc_gprs_obs(rs) >= 0.66 || return rs

    # Skip region if only one group
    rs.L > 1 || return rs

    ## Estimate ϕs

    # Estimate parameter vector ϕs
    get_ϕhat_wgbs!(rs, config)

    # Leave if we couldn't estimate ϕ
    length(rs.ϕhat) > 0 || return rs
    config.verbose && print_log("ϕhat: $(rs.ϕhat)")
        # readline()

    ## Statistical summaries
    get_stat_sums!(rs, config)

    # Set estimation region as processed
    rs.proc = true

    # Return
    return rs

end
"""
    `analyze_bam(BAM,FASTA,CONFIG)`

    Function that estimates ϕ for all analysis regions and stores MML, NME, as well as the parameter vector ϕ.

    # Examples
    ```julia-repl
    julia> CpelNano.analyze_bam(bam,fasta,config)
    ```
"""
function analyze_bam(bam::String, fasta::String, config::CpelNanoConfig)::Nothing

    # Find chromosomes
    chr_names, chr_sizes = get_genome_info(fasta)

    # Create output directory if not existant
    isdir(config.out_dir) || mkdir(config.out_dir)

    # Get output file names
    config.out_files = OutputFiles(config.out_dir, config.out_prefix)

    # Check if output file exists
    check_output_exists(config) && return nothing

    # If BED file provided, identify intervals of interest
    print_log("Checking for target regions ...")
    targ_regs = isempty(config.bed_reg) ? nothing : read_bed_reg(config)

    # Loop over chromosomes
    for i = 1:length(chr_names)

        # Get chromosome name and size
        chr = chr_names[i]
        chr_size = chr_sizes[i]

        # Print info
        print_log("Processing chr $(chr) of length $(chr_size) ...")

        # Get analysis region start positions
        if config.informme_mode

            # Partition the chromosome as informME does
            chr_part = get_informME_chr_part(chr_size, config)

        else

            # Take into account group when partitioning
            if isnothing(targ_regs)
                # If no BED file provided partition entire chromosome
                print_log("Partitioning chr $(chr) ...")
                chr_part = get_chr_part(chr, config, fasta)
            else
                # If BED file provided, partition regions of interest. NOTE: ensure no split CpG groups
                print_log("Partitioning target regions in chr $(chr) ...")
                haskey(targ_regs, chr) || continue
                chr_part = get_chr_part(chr, config, fasta, targ_regs[chr])
            end

        end

        # Get FASTA record for chr
        fa_rec = get_chr_fa_rec(chr, fasta)

        # Process each analysis region in chr
        out_pmap = pmap(x -> pmap_analyze_reg_bam(x, chr, bam, fa_rec, config), chr_part)

        # Write chromosome output
        write_output(out_pmap, config)

    end

    # Return nothing
    return nothing

end


text:
module RDF

export RDFFile, read_channel, spike_times, sample_times, RDFParams, trial_labels

const SAMPLE_RATE = 24000.0f0

# ============================================================================ #
mutable struct RDFParams
    d::Dict{String, Any}
end
# ---------------------------------------------------------------------------- #
RDFParams(ifile::AbstractString) = RDFParams(parse_ini_file(ifile))
# ---------------------------------------------------------------------------- #
function trial_labels(p::RDFParams, x::Vector{UInt8})
    nv = length(p.d["values"])
    lab = Vector{Float64}()
    for k in x
        if k <= nv
            push!(lab, p.d["values"][k])
        end
    end
    return lab
end
# ---------------------------------------------------------------------------- #
@inline Base.getindex(p::RDFParams, k::AbstractString) = p.d[k]
Base.show(io::IO, p::RDFParams) = Base.show(io, "text/plain", p.d)
# ============================================================================ #
mutable struct RDFFile
    path::String
    ids::Vector{UInt8}
    types::Vector{UInt8}
    type_sizes::Dict{UInt8,Int}
    npacket::Int
    first_packet::Int
    first_packet_time::Float32
end
# ---------------------------------------------------------------------------- #
function RDFFile(path::AbstractString)
    r = open(path, "r") do io

        nchan = read(io, UInt8)
        ids = Vector{UInt8}(undef, nchan)
        read!(io, ids)
        types = Vector{UInt8}(undef, nchan)
        read!(io, types)

        np = read(io, Int32)

        offset = position(io)

        seek(io, offset + sizeof(UInt8))
        first_packet_time = read(io, Float32)

        type_sizes = Dict{UInt8,Int}(x => 0 for x in ids)

        return RDFFile(path, ids, types, type_sizes, np, offset,
            first_packet_time)
    end

    for k in r.ids
        r.type_sizes[k] = sizeof(get_eltype(r, k))
    end

    return r
end
# ---------------------------------------------------------------------------- #
function Base.show(io::IO, r::RDFFile)
    print(io, "RDFFile @ \"$(r.path)\" with ")
    nchan = length(r.ids)
    println(io, "$(nchan) channel(s):")
    for k = 1:nchan
        println("  channel $(k): id = $(r.ids[k]), type = $(get_eltype(r, k))")
    end
end
# ---------------------------------------------------------------------------- #
function read_channel(r::RDFFile, id::Integer)
    0 < id < 256 || error("Invalid channel: must be 0 < id < 256")

    elt = get_eltype(r, id)
    data = Vector{elt}(undef, 0)
    ts = Vector{Float32}(undef, 0)
    #idx = Vector{Int}(undef, 0)
    open(r.path, "r") do io
        seekend(io)
        len = position(io)
        seek(io, r.first_packet)

        for k in 1:r.npacket
            idp = read(io, UInt8)

            # correct packet arrival times to be relative to the first packet
            # of the file, this way continuous channels and event channels
            # will be synchronized "automatically"
            time = read(io, Float32) - r.first_packet_time
            nel = read(io, Int32)

            if idp == id
                #push!(idx, length(data))
                tmp = Vector{elt}(undef, nel)
                read!(io, tmp)
                append!(data, tmp)
                push!(ts, time)
            else
                pos = position(io) + nel * r.type_sizes[idp]
                if 0 < pos <= len
                    seek(io, pos)
                end
            end
        end
    end
    return data, Vector{Float64}(ts)#, idx
end
# ---------------------------------------------------------------------------- #
function get_eltype(r::RDFFile, id::Integer)
    kchan = findfirst(isequal(id), r.ids)

    kchan == nothing && error("Invalid id: $(id)")

    typ = r.types[kchan]

    if typ & 0x01 == 0x00
        if typ == 0x02
            elt = UInt8
        elseif typ == 0x04
            elt = UInt16
        elseif typ == 0x06
            elt = UInt32
        else
            error("Invalid type: $(typ)")
        end
    else
        if typ == 0x03
            elt = Int8
        elseif typ == 0x05
            elt = Int16
        elseif typ == 0x07
            elt = Int32
        elseif typ == 0x0f
            elt = Float32
        else
            error("Invalid type: $(typ)")
        end
    end

    return elt
end
# ---------------------------------------------------------------------------- #
function spike_times(data::AbstractVector{<:Real})
    win = 64
    thr = 0.97f0
    ts = Vector{Float64}()
    k = 1
    si = (1.0f0 / SAMPLE_RATE)
    while k < length(data)
        if data[k] > thr
            ke = min(k+win-1, length(data))
            _, kmx = findmax(data[k:ke])
            if kmx != nothing
                kmx += k - 1
                push!(ts, (kmx - 1) * si)
                k += win
            else
                k = length(data) + 1
            end
        else
            k += 1
        end
    end
    return ts
end
# ---------------------------------------------------------------------------- #
@inline sample_times(data::AbstractVector{<:Real}) =
    range(0, length(data) * (1.0 / SAMPLE_RATE), length=length(data))
# ---------------------------------------------------------------------------- #
function parse_ini_file(pth::AbstractString)
    isfile(pth) || error("Cannot ind file \"$(pth)\"")
    d = Dict{String, Any}()
    for line in eachline(pth, keep=false)
        m = match(r"([^\=\s]+)\s*\=\s*([^\n]+)", line)
        if m != nothing
            if m[2][1] == '['
                ke = findfirst(isequal(']'), m[2])
                ke = ke == nothing ? length(m[2]) : ke - 1
                val = map(x -> parse(Float64, x), split(m[2][2:ke], r"[,\s]+"))
                d[m[1]] = val
            else
                val = tryparse(Float64, m[2])
                if val == nothing
                    d[m[1]] = m[2]
                else
                    d[m[1]] = val
                end
            end
        else
            @warn("Invalid statement: \"$(line)\"")
        end
    end
    return d
end
# ============================================================================ #
end


text:
const Meter = BaseUnit("m")
const Second = BaseUnit("s")
const KiloGram = BaseUnit("g", prefix=3)
const Kilogram = KiloGram
const Kelvin = BaseUnit("K")
const Mole = BaseUnit("mol")
const Candela = BaseUnit("cd")
const Radian = BaseUnit("rad")
const Ampere = BaseUnit("A")
const Photon = BaseUnit("γ")
const Entity = BaseUnit("entity")

const Gram = Milli*KiloGram
const Newton = DerivedUnit("N",KiloGram*Meter/Second^2)
const Joule = DerivedUnit("J", Newton*Meter)
const Coulomb = DerivedUnit("C", Ampere*Second)
const Volt = DerivedUnit("V", Joule/Coulomb)
const Ohm = DerivedUnit("Ω", Volt/Ampere)
const Hertz = DerivedUnit("Hz", 1/Second)
const Steradian = DerivedUnit("sr", Radian^2)
const Pascal = DerivedUnit("Pa", Newton/Meter^2)
const Watt = DerivedUnit("W",Joule/Second)
const Farad = DerivedUnit("F",Coulomb/Volt)
const Siemens = DerivedUnit("S", Ohm^-1)
const Weber = DerivedUnit("Wb", Joule/Ampere)
const Tesla = DerivedUnit("T", Newton/(Ampere*Meter))
const Henry = DerivedUnit("H", Weber/Ampere)
const Lumen = DerivedUnit("lm", Candela/Steradian)
const Lux = DerivedUnit("lx", Lumen/Meter^2)
const Becquerel = DerivedUnit("Bq",Photon/Second)
const Katal = DerivedUnit("kat",Mole/Second)

export Meter, Second, KiloGram, Kilogram, Kelvin, Mole, Candela, Radian, Ampere, Photon, Gram
export Entity, Newton, Joule, Coulomb, Volt, Ohm, Hertz, Steradian, Pascal, Watt, Farad
export Siemens, Weber, Tesla, Henry, Lumen, Lux, Becquerel, Katal


text:
# AdaptivePenalization penalty
# ------------------------------------------------------------
@with_kw mutable struct AdaptiveData{T <: AbstractFloat}
    r::T = 1.0
end

@with_kw mutable struct AdaptiveParameters{T <: AbstractFloat}
    ζ::T = 0.1
    γ₁::T = 1e-5
    γ₂::T = 0.01
    γ₃::T = 0.25
    σ::T = 1e-5
    α::T = 0.95
    θ::T = 1.1
    ν::T = 0.1
    β::T = 1.1
    η::T = 1.25
end

struct AdaptivePenalization{T <: AbstractFloat} <: AbstractPenalization
    data::AdaptiveData{T}
    parameters::AdaptiveParameters{T}

    function AdaptivePenalization(r::AbstractFloat; kw...)
        T = typeof(r)
        return new{T}(AdaptiveData{T}(; r = r), AdaptiveParameters{T}(;kw...))
    end
end
function penalty(::AbstractProgressiveHedgingSolver, penalty::AdaptivePenalization)
    return penalty.data.r
end
function init_penalty!(ph::AbstractProgressiveHedgingSolver, penalty::AdaptivePenalization)
    update_dual_gap!(ph)
    @unpack δ₂ = ph.data
    @unpack ζ = penalty.parameters
    penalty.data.r = max(1., 2*ζ*abs(calculate_objective_value(ph)))/max(1., δ₂)
end
function update_penalty!(ph::AbstractProgressiveHedgingSolver, penalty::AdaptivePenalization)
    @unpack δ₁, δ₂ = ph.data
    @unpack r = penalty.data
    @unpack γ₁, γ₂, γ₃, σ, α, θ, ν, β, η = penalty.parameters

    δ₂_prev = length(ph.dual_gaps) > 0 ? ph.dual_gaps[end] : Inf

    μ = if δ₁/norm(ph.ξ,2)^2 >= γ₁
        if (δ₁-δ₂)/(1e-10 + δ₂) > γ₂
            α
        elseif (δ₂-δ₁)/(1e-10 + δ₁) > γ₃
            θ
        else
            1.
        end
    elseif δ₂ > δ₂_prev
        if (δ₂-δ₂_prev)/δ₂_prev > ν
            β
        else
            1.
        end
    else
        η
    end
    penalty.data.r = μ*r
end

# API
# ------------------------------------------------------------
"""
    Adaptive

...
# Parameters

...
"""
struct Adaptive{T <: AbstractFloat} <: AbstractPenalizer
    r::T
    parameters::Dict{Symbol,Any}

    function Adaptive(r::AbstractFloat; kw...)
        T = typeof(r)
        return new{T}(r, Dict{Symbol,Any}(kw))
    end
end
Adaptive(; kw...) = Adaptive(1.0; kw...)

function (adaptive::Adaptive)()
    return AdaptivePenalization(adaptive.r; adaptive.parameters...)
end

function str(::Adaptive)
    return "adaptive penalty"
end


text:
module LennardJones

include("MD_Structure.jl")

using .MD_Structure

export MD

end


text:
@testset "NegativeBinomialLikelihood{<:NBParamFailure,<:LogisticLink}" begin
    test_auglik(NegativeBinomialLikelihood(NBParamFailure(10), LogisticLink()))
    test_auglik(NegativeBinomialLikelihood(NBParamFailure(5.5), LogisticLink()))
end


text:
# This example computes the space-charge field of a particle beam with single momentum.
# The particle beam moves along z-direction with γ=30.

using NChargedBodyTreecode

# number of particles
const N = 10000

# create position and momentum distribution of N particles
positions = rand(3, N)
momenta = zeros(3, N)
momenta[3, :] .= 29.9833287

# create a particle beam in which each particle's charge and mass are -1 and 1
beam = Particles(; pos=positions, mom=momenta, charge=-1.0, mass=1.0)

# create 4 copies of a particle beam
beam0 = deepcopy(beam)
beam1 = deepcopy(beam)
beam2 = deepcopy(beam)
beam3 = deepcopy(beam)

# update particle field by brute-force method (lambda is a characteristic length for the normalization of length quantity)
updateParticlesField!(beam0, BruteForce(); lambda=1.0)

# update particle field by different treecode methods with the following parameters:
#   n: degree of interpolation
#   N0: maximum number of particles in the leaf cluster
#   eta: admissibility parameter
#   lambda: a characteristic length for the normalization of length quantity
####
const n = 4
const N0 = 128
const eta = 0.5
updateParticlesField!(beam1, TreecodeStretch(eta=eta, N0=N0, n=n); lambda=1.0)
updateParticlesField!(beam2, TreecodeAvgRestFrame(eta=eta, N0=N0, n=n); lambda=1.0)
updateParticlesField!(beam3, TreecodeUniform(eta=eta, N0=N0, n=n); lambda=1.0)

# relative errors of different treecodes (comparing to brute-force)
error1 = relerror(beam1, beam0)
error2 = relerror(beam2, beam0)
error3 = relerror(beam3, beam0)
println("----------------")
println("relative errors")
println("----------------")
println("Treecode-Stretch: $(error1)")
println("Treecode-AVGRF: $(error2)")
println("Treecode-Uniform: $(error3)")


text:
# Regenerate results using the parameter array
println(now())
using Debug
# Initialize simulation functions
fpath="../src/new_SimulationFunctions.jl" # Loads up DataFile , SolveBalances , Balances, Kinetics, Control
include(fpath)

readdir,sim_select = ARGS[1],eval(parse(ARGS[2]))

# Establish simulation parameters
TSTART = 0.0;
Ts = 1.0; # 0.1
TSTOP = 250.0; # 100 ; 500

# Initialize intial data dictionary
data_dictionary = DataFile(TSTART,TSTOP,Ts)
len_cpa = length(data_dictionary["CONTROL_PARAMETER_ARRAY"])
len_rca = length(data_dictionary["RATE_CONSTANT_ARRAY"])
len_sca = length(data_dictionary["SATURATION_CONSTANT_ARRAY"])
sca_indices,sca_values = [],[]
for i in 1:len_sca
  val = data_dictionary["SATURATION_CONSTANT_ARRAY"][i]
if val !=0
  push!(sca_indices,i)
  push!(sca_values,val)
end
end

sca_pack = zip(sca_indices,sca_values)

 # CPA, let's give access to gain AND power
 cpa_gain = data_dictionary["CONTROL_PARAMETER_ARRAY"][:,1] # Gain
 cpa_order= data_dictionary["CONTROL_PARAMETER_ARRAY"][:,2] # Power

 # get lengths of new arrays
 len_cpa_gain = length(cpa_gain)
 len_cpa_order = length(cpa_order)
 len_sca_values = length(sca_values)

 len_sca_prime = length(sca_values)
   A = len_rca + len_cpa_gain
   B = A + len_cpa_order
   C = B + len_sca_prime

# read parameter arrays resulting from POETs

PC = readdlm(string(readdir,"/PC.dat")) #readdir

# PC is an array of parameter arrays
PC_results =[]
fraction = 1
popsize = Int(length(PC[1,:])*fraction)
println("Model Population Size ", popsize)

simulations = ["edit_file_none.jl","edit_file_tgfb.jl","edit_file_mapkin.jl","edit_file_ss.jl","edit_file_ssecad.jl","edit_file_ss_tgfb3ab.jl","edit_file_ss_dnlef1.jl","edit_file_tgfb3.jl","edit_file_dnsmad.jl","edit_file_tgfb3_dnlef1.jl","edit_file_tgfb_tgfb3ab.jl"]
# sim_select
@show sim_select
simulations = [simulations[x] for x in sim_select]
if sim_select == Any[]
popsize =1
end
@show simulations

for i in 1:popsize
  println("Regen ",i," of ",popsize)
  # Refactor data dictionary for model execution
  # need to specify implicit strucutre of parameter_array
  P = PC[:,i]
  data_dictionary["RATE_CONSTANT_ARRAY"] = P[1:len_rca]
  data_dictionary["CONTROL_PARAMETER_ARRAY"][:,1] = P[len_rca+1:A]
  data_dictionary["CONTROL_PARAMETER_ARRAY"][:,2] = P[A+1:B]
  # Rebuild SCA
  sca_pack = zip(sca_indices,P[B+1:C]) # new values
  for (i,val) in sca_pack
  data_dictionary["SATURATION_CONSTANT_ARRAY"][i]=val
  end

  # Run to Steady State
 #println("running to SS")
  XSS = EstimateSteadyState(data_dictionary)
  initial_condition_array = vec(XSS);
  data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
  xss_data_dict = deepcopy(data_dictionary)

  sim_data = [] # Keep data in RAM
  for sim in simulations
    ##### What data_dictionary should be loaded here??? clearly the active one
    active_data_dictionary = deepcopy(xss_data_dict)
    (TSIM,X) =  SolveModelGenericPOETS(TSTART,TSTOP,Ts,active_data_dictionary,string("../src/",sim))
    push!(sim_data,X)
    # now reset the data dictionary
  end
  push!(PC_results,sim_data)

end

# now get mean and std

# indices for simulation data array sim_data
# 1 "edit_file_none.jl"
# 2 "edit_file_tgfb.jl"
# 3 "edit_file_mapkin.jl"
# 4 "edit_file_ss.jl"
# 5 "edit_file_ssecad.jl"
# 6 "edit_file_ss_tgfb3ab.jl"
# 7 "edit_file_ss_dnlef1.jl"
# 8 "edit_file_tgfb3.jl"
# 9 "edit_file_dnsmad.jl"
# 10 "edit_file_tgfb3_dnlef1.jl"

# LEF 1 protein (technically this not not total, due to active amt)
index_lef = 58
# LEF 1 active
index_lef_act = 59
# Snail_slug total protein
index_ss = 62
# ecadherin protein (unbound)
index_ecad = 51
# ecad_bcat index
index_ecad_bcat = 53
# vimentin protein
index_vim = 63
# tgfb3 protein
index_tgfb3 = 92
# This is fit to ACTIVE_TCF4, data is IP from bcat and tcf-4
index_bcat_tcf = 57 # i.e. this index is for Active_TCF4
# Should probably be total protein, so sum index_bcat with index_tcf
index_tcf = 56

index_psmad = 40
index_lef_mrna = 84
index_ecad_mrna = 16

index_gsk3_p = 44
## Designate arrays for each figure

F2A = Array{Float64,1}([])
F2B = Array{Float64,1}([])
F2C = Array{Float64,1}([])
F2D = Array{Float64,1}([])
F2E = Array{Float64,1}([])
F2F = Array{Float64,1}([])
F2G = Array{Float64,1}([])
F2H = Array{Float64,1}([])
F2I = Array{Float64,1}([])
F2J = Array{Float64,1}([])
F2K = Array{Float64,1}([])
F2L = Array{Float64,1}([])

F3B = Array{Float64,1}([])

t0  = 1
t24 = 74
t48 = 149
t72 = 224
# loop over results and build each figure

function norm_dat(data)
  min = minimum(data)
  max = maximum(data)

#  if max == 0
  #  max = 1
#  end
SMALL = 1e-3
if sum(data) > SMALL
  normalized_data = (data-min)/(max-min)
else
  normalized_data = zeros(data)
end
end

#
#out_dir = "../poets/results_poets/"
#writedlm(string(out_dir,"ensemble_dynamics.dat"),PC_results)

# Panel experiments

# indices for simulation data array sim_data
# 1 "edit_file_none.jl"
# 2 "edit_file_tgfb.jl"
# 3 "edit_file_mapkin.jl"
# 4 "edit_file_ss.jl"
# 5 "edit_file_ssecad.jl"
# 6 "edit_file_ss_tgfb3ab.jl"
# 7 "edit_file_ss_dnlef1.jl"
# 8 "edit_file_tgfb3.jl"
# 9 "edit_file_dnsmad.jl"
# 10 "edit_file_tgfb3_dnlef1.jl"

pan_a = pan_b = [1,2,3]
pan_c = pan_d = [1,4,5]
pan_e = [1,4,6]
pan_f = [1,4,7]
pan_g = [1,8,9]
pan_h = pan_i = [1,8,9,10] #

pan3b = [1,8,11]

for state in PC_results

if sim_select == pan_a

# PANEL A
bar1 = state[1][t48,index_ss]
bar2 = state[2][t48,index_ss]
bar3 = state[3][t48,index_ss]
bar1,bar2,bar3=norm_dat([bar1,bar2,bar3])
push!(F2A,bar1)
push!(F2A,bar2)
push!(F2A,bar3)

# PANEL B
bar1 = state[1][t48,index_lef] + state[1][t48,index_lef_act]
bar2 = state[2][t48,index_lef] + state[2][t48,index_lef_act]
bar3 = state[3][t48,index_lef] + state[3][t48,index_lef_act]
bar1,bar2,bar3=norm_dat([bar1,bar2,bar3])
push!(F2B,bar1)
push!(F2B,bar2)
push!(F2B,bar3)

end

if sim_select == pan_c

# PANEL C
bar1 = state[1][t48,index_tgfb3]
bar2 = state[2][t48,index_tgfb3]
bar3 = state[3][t48,index_tgfb3]
bar1,bar2,bar3=norm_dat([bar1,bar2,bar3])
push!(F2C,bar1)
push!(F2C,bar2)
push!(F2C,bar3)
# PANEL D
bar1 = state[1][t48,index_vim]
bar2 = state[2][t48,index_vim]
bar3 = state[3][t48,index_vim]
bar1,bar2,bar3=norm_dat([bar1,bar2,bar3])
push!(F2D,bar1)
push!(F2D,bar2)
push!(F2D,bar3)

end

if sim_select == pan_e
# PANEL E
bar1 = state[1][t48,index_lef] + state[1][t48,index_lef_act]
bar2 = state[2][t48,index_lef] + state[2][t48,index_lef_act]
bar3 = state[3][t48,index_lef] + state[3][t48,index_lef_act]
bar1,bar2,bar3=norm_dat([bar1,bar2,bar3])
push!(F2E,bar1)
push!(F2E,bar2)
push!(F2E,bar3)

end

if sim_select == pan_f
# PANEL F # index_ecad_bcat
bar1 = state[1][t48,index_ecad] + state[1][t48,index_ecad_bcat]
bar2 = state[2][t48,index_ecad] + state[2][t48,index_ecad_bcat]
bar3 = state[3][t48,index_ecad] + state[3][t48,index_ecad_bcat]
bar1,bar2,bar3=norm_dat([bar1,bar2,bar3])
push!(F2F,bar1)
push!(F2F,bar2)
push!(F2F,bar3)
end

if sim_select == pan_g

# PANEL G
bar1 = state[1][t48,index_lef] + state[1][t48,index_lef_act]
bar2 = state[2][t48,index_lef] + state[2][t48,index_lef_act]
bar3 = state[3][t48,index_lef] + state[3][t48,index_lef_act]
bar1,bar2,bar3=norm_dat([bar1,bar2,bar3])
push!(F2G,bar1)
push!(F2G,bar2)
push!(F2G,bar3)

end

if sim_select == pan_h
# PANEL H
bar1 = state[1][t48,index_vim]
bar2 = state[2][t48,index_vim]
bar3 = state[3][t48,index_vim]
bar4 = state[4][t48,index_vim]
bar1,bar2,bar3,bar4=norm_dat([bar1,bar2,bar3,bar4])
push!(F2H,bar1)
push!(F2H,bar2)
push!(F2H,bar3)
push!(F2H,bar4)

# PANEL I
bar1 = state[1][t48,index_ecad] + state[1][t48,index_ecad_bcat]
bar2 = state[2][t48,index_ecad] + state[2][t48,index_ecad_bcat]
bar3 = state[3][t48,index_ecad] + state[3][t48,index_ecad_bcat]
bar4 = state[4][t48,index_ecad] + state[4][t48,index_ecad_bcat]
bar1,bar2,bar3,bar4=norm_dat([bar1,bar2,bar3,bar4])
push!(F2I,bar1)
push!(F2I,bar2)
push!(F2I,bar3)
push!(F2I,bar4)

end

if sim_select == pan3b

# PANEL ####
bar1 = state[1][t48,index_gsk3_p]
bar2 = state[2][t48,index_gsk3_p]
bar3 = state[3][t48,index_gsk3_p]
bar1,bar2,bar3=norm_dat([bar1,bar2,bar3])
push!(F3B,bar1)
push!(F3B,bar2)
push!(F3B,bar3)

end

# NEED TO PUT qPCR EXPERIMENTAL DATA HERE AND REFERENCE BELOW

# PANEL J

if sim_select == [2] # this set of panels doesn't work in the same way, there is no cross validation

bar1 = state[1][t0,index_ecad_mrna]
bar2 = state[1][t24,index_ecad_mrna]
bar3 = state[1][t48,index_ecad_mrna]
bar4 = state[1][t72,index_ecad_mrna]
bar1,bar2,bar3,bar4=norm_dat([bar1,bar2,bar3,bar4])
push!(F2J,bar1)
push!(F2J,bar2)
push!(F2J,bar3)
push!(F2J,bar4)

# PANEL K
bar1 = state[1][t0,index_psmad]
bar2 = state[1][t24,index_psmad]
bar3 = state[1][t48,index_psmad]
bar4 = state[1][t72,index_psmad]
bar1,bar2,bar3,bar4=norm_dat([bar1,bar2,bar3,bar4])
push!(F2K,bar1)
push!(F2K,bar2)
push!(F2K,bar3)
push!(F2K,bar4)

# PANEL L

bar1 = state[1][t0,index_lef_mrna]
bar2 = state[1][t24,index_lef_mrna]
bar3 = state[1][t48,index_lef_mrna]
bar4 = state[1][t72,index_lef_mrna]
bar1,bar2,bar3,bar4=norm_dat([bar1,bar2,bar3,bar4])
push!(F2L,bar1)
push!(F2L,bar2)
push!(F2L,bar3)
push!(F2L,bar4)
end
# indices for simulation data array sim_data
# 1 "edit_file_none.jl"
# 2 "edit_file_tgfb.jl"
# 3 "edit_file_mapkin.jl"
# 4 "edit_file_ss.jl"
# 5 "edit_file_ssecad.jl"
# 6 "edit_file_ss_tgfb3ab.jl"
# 7 "edit_file_ss_dnlef1.jl"
# 8 "edit_file_tgfb3.jl"
# 9 "edit_file_dnsmad.jl"
# 10 "edit_file_tgfb3_dnlef1.jl"
end

# get mean and standard deviation
N = length(PC_results)
# bar1=sum(F2A[1:3:end])/N
# bar2=sum(F2A[2:3:end])/N
# bar3=sum(F2A[3:3:end])/N

# println(F2A[3:3:end])
#std(Array{Float64,1}(F2A[3:3:end]))

#out_dir = "../poets/results_poets/cross_validation/"
out_dir = string(readdir,"/")

if sim_select == pan_a # AB
  F2A = [sum(F2A[1:3:end])/N,sum(F2A[2:3:end])/N,sum(F2A[3:3:end])/N,std(F2A[1:3:end]),std(F2A[2:3:end]),std(F2A[3:3:end])]
  F2B = [sum(F2B[1:3:end])/N,sum(F2B[2:3:end])/N,sum(F2B[3:3:end])/N,std(F2B[1:3:end]),std(F2B[2:3:end]),std(F2B[3:3:end])]
writedlm(string(out_dir,"F2A",".dat"),F2A)
writedlm(string(out_dir,"F2B",".dat"),F2B)
end

if sim_select == pan_c # CD
F2C = [sum(F2C[1:3:end])/N,sum(F2C[2:3:end])/N,sum(F2C[3:3:end])/N,std(F2C[1:3:end]),std(F2C[2:3:end]),std(F2C[3:3:end])]
F2D = [sum(F2D[1:3:end])/N,sum(F2D[2:3:end])/N,sum(F2D[3:3:end])/N,std(F2D[1:3:end]),std(F2D[2:3:end]),std(F2D[3:3:end])]
writedlm(string(out_dir,"F2C",".dat"),F2C)
writedlm(string(out_dir,"F2D",".dat"),F2D)
end

if sim_select == pan_e #E
F2E = [sum(F2E[1:3:end])/N,sum(F2E[2:3:end])/N,sum(F2E[3:3:end])/N,std(F2E[1:3:end]),std(F2E[2:3:end]),std(F2E[3:3:end])]
writedlm(string(out_dir,"F2E",".dat"),F2E)
end
if sim_select == pan_f #F
F2F = [sum(F2F[1:3:end])/N,sum(F2F[2:3:end])/N,sum(F2F[3:3:end])/N,std(F2F[1:3:end]),std(F2F[2:3:end]),std(F2F[3:3:end])]
writedlm(string(out_dir,"F2F",".dat"),F2F)
end
if sim_select == pan_g #G
F2G = [sum(F2G[1:3:end])/N,sum(F2G[2:3:end])/N,sum(F2G[3:3:end])/N,std(F2G[1:3:end]),std(F2G[2:3:end]),std(F2G[3:3:end])]
writedlm(string(out_dir,"F2G",".dat"),F2G)
end
if sim_select == pan_h #HI
F2H = [sum(F2H[1:4:end])/N,sum(F2H[2:4:end])/N,sum(F2H[3:4:end])/N,sum(F2H[4:4:end])/N,std(F2H[1:4:end]),std(F2H[2:4:end]),std(F2H[3:4:end]),std(F2H[4:4:end])]
F2I = [sum(F2I[1:4:end])/N,sum(F2I[2:4:end])/N,sum(F2I[3:4:end])/N,sum(F2I[4:4:end])/N,std(F2I[1:4:end]),std(F2I[2:4:end]),std(F2I[3:4:end]),std(F2I[4:4:end])]
writedlm(string(out_dir,"F2H",".dat"),F2H)
writedlm(string(out_dir,"F2I",".dat"),F2I)
end

if sim_select == pan3b # AB
F3B = [sum(F3B[1:3:end])/N,sum(F3B[2:3:end])/N,sum(F3B[3:3:end])/N,std(F3B[1:3:end]),std(F3B[2:3:end]),std(F3B[3:3:end])]
writedlm(string(out_dir,"F3B",".dat"),F3B)
end

if sim_select == []
#F2J = [sum(F2J[1:4:end])/N,sum(F2J[2:4:end])/N,sum(F2J[3:4:end])/N,sum(F2J[4:4:end])/N,std(F2J[1:4:end]),std(F2J[2:4:end]),std(F2J[3:4:end]),std(F2J[4:4:end])]
#F2K = [sum(F2K[1:4:end])/N,sum(F2K[2:4:end])/N,sum(F2K[3:4:end])/N,sum(F2K[4:4:end])/N,std(F2K[1:4:end]),std(F2K[2:4:end]),std(F2K[3:4:end]),std(F2K[4:4:end])]
#F2L = [sum(F2L[1:4:end])/N,sum(F2L[2:4:end])/N,sum(F2L[3:4:end])/N,sum(F2L[4:4:end])/N,std(F2L[1:4:end]),std(F2L[2:4:end]),std(F2L[3:4:end]),std(F2L[4:4:end])]

# J E-cad mRNA
ecad_mrna_0,sd1 = 1.0,0.04
ecad_mrna_24,sd2 = 0.6,0.05
ecad_mrna_48,sd3 = 0.3,0.03
ecad_mrna_72,sd4 = 0.0,0.001

# K pSmad
psmad_0,sd5 = 0.0,0.001
psmad_24,sd6 = 0.12,0.02
psmad_48,sd7 = 0.88,0.06
psmad_72,sd8 = 1.0,0.06

# L mRNA LEF1
lef1_mrna_0,sd9 = 0.0,0.001
lef1_mrna_24,sd10 = 0.05,0.01
lef1_mrna_48,sd11 = 0.62,0.06
lef1_mrna_72,sd12 = 1.0,0.05

F2J = [ecad_mrna_0,ecad_mrna_24,ecad_mrna_48,ecad_mrna_72,sd1,sd2,sd3,sd4]
F2K = [psmad_0,psmad_24,psmad_48,psmad_72,sd5,sd6,sd7,sd8]
F2L = [lef1_mrna_0,lef1_mrna_24,lef1_mrna_48,lef1_mrna_72,sd9,sd10,sd11,sd12]

writedlm(string(out_dir,"F2J",".dat"),F2J)
writedlm(string(out_dir,"F2K",".dat"),F2K)
writedlm(string(out_dir,"F2L",".dat"),F2L)
end
println(now())


text:
using GapFindFill, Test, Gurobi
include("loadGAMS.jl")

println("======Tesing Gap Filling===")
println("======loading data===")
#*************SETS********************************
compoundList = read_GAMS_format_file("gapfill/compounds.txt")
@assert(length(compoundList) == 1822)
reactionList = read_GAMS_format_file("gapfill/reactions.txt")
# reactionList = modelRxnList + databaseRxnList
@assert(length(reactionList) == 2888)
modelRxnList = read_GAMS_format_file("gapfill/model_reactions.txt")
@assert(length(modelRxnList) == 1578)
databaseRxnList = read_GAMS_format_file("gapfill/database_reactions.txt")
@assert(length(databaseRxnList) == 1310)
revRxnList = read_GAMS_format_file("gapfill/reversible_rxns.txt")
@assert(length(revRxnList) == 1959)
noProdCompList = read_GAMS_format_file("gapfill/no_production_metabolites.txt")
@assert(length(noProdCompList) <= 8)
cytCompoundList = read_GAMS_format_file("gapfill/cytosolic_metabolites.txt")
@assert(length(cytCompoundList) == 714)
extCompoundList = read_GAMS_format_file("gapfill/extracellular_metabolites.txt")
@assert(length(extCompoundList) == 678)

#*****************PARAMETERS*****************************
fluxUBList = read_GAMS_format_file("gapfill/upperbounds_on_fluxes.txt")
@assert(length(fluxUBList) == 2888)
fluxLBList = read_GAMS_format_file("gapfill/lowerbounds_on_fluxes.txt")
@assert(length(fluxLBList) == 2888)
stoiList = read_GAMS_format_file("gapfill/S_matrix.txt")
@assert(length(stoiList) > 8000)

#****************PARAMETERS PREPROCESSING***************************
# map rxn_name to index
rxnDict = Dict{String, Integer}()
for (id, name) in enumerate(reactionList)
  rxnDict[name] = id
end
# map compound_name to index
compoundDict = Dict{String, Integer}()
for (id, name) in enumerate(compoundList)
  compoundDict[name] = id
end
# number of reactions & compounds
noRxn = length(reactionList)
noCompound = length(compoundList)
# isRev -- isRev[i] is true if rxn in rxnDict(rxn, i) is reversible
isRev = zeros(Bool, noRxn)
for el in revRxnList
  isRev[rxnDict[el]] = true
end
@assert(length(revRxnList) == sum(isRev))
# isMd -- isMd[i] is true if rxn in rxnDict(rxn, i) is in model
isMd = zeros(Bool, noRxn)
for el in modelRxnList
  isMd[rxnDict[el]] = true
end
@assert(length(modelRxnList) == sum(isMd))
# isDb -- isDb[i] is true if rxn in rxnDict(rxn, i) is in database
isDb = zeros(Bool, noRxn)
for el in databaseRxnList
  isDb[rxnDict[el]] = true
end
@assert(length(databaseRxnList) == sum(isDb))
# isCyt -- isCyt[i] is true if cp in compoundDict(cp, i) is in cytosol
isCyt = zeros(Bool, noCompound)
for el in cytCompoundList
  isCyt[compoundDict[el]] = true
end
@assert(length(cytCompoundList) == sum(isCyt))
# isExt -- isExt[i] is true if cp in compoundDict(cp, i) is in extracellular
isExt = zeros(Bool, noCompound)
for el in extCompoundList
  isExt[compoundDict[el]] = true
end
@assert(length(extCompoundList) == sum(isExt))
# no production metabolites list
noProdID = Array{Integer, 1}()
for el in noProdCompList
  push!(noProdID, compoundDict[el])
end
# get S matrix from stoiList
stoiMatrix = zeros(Float64, noCompound, noRxn)
for el in stoiList
  stoiMatrix[compoundDict[el[1]],rxnDict[el[2]]] = el[3]
end
println("size of S matrix: ", size(stoiMatrix))
# get upperbound on fluxes
fluxUB = zeros(Float64, noRxn)
for (name, val) in fluxUBList
  fluxUB[rxnDict[name]] = val
end
println("big-M in upperbound: ", maximum(fluxUB))
# get lowerbound on fluxes
fluxLB = zeros(Float64, noRxn)
for (name, val) in fluxLBList
  fluxLB[rxnDict[name]] = val
end
println("big-M in lowerbousnd: ", minimum(fluxLB))

isExt = zeros(Bool, noCompound)

println("======calling fill_gaps_min===")
# test fill_gaps_min
results = fill_gaps_min(isMd, isDb, isRev, isCyt, isExt, noProdID,
    stoiMatrix, fluxLB, fluxUB)
println("======results report===")
for (key, val) in results
  println("To fix metabolite $(compoundList[key])")
  for el in val
    println(join([reactionList[el2] for el2 in el], " "))
  end
  println("\n")
end

println("======Gap Filling Test Passed===")


text:
# src/addresses.jl
# Address Lists for Wallets

export lists,
       list,
       enlist!,
       swaplist!,
       onlist,
       delist!,
       addresses,
       address,
       whitelist,
       whitelist!,
       onwhitelist,
       blacklist,
       blacklist!,
       onblacklist

"""```
StringLike === Union{Symbol, <: AbstractString}
```
`String`-like objects are either `Symbol`s or `AbstractString`s.
"""
const StringLike = Union{Symbol,<:AbstractString}

"""```
lists(wallet::AbstractWallet)
```
Return all the lists on the wallet.
"""
function lists(wallet::AbstractWallet)
    missing_api("lists", wallet)
end

"""```
list(wallet::AbstractWallet, tag::StringLike)
```
Return specific list from `wallet` corresponding to `tag`.
"""
function list(wallet::AbstractWallet, tag::StringLike)
    return getindex(lists(wallet), tag)
end

"""```
enlist!(wallet::AbstractWallet, target::AbstractWalletOrAddress, tag)
```
Add `target` to the `tag` list on the given `wallet`.
"""
function enlist!(wallet::AbstractWallet, target::AbstractWalletOrAddress, tag::StringLike)
    missing_api("enlist", wallet, target, tag)
end

"""```
swaplist!(wallet::AbstractWallet, target::AbstractWalletOrAddress, tag1, tag2)
```
Swap `target` from `tag1` list to `tag2` list in `wallet`.
"""
function swaplist!(
    wallet::AbstractWallet,
    target::AbstractWalletOrAddress,
    tag1::StringLike = :black,
    tag2::StringLike = :white,
)
    if onlist(wallet, target, tag1)
        delist!(wallet, target, tag1)
        enlist!(wallet, target, tag2)
    elseif onlist(wallet, target, tag2)
        delist!(wallet, target, tag2)
        enlist!(wallet, target, tag1)
    else
        throw(ArgumentError("`$target` was not on `$tag1` or `$tag2` list in `$wallet`"))
    end
end

"""```
onlist(wallet::AbstractWallet, target::AbstractWalletOrAddress, tag)
```
Check if `target` is on the `tag` list on the given `wallet`.
"""
function onlist(wallet::AbstractWallet, target::AbstractWalletOrAddress, tag::StringLike)::Bool
    return target in list(wallet, tag)
end

"""```
delist!(wallet::AbstractWallet, target::AbstractWalletOrAddress, tag::StringLike)
```
Remove `target` from `tag` list on the `wallet`.
"""
function delist!(wallet::AbstractWallet, target::AbstractWalletOrAddress, tag::StringLike)
    missing_api("delist!", wallet, target, tag)
end

"""```
delist!(wallet::AbstractWallet, target::AbstractWalletOrAddress)
```
Remove `target` from any list on the `wallet`.
"""
function delist!(wallet::AbstractWallet, target::AbstractWalletOrAddress)
    missing_api("delist!", wallet, target)
end

"""```
addresses(wallet::AbstractWallet)
```
Return `addresses` of `wallet` for depositing.
"""
function addresses(wallet::AbstractWallet)
    return list(wallet, :addresses)
end

"""```
address(wallet::AbstractWallet{K}, asset::K)
```
Return the `address` of `wallet` for depositing the given `asset`.
"""
function address(wallet::AbstractWallet{K}, asset::K) where {K}
    return getindex(addresses(wallet), asset)
end

"""```
whitelist(wallet::AbstractAddress)
```
Return current whitelist for the given `wallet`.
"""
function whitelist(wallet::AbstractWallet)
    return list(wallet, :white)
end

"""```
whitelist!(wallet::AbstractWallet, target::AbstractWalletOrAddress)
```
Add `target` to the `wallet` whitelist.
"""
function whitelist!(wallet::AbstractWallet, target::AbstractWalletOrAddress)
    return enlist!(wallet, target, :white)
end

"""```
onwhitelist(wallet::AbstractWallet, target::AbstractWalletOrAddress)
```
Check if `target` is on the `wallet` whitelist.
"""
function onwhitelist(wallet::AbstractWallet, target::AbstractWalletOrAddress)::Bool
    return onlist(wallet, target, :white)
end

"""```
blacklist(wallet::AbstractAddress)
```
Return current blacklist for the given `wallet`.
"""
function blacklist(wallet::AbstractWallet)
    return list(wallet, :black)
end

"""```
blacklist!(wallet::AbstractWallet, target::AbstractWalletOrAddress)
```
Add `target` to the `wallet` blacklist.
"""
function blacklist!(wallet::AbstractWallet, target::AbstractWalletOrAddress)
    return enlist!(wallet, target, :black)
end

"""```
onblacklist(wallet::AbstractWallet, target::AbstractWalletOrAddress)
```
Check if `target` is on the `wallet` blacklist.
"""
function onblacklist(wallet::AbstractWallet, target::AbstractWalletOrAddress)::Bool
    return onlist(wallet, target, :black)
end


text:
using Base: sign_mask
using ZigZagBoomerang, Distributions, ForwardDiff, LinearAlgebra, SparseArrays, StructArrays
using StatsBase, JLD2
const ZZB = ZigZagBoomerang
using Test, Random
using BenchmarkTools
using ForwardDiff: Dual, value

StatsBase.cov2cor(C) = cov2cor(C, sqrt.(diag(C)))
cd(@__DIR__)
outer(x) = x*x'
d = 5
Random.seed!(3)
μ1 = 0.7randn(d)
μ2 = 0.7randn(d)
μ = μ1
Σ1 = 0.1*outer(randn(d,d))
Σ2 = 0.2*outer(randn(d,d))
Γ = sparse(inv(Σ1))

T = 8000.0
ℓ_(x, P)  = log(0.5pdf(MvNormal(P.μ1, P.Σ1), x)+0.5pdf(MvNormal(P.μ2, P.Σ2), x)) #+ randn()
ℓ(x, P) = ℓ_(x, P)

if true
    T = 4000.0

    Γ = sparse(inv( [Σ1 0I 0I;0I Σ1 0I;0I 0I Σ1]))
    d = size(Γ, 1)
    ℓ(x, P) = ℓ_(x[1:5], P) + ℓ_(x[6:10], P) +  ℓ_(x[11:15], P)
    μ = [μ1; μ1; μ1]
end

function negpartiali(f, d)
    id = collect(I(d))
    ith = [id[:,i] for i in 1:d]
    function (x, i, args...)
        sa = StructArray{Dual{}}((x, ith[i]))
        δ = -f(sa, args...).partials[]
        return δ
    end
end

function gradientsi(f, d)
    id = collect(1.0*I(d))
    ith = [id[:,i] for i in 1:d]
    z = zeros(d)
    function (t, x, θ, i, t′, F, args...)
       # sa =  (reinterpret(Dual{:Sion, Dual{:Confu, Float64, 1}, 1}, StructVector{Tuple{Float64, Float64, Float64, Float64}}((x, θ, ith[i], z))))
        sa = reinterpret(Dual{:Sion, Dual{:Confu, Float64, 1}, 1}, collect(zip(x, θ, ith[i], z)))
        δ = -f(sa, args...)
        if isnan(δ.partials[].value)
            @show t, x, θ, i, t′
            error("nan found")
        end
        return δ.partials[].value, (δ.partials[].partials[])*θ[i]
    end
end
∇ϕi = negpartiali(ℓ, d)
∇ϕi2 = gradientsi(ℓ, d)
function gradients(ℓ, d)
    function (y, t, x, θ, args...)
        x_ = x + Dual{:Bz0iEasFfw}(0.0, 1.0)*θ
        y_ = ForwardDiff.gradient(x->-ℓ(x, args...), x_)
        y .= value.(y_)
        y, dot(θ, y_).partials[]
    end
end

∇ϕ2! = gradients(ℓ, d)

t0 = 0.0
x0 = 1.01μ
θ0 = 1.0*ones(d)
c = .1

P = (;μ1=μ1, μ2=μ2, Σ1=Σ1 , Σ2=Σ2)

if true
    Z = ZigZag(Γ, μ) # parameters Γ, μ ignored
    #Z = ZigZag(∅, ∅, ∅, 0.0, 0.0, 1.0)
    if d > 10
        trace, final, (acc, num), cs = @time spdmp(∇ϕi2, t0, x0, θ0, T, ZZB.LocalBound(c*ones(d)),
            Z, P; adapt=true, progress=true)
    else
        trace, final, (acc, num), cs = @time pdmp(∇ϕi2, t0, x0, θ0, T, ZZB.LocalBound(c*ones(d)),
            Z, P; adapt=true, progress=true)
    end
else
    Z = BouncyParticle(Γ, μ, 1.0, 0.8, nothing) # parameters Γ, μ ignored
    trace, final, (acc, num), cs = @time pdmp(∇ϕ2!, t0, x0, θ0, T, ZZB.LocalBound(c),
        Z, P ; adapt=true, progress=true)
end
cs isa Vector ? @show(summarystats(cs)) : @show(cs)
t, x = ZigZagBoomerang.sep(trace)
r = eachindex(t)[end-100:end-1]
r = eachindex(t)[2:end-1]
cx = cummean(trace)

using GLMakie
if false
fig2 = Figure(resolution=(2000,2000))
e = 4
for i in 1:min(d, e^2)
    u = CartesianIndices((e,e))[i]
    lines(fig2[u[1],u[2]], t[r], getindex.(x, i)[r])
    if cs isa Vector
        lines!(fig2[u[1],u[2]], (cx[i].t), (cx[i].y), linewidth=3.0, color=:orange)
    end
end
display(fig2)
end
if true
    fig3 = Figure(resolution=(2000,2000))
    e = min(d, 5)
    for i in 1:min(d^2, e^2)
        u = CartesianIndices((e,e))[i]
        if u[1] == u[2]
            lines(fig3[u[1],u[2]], t[r], getindex.(x, u[1])[r])
            lines!(fig3[u[1],u[2]], t[r], fill(μ1[u[1]], length(r)), color=:green)
            lines!(fig3[u[1],u[2]], t[r], fill(μ2[u[1]], length(r)), color=:green)

        else
            lines(fig3[u[1],u[2]], getindex.(x, u[1])[r],  getindex.(x, u[2])[r], color=(:black,0.02))
            scatter!(fig3[u[1],u[2]], [
                    μ1[u[1]]
                    μ2[u[1]]
                    ],[
                    μ1[u[2]]
                    μ2[u[2]]
                    ], markersize=5.0, color=:green)
        end
    end
    display(fig3)
end


text:
@testset "RSwM2" begin

using DiffEqNoiseProcess

W = WienerProcess(0.0,0.0,rswm=RSWM(adaptivealg=:RSwM2))

dt = 0.2
calculate_step!(W,dt,nothing,nothing)
W.dt,W.curt,W.curW,W.dW

reject_step!(W,0.1,nothing,nothing)
W.dt,W.curt,W.curW,W.dW
DiffEqNoiseProcess.first(W.S₁)

accept_step!(W,0.2,nothing,nothing)
W.dt,W.curt,W.curW,W.dW

reject_step!(W,0.1,nothing,nothing)
W.dt,W.curt,W.curW,W.dW
DiffEqNoiseProcess.first(W.S₁)

for i in 1:100
  reject_step!(W,0.05,nothing,nothing)
  accept_step!(W,0.1,nothing,nothing)
end

dt = dt/100
for t in dt:dt:W.t[end]-dt
  W(t)
end

end


text:
export
        unchecked_value,
		set_unchecked_value,
        go_to_current_position,
        go_to,
        link_to

Base.next( sysc::SysCoordinates ) =
    ccall( (:syscoord_Next, LIB_SMILE), Void, (Ptr{Void},), sysc.ptr )

Base.getindex(sysc::SysCoordinates, index::Integer) =
    ccall( (:syscoord_GetIndex, LIB_SMILE), Cint, (Ptr{Void},Cint), sysc.ptr, index )
Base.setindex!(sysc::SysCoordinates, value::Integer, index::Integer) =
    ccall( (:syscoord_SetIndex, LIB_SMILE), Void, (Ptr{Void},Cint,Cint), sysc.ptr, index, value )

unchecked_value(sysc::SysCoordinates) =
    ccall( (:syscoord_UncheckedValue, LIB_SMILE), Cdouble, (Ptr{Void},), sysc.ptr )

set_unchecked_value( sysc::SysCoordinates, val::Real ) =
	ccall( (:syscoord_SetUncheckedValue, LIB_SMILE), Void, (Ptr{Void},Float64), sysc.ptr, val )

go_to(sysc::SysCoordinates, theIndex::Integer) =
    ccall( (:syscoord_GoTo, LIB_SMILE), Cint, (Ptr{Void},Cint), sysc.ptr, theIndex)
go_to_current_position(sysc::SysCoordinates) =
    ccall( (:syscoord_GoToCurrentPosition, LIB_SMILE), Void, (Ptr{Void},), sysc.ptr)

link_to(sysc::SysCoordinates, dmat::DMatrix) =
    ccall( (:syscoord_LinkTo_DMatrix, LIB_SMILE), Void, (Ptr{Void},Ptr{Void}), sysc.ptr, dmat.ptr)
link_to(sysc::SysCoordinates, nodedef::NodeDefinition) =
    ccall( (:syscoord_LinkTo_NodeDefinition, LIB_SMILE), Void, (Ptr{Void},Ptr{Void}), sysc.ptr, nodedef.ptr)
link_to(sysc::SysCoordinates, nodeval::NodeValue) =
    ccall( (:syscoord_LinkTo_NodeValue, LIB_SMILE), Void, (Ptr{Void},Ptr{Void}), sysc.ptr, nodeval.ptr)