text:
module SubglacialPlumes

using Parameters, DifferentialEquations

export Params, solve!, parametrize!

#Abstract types
abstract type AbstractModel{T <: Real, N<: Integer} end

"""
Struct to hold model parameters.
#Format: fieldname::Type = default_value.
#T & N are type parameters, usually real numbers (e.g. Float64) and integers
#(e.g. Int64) respectively. Variable names and default values from Jenkins 2011.
Also store the specified ambient temperature and salinity and the initial geometry here.
"""
@with_kw struct Params{T <: Real}
E0::T = 3.6e-2 #entrainment coefficient
Cd::T = 2.5e-3 #drag coefficient
St::T = 5.9e-4 #Stanton number
λ1::T = -5.73e-2 #Seawater freezing point slope
λ2::T = 8.32e-2 #Seawater freezing point offset
λ3::T = 7.61e-4 #depth dependence freezing point
L::T  = 3.35e5 #Latent heat of fusion for ice
ci::T  = 2.009e3 #Specific heat capacity for ice
cw::T = 3.974e3 #Specific heat capacity for water
βs::T = 7.86e-4 #haline contraction coefficient
βt::T = 3.87e-5 #thermal expansion coefficnet
g::T  = 9.81    #gravitational acceleration
StT::T= 1.1e-3  #Thermal Stanton number
StS::T= 3.1e-5  #Haline Stanton number
secs_per_year:: T = 365.25*24*60*60
Si::T = 0.0       #salinity of ice
Ti::T = 0.0     #tempatature of ice
ρ0::T = 1000.0    #reference density of water
geometry::Array{T,2}  = linear_geometry(0.01, -1000) #default geometry is one with slope o.01 and grounding line depth -1000
ambient_grid::Array{T,1} = Array(range(geometry[2,1], stop = 0, length = 100)); @assert ~any(diff(diff(ambient_grid)) .> 1e-5) #default ambient grid is linearly space between grounding line depth and surface
#assert regular grid in the ambient salinity
Sa::Array{T,1} = 34.6 * ones(size(ambient_grid)); @assert size(Sa) == size(ambient_grid)#default initial ambient is 34.6 everywhere
Ta::Array{T,1}  = 0.5  *ones(size(ambient_grid)); @assert size(Ta) == size(ambient_grid)#default initial ambient is 0.5C everywhere
end

"""
Structure for holding variables defined on the (user-specified) grid and grid properties
"""
@with_kw struct Grid{R <: Real, N <: Integer}
n::N #number of grid points
geometry_length::R #arc length of the geometry
x::Array{R,1}; @assert size(x) == (n,) #x-coordinates of geometry
z::Array{R,1}; @assert size(z) == (n,) #z-coordinates of geometry
dz::Array{R,1}; @assert size(dz) == (n,) #slope of geometry (dZb/dx)
s::Array{R,1}; @assert size(s) == (n,) #arc length parameter at grid points
U::Array{R,1} = zeros(n,); @assert size(U) == (n,)
D::Array{R,1} = zeros(n,); @assert size(D) == (n,)
S::Array{R,1} = zeros(n,); @assert size(S) == (n,)
T::Array{R,1} = zeros(n,); @assert size(T) == (n,)
Δρ::Array{R,1} = zeros(n,); @assert size(Δρ) == (n,)
ΔT::Array{R,1} = zeros(n,); @assert size(ΔT) == (n,)
m::Array{R,1} = zeros(n,); @assert size(m) == (n,)
mparam::Array{R,1} = zeros(n,); @assert size(m) == (n,) #store a parametrization of the melt rate
end

"""
structure for storing useful quantities computed from the input parameters
"""
@with_kw struct Store{R <: Real}
    zgl::R #depth of the grounding line
    v0::Array{R,1} = zeros(4,) #store initial conditions
    tau::R #thermal driving at the grounding line
    l0::R #z -length scale associated with freezing point dependence
    dρa_dz::Array{R,1}; #ambient density gradient, defined on ambient grid
end

"""
Structure to hold each of the above structures when plume is instantiated
"""
@with_kw struct State{T <: Real, N<: Integer} <: AbstractModel{T,N}
params::Params{T}
grid::Grid{T,N}
store::Store{T}
end

#useful functions

#define for testing
f(x,y) = 2x + 3y

#local freezing point of water
temp_freezing(S,Z,params) = params.λ1 * S + params.λ2 + params.λ3 * Z

#ambient freezing point
temp_af(Sa, Zb,params) = temp_freezing(Sa, Zb,params)

#meltwater freezing point
temp_if(Zb,params) = temp_freezing(params.Si, Zb, params)

#effective meltwater temperature
temp_ief(Z, params) = temp_freezing(params.Si,Z,params) - params.L/params.cw - params.ci*(temp_freezing(params.Si,Z,params) - params.Ti)/params.cw

#density contrast (linear equation of state)
density_contrast(S,Sa, T, Ta, params) = params.ρ0 * (params.βs * (Sa - S) - params.βt*(Ta - T))

#effective meltwater density contrast
density_contrast_effective(Sa, Ta, Z,params) = density_contrast(params.Si, Sa,  temp_ief(Z,params), Ta, params)

#melt rate factor
m_naught(S,Z,params) = params.St / (temp_freezing(S,Z,params) - temp_ief(Z,params))
m_naught(params) =  params.cw * params.St / params.L #approximation using L/c >> Tf - Tfi

#return the length of the geometry using linear interpolation
get_geometry_length(geometry) = sum( sqrt.(diff(geometry[1,:]).^2 .+ diff(geometry[2,:]).^2))

linear_geometry(α, zgl; n = 1000) = vcat(Array(range(0, stop = abs(zgl)/α, length = n))',Array( zgl .+ α* range(0, stop = abs(zgl)/α, length = n))') #initial geoemetry of the configuration. Default is linear with GL depth -1000 and slope 0.01

get_tau(Ta0,Sa0, zgl, params) = Ta0 - temp_freezing(Sa0, zgl, params)

#create plume state from input parameters
function start(params)
    #accepts input Z = Z_b(X) in form geometry = [Xb; Zb(Xb)].
    # Xb(1) specifies grounding line position, Zb(1) specifies grounding line position
    grid = Grid(
        n = size(params.geometry)[2],
        geometry_length = get_geometry_length(params.geometry),
        x = params.geometry[1,:],
        z = params.geometry[2,:],
        s = get_arclength(params.geometry),
        dz = get_slope(params.geometry))

    store = Store(
        zgl = params.geometry[2,1],
        tau = get_tau(params.Ta[1], params.Sa[1], params.geometry[2,1], params),
        l0  = get_tau(params.Ta[1], params.Sa[1], params.geometry[2,1], params)/params.λ3,  #z lengthscale associated with freezing point dependence. Uses slope at first grid point as angle scale
        dρa_dz = get_ambient_density_gradient(params.Sa, params.Ta, params)
    )
    plume=State(params, grid, store)
    return plume
end

"""
Returns the gradient in ambient density at grid points. Computed using centered finite differences, except for the endpoint, which use one sided FD.
"""
function get_ambient_density_gradient(Sa, Ta, params)
    #compute the associated ambient density using linear equation of state
    dz = diff(params.ambient_grid)[1] #we use a regular grid
    ρa = zeros(size(Sa))
    dρa_dz = zeros(size(Sa))
    @. ρa = params.ρ0 * (params.βs * Sa - params.βt * Ta)
    dρa_dz[2:end-1] = (-ρa[1:end-2] + ρa[3:end])/2/dz
    dρa_dz[1] = (-3/2 * ρa[1] + 2*ρa[2] -1/2 *ρa[3])/dz
    dρa_dz[end] = (1/2* ρa[end-2] - 2*ρa[end-1] + 3/2* ρa[end])/dz
    return dρa_dz
end

"""
returns the slope Z_b'(X) at grid points X. Computed using centered finite differences, except for the end points which use one-sided differences
"""
function get_slope(geometry)
    dx = diff(geometry[1,:])[1] #we use a regular grid
    slope = zeros(1,size(geometry)[2])
    slope[2:end-1] = (-geometry[2,1:end-2] + geometry[2,3:end])/2 /dx

    #one sided for all grid points before final
    slope[1:end-1] = (geometry[2,2:end] - geometry[2,1:end-1])/ dx

    #one sided fd for slope at both ends using 23 point stencil
    #slope[1] = (-3/2 * geometry[2,1] + 2*geometry[2,2] -1/2 *geometry[2,3])/dx
    slope[end] =(1/2* geometry[2,end-2] - 2*geometry[2,end-1] + 3/2* geometry[2,end])/dx
    return slope[1,:]
end

"""
Returns the arc length of the curve specified by the geometry
"""
function get_arclength(geometry)
    #returns the arc length at grid points
    s = zeros(1,size(geometry)[2])
    section_lengths = sqrt.(diff(geometry[1,:]).^2 .+ diff(geometry[2,:]).^2)
    s[2:end] =  [sum(section_lengths[1:i]) for i in 1:length(s)-1] #take sum over all elements with index below in section_lengths
    return s[1,:]
end

"""
Returns the value of variable Var (defined on the same grid as S) at the local position s
"""
function get_local_var(s_local, s, var)

    mxval, mxindx= findmax(-broadcast(abs, s .-s_local); dims = 1); #returns index of grid point with haf closest to zeros
    mxindx = mxindx[] #remove array packaging
    if mxindx > length(s) - 1 #we may encounter this case close to the front
        local_var = var[end]
        return local_var
    end

    if (mxindx == 1) && s_local < s[mxindx] #if local point is smaller than all grid points in s
        local_var = var[1]
    end

    if s[mxindx] > s_local #s_local is smaller than nearest grid point, so s_local ∈ (s[mxindx -1], s[mxindx])
        local_var = ((s[mxindx] - s_local)*var[mxindx-1] + (s_local - s[mxindx .- 1])*var[mxindx])/(s[mxindx] - s[mxindx - 1])
    else #s_local is largest than s at nearest grid point so s_local ∈ (s[mxindx], s[mxindx +1])
        local_var = ((s[mxindx+1] - s_local)*var[mxindx] + (s_local - s[mxindx])*var[mxindx+1])/(s[mxindx+1] - s[mxindx])
    end
    return local_var
end

"""
Converts a value of Δρ, ΔT, Z, Sa(Z), Ta(Z) to the corresponding values of salinity S and temperature T
"""
function convert_to_st(Δρ, ΔT, Z, Sa, Ta, params)
    #convert a value of Δρ and ΔT to T, S
    S = params.βs * Sa - Δρ./params.ρ0 - params.βt * (Ta - ΔT - params.λ2 - params.λ3 * Z)
    S = S/(params.βs - params.λ1 * params.βt)
    T = ΔT + temp_freezing(S, Z, params)
    return T, S
end

################### ODE Functions #############################################
"""
Returns the right hand side of the plume equations (equations 14-17 in Jenkins 2011)
"""
function get_rhs(du, v, plume, s)
    #unpack parameters
    D = v[1]
    U = v[2]
    Δρ = v[3]
    ΔT = v[4]
    @unpack params, grid, store = plume

    #recover salinity and temperature
    #M0 = ..
    α = get_local_var(s, grid.s, grid.dz)
    Z = get_local_var(s, grid.s, grid.z)
    Sa = get_local_var(Z, params.ambient_grid, params.Sa)
    Ta = get_local_var(Z, params.ambient_grid, params.Ta)
    dρa_dz = get_local_var(Z, params.ambient_grid, store.dρa_dz)
    M0 = m_naught(params) #approximation using L/c >> Tf - Tfi

    du[1] = params.E0*U*α + M0*U*ΔT
    du[2] = params.g*D*Δρ*α/params.ρ0 - params.Cd*U^2
    du[3] = (M0*U*ΔT*density_contrast_effective(Sa, Ta, Z,params)  +
            α * dρa_dz * D * U);
    du[4] = ((Ta - temp_af(Sa, Z, params))*params.E0*α*U +
            (temp_ief(Z, params) - temp_freezing(params.Si, Z, params))*M0*U*ΔT -
            params.λ3*α*D*U)
end

"""
Specifies the mass matrix A, where plume equations are specified by A * dv/dt = f (v = D, U, Δρ, ΔT)
"""
function update_func(A, v, plume, s)
    #unpack parameters
    D = v[1]
    U = v[2]
    Δρ = v[3]
    ΔT = v[4]

    A[1,1] = U
    A[1,2] = D
    A[1,3] = 0
    A[1,4] = 0
    A[2,1] = U^2
    A[2,2] = 2*D*U
    A[2,3] = 0
    A[2,4] = 0
    A[3,1] = U*Δρ
    A[3,2] = D*Δρ
    A[3,3] = D*U
    A[3,4] = 0
    A[4,1] = U*ΔT
    A[4,2] = D*ΔT
    A[4,3] = 0
    A[4,4] = D*U
end

"""
    Update the initial conditions stored in grid.
    Overload this method to use initial conditions other than zero flux initial conditions.
"""
function update_ic!(plume::AbstractModel; xc = 10)
    @unpack params, grid, store = plume
    αlocal = grid.dz[1]
    E  = params.E0 * αlocal
    Ltilde = params.L + params.ci*(temp_freezing(params.Si,store.zgl,params) - params.Ti)

    #get the ambient temperature and salnity at the grounding line
    Sa = get_local_var(store.zgl, params.ambient_grid, params.Sa)
    Ta = get_local_var(store.zgl, params.ambient_grid, params.Ta)

    v0 = zeros(4,) #explicitly define to ensure size compatibility
    v0[4] = E/(E + params.St) * store.tau;
    v0[3] = params.St * params.cw  * v0[4] * density_contrast_effective(Sa, Ta, store.zgl,params) / E /Ltilde
    v0[2] = (E * v0[3] * αlocal * params.g / params.ρ0 /(2*E + (3/2*params.Cd)))^(1/2) * xc^(1/2)
    v0[1] = 2/3 * E * xc
    store.v0 .= v0
    return plume
end

"""
    Solve the plume equations. Variables stored on grid are modified. solution span set span sspan: sovles on the whole s domain if there are no points at which dz goes negative, else solves until the first point where dz goes negative.
"""
function solve!(plume; sspan = (0, (any(plume.grid.dz .< 1e-5) ? plume.grid.s[findfirst(plume.grid.dz .< 1e-5)] : plume.grid.s[end])))
    @unpack params, grid, store = plume

    #endow plume with an initial condition (overload to change from no-discharge ic)
    update_ic!(plume)

    #set up differential equation
    M = DiffEqArrayOperator(ones(4,4), update_func = update_func)
    prob = ODEProblem(ODEFunction(get_rhs, mass_matrix = M), store.v0, sspan, plume)

    #set condition to trigger integration termination when velocity small
    condition(v,t, integrator) = v[2] - eps()
    affect!(integrator) = terminate!(integrator)
    cb = ContinuousCallback(condition, affect!)

    #note that the below does not account for the case where the plume terminates before edge of ice shelf
    sol = solve(prob, Rodas4(), callback = cb, saveat = plume.grid.s)
    nout = length(sol[1,:]) #number of cells solution fills
    grid.D[1:nout] .= sol[1,:]
    grid.U[1:nout] .= sol[2,:]
    grid.Δρ[1:nout] .= sol[3,:]
    grid.ΔT[1:nout] .= sol[4,:]

    #find corresponding salinity and temperature (no vectorization at params not vectorized)
    #for i = 1:length(grid.n)
    #    grid.T[i], grid.S[i] = SubglacialPlumes.convert_to_st(plume.grid.Δρ[i], plume.grid.ΔT[i], plume.grid.z[i], plume.grid.Sa[i], plume.grid.Ta[i], params)
    #end

    update_melt_rates!(plume)
    return plume
end

"""
Updates the melt rates stored in the grid structure to correspond to the values of U, ΔT stored in same structure
"""
function update_melt_rates!(plume)
    @unpack params, grid = plume
    M0 = m_naught(params) #approximation using L/c >> Tf - Tfi
    grid.m .= params.secs_per_year * M0.* grid.U .* grid.ΔT
    return plume
end

"""
    Update the mparam field to include a parametrization of the melt rate.
    Second argument specifies the method of paramtrization (default is Lazeroms w/out ad hoc geometric dependence)
"""
function parametrize!(plume; method = "Lazeroms")
    if method == "Lazeroms"
        parametrization_lazeroms!(plume)
    elseif method == "Lazeroms_AHG"
        parametrization_lazeroms_AHG!(plume)
    else
        error("Parametrization method must be one of following: Lazeroms, Lazeroms_AHG")
    end

end

"""
    Modify the melt rate parametrization field to the standard Lazeroms parametrization
"""
function parametrization_lazeroms!(plume)
    @unpack params, grid, store = plume
    ΔTscale   = params.E0 *grid.dz[1] * store.tau/params.St;
    Uscale = sqrt(params.βs * grid.Sa[1] * params.g * store.l0 * store.tau * params.E0 * grid.dz[1]/(params.L/params.cw) / params.Cd);
    snd = grid.s./(store.l0 / grid.dz[1]) #dimensionless version of arclength parameter s

    #need to add control for snd > 1
    M0 = m_naught(params) #approximation using L/c >> Tf - Tf
    @. grid.mparam = params.secs_per_year * M0 * Uscale *  ΔTscale * (3*(1 - snd)^(4/3) - 1).*(1 - (1 - snd)^(4/3))^(1/2) /2 /sqrt(2)
    return plume
end

"""
    Modify the melt rate parametrization field to the Lazeroms ad-hoc geometry parametrization
"""
function parametrization_lazeroms_AHG!(plume)
    @unpack params, grid, store = plume
    ΔTscale = zeros(grid.n,)
    Uscale  = zeros(grid.n,)
    @. ΔTscale = params.E0 *grid.dz * store.tau/params.St; #vector: different angle in scale at each grid point
    @. Uscale = sqrt(params.βs * grid.Sa[1] * params.g * store.l0 * store.tau * params.E0 * grid.dz/(params.L/params.cw) / params.Cd);
    snd = grid.s./(store.l0 / grid.dz[1]) #dimensionless version of arclength parameter s
    #need to add control for snd > 1
    M0 = m_naught(params) #approximation using L/c >> Tf - Tf
    @. grid.mparam = params.secs_per_year * M0 * Uscale *  ΔTscale * (3*(1 - snd)^(4/3) - 1).*(1 - (1 - snd)^(4/3))^(1/2) /2 /sqrt(2)
    return plume
end

end # module

max_stars_repo_name: alextbradley/SubglacialPlumes

max_issues_repo_name: alextbradley/SubglacialPlumes

max_forks_repo_name: alextbradley/SubglacialPlumes


text:
Base.@kwdef struct List
    entries::Observable{SimpleList}
    value::Observable{String}
    selected::Observable{Union{Int, Nothing}}=Observable{Union{Int, Nothing}}(nothing) # 0-based indexing
    hidden::Observable{Bool}=Observable(true)
    keydown::Observable{String}=Observable("")
end

function List(keys::Observable{Vector{String}}, value::Observable{String}; kwargs...)
    entries = lift(to_entries, keys)
    return List(; entries, value, kwargs...)
end

function List(entries::Observable{SimpleList}, value::Observable{String}; kwargs...)
    return List(; entries, value, kwargs...)
end

function to_entries(keys::AbstractVector{<:AbstractString})
    return Any[SimpleDict("key" => String(key), "value" => String(key)) for key in keys]
end

function jsrender(session::Session, l::List)

    activeClasses = ["text-gray-900", "bg-gray-200"]
    inactiveClasses = ["text-gray-700"]
    fixedClasses = ["cursor-pointer", "px-4", "py-2", "bg-white"]
    itemClasses = vcat(fixedClasses, inactiveClasses, "hover:" .* activeClasses)

    hidden, selected = l.hidden, l.selected

    list = DOM.ul(;
        class="border-2 border-gray-200 max-h-64",
        role="menu",
        style="position: absolute; left:0; right:0; top:0.5rem; overflow-y: scroll;"
    )

    onjs(session, l.entries, js"""
        function (entries) {
            const keys = entries.map(e => e.key);
            const values = entries.map(e => e.value);
            const list = $(list);
            const classes = $(itemClasses);
            while (list.childNodes.length > keys.length) {
                list.removeChild(list.lastChild);
            }
            for (let i = 0; i < keys.length; i++) {
                if (i >= list.children.length) {
                    const node = document.createElement("li");
                    node.classList.add("cursor-pointer");
                    node.role = "menuitem";
                    node.tabIndex = -1;
                    node.style.display = "block";
                    node.onclick = function (event) {
                        JSServe.update_obs($(l.value), event.target.dataset.value);
                    };
                    node.classList.add(...classes);
                    list.appendChild(node);
                }
                const child = list.children[i];
                child.innerText = keys[i];
                child.dataset.value = values[i];
            }
            JSServe.update_obs($(selected), null);
        }
    """)

    # Add behavior when user presses key
    onjs(session, l.keydown, js"""
        function(key) {
            const selected = JSServe.get_observable($(selected));
            const children = $(list).children;
            const len = children.length;
            if (key == "ArrowDown") {
                JSServe.update_obs($(selected), $(UtilitiesJS).cycle(selected, len, 1))
            } else if (key == "ArrowUp") {
                JSServe.update_obs($(selected), $(UtilitiesJS).cycle(selected, len, -1))
            } else if (key == "Enter" || key == "Tab") {
                if (selected !== null && len > 0) {
                    const child = children[selected];
                    JSServe.update_obs($(l.value), child.dataset.value);
                }
            }
        }
    """)

    # Style selected options
    onjs(session, selected, js"idx => $(UtilitiesJS).styleSelected($(list).children, idx, $activeClasses, $inactiveClasses)")
    # When out of focus, unselect
    onjs(session, hidden, js"hidden => hidden && JSServe.update_obs($(selected), null)")

    notify!(l.entries)

    ui = DOM.div(list; style="position: relative; z-index: 1;", hidden)
    return jsrender(session, ui)
end

ext: jl


text:
using Test
using ConvDiffMIPDECO
using jInv.Mesh
using LinearAlgebra

domain = [0 1. 0 1 0 1]
n  = [7 8 9]
M2D  = getRegularMesh(domain[1:4],n[1:2])

idg,idgn,idint = getGhostIndices(n[1:2])

xn = getNodalGrid(getPaddedMesh(M2D))

@test norm(xn[idg[1],1].- (-M2D.h[1]/2))<1e-15
@test norm(xn[idg[2],2].- (-M2D.h[2]/2))<1e-15
@test norm(xn[idg[3],1].- (1+M2D.h[1]/2))<1e-15
@test norm(xn[idg[4],2].- (1+M2D.h[2]/2))<1e-15
@test norm(xn[idgn[1],1].- (M2D.h[1]/2))<1e-15
@test norm(xn[idgn[2],2].- (M2D.h[2]/2))<1e-15
@test norm(xn[idgn[3],1].- (1-M2D.h[1]/2))<1e-15
@test norm(xn[idgn[4],2].- (1-M2D.h[2]/2))<1e-15
@test length(idgn[1])==n[2]
@test length(idgn[2])==n[1]
@test length(idgn[3])==n[2]
@test length(idgn[4])==n[1]
@test length(idg[1])==n[2]
@test length(idg[2])==n[1]
@test length(idg[3])==n[2]
@test length(idg[4])==n[1]
@test length(idint)==prod(n[1:2])

M3D  = getRegularMesh(domain,n)

idg,idgn,idint = getGhostIndices(n)

xn = getNodalGrid(getPaddedMesh(M3D))

@test norm(xn[idg[1],1].- (-M3D.h[1]/2))  < 1e-15
@test norm(xn[idg[2],2].- (-M3D.h[2]/2))  < 1e-15
@test norm(xn[idg[3],3].- (-M3D.h[3]/2))  < 1e-15
@test norm(xn[idg[4],1].- (1+M3D.h[1]/2)) < 1e-15
@test norm(xn[idg[5],2].- (1+M3D.h[2]/2)) < 1e-15
@test norm(xn[idg[6],3].- (1+M3D.h[3]/2)) < 1e-15
@test norm(xn[idgn[1],1].- (M3D.h[1]/2))  < 1e-15
@test norm(xn[idgn[2],2].- (M3D.h[2]/2))  < 1e-15
@test norm(xn[idgn[3],3].- (M3D.h[3]/2))  < 1e-15
@test norm(xn[idgn[4],1].- (1-M3D.h[1]/2))  < 1e-15
@test norm(xn[idgn[5],2].- (1-M3D.h[2]/2))  < 1e-15
@test norm(xn[idgn[6],3].- (1-M3D.h[3]/2))  < 1e-15

@test length(idgn[1])==prod(n[2:3])
@test length(idgn[2])==prod(n[[1,3]])
@test length(idgn[3])==prod(n[1:2])
@test length(idgn[4])==prod(n[2:3])
@test length(idgn[5])==prod(n[[1,3]])
@test length(idgn[6])==prod(n[1:2])

@test length(idg[1])==prod(n[2:3])
@test length(idg[2])==prod(n[[1,3]])
@test length(idg[3])==prod(n[1:2])
@test length(idg[4])==prod(n[2:3])
@test length(idg[5])==prod(n[[1,3]])
@test length(idg[6])==prod(n[1:2])

@test length(idint)==prod(n)

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
module ETLCDBReader

using StringEncodings
using Images

## Base of ETL img
abstract type ETLImg end

struct ETLMImg <: ETLImg
    etlnum :: UInt8
    sheetnum :: UInt8
    writernum :: UInt16
    serialnum :: UInt32
    charcode :: Char
    malefemale :: UInt8
    ageofwriter :: UInt8
    industryclass :: UInt8
    occupationclass :: UInt8
    data :: Array{Gray{N0f8},2}
end

struct ETLKImg <: ETLImg
    etlnum :: UInt8
    serialnum :: UInt64
    contents :: UInt8
    style :: UInt64
    charcode :: Char
    data :: Array{Gray{N0f8},2}
end

struct ETLCImg <: ETLImg
    etlnum :: UInt8
    sheetnum :: UInt8
    writernum :: UInt16
    serialnum :: UInt32
    charcode :: Char
    malefemale :: UInt8
    ageofwriter :: UInt8
    industryclass :: UInt8
    occupationclass :: UInt8
    data :: Array{Gray{N0f8},2}
end

struct ETLGImg <: ETLImg
    etlnum :: UInt8
    sheetnum :: UInt8
    writernum :: UInt16
    serialnum :: UInt32
    charcode :: Char
    malefemale :: UInt8
    ageofwriter :: UInt8
    industryclass :: UInt8
    occupationclass :: UInt8
    data :: Array{Gray{N0f8},2}
end

struct ETLBImg <: ETLImg
    etlnum :: UInt8
    serialnum :: UInt32
    charcode :: Char
    writernum :: UInt16
    data :: Array{Gray{N0f8},2}
end

const ETL = Array{<:ETLImg,1}

const ETLM = Array{ETLMImg,1}
const ETLK = Array{ETLKImg,1}
const ETLC = Array{ETLCImg,1}
const ETLG = Array{ETLGImg,1}
const ETLB = Array{ETLBImg,1}

etl9bwh = (64,63)

function bytes2uint(bytes::Array{UInt8}, T::Type{<:Unsigned})::T
    tmp::T = 0
    for b in bytes
        tmp += b
        tmp <<= 8
    end
    return tmp
end

getetlnum(img::ETLImg)::UInt8 = img.etlnum
getetlnum(imgs::ETL)::UInt8 = imgs[1].etlnum

function getetl9g(dirpath::String)::ETLG
end

function bytes2bit(bytes)
    data = Vector{UInt8}(undef, length(bytes)*8)
    i = 0
    for b in bytes
        for j in 1:8
            data[8*i+j] = (bytes[i+1] & (1 << (8-j))) == 0 ? 0 : 1
        end
        i += 1
    end
    return data
end

function showterm(img::ETLImg)
    s = size(img.data)
    for i in 1:s[2]
        for j in 1:s[1]
            print(img.data[j, i] == 0 ? "00" : "11")
        end
        print("\n")
    end
end

function getetl9b(dirpath::String, getrange = 1:0xffffffffffffffff)::ETLB
    if dirpath[end] != '/'
        dirpath = dirpath * '/'
    end
    etls = Vector{ETLBImg}()
    etlnum = 9
    etlfiles = filter(x -> occursin(r"ETL9B_[0-9]+", x), readdir(dirpath))
    writernum = 1
    serialnum = 1
    for filename in etlfiles
        open(dirpath * filename, "r") do file
            read(file, 576)
            while !eof(file)
                sheetnum = bytes2uint(read(file, 2), UInt16)
                charcode = decode(read(file, 2), "JIS_X0208")[1]
                read(file, 4)
                data = Gray.(N0f8.(reshape(bytes2bit(read(file, 504)), etl9bwh)))'
                read(file, 64)
                if minimum(getrange) <= serialnum <= maximum(getrange)
                    etlimg = ETLBImg(etlnum, serialnum, charcode, writernum, data)
                    push!(etls, etlimg)
                end

                writernum += (serialnum % 3036 == 0) ? 1 : 0
                serialnum += 1
            end
        end
    end
    return etls
end

end # module


text:
# Allocates a new output array and run [`downsample!`](@ref).
function downsample(a::AbstractMatrix, aggregator, scale)
    downsample!(initdownsample(a, scale), a, aggregator, scale)
end

# Downsamples matrix `a` to another matrix `out` of the correct size.
# - `aggregator` is a function such as `mean` or `sum` that can combine the
    # value of multiple cells to generate the downsampled cell.
# - `scale` is the downsampling factor.
function downsample!(out::AbstractMatrix, a::AbstractMatrix, aggregator, scale)
    scale == 1 && return out .= a

    h, w = size(out)
    for i = 1:h, j = 1:w
        i1, j1 = upsample_index((i, j), scale)
        i2, j2 = min.(size(a), (i1, j1) .+ scale .- 1)
        cells = collect(skipmissing(a[i1:i2, j1:j2]))
        out[i, j] = length(cells) > 0 ? aggregator(cells) : missing
    end
    return out
end

# Generate an array for downsampling array `a` by `scale`.
initdownsample(a, scale) = similar(a, downsample_index(size(a), scale))

# Convert indicies from the downsampled array to the larger original array.
upsample_index(index, scale) = ((index .- 1) .* scale) .+ 1

# Convert indicies from the original array to the downsampled array.
downsample_index(index, scale) = ((index .- 1) .÷ scale) .+ 1


text:
using Documenter
using ClusterValidityIndices
using DemoCards

# Generate the demo files
# this is the relative path to docs/
demopage, postprocess_cb, demo_assets = makedemos("examples")

assets = [
    joinpath("assets", "favicon.ico")
]

# if there are generated css assets, pass it to Documenter.HTML
isnothing(demo_assets) || (push!(assets, demo_assets))

# Make the documentation
makedocs(
    modules=[ClusterValidityIndices],
    format=Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets = assets,
    ),
    # format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Tutorial" => [
            "Guide" => "man/guide.md",
            # "Examples" => "man/examples.md",
            demopage,
            "Contributing" => "man/contributing.md",
            "Index" => "man/full-index.md"
        ]
    ],
    repo="https://github.com/AP6YC/ClusterValidityIndices.jl/blob/{commit}{path}#L{line}",
    sitename="ClusterValidityIndices.jl",
    authors="Sasha Petrenko",
    # assets=String[],
)

# 3. postprocess after makedocs
postprocess_cb()

deploydocs(
    repo="github.com/AP6YC/ClusterValidityIndices.jl.git",
    devbranch="develop",
)


text: precompile(Tuple{typeof(CG.julia_cg), Ptr{Nothing}, Ptr{Float64}, Ptr{Float64}, UInt64})

max_stars_repo_name: kmsquire/libcg

max_issues_repo_name: kmsquire/libcg

max_forks_repo_name: kmsquire/libcg


text:
using Pkg

# official registry https://pkg.julialang.org/

# activate current project
Pkg.activate(pwd())

# install missing packages
Pkg.instantiate()

# current status
Pkg.status()

# restore to default
Pkg.activate()

# managing projects
#
# I should not use this, as it seems to decide what to update based on
# #update-#downgrade, which downgrades Flux. But I can pin Flux.
#
Pkg.update()
# precompile command is availabe at PKG repl, the usage here is not very straightforward
# Pkg.precompile

LOAD_PATH
DEPOT_PATH
# ENV["JULIA_LOAD_PATH"]

Pkg.add("BSON")
Pkg.add("Images")
Pkg.add("MLDatasets")
Pkg.add("TensorBoardLogger")
Pkg.add("LoggingExtras")
Pkg.add("ProgressMeter")
Pkg.add("Plots")
Pkg.add("FileIO")
# for showing images
Pkg.add("ImageMagick")

# trying to use Flux master. The old Flux has BN gradient bug.
#
# Pkg.add("Flux")
Pkg.add(PackageSpec(name="Flux", version="0.9"))
Pkg.pin(PackageSpec(name="Flux", version="0.9"))
# Pkg.add(PackageSpec(name="Flux", rev="master"))

# FIXME has_cudnn instead of libcudnn
#
# develop Flux, and switch to 0.9 branch, do some fix
Pkg.develop("Flux")
# then "resolve" so that Tracker is added to current project's .toml. Then using
# Flux would work.
Pkg.resolve()

Pkg.add("CuArrays")

# used packages
# Pkg.add(PackageSpec(url="https://github.com/jaypmorgan/Adversarial.jl.git"))

# Pkg.add("Plots")

# FIXME Augmentor not on Julia 1.0
# Pkg.add("Augmentor")
# Pkg.add(PackageSpec(url="https://github.com/Evizero/Augmentor.jl"))

# CUDA support
#
# Pkg.add("CUDAapi")
# Pkg.add("NNlib")
# Pkg.add("CUDAnative")
# Pkg.add("CUDAdrv")

# For cutting edge Zygote, I'm not sure whether these masters are needed
#
# Pkg.add(PackageSpec(name="CUDAnative", rev="master"))
# Pkg.add(PackageSpec(name="CUDAdrv", rev="master"))
# Pkg.add(PackageSpec(name="CuArrays", rev="master"))

# Zygote is unstable
#
# Pkg.add("ZygoteRules")
# Pkg.add("Zygote")
# only the master export pullback
# Pkg.add(PackageSpec(name="ZygoteRules", rev="master"))
# Pkg.add(PackageSpec(name="Zygote", rev="master"))

# Pkg.add("Metalhead")
# Pkg.add("FixedPointNumbers")
# Pkg.add("Distances")
# Pkg.add("Distributions")
# Pkg.add("FileIO")

max_stars_repo_name: lihebi/AdvAE

max_issues_repo_name: lihebi/AdvAE

max_forks_repo_name: lihebi/AdvAE


text:
module thrusters

export output, feedback_loop

using aoc2019.computer

function output(program, phases, initial)
    ## setup computer io
    io = [Channel{Int64}(256) for phase in phases]
    ## add a channel for the thruster
    thruster = Channel{Int64}(256)
    append!(io, [thruster])
    for i in 1:length(phases)
        put!(io[i], phases[i])
    end
    ## start computers
    for i in 1:length(phases)
        @async computer.run(copy(program), io[i], io[i + 1])
    end
    ## provide initial signal
    put!(io[1], initial)
    ## read input to thruster
    return take!(thruster)
end

function feedback_loop(program, phases, initial)
    ## setup computer io
    io = [Channel{Int64}(256) for phase in phases]
    for i in 1:length(phases)
        put!(io[i], phases[i])
    end
    ## start amplifiers
    amps = []
    for i in 1:length(phases)
        amp = @async computer.run(copy(program), io[i], io[mod(i + 1, 1:length(phases))])
        push!(amps, amp)
    end
    ## provide initial signal
    put!(io[1], initial)
    ## wait for all computers to halt
    while !all([istaskdone(amp) for amp in amps]); yield(); end
    ## get any outputs -- there should only be one
    for channel in io
        if isready(channel)
            return take!(channel)
        end
    end
end

end


text:
using Plots, SNN

SNN.IF_curve(SNN.IF, 0:0.1:20)
SNN.IF_curve(SNN.IZ, 0:0.1:20)
SNN.IF_curve(SNN.HH, 1e-5*(0:0.01:4))


text:
using JLD2, Attitude, MATLAB, SatelliteDynamics, Dierckx
# include ksfunctions.jl
include(joinpath(dirname(dirname(@__FILE__)),"ks_functions.jl"))
@load "rate_control_transfers/30_v3_day_transfer.jld2" X U
dscale = 1e7 # m
tscale = 20000 # s
uscale = 10000.0
dt = 3e-2

t_days,t_hist = get_time_transfer(X,dt,tscale)

r_eci_hist = dscale*mat_from_vec([x_from_u(X[i][1:4]) for i = 1:length(X)])
v_eci_hist = (dscale/tscale)*mat_from_vec([xdot_from_u(X[i][1:4],X[i][5:8]) for i = 1:length(X)])

# function discretize_t(t_hist,n)
#     # new time vec
#     new_t = [0.0]
#     for i = 1:(length(t_hist)-1)
#
#         # time spacing betwen two knot points
#         δt = t_hist[i+1] - t_hist[i]
#
#         # add intermediate points between the knot points
#         for j = 1:n
#             push!(new_t,δt*j/n + t_hist[i])
#         end
#     end
#     return new_t
# end
new_t = discretize_t(t_hist,4)
mat"
$rnew = spline($t_hist,$r_eci_hist,$new_t)/1000;
"
mat"
rnew = spline($t_hist,$r_eci_hist,$new_t)/1000;
rnew = round(rnew,3,'Significant')
figure
hold on
plot3(rnew(1,:),rnew(2,:),rnew(3,:),'k')
[X,Y,Z] = sphere(10)
X = round(X,3,'Significant')
Y = round(Y,3,'Significant')
Z = round(Z,3,'Significant')
Re = 6010
h = surf(Re*X,Re*Y,Re*Z)
axis equal
view(-69,14)
hold off
h.EdgeAlpha = 0
h.FaceColor = 'blue'
view([58 2])
%axis equal
hold off
addpath('/Users/kevintracy/devel/Low-Thrust-TrajOpt/matlab2tikz/src')
%matlab2tikz('full_30_traj.tex')
"
# TODO: Spline this in a data efficient way
mat"
figure
hold on
plot3($r_eci_hist(1,:),$r_eci_hist(2,:),$r_eci_hist(3,:))
axis equal
hold off
"

# ANIMATION STUFF
# output_to_matlab = (new_time_vec = new_t, x_hist_t_correct = rnew)
# file = matopen("matlab_video_data.mat", "w")
# write(file, "output_to_matlab", output_to_matlab)
# close(file)
# mat"""
# [x,y,z] = sphere(20);
# new_time_vec = $new_t;
# x_hist_t_correct = $rnew;
# a_el_1 = [0 90];
# a_el_2 = [-90 0];
#
#
# f1 = figure('units','normalized','outerposition',[0 0 1 1]);
#
# iter = 0;
#
# len = length(1:200:length(x_hist_t_correct));
#
# %v = VideoWriter('newfile.avi');
# %open(v);
#
# [x,y,z] = sphere(20);
#
# $for i = 1:200:length(x_hist_t_correct)
# for i = 1:5
#     i
#     iter = iter+1;
#     figure(f1)
#     %title('delta T' = ',num2str(floor(new_time_vec(i))),' Days'])
#     hold on
#     cla
#     hold on
#
#     plot3(x_hist_t_correct(1,1:i),x_hist_t_correct(2,1:i),x_hist_t_correct(3,1:i))
#     surf(6378*x,6378*y,6378*z,'FaceColor','blue','EdgeColor', 'black')
#     plot3(x_hist_t_correct(1,i),x_hist_t_correct(2,i),x_hist_t_correct(3,i),'o')
#     axis equal
#     axis(dscale*[-6.5 4.5 -4.5 4.5 -1.8 1.8])
#     %view((i/500),10)
#     view(-60,5)
#     title('GTO-GEO Trajectory Optimization')
#     xlabel('ECI I (m)')
#     ylabel('ECI J (m)')
#     zlabel('ECI K (m)')
#     %view(0 - 90*i/(50*len),90 - 90*i/(50*len))
#
#     drawnow
#      %frame= getframe(gcf);
#      %writeVideo(v,frame)
#
#
#
# end
#
# %close(v);
#
# """

Xm = mat_from_vec(X)
U_real = Xm[12:14,:]

mat"
figure
hold on
title('Thrust Vector')
plot($t_hist,$U_real')
ylabel('Thrust Acceleration (0.0001m/s^2)')
xlabel('Time (days)')
hold off
xlim([0,$t_hist(end)])
legend('u_x','u_y','u_z')
%saveas(gcf,'long_traj.png')
"

U_real_vec = [X[i][12:14] for i = 1:length(X)]
Unorm = [norm(U_real_vec[i]) for i = 1:length(X)]

mat"
figure
hold on
plot($t_hist,$Unorm/$uscale)
xlim([$t_hist(4),$t_hist(end)])
xlabel('Time (days)')
ylabel('Thrust Magnitude (m/s^2)')
hold off
"

angle_hist = get_rtn_in_out(U_real,r_eci_hist,v_eci_hist)

angles = mat_from_vec(angle_hist)

mat"
figure
hold on
plot(round($t_hist,4,'Significant'),round(rad2deg($angles'),4,'Significant'),'linewidth',1.2)
xlabel('Time (days)')
xlim([0 $t_hist(end)])
ylabel('Angle (degrees)')
legend('In Plane','Out of Plane','Location','southwest')
hold off
addpath('/Users/kevintracy/devel/Low-Thrust-TrajOpt/matlab2tikz/src')
%matlab2tikz('thrust_angles_30.tex')
"

# picture in picture
mat"
figure
hold on
plot(round($t_hist,4,'Significant'),round(rad2deg($angles'),4,'Significant'),'linewidth',1.2)
xlabel('Time (days)')
xlim([0 $t_hist(end)])
ylabel('Angle (degrees)')
legend('In Plane','Out of Plane','Location','southwest')
a2 = axes()
a2.Position = [0.1200 0.2600 0.2 0.2]; % xlocation, ylocation, xsize, ysize
plot(a2,$t_hist(1:500),$angles(:,1:500)'); axis tight
hold off
addpath('/Users/kevintracy/devel/Low-Thrust-TrajOpt/matlab2tikz/src')
%matlab2tikz('thrust_angles_30.tex')
"

# now let's get angle from velocity vector
angle_from_v = [acos(dot(normalize(v_eci_hist[:,i]),normalize(U_real[:,i]))) for i = 1:size(U_real,2)]

mat"
figure
hold on
plot(rad2deg($angle_from_v))
hold off
"
println("Time in days:  $t_days")

# now let's get orbital elements
const SD = SatelliteDynamics
oe_hist = [SD.sCARTtoOSC([r_eci_hist[:,i];v_eci_hist[:,i]]) for i = 1:size(r_eci_hist,2)]

a_hist = [oe_hist[i][1] for i = 1:length(oe_hist)]
e_hist = [oe_hist[i][2] for i = 1:length(oe_hist)]
i_hist = [rad2deg(oe_hist[i][3]) for i = 1:length(oe_hist)]

mat"
figure
hold on
subplot(3,1,1)
hold on
plot($t_hist,$a_hist/1000)
h = plot($t_hist, 42e3*ones(length($t_hist),1),'r');
legend([h(1)],'GEO SMA','Location','SouthEast')
ylabel('Semi-Major Axis (km)')
xlabel('Time (days)')
xlim([0,$t_hist(end)])
ylim([20e3,46e3])
subplot(3,1,2)
hold on
plot($t_hist,$e_hist)
ylabel('Eccentricity')
xlabel('Time (days)')
xlim([0,$t_hist(end)])
subplot(3,1,3)
hold on
plot($t_hist,$i_hist)
ylabel('Inclination (deg)')
xlabel('Time (days)')
xlim([0,$t_hist(end)])
hold off
"

angle_hist = get_rtn_in_out(U_real,r_eci_hist,v_eci_hist)
angles = mat_from_vec(angle_hist)

thirty = (t_hist = t_hist, a_hist = a_hist, e_hist = e_hist, i_hist = i_hist, Unorm = Unorm,angles = angles)

using JLD2
@save "plotting_data/plotting_data_30.jld2" thirty

max_stars_repo_name: RoboticExplorationLab/KSLowThrust

max_issues_repo_name: RoboticExplorationLab/KSLowThrust

max_forks_repo_name: RoboticExplorationLab/KSLowThrust


text:
module PkgDev

using Pkg, LibGit2, RegistryTools, URIs
import GitHub
import PkgButlerEngine
import Base64
using UUIDs
import DocumentFormat
using FilePathsBase
import TOML

include("tag.jl")
include("pkgbutler.jl")
include("format.jl")

# remove extension .jl
const PKGEXT = ".jl"
splitjl(pkg::AbstractString) = endswith(pkg, PKGEXT) ? pkg[1:end-length(PKGEXT)] : pkg

function get_repo_onwer_from_url(pkg_url)
    startswith(pkg_url, "git@") &&  error("Only packages that use https as the git transport protocol are supported.")

    pkg_owner_repo_name = if lowercase(splitext(URI(pkg_url).path)[2])==".git"
        splitext(URI(pkg_url).path)[1][2:end]
    else
        URI(pkg_url).path[2:end]
    end

    return pkg_owner_repo_name
end

"""
    config()

Interactive configuration of the development environment.

PkgDev.jl operations require `git` minimum configuration that keeps user signature
(user.name & user.email).
"""
function config(force::Bool=false)
    # setup global git configuration
    cfg = LibGit2.GitConfig(LibGit2.Consts.CONFIG_LEVEL_GLOBAL)
    try
        println("PkgDev.jl configuration:")

        username = LibGit2.get(cfg, "user.name", "")
        if isempty(username) || force
            username = LibGit2.prompt("Enter user name", default=username)
            LibGit2.set!(cfg, "user.name", username)
        else
            println("User name: $username")
        end

        useremail = LibGit2.get(cfg, "user.email", "")
        if isempty(useremail) || force
            useremail = LibGit2.prompt("Enter user email", default=useremail)
            LibGit2.set!(cfg, "user.email", useremail)
        else
            println("User email: $useremail")
        end

        # setup github account
        ghuser = LibGit2.get(cfg, "github.user", "")
        if isempty(ghuser) || force
            ghuser = LibGit2.prompt("Enter GitHub user", default=(isempty(ghuser) ? username : ghuser))
            LibGit2.set!(cfg, "github.user", ghuser)
        else
            println("GitHub user: $ghuser")
        end
    finally
        finalize(cfg)
    end
    lowercase(LibGit2.prompt("Do you want to change this configuration?", default="N")) == "y" && config(true)
    return
end

end # module


text:
using Winston

# Each example will save a file. Put the file in a temporary directory,
# so they don't clutter up whatever directory this is being run from.
cwd = pwd()
# First we have to get the absolute path to the examples...
tls = task_local_storage()
dir = dirname(tls[:SOURCE_PATH])
direxample = normpath(joinpath(dir, "..", "examples"))
# Now switch to a temporary directory
dirdump = mktempdir()
cd(dirdump)

for i = 1:6
    include(joinpath(direxample, "example"*string(i)*".jl"))
    Winston.display(i == 3 ? t2 : p)
end

cd(cwd)


text:
"""
    abstract type PopType
Supertype for types used in PopGen.jl
"""
abstract type PopType end       # may be removed in the future

"""
    PopObj(samples::DataFrame, loci::DataFrame)
The data struct used for the PopGen population genetics ecosystem. You are
STRONGLY discouraged from manually creating dataframes to pass into a PopObj,
and instead should use the provided genepop, csv, or vcf file importers.

- `samples` ::DataFrame individual/sample data with the columns:
    - `name` ::String the individual/sample names
    - `population` ::String population names/numbers
    - `ploidy` ::Int8 ploidy in order of `ind`
    - `longitude` ::Float64 longitude values
    - `latitude` ::Float64 latitude values
- `loci` ::DataFrame loci and their genotypes
    - columns are named by loci
    - genotypes are Tuples of Int16 or Int16, arraged in order of `.samples.name`
"""
mutable struct PopObj <: PopType
    samples::DataFrame
    loci::DataFrame
    function PopObj(x::DataFrame, y::DataFrame)
        if sort(names(x)) != [:latitude, :longitude, :name, :ploidy, :population]
            error("Incorrect column names in samples dataframe.
            Columns should be: name, population, ploidy, longitude, latitude")
        end
        size(x,1) != size(y,1) && error("length mismatch of dataframes. samples: $(size(x,1)) | loci: $(size(y, 1))")
        typeof(x.name) != Array{String,1} && @error ":name values must be of type String"
        new(x,y)
    end
end

"""
    PopOpt(samples::DataFrame, loci::DataFrame)
An immutable version of a `PopObj`, used for under-the-hood processes in PopGen.jl.
"""
struct PopOpt <: PopType
    samples::DataFrame
    loci::DataFrame
    PopOpt(x::PopObj) = new(x.samples, x.loci)
end


text:
using WGANGP
using Test
using Random

include("wgangp.jl")


text:
# Multiple single neuron examples derived from the Prinz et al STG model.
# To try some of the other neuron variants, redefine the `channels` vector, using any of the
# additional commented out parameters listed below. Then re-run the neuron constructor,
# simulation, solve and plot as before.

using Conductor, OrdinaryDiffEq, Plots

include(joinpath(@__DIR__, "prinz_kinetics.jl"));

############################################################################################
# Prinz 2003 J. Neurophysiol
############################################################################################

# Figure 2C Bursting neuron
channels = [NaV(100mS/cm^2),
            CaT(0mS/cm^2),
            CaS(4mS/cm^2),
            KA(0mS/cm^2),
            KCa(15mS/cm^2),
            Kdr(50mS/cm^2),
            H(.02mS/cm^2),
            leak(.03mS/cm^2)]

@named neuron = Soma(channels, gradients, area = area, V0 = -50mV, aux = [calcium_conversion]);

t = 2500
sim = Simulation(neuron, time = t*ms)
solution = solve(sim, Rosenbrock23());

# Plot at 5kHz sampling
fig = plot(solution; plotdensity=Int(t*4), size=(1200,800), vars=[neuron.sys.Vₘ]);
gui(fig)

# Uncomment and eval `png(...)` to save as PNG
# png(fig, "figure_simulated")

# Evaluate the alternate channel vectors below to experiment with different styles of
# neurons, derived directly from the original publications...

#=
# Figure 3C Irregular bursting
channels = [NaV(400mS/cm^2),
            CaT(0mS/cm^2),
            CaS(8mS/cm^2),
            KA(50mS/cm^2),
            KCa(20mS/cm^2),
            Kdr(50mS/cm^2),
            H(.04mS/cm^2),
            leak(0mS/cm^2)]
=#

#=
# Figure 5A broad APs - qualitatively similar
channels = [NaV(200mS/cm^2),
            CaT(0mS/cm^2),
            CaS(2mS/cm^2),
            KA(0mS/cm^2),
            KCa(15mS/cm^2),
            Kdr(0mS/cm^2),
            H(.03mS/cm^2),
            leak(.04mS/cm^2)]
=#

#=
# Figure 5B "spike triplets" - doesn't reproduce; spikes repetitively
channels = [NaV(100mS/cm^2),
            CaT(0mS/cm^2),
            CaS(10mS/cm^2),
            KA(50mS/cm^2),
            KCa(10mS/cm^2),
            Kdr(50mS/cm^2),
            H(.03mS/cm^2),
            leak(.05mS/cm^2)]
=#

#=
# Figure 5C Plateau burster - doesnt reproduce,
channels = [NaV(400mS/cm^2),
            CaT(2.5mS/cm^2),
            CaS(10mS/cm^2),
            KA(20mS/cm^2),
            KCa(5mS/cm^2),
            Kdr(25mS/cm^2),
            H(.04mS/cm^2),
            leak(.03mS/cm^2)]
=#

#=
# Figure 5D (Parabolic burster) - bursts but not sure if "parabolic"
channels = [NaV(100mS/cm^2),
            CaT(0mS/cm^2),
            CaS(6mS/cm^2),
            KA(10mS/cm^2),
            KCa(10mS/cm^2),
            Kdr(50mS/cm^2),
            H(.03mS/cm^2),
            leak(.05mS/cm^2)]
=#

#=
# figure 5G - Bursting with 17-15-17-15... APs - not sure if pattern matches
channels = [NaV(500mS/cm^2),
            CaT(2.5mS/cm^2),
            CaS(8mS/cm^2),
            KA(0mS/cm^2),
            KCa(15mS/cm^2),
            Kdr(75mS/cm^2),
            H(.05mS/cm^2),
            leak(0mS/cm^2)]
=#

############################################################################################
# Prinz 2003 - J. Neuroscience
############################################################################################

#=
# bursting neuron
channels = [NaV(200mS/cm^2),
            CaT(2.5mS/cm^2),
            CaS(4mS/cm^2),
            KA(50mS/cm^2),
            KCa(5mS/cm^2),
            Kdr(100mS/cm^2),
            H(.01mS/cm^2),
            leak(.01mS/cm^2)]
=#

#=
# tonic spiking neuron
channels = [NaV(200mS/cm^2),
            CaT(0mS/cm^2),
            CaS(4mS/cm^2),
            KA(10mS/cm^2),
            KCa(10mS/cm^2),
            Kdr(125mS/cm^2),
            H(.05mS/cm^2),
            leak(.04mS/cm^2)]
=#

############################################################################################
# Prinz 2004 - Nature Neurosci. -- (all below are intrinsically bursting ABPD)
############################################################################################

#=
# AB/PD 1
channels = [NaV(400mS/cm^2),
            CaT(2.5mS/cm^2),
            CaS(6mS/cm^2),
            KA(50mS/cm^2),
            KCa(10mS/cm^2),
            Kdr(100mS/cm^2),
            H(.01mS/cm^2),
            leak(0mS/cm^2)]
=#

#=
# AB/PD 2
channels = [NaV(100mS/cm^2),
            CaT(2.5mS/cm^2),
            CaS(6mS/cm^2),
            KA(50mS/cm^2),
            KCa(5mS/cm^2),
            Kdr(100mS/cm^2),
            H(.01mS/cm^2),
            leak(0mS/cm^2)]
=#

#=
# AB/PD 3
channels = [NaV(200mS/cm^2),
            CaT(2.5mS/cm^2),
            CaS(4mS/cm^2),
            KA(50mS/cm^2),
            KCa(5mS/cm^2),
            Kdr(50mS/cm^2),
            H(.01mS/cm^2),
            leak(0mS/cm^2)]
=#

#=
# AB/PD 4
channels = [NaV(200mS/cm^2),
            CaT(5mS/cm^2),
            CaS(4mS/cm^2),
            KA(40mS/cm^2),
            KCa(5mS/cm^2),
            Kdr(125mS/cm^2),
            H(.01mS/cm^2),
            leak(0mS/cm^2)]
=#

#=
# AB/PD 5
channels = [NaV(300mS/cm^2),
            CaT(2.5mS/cm^2),
            CaS(2mS/cm^2),
            KA(10mS/cm^2),
            KCa(5mS/cm^2),
            Kdr(125mS/cm^2),
            H(.01mS/cm^2),
            leak(0mS/cm^2)]
=#


text:
function f!(resid, x)
    for i in eachindex(x)
        resid[i] = sin(x[i]) + x[i]^3
    end
end
x = ones(5)
@test Sundials.kinsol(f!, x, linear_solver = :Band, jac_upper = 0, jac_lower = 0) ==
      Sundials.kinsol(f!, x)

function f_iip(du, u, p, t)
    @. du = sin(u) + u^3
end
function f_iip(du, u, p)
    @. du = sin(u) + u^3
end
function f_oop(u, p, t)
    return @. sin(u) + u^3
end
function f_oop(u, p)
    return @. sin(u) + u^3
end

#Testing KINSOL()
u0 = ones(5)
prob_iip = SteadyStateProblem(f_iip, u0)
prob_oop = SteadyStateProblem(f_oop, u0)
@test solve(prob_iip, KINSOL()) == solve(prob_iip, KINSOL(linear_solver=:Band, jac_upper=0, jac_lower=0))
@test solve(prob_oop, KINSOL()) == solve(prob_oop, KINSOL(linear_solver=:Band, jac_upper=0, jac_lower=0))
prob_iip = NonlinearProblem(f_iip, u0)
prob_oop = NonlinearProblem(f_oop, u0)
@test solve(prob_iip, KINSOL()) == solve(prob_iip, KINSOL(linear_solver=:Band, jac_upper=0, jac_lower=0))
@test solve(prob_oop, KINSOL()) == solve(prob_oop, KINSOL(linear_solver=:Band, jac_upper=0, jac_lower=0))


text:
test = (5764801, 17807724)
puzzle = (8987316, 14681524)
modulus = 20201227
seed = 7

mods = let mods = Vector(undef, modulus - 1)
    mods[1] = seed
    for i in 2:modulus - 1
        mods[i] = (7 * mods[i - 1]) % modulus
    end
    mods
end
@assert mods[end] == 1

function q1(card, door)::Int
    loops = (findfirst(==(card), mods), findfirst(==(door), mods))
    return mods[(*(loops...) - 1) % (modulus - 1) + 1]
end

@assert q1(test...) == 14897079

println("Q1: ", q1(puzzle...))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
name = "libvorbis"
version = v"1.3.6"
# Collection of sources required to build imagemagick

sources = [
    "https://downloads.xiph.org/releases/vorbis/libvorbis-1.3.6.tar.xz" =>
    "af00bb5a784e7c9e69f56823de4637c350643deedaf333d0fa86ecdba6fcb415",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd libvorbis-1.3.6/
./configure --prefix=$prefix --host=$target
make -j${ncore}
make install
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    # Windows
    Windows(:i686),
    Windows(:x86_64),

    # linux
    Linux(:i686, :glibc),
    Linux(:x86_64, :glibc),
    Linux(:aarch64, :glibc),
    Linux(:armv7l, :glibc),
    Linux(:powerpc64le, :glibc),

    # musl
    Linux(:i686, :musl),
    Linux(:x86_64, :musl),
    Linux(:aarch64, :musl),
    Linux(:armv7l, :musl),

    # The BSD's
    FreeBSD(:x86_64),
    MacOS(:x86_64),
]

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

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaIO/OggBuilder/releases/download/v1.3.3-7/build_Ogg.v1.3.3.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: SimonDanisch/LibVorbisBuilder

max_issues_repo_name: SimonDanisch/LibVorbisBuilder

max_forks_repo_name: SimonDanisch/LibVorbisBuilder


text:
"""
Find index of location of coordinate or dipole in leadfield
"""
find_location(l, d::Union{Dipole, Coordinate}) = find_location(l::Leadfield, d.x, d.y, d.z)

function find_location(l::Leadfield, x::Number, y::Number, z::Number)

    valid_x = l.x .== x
    valid_y = l.y .== y
    valid_z = l.z .== z

    idx = find(valid_x & valid_y & valid_z)

    if isempty(idx)

        dists = [euclidean([l.x[i], l.y[i], l.z[i]], [x, y, z]) for i = 1:length(l.x)]
        idx = findfirst(dists, minimum(dists))
    else
        idx = idx[1]
    end

    return idx
end


text:
module MyTypes

import Base: ==
export MyType, MyContainer

## Objects we want to save
# data in MyType is always of length 5, and that is the basis for a more efficient serialization
struct MyType{T}
    data::Vector{T}
    id::Int

    function MyType{T}(v::Vector{T}, id::Integer) where T
        length(v) == 5 || error("All vectors must be of length 5")
        new{T}(v, id)
    end
end
MyType(v::Vector{T}, id::Integer) where {T} = MyType{T}(v, id)
Base.eltype(::Type{MyType{T}}) where {T} = T
==(a::MyType, b::MyType) = a.data == b.data && a.id == b.id

struct MyContainer{T}
    objs::Vector{MyType{T}}
end
Base.eltype(::Type{MyContainer{T}}) where {T} = T
==(a::MyContainer, b::MyContainer) = length(a.objs) == length(b.objs) && all(i->a.objs[i]==b.objs[i], 1:length(a.objs))

end # MyTypes

### Here are the definitions needed to implement the custom serialization
# If you prefer, you could include these definitions in the MyTypes module
module MySerializer

using HDF5, JLD, ..MyTypes
using Compat

## Defining the serialization format
mutable struct MyContainerSerializer{T}
    data::Matrix{T}
    ids::Vector{Int}
end
MyContainerSerializer(data::Matrix{T},ids) where {T} = MyContainerSerializer{T}(data, ids)
Base.eltype(::Type{MyContainerSerializer{T}}) where {T} = T
Base.eltype(::MyContainerSerializer{T}) where {T} = T

JLD.readas(serdata::MyContainerSerializer) =
    MyContainer([MyType(serdata.data[:,i], serdata.ids[i]) for i = 1:length(serdata.ids)])
function JLD.writeas(data::MyContainer{T}) where T
    ids = [obj.id for obj in data.objs]
    n = length(data.objs)
    vectors = Matrix{T}(undef, 5, n)
    for i = 1:n
        vectors[:,i] = data.objs[i].data
    end
    MyContainerSerializer(vectors, ids)
end

end # MySerializer

using ..MyTypes, JLD, Compat.Test

obj1 = MyType(rand(5), 2)
obj2 = MyType(rand(5), 17)
container = MyContainer([obj1,obj2])
filename = joinpath(tempdir(), "customserializer-$(randstring()).jld")
jldopen(filename, "w") do file
    write(file, "mydata", container)
end

container_r = jldopen(filename) do file
    obj = file["mydata"]
    dtype = JLD.datatype(obj.plain)
    @test JLD.jldatatype(JLD.file(obj), dtype) === MySerializer.MyContainerSerializer{Float64}
    read(file, "mydata")
end

@test container_r == container


text:
abstract type Ordering end

Base.iterate(ordering::Ordering, state = 0) = state > 0 ? nothing : (ordering, state + 1)
Base.length(ordering::Ordering) = 1

Base.show(io::IO, ordering::Ordering) = print(io, string(ordering))
Base.show(io::IO, ::MIME"application/prs.juno.inline", ordering::Ordering) = print(io, string(ordering))

struct EstimatorOrdering{TE<:Estimator} <: Ordering
    estimator::TE
    orientation::Symbol
end

function EstimatorOrdering(estimator::Estimator; orientation::Symbol = :superiority)
    EstimatorOrdering{typeof(estimator)}(estimator, orientation)
end

string(ordering::EstimatorOrdering) = "EstimatorOrdering"

function more_extreme(x1::TI, x2::TI, x1_::TI, x2_::TI,
    ordering::EstimatorOrdering, design::TD; orientation = nothing)  where {TI<:Integer,TD<:AbstractDesign}

    orientation = (orientation == nothing) ? ordering.orientation : orientation
    if orientation == :superiority
        return ordering.estimator(x1, x2, design) >= ordering.estimator(x1_, x2_, design)
    end
    if orientation == :inferiority
        return ordering.estimator(x1, x2, design) <= ordering.estimator(x1_, x2_, design)
    end
    error("orientation must be :superiororit or :inferiority")
end

function more_extreme(x1::TI, x2::TI, x1_::TI, x2_::TI, estimator::TE,
    design::TD; orientation = nothing)  where {TI<:Integer,TE<:Estimator,TD<:AbstractDesign}

    return more_extreme(x1, x2, x1_, x2_, EstimatorOrdering(estimator), design,
        orientation = oreintation)
end

function p_value(x1::TI, x2::TI, p0::TR, ordering::TO, design::TD; orientation = nothing) where {TI<:Integer,TR<:Real,TO<:Ordering,TD<:AbstractDesign}

    XX   = sample_space(design)
    inds = more_extreme.(XX[:,1], XX[:,2], x1, x2, ordering, design, orientation = orientation)
    return min(1, max(0, sum(
        pmf.(XX[inds,2], n2.(design, XX[inds,1]), p0) .*
        pmf.(XX[inds,1], n1(design), p0)
    )  ) )
end

function p_value(x1::TI, x2::TI, p0::TR, estimator::TE, design::TD; orientation::Symbol = :superiority) where {TI<:Integer,TR<:Real,TE<:Estimator,TD<:AbstractDesign}

    return p_value(x1, x2, p0, EstimatorOrdering(estimator, orientation), design)
end

function compatible(ordering::Ordering, design::AbstractDesign, pnull::Real, α::Real)

    @assert (0 <= pnull <= 1) "pnull must be in [0,1]"

    XX                      = sample_space(design)
    design_rejects          = reject.(XX[:,1], XX[:,2], design)
    pvals                   = p_value.(XX[:,1], XX[:,2], pnull, ordering, design)
    ordering_rejects        = pvals .< α
    incompatibility_degree  = sum(design_rejects .& .!ordering_rejects)

    df = DataFrames.DataFrame(
        x1 = XX[:,1],
        x2 = XX[:,2],
        design_rejects = design_rejects,
        ordering_rejects = ordering_rejects,
        p_value = pvals
    )

    return Dict{String,Any}(
        "compatible" => incompatibility_degree == 0,
        "incompatibility degree" => incompatibility_degree,
        "details" => df
    )
end

function mlecompatible(design::AbstractDesign, pnull::Real, α::Real)

    return compatible(EstimatorOrdering(MaximumLikelihoodEstimator()), design, pnull, α)
end


text:
# First generate some random data and test functions in Base on it
const NVALS = 1000

const input = Dict(
    t=>[
        [ (randindomain(t, NVALS, domain),) for (_, _, domain) in base_unary_real ];
        [ (randindomain(t, NVALS, domain1), randindomain(t, NVALS, domain2))
            for (_, _, domain1, domain2) in base_binary_real ]
    ]
    for t in (Float32, Float64)
)

const fns = [[x[1:2] for x in base_unary_real]; [x[1:2] for x in base_binary_real]]

# output = Dict(t=>[fns[i](input[t][i]...) for i = 1:length(fns)] for t in (Float32, Float64))

@testset "Definitions and Comparison with Base for Reals" begin

  for t in (Float32, Float64), i = 1:length(fns)
    inp = input[t][i]
    mod, fn = fns[i]
    base_fn = getproperty(mod, fn)
    vml_fn = getproperty(IntelVectorMath, fn)
    vml_fn! = getproperty(IntelVectorMath, Symbol(fn, "!"))

    Test.@test parentmodule(vml_fn) == IntelVectorMath

    # Test.test_approx_eq(output[t][i], fn(input[t][i]...), "Base $t $fn", "IntelVectorMath $t $fn")
    baseres = base_fn.(inp...)
    Test.@test vml_fn(inp...) ≈ base_fn.(inp...)

    # cis changes type (float to complex, does not have mutating function)
    if length(inp) == 1
      if fn != :cis
        vml_fn!(inp[1])
        Test.@test inp[1] ≈ baseres
      end
    elseif length(inp) == 2
      out = similar(inp[1])
      vml_fn!(out, inp...)
      Test.@test out ≈ baseres
    end

  end

end

@testset "Error Handling and Settings" begin

  # Verify that we still throw DomainErrors
  Test.@test_throws DomainError IntelVectorMath.sqrt([-1.0])
  Test.@test_throws DomainError IntelVectorMath.log([-1.0])

  # Setting accuracy
  vml_set_accuracy(VML_LA)
  Test.@test vml_get_accuracy() == VML_LA

  vml_set_accuracy(VML_EP)
  Test.@test vml_get_accuracy() == VML_EP

end


text:
immutable MvBrownianMotion{Cov<:AbstractPDMat} <: ContinuousMultivariateStochasticProcess
  timepoints::Vector{Float64}
  npoints::Int64
  Σ::Cov

  function MvBrownianMotion{Cov<:AbstractPDMat}(t::Vector{Float64}, n::Int64, Σ::Cov)
    t[1] > 0.0 || error("First provided time point must be positive for Brownian motion.")
    issorted(t, lt=<=) || error("The time points must be strictly sorted.")
    int64(length(t)) == n || error("Number of time points must be equal to the vector holding the time points.")
    MvBrownianMotion{Cov<:AbstractPDMat}(t, n, Σ)
  end
end


text:
mutable struct Blob <: GitHubType
    content::Union{String, Nothing}
    encoding::Union{String, Nothing}
    url::Union{HTTP.URI, Nothing}
    sha::Union{String, Nothing}
    size::Union{Int, Nothing}
end

Blob(data::Dict) = json2github(Blob, data)
Blob(sha::AbstractString) = Blob(Dict("sha" => sha))

namefield(blob::Blob) = blob.sha

@api_default function blob(api::GitHubAPI, repo, blob_obj; options...)
    result = gh_get_json(api, "/repos/$(name(repo))/git/blobs/$(name(blob_obj))"; options...)
    return Blob(result)
end

@api_default function create_blob(api::GitHubAPI, repo; options...)
    result = gh_post_json(api, "/repos/$(name(repo))/git/blobs"; options...)
    return Blob(result)
end


text:
using Documenter

"""
Wrapper for src/simulation code. This is duplicitous since this should really
be inside of src/simulation/load.jl but I had issues with threading and module
loading in the past so this approach should work for now even if it is clunky.
"""
module Simulation
    packages = [:StatsBase,
                :Distributions,
                :DataFrames,
                :HypothesisTests,
                :IterTools,
                :DocStringExtensions,
                :Combinatorics,
                :DataStructures]

    for package in packages
        eval(:(using $package))
    end

    filenames = ["common.jl", "utils.jl", "library.jl", "transfection.jl",
                 "selection.jl", "sequencing.jl", "processing.jl",
                 "designs.jl"]
    for filename in filenames
        include(joinpath(Base.source_dir(), "..", "src", "simulation", filename))
    end

    export FacsScreen, GrowthScreen, Library, Sampleable, Delta, CRISPRi, CRISPRn,
    construct_library, transfect, select, counts_to_freqs, sequencing,
    differences_between_bins, auprc
end
using Simulation

makedocs(
    modules = [Simulation],
    clean = false,
    format = :html,
    sitename = "Crispulator.jl",
    pages = Any[
        "Home" => "index.md",
        "Tutorial" => "tutorial.md",
        "Custom Simulations" => "custom.md",
        "Simulation Internals" => "internals.md"
    ]
)

deploydocs(
    repo   = "github.com/tlnagy/Crispulator.jl.git",
    julia  = "release",
    osname = "linux",
    deps = nothing,
    make = nothing,
    target = "build",
)


text:
function GzipFilter()
    return GzipFilter(())
end

function filter(obj::GzipFilter, arg0::NamedData)
    return jcall(obj, "filter", NamedData, (NamedData,), arg0)
end


text:
# Adapted from https://github.com/bensadeghi/DecisionTree.jl
__precompile__()

module Estimators
import Base: length, convert, promote_rule, show, start, next, done
export Estimator, Leaf, Node,
       depth,
       fit_regression_tree,
       predict, assign_leaves

float(x) = map(Float64, x)
neg(arr) = map(!, arr) # `.!arr` is invalid in 0.5, and `!arr` triggers a warning in 0.6.
const NO_BEST=(0,0)

# converts a boolean vector into an integer by summing the appropriate powers of 2
boolvec2int(v::Vector{Bool}) = sum(b*v for (b,v) in zip(2.^(0:(length(v)-1)), v))

abstract type Estimator
end

immutable Leaf <: Estimator
    id::Integer
    value::Real
end
    length(leaf::Leaf) = 1
    depth(leaf::Leaf) = 0

Leaf(id::Vector{Bool}, value::Real) = Leaf(boolvec2int(id), value)
Leaf(value::Real) = Leaf(0,value)

immutable Node <: Estimator
    node_id::Integer
    featid::Integer
    featval::Any
    left::Union{Leaf,Node}
    right::Union{Leaf,Node}
end
    length(tree::Node) = length(tree.left) + length(tree.right)
    depth(tree::Node) = 1 + max(depth(tree.left), depth(tree.right))

function Node(id::Vector{Bool}, featid::Integer, featval::Any, left::Union{Leaf,Node}, right::Union{Leaf,Node})
    return Node(boolvec2int(id), featid, featval, left, right)
end

convert(::Type{Node}, x::Leaf) = Node(0, nothing, x, Leaf(nothing,[nothing])) #makes a node with that leaf as it's left node
promote_rule(::Type{Node}, ::Type{Leaf}) = Node # this tells julia to call convert(::node, leaf) on the leaf when this happens
promote_rule(::Type{Leaf}, ::Type{Node}) = Node

function fit_regression_tree{T<:Float64, U<:Real}(X::Matrix{U}, Y::Vector{T}; min_samples_leaf=5, max_depth=-1, node_id::Vector{Bool}=[true])
    if max_depth < -1
        error("Unexpected value for max_depth: $(max_depth) (expected: max_depth >= 0, or max_depth = -1 for infinite depth)")
    end
    if length(Y) <= min_samples_leaf || max_depth==0
        return Leaf(node_id, mean(Y))

    end
    S = _split_mse(X, Y)
    if S == NO_BEST
        return Leaf(node_id, mean(Y))
    end
    feat_id, thresh = S
    split = X[:,feat_id] .< thresh
    return Node(node_id, feat_id, thresh,
                fit_regression_tree(X[split,:], Y[split],
                    min_samples_leaf=min_samples_leaf, max_depth=max(max_depth-1, -1), node_id=vcat(true,node_id)),
                fit_regression_tree(X[neg(split),:], Y[neg(split)],
                    min_samples_leaf=min_samples_leaf, max_depth=max(max_depth-1, -1), node_id=vcat(false,node_id))
                )
end

function _split_mse{T<:Float64, U<:Real}(X::Matrix{U}, Y::Vector{T})
    N, p = size(X)

    best = NO_BEST
    best_val = -Inf

    for i in 1:p
        ord = sortperm(X[:,i])
        X_i = X[ord,i]
        Y_i = Y[ord]
        if N > 100
            if VERSION >= v"0.4.0-dev"
                domain_i = quantile(X_i, linspace(0.01, 0.99, 99);
                                    sorted=true)
            else  # sorted=true isn't supported on StatsBase's Julia 0.3 version
                domain_i = quantile(X_i, linspace(0.01, 0.99, 99))
            end
        else
            domain_i = X_i
        end
        value, thresh = _best_mse_loss(X_i, Y_i, domain_i)
        if value > best_val
            best_val = value
            best = (i, thresh)
        end
    end

    return best
end

function _best_mse_loss{T<:Float64, U<:Real}(X::Vector{U}, Y::Vector{T}, domain)
    # True, but costly assert. However, see
    # https://github.com/JuliaStats/StatsBase.jl/issues/164
    # @assert issorted(X) && issorted(domain)
    best_val = -Inf
    best_thresh = 0.0
    s_l = s2_l = zero(T)
    su = sum(Y)::T
    su2 = zero(T); for l in Y su2 += l*l end  # sum of squares
    nl = 0
    n = length(Y)
    i = 1
    # Because the `X` are sorted, below is an O(N) algorithm for finding
    # the optimal threshold amongst `domain`. We simply iterate through the
    # array and update s_l and s_r (= sum(Y) - s_l) as we go. - @cstjean
    @inbounds for thresh in domain
        while i <= length(Y) && X[i] < thresh
            l = Y[i]

            s_l += l
            s2_l += l*l
            nl += 1

            i += 1
        end
        s_r = su - s_l
        s2_r = su2 - s2_l
        nr = n - nl
        # This check is necessary I think because in theory all Y could
        # be the same, then either nl or nr would be 0. - @cstjean
        if nr > 0 && nl > 0
            loss = s2_l - s_l^2/nl + s2_r - s_r^2/nr
            if -loss > best_val
                best_val = -loss
                best_thresh = thresh
            end
        end
    end
    return best_val, best_thresh
end

predict(leaf::Leaf, feature::Vector) = leaf.value

function predict(tree::Node, X::Vector)
    if tree.featval == nothing
        return predict(tree.left, X)
    elseif X[tree.featid] < tree.featval
        return predict(tree.left, X)
    else
        return predict(tree.right, X)
    end
end

function predict(tree::Union{Leaf,Node}, X::Matrix)
    N = size(X,1)
    predictions = Array{Any}(N)
    for i in 1:N
        predictions[i] = predict(tree, X[i,:])
    end
    println("predictions ->")
    println(predictions)
    if typeof(predictions[1]) <: Float64
        return float(predictions)
    else
        return predictions
    end
end

function predict(estimator_dict::Dict, X::Matrix)
    return Dict(k=>predict(v,X) for (k,v) in estimator_dict)
end

assign_leaves(leaf::Leaf, feature::Vector) = leaf.id

function assign_leaves(tree::Node, X::Vector)
    if tree.featval == nothing
        return assign_leaves(tree.left, X)
    elseif X[tree.featid] < tree.featval
        return assign_leaves(tree.left, X)
    else
        return assign_leaves(tree.right, X)
    end
end

function assign_leaves(tree::Union{Leaf,Node}, X::Matrix)
    N = size(X,1)
    assignments = Array{Int}(N)
    for i in 1:N
        assignments[i] = assign_leaves(tree, X[i,:])
    end
    return assignments
end

function assign_leaves(estimator_dict::Dict, X::Matrix)
    return Dict(k=>assign_leaves(v,X) for (k,v) in estimator_dict)
end

end #module

max_stars_repo_name: naskoD/bachelorThesis

max_issues_repo_name: naskoD/bachelorThesis

max_forks_repo_name: naskoD/bachelorThesis


text:
################################################################################
#
#  Constructors
#
################################################################################

@doc Markdown.doc"""
    hermitian_space(K::NumField, n::Int) -> HermSpace

Create the Hermitian space over `K` with dimension `n` and Gram matrix equal to
the identity matrix. The number field `K` must be a quadratic extension, that
is, `degree(K) == 2` must hold.
"""
function hermitian_space(K::NumField, n::Int)
  G = identity_matrix(K, n)
  return HermSpace(K, G)
end

@doc Markdown.doc"""
    hermitian_space(K::NumField, G::MatElem) -> HermSpace

Create the Hermitian space over `K` with Gram matrix equal to the identity
matrix. The matrix `G` must be square and Hermitian with respect to the non-trivial
automorphism of `K`. The number field `K` must be a quadratic extension, that
is, `degree(K) == 2` must hold.
"""
function hermitian_space(K::NumField, G::MatElem)
  return HermSpace(K, G)
end

################################################################################
#
#  String I/O
#
################################################################################

function Base.show(io::IO, V::HermSpace)
  print(io, "Hermitian space over\n")
  println(io, base_ring(V))
  println(io, "with Gram matrix")
  print(io, gram_matrix(V))
end

################################################################################
#
#  Absolute field
#
################################################################################

# this is used internally to accelerate computations by passing to an absolute
# field
function absolute_field(V::HermSpace)
  c = get_special(V, :absolute_field)
  if c === nothing
    Eabs, EabsToE, KtoE = absolute_field(base_ring(V))
    set_special(V, :absolute_field => (Eabs, EabsToE, KtoE))
    return Eabs, EabsToE, KtoE
  else
    return c::Tuple{AnticNumberField,Hecke.NfToNfRel,NfToNfMor}
  end
end

################################################################################
#
#  Predicates
#
################################################################################

isquadratic(V::HermSpace) = false

ishermitian(V::HermSpace) = false

_base_algebra(V::HermSpace) = V.E

################################################################################
#
#  Properties
#
################################################################################

involution(V::HermSpace) = V.involution

fixed_field(V::HermSpace) = V.K

################################################################################
#
#  Inner product
#
################################################################################

function _inner_product(G, v, w, involution)
  return _inner_product(G, v, [involution(x) for x in w])
end

function inner_product(V::HermSpace, v::Vector, w::Vector)
  _inner_product(gram_matrix(V), v, w, involution(V))
end

################################################################################
#
#  Diagonalization
#
################################################################################

function diagonal(V::HermSpace)
  D, _ = _gram_schmidt(gram_matrix(V), involution(V))
  return map(fixed_field(V), diagonal(D))
end

################################################################################
#
#  Hasse and Witt invariant
#
################################################################################

function hasse_invariant(L::HermSpace, p)
  throw(error("The space must be quadratic"))
end

function witt_invariant(L::HermSpace, p)
  throw(error("The space must be quadratic"))
end

################################################################################
#
#  Local equivalence
#
################################################################################

function isequivalent(L::HermSpace{AnticNumberField}, M::HermSpace{AnticNumberField}, p::fmpz)
  return _isequivalent(L, M, p)
end

function isequivalent(L::HermSpace, M::HermSpace, p::NfOrdIdl)
  return _isequivalent(L, M, p)
end

function _isequivalent(L::HermSpace, M::HermSpace, p)
  base_ring(L) != base_ring(M) && error("Both spaces must have the same base field")
  A = gram_matrix(L)
  B = gram_matrix(M)
  if A == B
    return true
  end

  if rank(L) != rank(M)
    return false
  end

  return islocal_norm(base_ring(L), det(L) * det(M), p)[1]
end

function isequivalent(L::HermSpace, M::HermSpace, P::InfPlc)
  if L == M
    return true
  end

  if iscomplex(P)
    return true
  end

  DL = diagonal(L)
  DM = diagonal(M)
  iL = count(d -> isnegative(d, P), DL)
  iM = count(d -> isnegative(d, P), DM)
  return iL == iM
end

################################################################################
#
#  Global equivalence
#
################################################################################

function isequivalent(M::HermSpace, L::HermSpace)
  if gram_matrix(M) == gram_matrix(L)
    return true
  end

  if rank(M) != rank(L)
    return false
  end

  E = base_ring(M)
  # I could replace this with a islocal_norm at the ramified primes + primes
  # dividing right hand side
  return isnorm(E, det(M) * det(L))[1]
end

################################################################################
#
#  Isotropic
#
################################################################################

isisotropic(V::HermSpace, p::InfPlc) = _isisotropic(V, p)

function isisotropic(V::HermSpace, q)
  if nf(order(q)) == base_ring(V)
    p = minimum(q)
  else
    p = q
  end
  @assert fixed_field(V) == nf(order(p))
  r = rank(V)
  if r >= 3
    return true
  elseif r == 0
    return false
  end
  d = det(V)
  if r == 1
    return d == 0
  end
  return islocal_norm(base_ring(V), -d, p)
end

################################################################################
#
#  Hyperbolic spaces
#
################################################################################

function _islocally_hyperbolic_hermitian_detclass(rk, d, E, K, p)
  if isodd(rk)
    return false
  end
  if d == 1
    if iseven(div(rk, 2))
      return true
    else
      return islocal_norm(E, K(-1), p)
    end
  else
    if iseven(div(rk, 2))
      return false
    else
      return !islocal_norm(E, K(-1), p)
    end
  end
end

function islocally_hyperbolic(V::HermSpace, p)
  rk = rank(V)
  if isodd(rk)
    return false
  end
  return islocal_norm(base_ring(V), det(V) * (-1)^(div(rk, 2)), p)
end

################################################################################
#
#  Embeddings
#
################################################################################

# Hermitian forms are uniquely determined by their rank and determinant class
# Thus there is no restriction to embeddings.

@doc Markdown.doc"""
    islocally_represented_by(U::HermSpace, V::HermSpace, p)

Return whether $U$ is represented by $V$ locally at $\mathfrak p$.
"""
function islocally_represented_by(U::HermSpace, V::HermSpace, p)
  if rank(U) < rank(V)
    return false
  elseif rank(U) == rank(V)
    return isequivalent(U, V, p)
  else
    return true
  end
end

# There are no restrictions, since spaces are uniquely determined by their
# rank and determinant.

@doc Markdown.doc"""
    isrepresented_by(U::HermSpace, V::HermSpace)

Return whether $U$ is represented by $V$, that is, whether $U$ embeds into $V$.
"""
function isrepresented_by(U::HermSpace, V::HermSpace)
  v = rank(V) - rank(U)
  if v < 0
    return false
  elseif v == 0
    return isequivalent(U, V)
  else
    return true
  end
end


text:
#
#= = # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Description
#
#   This file contains functions to retrieve space indices for many models.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # = =#

export get_space_index, init_space_indices

# Interpolation types used in space indices.
_space_indices_itp_constant{T} =
    Interpolations.GriddedInterpolation{T,1,T,Gridded{Interpolations.Constant{Interpolations.Nearest}},Tuple{Array{Float64,1}}}

_space_indices_itp_linear{T} =
    Interpolations.GriddedInterpolation{T,1,T,Gridded{Linear},Tuple{Array{Float64,1}}}

include("./dtcfile.jl")
include("./fluxtable.jl")
include("./solfsmy.jl")
include("./wdcfiles.jl")
include("./dstfiles.jl")
include("./space_indices_helpers.jl")

################################################################################
#                                Private Macros
################################################################################

macro _check_data(itp, DT)
    quote
        if ($(esc(DT)) < $(esc(itp)).knots[1][1]) ||
           ($(esc(DT)) > $(esc(itp)).knots[1][end])
            error("The data for the requested Julian Day is not available!")
        end
    end
end

################################################################################
#                               Public Functions
################################################################################

"""
    get_space_index(T, DT::Number; ...)

Return the space index `T` at the day `DT` [Julian Day]. `T` can be:

## Daily 10.7-cm solar flux

The daily 10.7-cm solar flux can be obtained using:

* `F10()`: 10.7-cm adjusted solar flux \\[10⁻²² W/(M² Hz)].
* `F10adj()`: 10.7-cm adjusted solar flux \\[10⁻²² W/(M² Hz)].
* `F10obs()`: 10.7-cm observed solar flux \\[10⁻²² W/(M² Hz)].

These indices require `fluxtable` (see `init_space_indices`).

## Daily average 10.7-cm solar flux

The daily average 10.7-cm solar flux, centered at `DT`, can be obtained using:

* `F10M()`: 10.7-cm adjusted solar flux \\[10⁻²² W/(M² Hz)].
* `F10Madj()`: 10.7-cm adjusted solar flux \\[10⁻²² W/(M² Hz)].
* `F10Mobs()`: 10.7-cm observed solar flux \\[10⁻²² W/(M² Hz)].

In this case, the keyword `window::Int` can be passed to select the size of the
window. By default, it is selected as 81.

These indices require `fluxtable` (see `init_space_indices`).

# Daily Kp and Ap

* `Kp()`: Kp index (daily mean).
* `Kp_vect()`: A vector containing the Kp index for the following hours of the
               day: 0-3h, 3-6h, 6-9h, 9-12h, 12-15h, 15-18h, 18-20h, 20-23h.
* `Ap()`: Ap index (daily mean).
* `Ap_vect()`: A vector containing the Ap index for the following hours of the
               day: 0-3h, 3-6h, 6-9h, 9-12h, 12-15h, 15-18h, 18-20h, 20-23h.

These indices require `wdcfiles` (see `init_space_indices`).

# Daily S10, M10, and Y10

* `S10()`: EUV index (26-34 nm) scaled to F10.7.
* `M10()`: MG2 index scaled to F10.7.
* `Y10()`: Solar X-ray & Lya index scaled to F10.7.

These indices require `solfsmy` (see `init_space_indices`).

# 81-day centered average of S10, M10, and Y10.

* `S81a`: EUV 81-day averaged centered index.
* `M81a`: MG2 81-day averaged centered index.
* `Y81a`: Solar X-ray & Lya 81-day averaged centered index.

These indices require `solfsmy` (see `init_space_indices`).

# Exospheric temperature variation due to Dst

* `DstΔTc`: Exospheric temperature variation due to `Dst` [K].

This index requires `dtcfile` (see `init_space_indices`).

"""
@inline get_space_index(::SPI{:F10}, DT::Number) =
    get_space_index(Val(:F10adj), DT)

@inline function get_space_index(::SPI{:F10obs}, DT::Number)
    @_check_data(get(_fluxtable_data).F10obs, DT)
get(_fluxtable_data).F10obs(DT)
end

@inline function get_space_index(::SPI{:F10adj}, DT::Number)
    @_check_data(get(_fluxtable_data).F10adj, DT)
get(_fluxtable_data).F10adj(DT)
end

@inline get_space_index(::SPI{:F10M}, DT::Number; window::Int=81) =
    get_space_index(Val(:F10Madj), DT; window=window)

@inline function get_space_index(::SPI{:F10Mobs}, DT::Number; window::Int=81)
    Δ = floor(Int, (window - 1) / 2)
    @_check_data(get(_fluxtable_data).F10obs, DT - Δ)
    @_check_data(get(_fluxtable_data).F10obs, DT + Δ)
    mean(get(_fluxtable_data).F10obs(DT - Δ:1:DT + Δ))
end

@inline function get_space_index(::SPI{:F10Madj}, DT::Number; window::Int=81)
    Δ = floor(Int, (window - 1) / 2)
    @_check_data(get(_fluxtable_data).F10adj, DT - Δ)
    @_check_data(get(_fluxtable_data).F10adj, DT + Δ)
    mean(get(_fluxtable_data).F10adj(DT - Δ:1:DT + Δ))
end

@inline function get_space_index(::SPI{:Kp_vect}, DT::Number)
    @_check_data(get(_wdc_data).Kp, DT)
    get(_wdc_data).Kp(DT)
end

@inline function get_space_index(::SPI{:Ap_vect}, DT::Number)
    @_check_data(get(_wdc_data).Ap, DT)
    get(_wdc_data).Ap(DT)
end

@inline get_space_index(::SPI{:Kp}, DT::Number) =
    mean(get_space_index(Val(:Kp_vect), DT))

@inline get_space_index(::SPI{:Ap}, DT::Number) =
    mean(get_space_index(Val(:Ap_vect), DT))

@inline function get_space_index(::SPI{:Dst_vect}, DT::Number)
    @_check_data(get(_dst_data).Dst, DT)
    get(_dst_data).Dst(DT)
end

@inline get_space_index(::SPI{:Dst}, DT::Number) =
    mean(get_space_index(Val(:Dst_vect), DT))

@inline function get_space_index(::SPI{:S10}, DT::Number)
    @_check_data(get(_solfsmy_data).S10, DT)
    get(_solfsmy_data).S10(DT)
end

@inline function get_space_index(::SPI{:S81a}, DT::Number)
    @_check_data(get(_solfsmy_data).S81a, DT)
    get(_solfsmy_data).S81a(DT)
end

@inline function get_space_index(::SPI{:M10}, DT::Number)
    @_check_data(get(_solfsmy_data).M10, DT)
    get(_solfsmy_data).M10(DT)
end

@inline function get_space_index(::SPI{:M81a}, DT::Number)
    @_check_data(get(_solfsmy_data).M81a, DT)
    get(_solfsmy_data).M81a(DT)
end

@inline function get_space_index(::SPI{:Y10}, DT::Number)
    @_check_data(get(_solfsmy_data).Y10, DT)
    get(_solfsmy_data).Y10(DT)
end

@inline function get_space_index(::SPI{:Y81a}, DT::Number)
    @_check_data(get(_solfsmy_data).Y81a, DT)
    get(_solfsmy_data).Y81a(DT)
end

@inline function get_space_index(::SPI{:DstΔTc}, DT::Number)
    @_check_data(get(_dtcfile_data).DstΔTc, DT)
    get(_dtcfile_data).DstΔTc(DT)
end

"""
    init_space_indices(...)

Initialize all space indices. The files that will be initialized must be
indicated by the array of symbols passed to the keyword argument
`enabled_files`. If this is `nothing`, which is the default, then all files will
be initialized. The symbol related to each file is described next.

Notice that the initialization process can be changed by a set of keywords as
described next.

## DTCFILE

**Symbol**: `:dtcfile`

This file contains the exospheric temperature variation caused by the Dst index.
This is used for the JB2008 atmospheric model.

### Keywords

* `dtcfile_path`: Path for the file `DTCFILE.TXT`. If `nothing`, then it will be
                  downloaded. (**Default** = `nothing`)
* `dtcfile_force_download`: If `true`, then the file will always be downloaded
                            if the path is not specified. (**Default** =
                            `false`).

## fluxtable

**Symbol**: `:fluxtable`

This file contains the F10.7 flux data in different formats.

### Keywords

* `fluxtable_path`: Path for the file `fluxtable.txt`. If `nothing`, then it
                    will be downloaded. (**Default** = `nothing`)
* `fluxtable_force_download`: If `true`, then the file will always be downloaded
                              if the path is not specified.
                              (**Default** = `false`).

## SOLFSMY

**Symbol**: `:solfsmy`

This files contains the indices necessary for the JB2008 atmospheric model.

### Keywords

* `solfsmy_path`: Path for the file `SOLFSMY.TXT`. If `nothing`, then it will be
                  downloaded. (**Default** = `nothing`)
* `solfsmy_force_download`: If `true`, then the file will always be downloaded
                            if the path is not specified. (**Default** =
                            `false`).

## WDC Files

**Symbol**: `:wdcfiles`

This set of files contain the Kp and Ap indices.

### Keywords

* `wdcfiles_path`: Path for the directory with the WDC files. If `nothing`, then
                   they will be downloaded. (**Default** = `nothing`)
* `wdcfiles_force_download`: If `true`, then the files will always be downloaded
                            if the path is not specified. (**Default** =
                            `false`).
* `wdcfiles_oldest_year`: Oldest year in which the WDC file will be obtained.
                          (**Default** = past 3 years).
* `wdcfiles_newest_year`: Newest year in which the WDC file will be obtained. If
                          it is `nothing`, then it defaults to the current year.
                          (**Default** = `nothing`).

"""
function init_space_indices(;
            enabled_files=nothing,
            dtcfile_path=nothing,
            dtcfile_force_download=false,

            fluxtable_path=nothing,
            fluxtable_force_download=false,

            solfsmy_path=nothing,
            solfsmy_force_download=false,

            wdcfiles_dir=nothing,
            wdcfiles_force_download=false,
            wdcfiles_oldest_year=year(now()) - 3,
            wdcfiles_newest_year=nothing,

            dstfiles_dir=nothing,
            dstfiles_force_download=false,
            dstfiles_oldest_year=year(now()) - 3,
            dstfiles_newest_year=nothing
        )
    dtcfile   = (enabled_files == nothing) || (:dtcfile in enabled_files)
    fluxtable = (enabled_files == nothing) || (:fluxtable in enabled_files)
    solfsmy   = (enabled_files == nothing) || (:solfsmy in enabled_files)
    wdcfiles  = (enabled_files == nothing) || (:wdcfiles in enabled_files)
    dstfiles  = (enabled_files == nothing) || (:dstfiles in enabled_files)

    dtcfile && _init_dtcfile(local_path=dtcfile_path,
                            force_download=dtcfile_force_download)

    fluxtable && _init_fluxtable(local_path=fluxtable_path,
                                force_download=fluxtable_force_download)

    solfsmy && _init_solfsmy(local_path=solfsmy_path,
                            force_download=solfsmy_force_download)

    wdcfiles && _init_wdcfiles(local_dir=wdcfiles_dir,
                            force_download=wdcfiles_force_download,
                            wdcfiles_oldest_year=wdcfiles_oldest_year,
                            wdcfiles_newest_year=wdcfiles_newest_year)

    dstfiles && _init_dstfiles(local_dir=dstfiles_dir,
                            force_download=dstfiles_force_download,
                            dstfiles_oldest_year=dstfiles_oldest_year,
                            dstfiles_newest_year=dstfiles_newest_year)

    nothing
end


text:
using Pkg.Artifacts
using JSON
using NodeJS
using URIs

function build()
    package_dict = JSON.parsefile(joinpath(@__DIR__, "package.json"))
    pkgname = package_dict["name"]
    version = VersionNumber(package_dict["version"])
    host = "https://github.com/MichaelHatherly/NomnomlJS.jl/releases/download"

    build_path = joinpath(@__DIR__, "build")

    if ispath(build_path)
        rm(build_path, force=true, recursive=true)
    end

    mkpath(build_path)

    artifact_toml = joinpath(build_path, "..", "Artifacts.toml")

    npm = NodeJS.npm_cmd()

    product_hash = create_artifact() do artifact_dir
        cp(joinpath(@__DIR__, "index.js"), joinpath(artifact_dir, "index.js"))
        cp(joinpath(@__DIR__, "package.json"), joinpath(artifact_dir, "package.json"))
        cp(joinpath(@__DIR__, "package-lock.json"), joinpath(artifact_dir, "package-lock.json"))

        run(Cmd(`$npm install`, dir = artifact_dir))
    end

    archive_filename = "$pkgname-$version.tar.gz"

    download_hash = archive_artifact(product_hash, joinpath(build_path, archive_filename))

    bind_artifact!(
        artifact_toml,
        "nomnoml",
        product_hash,
        force = true,
        download_info = Tuple[
            (
                "$host/$(escapeuri(string(version)))/$archive_filename",
                download_hash,
            ),
        ],
    )
end

build()


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*-----------------------------------------------------------------------
#* Author: Tiago Os򱨯rio
#* $Id$
#*---------------------------------------------------------------------
type Sum
	Sum()=begin
		new(
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"input signal 1"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"input signal 2"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"output signal"
			)),
			[
				:(output = input1+input2),
			],
			[
				"Calculate output",
			],
			[:input1,:input2,:output,]
		)
	end
	input1::DanaReal
	input2::DanaReal
	output::DanaReal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Sum
function setEquationFlow(in::Sum)
	addEquation(1)
end
function atributes(in::Sum,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Icon]="icon/Sum"
	fields[:Brief]="Model Sum."
	fields[:Info]="== Inputs ==
* Two input signals.

== Outputs ==
* One output signal.
"
	drive!(fields,_)
	return fields
end
Sum(_::Dict{Symbol,Any})=begin
	newModel=Sum()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using DiffEqSensitivity
using Zygote
#using ReverseDiff
using Random
using LinearAlgebra
using DifferentialEquations#,DiffEqFlux
#using ForwardDiff
using GalacticOptim,Optim,NLopt
using Plots
#
#using GLOQ
include("../src/GLOQ.jl")
pyplot()
# BlackBoxOptim somehow downgrade some packages and as a result breaks the auto-differentiation with Zygote
# we should avoid it.
# NLopt seems to be really fast. Their LBFGS with box constrains is faster than Optim's FMinbox(LBFGS)

N_states = 4;
freqs = [4.10817777; 3.88303940; 3.61937188]
omegas = 2.0*pi.*freqs
gamma1   = [2.59451982e-05; 4.54296537e-05; 0.0]
gamma2   = [2.58005604e-05; 9.00000000e-05; 0.0]
omr = 2.0*pi*(3.8830355 - 2.0e-3)
HK_free = GLOQ.RotationFrameDiagonal(omegas,omr)
L1,L2 = GLOQ.RotationFrameLindblad(gamma1,gamma2)

width = 17.00
TC = 2.5*width
Ω = 0.5*pi/(TC*sqrt(2.0))
θ = 0.0
HK_control,HS_control = GLOQ.RotationFrameRamseyControl(Ω,θ,N_states)

fromState = 1
rho_u0 = [0.0;0.0;0.0;0.0]
rho_v0 = [0.0;0.0;0.0;0.0]
rho_u0[fromState+1] = 1.0

T_Ramsey = 1.0*GLOQ.GLOQ_MICRO_SEC
N_dark_times = 21#101
t_dark_times = zeros(Float64,N_dark_times)
dt = T_Ramsey/(N_dark_times-1)
for i = 1:N_dark_times
	dark_time = (i-1)*dt
	t_dark_times[i] = dark_time
end
rho_Ramsey_u,rho_Ramsey_v = GLOQ.RamseyForwardSolve(rho_u0,rho_v0,
			     omegas,omr,
				 gamma1,gamma2,
				 1, # initial state
				 TC,t_dark_times,N_states)
population_Ramsey = GLOQ.get_population(rho_Ramsey_u)

fig=plot(t_dark_times./GLOQ.GLOQ_MICRO_SEC,population_Ramsey)
display(fig)
#################################################
# Learn parameters
#################################################
p_true = [freqs;gamma1[1:2];gamma2[1:2]]
p_initial = [freqs.-1e-4;0.9.*gamma1[1:2];0.9.*gamma2[1:2]]
lower_bound = (0.85).*p_true
upper_bound = (1.15).*p_true

rho_Ramsey_u_guess,rho_Ramsey_v_guess = GLOQ.RamseyForwardSolve(rho_u0,rho_v0,
				 (2*pi).*p_initial[1:3],omr,
				 [p_initial[4:5];0.0],[p_initial[6:7];0.0],
				 1, # initial state
				 TC,t_dark_times,N_states;initial_type = "states")
population_guess= GLOQ.get_population(rho_Ramsey_u_guess)

plot!(fig,t_dark_times./GLOQ.GLOQ_MICRO_SEC,population_guess,line=(:dash),legend=:outerright,)
display(fig)
fig2 = plot(t_dark_times./GLOQ.GLOQ_MICRO_SEC,population_guess-population_Ramsey,legend=:outerright,)
display(fig2)
##################################################################################################
global function_call
function loss(p)
	_rho_Ramsey_u,_rho_Ramsey_v = GLOQ.RamseyForwardSolve(rho_u0,rho_v0,
				     (2*pi).*p[1:3],omr,
					 [p[4:5];0.0],[p[6:7];0.0],#gamma1,gamma2,
					 1, # initial state
					 TC,t_dark_times,N_states;
					 method="DiffEqDefault",
					 sensealg=BacksolveAdjoint()) # 144 secs for fminbox
	_population_Ramsey = GLOQ.get_population(_rho_Ramsey_u)
	_loss = sum(abs2,_population_Ramsey-population_Ramsey)/N_dark_times
	return _loss
end
println("Initial loss: ",loss(p_initial))
@time Zygote.gradient(loss,p_initial)

global function_call_num
function_call_num = 0
function loss_gala(p,pp)
	global function_call_num
	function_call_num += 1
	_rho_Ramsey_u,_rho_Ramsey_v = GLOQ.RamseyForwardSolve(rho_u0,rho_v0,
				     (2*pi).*p[1:3],omr,
					 [p[4:5];0.0],[p[6:7];0.0],#gamma1,gamma2,
					 1, # initial state
					 TC,t_dark_times,N_states;
					 method="DiffEqDefault",
					 sensealg=BacksolveAdjoint())
	_population_Ramsey = GLOQ.get_population(_rho_Ramsey_u)
	_loss = sum(abs2,_population_Ramsey-population_Ramsey)/N_dark_times
	#println("Function call ",function_call_num," Loss: ",_loss)
	return _loss
end

#=
struct OptimizationSolution{T, N, uType, P, A, Tf, O} <: AbstractOptimizationSolution{T, N}
    u::uType # minimizer
    prob::P # optimization problem
    alg::A # algorithm
    minimum::Tf
    retcode::Symbol
    original::O # original output of the optimizer
end
=#

optimization_function_fd = OptimizationFunction(loss_gala, GalacticOptim.AutoFiniteDiff())
prob_fd = GalacticOptim.OptimizationProblem(optimization_function_fd, p_initial,
										 lb = lower_bound, ub = upper_bound)

optimization_function_autozygote = OptimizationFunction(loss_gala, GalacticOptim.AutoZygote())
prob_autozygote = GalacticOptim.OptimizationProblem(optimization_function_autozygote, p_initial,
										 lb = lower_bound, ub = upper_bound)

function_call_num = 0

println("Optim Fminbox(LBFGS) Optimization starts")
@time sol_optim_minbox_lbfgs = GalacticOptim.solve(prob_autozygote ,Fminbox(LBFGS()),
						        outer_iterations = 20,
								iterations = 10,
								show_trace=true,
								f_tol = 1e-3,
								outer_f_tol = 1e-3)
println("Optim Fminbox(LBFGS) Optimization done")

function_call_num = 0
println("NLopt Nelder-Mead Optimization starts")
@time sol_nlopt_NedlerMead = solve(prob_autozygote , Opt(:LN_NELDERMEAD,length(p_initial)),
								   maxiters=200,
								   ftol_rel=1e-3)
println("NLopt Nelder-Mead Optimization done")

function_call_num = 0
println("NLopt LBFGS Optimization with Zygote gradient starts")
@time sol_nlopt_LBFGS = solve(prob_autozygote , Opt(:LD_LBFGS,length(p_initial)),
								   maxiters=200,
								   ftol_rel=1e-3)
println("NLopt LBFGS Optimization with Zygote gradient done")

function_call_num = 0
println("NLopt LBFGS Optimization with FD gradient starts")
@time sol_nlopt_LBFGS_fd = solve(prob_fd , Opt(:LD_LBFGS,length(p_initial)),
								   maxiters=200,
								   ftol_rel=1e-3)
println("NLopt LBFGS Optimization with FD gradient done")

#=
@time p_DEFlux = DiffEqFlux.sciml_train(loss, p_initial, Fminbox(LBFGS()),
										lower_bounds=lower_bound,
										upper_bounds=upper_bound,
										show_trace=true,
										iterations=10,
										outer_iterations=10,
										f_tol=1e-4,
										outer_f_tol=5e-4,
										outer_x_tol=1e-6)
=#
println("True p: ",p_true)
println("Initial p: ", p_initial," Initial error: ",p_initial-p_true)
println("\nOptim Fminbox-LBFGS: ",sol_optim_minbox_lbfgs.u," \nLoss:",sol_optim_minbox_lbfgs.minimum," \nError: ",sol_optim_minbox_lbfgs.u-p_true)
println("\nNLopt Nelder-Mead: ",sol_nlopt_NedlerMead.u," \nLoss:",sol_nlopt_NedlerMead.minimum," \nError: ",sol_nlopt_NedlerMead.u-p_true)
println("\nNLopt LBFGS with Zygote gradient: ",sol_nlopt_LBFGS.u," \nLoss:",sol_nlopt_LBFGS.minimum," \nError: ",sol_nlopt_LBFGS.u-p_true)
println("\nNLopt LBFGS with FD gradient: ",sol_nlopt_LBFGS_fd.u," \nLoss:",sol_nlopt_LBFGS_fd.minimum," \nError: ",sol_nlopt_LBFGS_fd.u-p_true)
#println("\nNLopt LBFGS with ReverseDiff gradient: ",sol_nlopt_LBFGS_rd.u," \nLoss:",sol_nlopt_LBFGS_rd.minimum," \nError: ",sol_nlopt_LBFGS_rd.u-p_true)

#######################################################################

rho_Ramsey_u_optim,rho_Ramsey_v_optim = GLOQ.RamseyForwardSolve(rho_u0,rho_v0,
				 (2*pi).*sol_nlopt_LBFGS.u[1:3],omr,
				 [sol_nlopt_LBFGS.u[4:5];0.0],[sol_nlopt_LBFGS.u[6:7];0.0],
				 1,# initial state
				 TC,t_dark_times,N_states;initial_type = "states")
population_optim = GLOQ.get_population(rho_Ramsey_u_optim)

plot!(fig2,t_dark_times./GLOQ.GLOQ_MICRO_SEC,population_optim-population_Ramsey,
	  line=(:dash),lab=:false,legend=:outerright,)
display(fig2)

fig3=plot(t_dark_times./GLOQ.GLOQ_MICRO_SEC,population_optim-population_Ramsey,
	  line=(:dash),lab=:false,legend=:outerright,)
display(fig3)

fig4=plot(t_dark_times./GLOQ.GLOQ_MICRO_SEC,population_Ramsey)
plot!(fig4,t_dark_times./GLOQ.GLOQ_MICRO_SEC,population_optim,line=(:dash))


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

using Markdown
using InteractiveUtils

# ╔═╡ d5257f82-11de-11eb-0176-b1858e7894ae
begin
	using DataFrames
	using Pipe
	using Gadfly
	using CSV
	using BenchmarkTools
end

# ╔═╡ a612e268-11e9-11eb-2305-0175760c625d
HTML("""
		<style>
			main {
				max-width: 90%;
			}
		</style>
	""")

# ╔═╡ e401509a-11de-11eb-2139-f90fc50ae6cf
begin
	datadirectory = "../data/titanic/"
	train_df = CSV.read(datadirectory * "train.csv")
	describe(train_df)
end

# ╔═╡ a6b3a7a4-11ea-11eb-0d7c-0797092fec7e
ENV["COLUMNS"] = 10

# ╔═╡ Cell order:
# ╠═a612e268-11e9-11eb-2305-0175760c625d
# ╠═d5257f82-11de-11eb-0176-b1858e7894ae
# ╠═e401509a-11de-11eb-2139-f90fc50ae6cf
# ╠═a6b3a7a4-11ea-11eb-0d7c-0797092fec7e

max_stars_repo_name: StatisticalMice/experiments

max_issues_repo_name: StatisticalMice/experiments

max_forks_repo_name: StatisticalMice/experiments