text:
"""
    LShapedSolver(lpsolver::AbstractMathProgSolver; <keyword arguments>)

Return an L-shaped algorithm object that can optimize a two-stage `StochasticPrograms`. Supply `lpsolver`, a MathProgBase solver capable of solving linear-quadratic problems.

The following L-shaped regularizations are available
- [`NoRegularization`](@ref):  L-shaped algorithm (default)
- [`RegularizedDecomposition`](@ref):  Regularized decomposition ?RegularizedDecomposition for parameter descriptions.
- [`TrustRegion`](@ref):  Trust-region ?TrustRegion for parameter descriptions.
- [`LevelSet`](@ref):  Level-set ?LevelSet for parameter descriptions.

The following aggregation schemes are available
- [`NoAggregation`](@ref):  Multi-cut L-shaped algorithm (default)
- [`PartialAggregation`](@ref):  ?PartialAggregation for parameter descriptions.
- [`FullAggregation`](@ref):  ?FullAggregation for parameter descriptions.
- [`DynamicAggregation`](@ref):  ?DynamicAggregation for parameter descriptions.
- [`ClusterAggregation`](@ref):  ?ClusterAggregation for parameter descriptions.
- [`HybridAggregation`](@ref):  ?HybridAggregation for parameter descriptions.

The following consolidation schemes are available
- [`NoConsolidation`](@ref)
- [`Consolidation`](@ref)

The following execution policies are available
- [`Serial`](@ref):  Classical L-shaped (default)
- [`Synchronous`](@ref): Classical L-shaped run in parallel
- [`Asynchronous`](@ref): Asynchronous L-shaped ?Asynchronous for parameter descriptions.

...
# Arguments
- `lpsolver::AbstractMathProgSolver`: MathProgBase solver capable of solving linear (and possibly quadratic) programs.
- `feasibility_cuts::Bool = false`: Specify if feasibility cuts should be used
- `subsolver::AbstractMathProgSolver = lpsolver`: Optionally specify a different solver for the subproblems.
- `regularize::AbstractRegularizer = DontRegularize()`: Specify regularization procedure (DontRegularize, RegularizedDecomposition/RD/WithRegularizedDecomposition, TrustRegion/TR/WithTrustRegion, LevelSet/LV/WithLevelSets).
- `aggregate::AbstractAggregator = DontAggregate()`: Specify aggregation procedure (DontAggregate, Aggregate, PartialAggregate, DynamicAggregate)
- `consolidate::AbstractConsolidator = DontConsolidate()`: Specify consolidation procedure (DontConsolidate, Consolidate)
- `execution::Execution = Serial`: Specify how algorithm should be executed (Serial, Synchronous, Asynchronous). Distributed variants requires worker cores.
- `crash::CrashMethod = Crash.None`: Crash method used to generate an initial decision. See ?Crash for alternatives.
- <keyword arguments>: Algorithm specific parameters, See `?LShaped` for list of possible arguments and default values.
...

## Examples

The following solves a stochastic program `sp` created in `StochasticPrograms.jl` using the L-shaped algorithm with GLPK as an `lpsolver`.

```jldoctest
julia> optimize!(sp, solver = LShapedSolver(GLPKSolverLP()))
L-Shaped Gap  Time: 0:00:00 (6 iterations)
  Objective:       -855.8333333333339
  Gap:             0.0
  Number of cuts:  8
  Iterations:      6
:Optimal
```
"""
mutable struct LShapedSolver{S <: SubSolver, E <: Execution, R <: AbstractRegularizer, A <: AbstractAggregator, C <: AbstractConsolidator} <: AbstractStructuredSolver
    lpsolver::MPB.AbstractMathProgSolver
    subsolver::S
    feasibility_cuts::Bool
    execution::E
    regularize::R
    aggregate::A
    consolidate::C
    crash::CrashMethod
    parameters::Dict{Symbol,Any}

    function LShapedSolver(lpsolver::MPB.AbstractMathProgSolver;
                           execution::Execution = Serial(),
                           feasibility_cuts::Bool = false,
                           regularize::AbstractRegularizer = DontRegularize(),
                           aggregate::AbstractAggregator = DontAggregate(),
                           consolidate::AbstractConsolidator = DontConsolidate(),
                           crash::CrashMethod = Crash.None(),
                           subsolver::SubSolver = lpsolver, kwargs...)
        S = typeof(subsolver)
        E = typeof(execution)
        R = typeof(regularize)
        A = typeof(aggregate)
        C = typeof(consolidate)
        return new{S, E, R, A, C}(lpsolver,
                                  subsolver,
                                  feasibility_cuts,
                                  execution,
                                  regularize,
                                  aggregate,
                                  consolidate,
                                  crash,
                                  Dict{Symbol,Any}(kwargs))
    end
end

function StructuredModel(stochasticprogram::StochasticProgram, solver::LShapedSolver)
    x₀ = solver.crash(stochasticprogram, solver.lpsolver)
    return LShaped(stochasticprogram, x₀, solver.lpsolver, get_solver(solver.subsolver), solver.feasibility_cuts, solver.execution, solver.regularize, solver.aggregate, solver.consolidate; solver.parameters...)
end

function add_params!(solver::LShapedSolver; kwargs...)
    push!(solver.parameters, kwargs...)
    for (k,v) in kwargs
        if k ∈ [:variant, :lpsolver, :subsolver, :feasibility_cuts, :execution, :regularize, :aggregate, :crash]
            setfield!(solver, k, v)
            delete!(solver.parameters, k)
        end
    end
    return nothing
end

function add_regularization_params!(solver::LShapedSolver; kwargs...)
    add_regularization_params!(solver.regularize; kwargs...)
end

function internal_solver(solver::LShapedSolver)
    return solver.lpsolver
end

function optimize_structured!(lshaped::AbstractLShapedSolver)
    return lshaped()
end

function fill_solution!(stochasticprogram::StochasticProgram, lshaped::AbstractLShapedSolver)
    # First stage
    first_stage = StochasticPrograms.get_stage_one(stochasticprogram)
    nrows, ncols = first_stage_dims(stochasticprogram)
    StochasticPrograms.set_decision!(stochasticprogram, decision(lshaped))
    μ = try
        MPB.getreducedcosts(lshaped.mastersolver.lqmodel)[1:ncols]
    catch
        fill(NaN, ncols)
    end
    StochasticPrograms.set_first_stage_redcosts!(stochasticprogram, μ)
    λ = try
        MPB.getconstrduals(lshaped.mastersolver.lqmodel)[1:nrows]
    catch
        fill(NaN, nrows)
    end
    StochasticPrograms.set_first_stage_duals!(stochasticprogram, λ)
    # Second stage
    fill_submodels!(lshaped, scenarioproblems(stochasticprogram))
    return nothing
end

function solverstr(solver::LShapedSolver)
    solver_str = "$(str(solver.execution))$(str(solver.regularize))"
    aggregate_str = str(solver.aggregate)
    if aggregate_str != ""
        return string(solver_str, " with ", aggregate_str)
    else
        return solver_str
    end
end

function solver_complexity(sp, ls)
    optimize!(sp, solver = ls)
    sm = spsolver_model(sp)
    return ncuts(sm)*niterations(sm)
end

function default_choice(given, default, null)
    if default isa null
        return given
    end
    return default
end


text:
struct Scene <: AbstractScene
    _vis::MeshCat.Visualizer
    _windows::Union{Blink.Window, Nothing}
    _root::AbstractSceneObject
    _root_transform::AbstractSceneObject

    function Scene(; openWindow::Bool = false)
        vis = MeshCat.Visualizer()

        win = nothing
        if (openWindow)
            window_defaults = Blink.@d(
                :title => "Glimmer Viewer",
                :width=>1200,
                :height=>800,
            )
            win = Blink.Window(window_defaults)

            Blink.AtomShell.opentools(win)
            open(vis, win)
        end

        root, root_transform = prepare_scene!(vis)

        return new(vis, win, root, root_transform)
    end

end

# returns the URL for the local server
url(s::Scene) = MeshCat.url(s._vis.core)

vis(s::Scene) = s._vis
win(s::Scene) = s._window
root(s::Scene) = s._root
root_transform(s::Scene) = s._root_transform

function clear(s::Scene)
    root_so = root(s)
    root_children = prop(root_so, :children)
    while(length(root_children) > 0)
        delete!(root_children[1])
        root_children = prop(root_so, :children)
    end

    # index = findfirst(x -> x===so, parent_children)
    # if (index === nothing)
    #     @error "Can't find the child $(name(so)) to delete"
    # end
    # deleteat!(parent_children, index)

    # v = vis(so)
    # delete!(v)
end

function prepare_scene!(vis::MeshCat.Visualizer)
    root_transform = EmptySceneObject(name="GlimmerTransform")
    root_transform_vis = vis["/Glimmer"]
    vis!(root_transform, root_transform_vis)

    root = EmptySceneObject(name="GlimmerRoot")
    root_vis = root_transform_vis["Root"]
    vis!(root, root_vis)

    # flip Y and Z axes
    # MeshCat.settransform!(root_transform_vis, tr2affine(Transform(unitX3(), unitZ3(), unitY3())))
    # rotation = SMatrix{3,3}(1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0) # X, Z, Y
    # MeshCat.settransform!(root_transform_vis, AffineMap(rotation, zero3()))
    set_Y_up!(root_transform_vis)

    return root, root_transform
end

function cameraTransform!(s::Scene, tr::AffineMap)
    cam_vis = vis(s)["/Cameras/default/rotated/<object>"]

    MeshCat.settransform!(cam_vis, tr)
end

function set_Y_up!(s::Scene)
    set_Y_up!(vis(s))
end

function set_Y_up!(vis::MeshCat.Visualizer)
    rotation = SMatrix{3,3}(1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0) # X, Z, Y
    MeshCat.settransform!(vis, AffineMap(rotation, zero3()))
end

function set_Z_up!(s::Scene)
    # MeshCat.settransform!(vis(s), tr2affine(Transform(unitX3(), unitY3(), unitZ3())))
    rotation = SMatrix{3,3}(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0) # X, Y, Z
    MeshCat.settransform!(vis(s), AffineMap(rotation, zero3()))
end

function grid!(s::Scene, visible::Bool)
    MeshCat.setprop!(vis(s)["/Grid"], "visible", visible)
end

function axes!(s::Scene, visible::Bool)
    MeshCat.setprop!(vis(s)["/Axes"], "visible", visible)
end

function cameraPlanes!(s::Scene, near::Float64 = 0.1, far::Float64 = 1000.0)
    cam_vis = vis(s)["/Cameras/default/rotated/<object>"]

    # MeshCat.setprop!(cam_vis, "zoom", 1.0)
    MeshCat.setprop!(cam_vis, "near", near)
    MeshCat.setprop!(cam_vis, "far", far)
end

ext: jl


text:
#=
    There are nine general cases for processing a*b
                   b <= 0    b <0< b    b >= 0
       a <=  0    LteLte      LteZer    LteGte
       a <0< a    ZerLte      ZerZer    ZerGte
       a >=  0    GteLte      GteZer    GteGte
=#

function (*){G1<:Grasp,G2<:Grasp,R<:Real}(a::Rvl{G1,R}, b::Rvl{G2,R})
    z = zero(R)
    if     a.hi <= z
       if     b.hi <= z
           mulLteLte(a,b)
       elseif b.lo >= z
           mulLteGte(a,b)
       else   # b straddles 0
           mulLteZer(a,b)
       end
    elseif a.lo >= z
       if     b.hi <= z
           mulGteLte(a,b)
       elseif b.lo >= z
           mulGteGte(a,b)
       else   # b straddles 0
           mulGteZer(a,b)
       end
    else   # a straddles 0
       if     b.hi <= z
           mulZerLte(a,b)
       elseif b.lo >= z
           mulZerGte(a,b)
       else   # b straddles 0
           mulZerZer(a,b)
       end
    end
end

for (fn,loa,lob,hia,hib) in [ (:mulLteLte, :(a.hi), :(b.hi), :(a.lo), :(b.lo)),
                              (:mulLteGte, :(a.lo), :(b.hi), :(a.hi), :(b.lo)),
                              (:mulLteZer, :(a.lo), :(b.hi), :(a.lo), :(b.lo)),
                              (:mulGteLte, :(a.hi), :(b.lo), :(a.lo), :(b.hi)),
                              (:mulGteGte, :(a.lo), :(b.lo), :(a.hi), :(b.hi)),
                              (:mulGteZer, :(a.hi), :(b.lo), :(a.hi), :(b.hi)),
                              (:mulZerLte, :(a.hi), :(b.lo), :(a.lo), :(b.lo)),
                              (:mulZerGte, :(a.lo), :(b.hi), :(a.hi), :(b.hi)),
                            ]
  @eval begin
    function ($fn){G1<:Grasp,G2<:Grasp,R<:Real}(a::Rvl{G1,R}, b::Rvl{G2,R})
        aLoIsOpen, aHiIsOpen = boundaries(G1)
        bLoIsOpen, bHiIsOpen = boundaries(G2)
        abGrasp = boundaries( (aLoIsOpen|bLoIsOpen), (aHiIsOpen|bHiIsOpen) )

        lo = (*)(($loa), ($lob), RoundDown)
        hi = (*)(($hia), ($hib), RoundUp)

        Rvl{abGrasp,R}(lo, hi)
    end
  end
end

(*){G<:Grasp,F<:AbstractFloat}(a::Rvl{G,F}, b::F) = (*)(a, Rvl{G,F}(b))
(*){G<:Grasp,F<:AbstractFloat}(a::F, b::Rvl{G,F}) = (*)(Rvl{G,F}(a), b)

(*){G<:Grasp,R<:Real}(a::Rvl{G,R}, b::R) = (*)(a, Rvl{G,R}(b))
# ambiguity with Bool
#(*){G<:Grasp,R<:Real}(a::R, b::Rvl{G,R}) = (*)(Rvl{G,R}(a), b)

(*){G<:Grasp,I<:WorkInt}(a::Rvl{G,Float64}, b::I) = (*)(a, Rvl{G,Float64}(convert(Float64,b)))
(*){G<:Grasp,I<:WorkInt}(a::I, b::Rvl{G,Float64}) = (*)(Rvl{G,Float64}(convert(Float64,a)), b)
(*){G<:Grasp,I<:WorkInt}(a::Rvl{G,Float32}, b::I) = (*)(a, Rvl{G,Float32}(convert(Float32,b)))
(*){G<:Grasp,I<:WorkInt}(a::I, b::Rvl{G,Float32}) = (*)(Rvl{G,Float32}(convert(Float32,a)), b)


text:
densityest_n_ds = [
    [
    (1, 3), # compile run
    (1, 125),
    (1, 250),
    (1, 500),
    (1, 1000),
    (1, 2000),
    (1, 3000),
    ],
    [
    (2, 2), # compile run
    (2, 10),
    (2, 20),
    (2, 30),
    (2, 40),
    (2, 50),
    ],
    [
    (3, 2), # compile run
    (3, 6),
    (3, 9),
    (3, 12),
    (3, 15),
    ],
    [
    (3, 2), # compile run
    (4, 4),
    (4, 6),
    (4, 8),
    ],
    [
    (3, 2), # compile run
    (8, 2),
    (8, 3),
    ],
    ]
densityest_insts(use_nat::Bool) = [
    [(500, n, 2d, false, use_nat, use_nat) for (n, d) in nds]
    for nds in densityest_n_ds
    ]

insts = OrderedDict()
insts["nat"] = (nothing, densityest_insts(true))
insts["ext"] = (:SOCExpPSD, densityest_insts(false))
return (DensityEstJuMP, insts)


text:
# Implement PatternToNumber

const NUCLEOTIDE_ORDER = Dict{Char, Int64}('A' => 0, 'C' => 1, 'G' => 2, 'T' => 3)

pattern_to_number(pattern) = foldl((acc, nucleotide) -> 4 * acc + NUCLEOTIDE_ORDER[nucleotide], pattern; init=0)

function main()
    dna = readline()
    result = pattern_to_number(dna)
    println(result)
end

main()


text:
using Test
using LinearAlgebra, LinearAlgebraX, Mods

@testset "Basics" begin
    A = big.(rand(Int, 10, 10) .% 100)

    include("../extras/hilbert.jl")

    H = hilbert(12)
    @test rankx(H) == 12

    @test detx(A) == det(Rational.(A))
    @test detx(A) == detx(A')
    @test A * invx(A) == eye(BigInt, 10)

    A = ones(Int, 3, 5)
    N = nullspacex(A)
    @test all(0 .== A * N)
end

@testset "Polynomials" begin
    using SimplePolynomials

    A = triu(ones(Int, 5, 5))
    p = char_poly(A)
    x = getx()
    @test p == (x - 1)^5

    A = ones(Int, 5, 5)
    p = char_poly(A)
    @test p == x^4 * (x - 5)
end

@testset "Modular" begin
    T = Mod{17}
    A = rand(T, 6, 6)
    p = char_poly(A)
    @test detx(A) == p(0)

    A = rand(Int, 20, 20) .% 100
    @test rankx(A) <= 20
end

@testset "Homogeneous" begin
    v = HVector(2, 3, 2)
    @test v[2] == 3 // 2
    @test dot(v, v) == 1
    w = HVector(20, 30, 20)
    @test v == w

    A = rand(Int, 3, 3) .% 20
    w = A * v
    @test w == HVector(A * v.data)
end

@testset "Permanent" begin
    A = ones(Int, 6, 6)
    @test permanent(A) == factorial(6)

    A = A - eye(Int, 6)
    @test permanent(A) == 265  # number of derangements of [6]

    A = triu(ones(5, 5))
    @test permanent(A) == 1
end


text:
"""
    Module Basis3D
Includes DG core functions for evaluating 3D bases on hexahedral elements
"""
module Basis3D

export vandermonde3D, grad_vandermonde3D
export chebyshev3D, equidistant_nodes3D, gauss_quad3D
export dPhi3D

import Basis1D
using LinearAlgebra, SparseArrays

"""
    vandermonde3D(p, x, y, z)
Return the 3D vandermonde matrix - the order 0:p Legendre polynomials
evaluated at 3D nodes (x,y,z)
"""
function vandermonde3D(p, x, y, z)
    V = Array{Float64,6}(undef,length(x),length(y),length(z),p+1,p+1,p+1)
    lx = Basis1D.vandermonde1D(p, x)
    ly = Basis1D.vandermonde1D(p, y)
    lz = Basis1D.vandermonde1D(p, z)
    for ipz = 1:p+1
        for ipy = 1:p+1
            for ipx = 1:p+1
                for iz = 1:length(z)
                    for iy = 1:length(y)
                        for ix = 1:length(x)
                            V[ix,iy,iz,ipx,ipy,ipz] = lx[ix,ipx]*ly[iy,ipy]*lz[iz,ipz]
                        end
                    end
                end
            end
        end
    end
    return V
end

"""
    grad_vandermonde3D(p, x, y, z)
Return the gradient of the order 0:p Legendre polynomials evaluated
at 3D nodes (x,y,z)
"""
function grad_vandermonde3D(p, x, y, z)
    dV = Array{Float64,7}(undef, length(x),length(y),length(z), p+1,p+1,p+1, 3)
    lx = Basis1D.vandermonde1D(p, x)
    ly = Basis1D.vandermonde1D(p, y)
    lz = Basis1D.vandermonde1D(p, z)
    dlx = Basis1D.grad_vandermonde1D(p, x)
    dly = Basis1D.grad_vandermonde1D(p, y)
    dlz = Basis1D.grad_vandermonde1D(p, z)
    # Derivative wrt x
    for ipz = 1:p+1
        for ipy = 1:p+1
            for ipx = 1:p+1
                for iz = 1:length(z)
                    for iy = 1:length(y)
                        for ix = 1:length(x)
                            dV[ix,iy,iz,ipx,ipy,ipz,1] = dlx[ix,ipx]*ly[iy,ipy]*lz[iz,ipz]
                        end
                    end
                end
            end
        end
    end
    # Derivative wrt y
    for ipz = 1:p+1
        for ipy = 1:p+1
            for ipx = 1:p+1
                for iz = 1:length(z)
                    for iy = 1:length(y)
                        for ix = 1:length(x)
                            dV[ix,iy,iz,ipx,ipy,ipz,2] = lx[ix,ipx]*dly[iy,ipy]*lz[iz,ipz]
                        end
                    end
                end
            end
        end
    end
    # Derivative wrt z
    for ipz = 1:p+1
        for ipy = 1:p+1
            for ipx = 1:p+1
                for iz = 1:length(z)
                    for iy = 1:length(y)
                        for ix = 1:length(x)
                            dV[ix,iy,iz,ipx,ipy,ipz,3] = lx[ix,ipx]*ly[iy,ipy]*dlz[iz,ipz]
                        end
                    end
                end
            end
        end
    end
    return dV
end

"""
    gauss_quad3D(p)
Return the 3D Gaussian quadrature points and weights (x,w) of order p on the domain [-1,1]^3
"""
function gauss_quad3D(p)
    (x1D,w1D) = Basis1D.gauss_quad(p)
    n = length(x1D)
    x3D = Array{Float64,4}(undef, 3, n,n,n)
    w3D = Array{Float64,3}(undef, n,n,n)
    for iz = 1:n
        for iy = 1:n
            for ix = 1:n
                x3D[1,ix,iy,iz] = x1D[ix]
                x3D[2,ix,iy,iz] = x1D[iy]
                x3D[3,ix,iy,iz] = x1D[iz]
                w3D[ix,iy,iz] = w1D[ix]*w1D[iy]*w1D[iz]
            end
        end
    end
    return (x3D, w3D)
end

"""
    chebyshev3D(p)
Return the 3D chebyshev points of order p on the domain [-1,1]^3
"""
function chebyshev3D(p)
    x1D = Basis1D.chebyshev(p)
    n = length(x1D)
    x3D = Array{Float64,4}(undef, 3, n,n,n)
    for iz = 1:n
        for iy = 1:n
            for ix = 1:n
                x3D[1,ix,iy,iz] = x1D[ix]
                x3D[2,ix,iy,iz] = x1D[iy]
                x3D[3,ix,iy,iz] = x1D[iz]
            end
        end
    end
    return x3D
end

"""
    equidistant_nodes3D(p)
Return the "order p" 3D points on the domain [-1,1]^3
"""
function equidistant_nodes3D(p)
    x1D = -1:2/p:1
    n = length(x1D)
    x3D = Array{Float64,4}(undef, 3, n,n,n)
    for iz = 1:n
        for iy = 1:n
            for ix = 1:n
                x3D[1,ix,iy,iz] = x1D[ix]
                x3D[2,ix,iy,iz] = x1D[iy]
                x3D[3,ix,iy,iz] = x1D[iz]
            end
        end
    end
    return x3D
end

"""
    interpolation_matrix3D(xyz_from, xyz_to)
Compute an interpolation matrix from a set of 3D points to another set of 3D points.
Assumes that points interpolating from provide enough accuracy (aka - they are
well spaced out and of high enough order), and define a cube. Points
interpolating onto can be of any size, but must be defined on this same cube.
Interpolation matrix ∈ ℜ^(3D size of xyz_to x 3D size of xyz_from)
"""
function interpolation_matrix3D(xyz_from, xyz_to)
    # Create nodal representation of reference bases
    (x_from, y_from, z_from) = xyz_from
    order = size(x_from,1) - 1 # Assumes order = (size of x_from) - 1
    @assert size(x_from,1) == size(y_from,1) == size(z_from,1)
    n_from = size(x_from,1)*size(y_from,1)*size(z_from,1)
    l_from = vandermonde3D(order, x_from, y_from, z_from)

    eye = diagm(0=>ones(n_from))
    V = reshape(l_from, n_from,n_from)
    coeffs_phi = V \ eye

    # Compute reference bases on the output points
    l_to = vandermonde3D(order, xyz_to...)
    n_to = prod(size.(xyz_to,1))
    V_to = reshape(l_to, n_to,n_from)

    # Construct interpolation matrix
    Interp3D = V_to*coeffs_phi
end

"""
    dPhi3D(xyz_from, xyz_to)
Compute the gradient of the basis functions defined by 3D points on another
set of 3D points.

Assumes that points interpolating from provide enough accuracy (aka - they are
well spaced out and of high enough order), and define a cube. Points
interpolating onto can be of any size, but must be defined on this same cube.
dPhi ∈ ℜ^(size of xyz_to × size of xyz_from × 3)
"""
function dPhi3D(xyz_from, xyz_to)
    dim = 3
    # Create nodal representation of reference bases
    (x_from, y_from, z_from) = xyz_from
    order = size(x_from,1) - 1 # Assumes order = (size of x_from) - 1
    @assert size(x_from,1) == size(y_from,1) == size(z_from,1)
    n_from = size(x_from,1)*size(y_from,1)*size(z_from,1)
    l_from = vandermonde3D(order, x_from, y_from, z_from)

    eye = diagm(0=>ones(n_from))
    V = reshape(l_from, n_from,n_from)
    coeffs_phi = V \ eye

    # Compute derivative of reference bases on the output points
    dl_to = grad_vandermonde3D(order, xyz_to...)
    n_to = prod(size.(xyz_to,1))
    dV = reshape(dl_to, n_to,n_from,dim)
    # Construct gradient of phi = dV*coeffs_phi
    dPhi_to = Array{Float64,3}(undef, n_to,n_from,dim)
    for l = 1:dim
        dPhi_to[:,:,l] = dV[:,:,l]*coeffs_phi
    end
    return dPhi_to
end

end

max_stars_repo_name: NoseKnowsAll/DGToolkit

max_issues_repo_name: NoseKnowsAll/DGToolkit

max_forks_repo_name: NoseKnowsAll/DGToolkit


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    UniformPartition(k, [shuffle])

A method for partitioning spatial objects uniformly into `k` subsets
of approximately equal size. Optionally `shuffle` the data (default
to `true`).
"""
struct UniformPartition <: PartitionMethod
  k::Int
  shuffle::Bool
end

UniformPartition(k::Int) = UniformPartition(k, true)

function partition(object, method::UniformPartition)
  n = nelements(object)
  k = method.k

  @assert k ≤ n "number of subsets must be smaller than number of points"

  inds = method.shuffle ? shuffle(1:n) : collect(1:n)
  subsets = collect(Iterators.partition(inds, n ÷ k))

  Partition(object, subsets)
end


text:
@everywhere begin
    using DrWatson
    import JSON
    using Logging
    using HDF5
    using PWS
end

f = ARGS[1]
dict = JSON.parsefile(projectdir("_research", "tmp", f))["params"]
@info "Read file" file = projectdir("_research", "tmp", f)

duration = dict["duration"]

num_responses = dict["num_responses"]
run_name = dict["run_name"]

mean_s = dict["mean_s"]
corr_time_s = dict["corr_time_s"]
corr_time_ratio = dict["corr_time_ratio"]
dtimes = collect(0.0:0.1:duration)

save_dict = Dict("Alg" => dict["algorithm"], "Duration" => duration, "Smean" => mean_s, "Xtimescale" => corr_time_s / corr_time_ratio)

if dict["algorithm"] == "ti"
    algorithm = TIEstimate(0, 16, dict["ti_samples"])
    save_dict["M"] = 16 * dict["ti_samples"]
elseif dict["algorithm"] == "annealing"
    algorithm = AnnealingEstimate(5, 100, 100)
elseif dict["algorithm"] == "directmc"
    algorithm = DirectMCEstimate(dict["directmc_samples"])
    save_dict["M"] = dict["directmc_samples"]
elseif dict["algorithm"] == "smc"
    algorithm = SMCEstimate(dict["smc_samples"])
    save_dict["M"] = dict["smc_samples"]
else
    error("Unsupported algorithm " * dict["algorithm"])
end

@info "Parameters" run_name duration num_responses algorithm mean_s corr_time_s corr_time_ratio

system_fn = () -> PWS.gene_expression_system(
    mean_s=mean_s,
    corr_time_s = corr_time_s,
    corr_time_x = corr_time_s / corr_time_ratio,
    dtimes=dtimes
)

mi = PWS.run_parallel(system_fn, algorithm, num_responses)
result = Dict("Samples" => mi, "DiscreteTimes" => dtimes)

function DrWatson._wsave(filename, result::Dict)
    h5open(filename, "w") do file
        PWS.write_hdf5!(file, result)
    end
end

filename = savename(save_dict, "hdf5")
local_path = datadir("gene_expression", run_name, filename)
tagsave(local_path, merge(dict, result), storepatch=false)
@info "Saved to" filename

# upload to SUN storage

include("smbclient.jl")
sun_path = joinpath(sun_home, "data", dict["algorithm"], run_name, filename)
mkpath(dirname(sun_path))
cp(local_path, sun_path)


text:
using Zygote
using LinearAlgebra
using GalacticOptim,NLopt,Optim
using Plots
using GLOQ

# System parameters for a simple two level open quantum system
N_states = 2; # number of states
freqs = [4.1] # transition frequency in GHz
omegas = 2.0*pi.*freqs # change to angular frequency
gamma1   = [1.0/(45.0*GLOQ.GLOQ_MICRO_SEC)] # Reciprocal of relaxation time - T1 (in units of ns)
gamma2   = [1.0/(24.0*GLOQ.GLOQ_MICRO_SEC)] # Reciprocal of dephasing time - T2 (in units of ns)
omr_ramsey = 2.0*pi*(4.1 - 5.0e-4) # drive frequency for the Ramsey experiment
omr_echo = 2.0*pi*4.1  # drive frequency for the Echo experiment
omr_t1 = 2.0*pi*4.1  # drive frequency for the T1 experiment
TC = 2.5*17.0 # total control time

# Initial state
initial_state = 0
state_u0 = [0.0;0.0]
state_v0 = [0.0;0.0]
state_u0[initial_state+1] = 1.0

# Forward solve to generate synthetic data for the Ramsey experiment
# Duration of the Ramsey experiment, largest dark time
T_Ramsey = 10.0*GLOQ.GLOQ_MICRO_SEC # convert micro-sec to nano-sec
# total number of dark time samples
N_dark_ramsey = 201
dt_ramsey = T_Ramsey/N_dark_ramsey
t_dark_ramsey = collect(range(0.0, T_Ramsey, length=N_dark_ramsey))
# Forward solve
rho_synthetic_ramsey_u,rho_synthetic_ramsey_v = GLOQ.RamseyForwardSolve(
				 state_u0,state_v0, # initial values, u for the real part, v for the imaginary part
			     omegas,omr_ramsey, # transition frequencies, drive frequency
				 gamma1,gamma2, # decay and dephasing
				 initial_state, # initial state
				 TC,t_dark_ramsey,N_states) # control time, dark time, total number of states
population_ramsey_synthetic = GLOQ.get_population(rho_synthetic_ramsey_u)

# Add noise to the synthetic Echo data
noisy_ramsey = copy(population_ramsey_synthetic)
additive_noise_ramsey = 0.025*randn(N_dark_ramsey)
noisy_ramsey .+= additive_noise_ramsey

# Physically, noisy data of population must be a probability.
# Shift and rescale the data so it is between [0,1]
for j = 1:N_states
	shift = minimum(noisy_ramsey[:,j])
	if(shift<0.0)
		noisy_ramsey[:,j] .-= shift
	end
end

for i = 1:N_dark_ramsey
	noisy_ramsey[i,:] ./= sum(noisy_ramsey[i,:])
end

# Echo experiment
# Duration of the Echo experiment, largest dark time
T_Echo = 30.0*GLOQ.GLOQ_MICRO_SEC # convert micro-sec to nano-sec
# total number of dark time samples
N_dark_echo = 601
# Step size of the Echo experiment
dt_echo = T_Echo/N_dark_echo
t_dark_echo = collect(range(0.0, T_Echo, length=N_dark_echo))
rho_synthetic_echo_u,rho_synthetic_echo_v = GLOQ.EchoForwardSolve(
				 state_u0,state_v0, # initial values, u for the real part, v for the imaginary part
			     omegas,omr_echo, # transition frequencies, drive frequency
				 gamma1,gamma2, # decay and dephasing
				 initial_state, # initial state
				 TC,t_dark_echo,N_states) # control time, dark time, total number of states
population_echo_synthetic = GLOQ.get_population(rho_synthetic_echo_u)

# Add noise to the synthetic Echo data
noisy_echo = copy(population_echo_synthetic)
additive_noise_echo = 0.025*randn(N_dark_echo)
noisy_echo .+= additive_noise_echo

# Physically, noisy data of population must be a probability.
# Shift and rescale the data so it is between [0,1]
for j = 1:N_states
	shift = minimum(noisy_echo[:,j])
	if(shift<0.0)
		noisy_echo[:,j] .-= shift
	end
end

for i = 1:N_dark_echo
	noisy_echo[i,:] ./= sum(noisy_echo[i,:])
end

# T1
# Duration of the T1-decay experiment, largest dark time
T_t1 = 40.0*GLOQ.GLOQ_MICRO_SEC # convert micro-sec to nano-sec
# total number of dark time samples
N_dark_t1 = 401
# Step size of the T1 experiment
dt_t1 = T_t1/N_dark_t1
t_dark_t1 = collect(range(0.0, T_t1, length=N_dark_t1))
rho_synthetic_t1_u,rho_synthetic_t1_v = GLOQ.T1ForwardSolve(
				 state_u0,state_v0, # initial values, u for the real part, v for the imaginary part
			     omegas,omr_t1, # transition frequencies, drive frequency
				 gamma1,gamma2, # decay and dephasing
				 initial_state, # initial state
				 TC,t_dark_t1,N_states) # control time, dark time, total number of states
population_t1_synthetic = GLOQ.get_population(rho_synthetic_t1_u)

# Add noise to the synthetic T1 data
noisy_t1 = copy(population_t1_synthetic)
additive_noise_t1 = 0.025*randn(N_dark_t1)
noisy_t1 .+= additive_noise_t1

# Physically, noisy data of population must be a probability.
# Shift and rescale the data so it is between [0,1]
for j = 1:N_states
	shift = minimum(noisy_t1[:,j])
	if(shift<0.0)
		noisy_t1[:,j] .-= shift
	end
end

for i = 1:N_dark_t1
	noisy_t1[i,:] ./= sum(noisy_t1[i,:])
end

# Define the loss function for the GalacticOptim
# p: phyiscal parameters:
#	 p[1] = transition frequency in GHz
#    p[2] = gamma1
#    p[3] = gamma2
# dummy_parameter: needed by GalacticOptim, one can just put [] here
function loss(p,dummy_parameter)
	# Ramsey
	_rho_ramsey_u,_rho_ramsey_v = GLOQ.RamseyForwardSolve(state_u0,state_v0,
				     (2*pi).*[p[1]],omr_ramsey,
					 [p[2]],[p[3]],#gamma1,gamma2,
					 initial_state, # initial state
					 TC,t_dark_ramsey,N_states)
	_population_ramsey = GLOQ.get_population(_rho_ramsey_u)
	# Echo
    _rho_echo_u,_rho_echo_v = GLOQ.EchoForwardSolve(state_u0,state_v0,
				     (2*pi).*[p[1]],omr_echo,
					 [p[2]],[p[3]],#gamma1,gamma2,
					 initial_state, # initial state
					 TC,t_dark_echo,N_states)
	_population_echo = GLOQ.get_population(_rho_echo_u)
	# T1
    _rho_t1_u,_rho_t1_v = GLOQ.T1ForwardSolve(state_u0,state_v0,
				     (2*pi).*[p[1]],omr_t1,
					 [p[2]],[p[3]],#gamma1,gamma2,
					 initial_state, # initial state
					 TC,t_dark_t1,N_states)
	_population_t1 = GLOQ.get_population(_rho_t1_u)

	_loss = sum(abs2,_population_ramsey-noisy_ramsey)*dt_ramsey+
			sum(abs2,_population_echo-noisy_echo)*dt_echo+
			sum(abs2,_population_t1-noisy_t1)*dt_t1

	return _loss
end

plot_callback = function(p,other_args)
	# Ramsey
	rho_ramsey_u,rho_ramsey_v = GLOQ.RamseyForwardSolve(state_u0,state_v0,
					 (2*pi).*[p[1]],omr_ramsey,
					 [p[2]],[p[3]],#gamma1,gamma2,
					 initial_state, # initial state
					 TC,t_dark_ramsey,N_states)
	population_ramsey = GLOQ.get_population(rho_ramsey_u)
	# Echo
    rho_echo_u,rho_echo_v = GLOQ.EchoForwardSolve(state_u0,state_v0,
				     (2*pi).*[p[1]],omr_echo,
					 [p[2]],[p[3]],#gamma1,gamma2,
					 initial_state, # initial state
					 TC,t_dark_echo,N_states)
	population_echo = GLOQ.get_population(rho_echo_u)

	# T1
    rho_t1_u,rho_t1_v = GLOQ.T1ForwardSolve(state_u0,state_v0,
				     (2*pi).*[p[1]],omr_t1,
					 [p[2]],[p[3]],#gamma1,gamma2,
					 initial_state, # initial state
					 TC,t_dark_t1,N_states)
	population_t1 = GLOQ.get_population(rho_t1_u)

	# Plot
	# Ramsey
	fig_ramsey = plot(t_dark_ramsey./GLOQ.GLOQ_MICRO_SEC,noisy_ramsey,label=["Noisy-0" "Noisy-1"],
			          line = (:dash,0.0), marker = ([:hex :hex], 5, 0.5),legend=:outerright,
					  title="Ramsey");
	plot!(fig_ramsey,t_dark_ramsey./GLOQ.GLOQ_MICRO_SEC,population_ramsey,label=["Noisy-0" "Opt-1"],legend=:outerright);
	# Echo
	fig_echo = plot(t_dark_echo./GLOQ.GLOQ_MICRO_SEC,noisy_echo,label=["Noisy-0" "Noisy-1"],
		  		     line = (:dash,0.0), marker = ([:hex :hex], 5, 0.5),legend=:outerright,
					 title="Echo");
	plot!(fig_echo,t_dark_echo./GLOQ.GLOQ_MICRO_SEC,population_echo,label=["Opt-0" "Opt-1"],legend=:outerright);
	# T1
	fig_t1=plot(t_dark_t1./GLOQ.GLOQ_MICRO_SEC,noisy_t1,label=["Noisy-0" "Noisy-1"],
			    line = (:dash,0.0), marker = ([:hex :hex], 5, 0.5),legend=:outerright,
				title="T1");
	plot!(fig_t1,t_dark_t1./GLOQ.GLOQ_MICRO_SEC,population_t1,label=["Opt-0" "Opt-1"],legend=:outerright);

	display( plot(fig_ramsey,fig_echo,fig_t1,layout=grid(3,1),size=[1000,1500],
				  legendfontsize=15,xtickfontsize=15,ytickfontsize=15,titlefontsize=18) )
	return false
end

p_true = [freqs;gamma1;gamma2] # values to generate synthetic data
# initial guess for the optimization
p_initial = [freqs.-0.5e-4;0.75.*gamma1;1.25.*gamma2]
# bounds for the optimization
lower_bound = (0.5).*p_true
upper_bound = (1.5).*p_true

# construct optimization object, use Zygote auto-differentiation to compute the gradient
loss_gradient = GalacticOptim.OptimizationFunction(loss, GalacticOptim.AutoZygote())
opt_prob = GalacticOptim.OptimizationProblem(loss_gradient, p_initial,
										 lb = lower_bound, ub = upper_bound)

println("Optim Fminbox(LBFGS) Optimization starts")
@time sol = GalacticOptim.solve(opt_prob ,Fminbox(LBFGS()),
								cb = plot_callback,
								outer_iterations = 20,
								iterations = 10,
								show_trace=true,
								f_tol = 1e-5,
								outer_f_tol = 1e-5)
println("Optim Fminbox(LBFGS) Optimization done")

# present the solutions
println("\nOptimized results: ",sol.u,
        "\nLoss: ",sol.minimum,
		"\nError: ",sol.u-p_true)


text:
# test not exported
@test_throws UndefVarError Pipeline
@test_throws UndefVarError AbstractPipeline
@test_throws UndefVarError ImmutablePipeline
@test AugmentorBase.AbstractPipeline <: Any
@test AugmentorBase.Pipeline <: AugmentorBase.AbstractPipeline
@test AugmentorBase.ImmutablePipeline <: AugmentorBase.Pipeline

@test typeof((FlipX(),FlipY())) <: AugmentorBase.AbstractPipeline
@test AugmentorBase.operations((FlipX(),FlipY())) === (FlipX(),FlipY())

@testset "ImmutablePipeline" begin
    @test_throws MethodError AugmentorBase.ImmutablePipeline()
    @test_throws MethodError AugmentorBase.ImmutablePipeline(())
    @test_throws MethodError AugmentorBase.ImmutablePipeline(1)
    @test_throws MethodError AugmentorBase.ImmutablePipeline((1,))
    @test_throws MethodError AugmentorBase.ImmutablePipeline{1}((1,))
    @test_throws MethodError AugmentorBase.ImmutablePipeline{2}((FlipX(),))
    @test_throws MethodError AugmentorBase.ImmutablePipeline{1}(FlipX())

    @test @inferred(AugmentorBase.ImmutablePipeline(FlipX())) === @inferred(AugmentorBase.ImmutablePipeline{1}((FlipX(),)))

    p = @inferred AugmentorBase.ImmutablePipeline(FlipX())
    @test p === @inferred(AugmentorBase.ImmutablePipeline((FlipX(),)))
    @test typeof(p) <: AugmentorBase.ImmutablePipeline{1}
    @test @inferred(length(p)) === 1
    @test @inferred(AugmentorBase.operations(p)) === (FlipX(),)

    p = @inferred AugmentorBase.ImmutablePipeline(FlipX(), FlipY())
    @test p === @inferred(AugmentorBase.ImmutablePipeline((FlipX(),FlipY())))
    @test typeof(p) <: AugmentorBase.ImmutablePipeline{2}
    @test @inferred(length(p)) === 2
    @test @inferred(AugmentorBase.operations(p)) === (FlipX(),FlipY())

    p = @inferred AugmentorBase.ImmutablePipeline(FlipX(), FlipY(), Rotate90())
    @test p === @inferred(AugmentorBase.ImmutablePipeline((FlipX(),FlipY(),Rotate90())))
    @test typeof(p) <: AugmentorBase.ImmutablePipeline{3}
    @test @inferred(length(p)) === 3
    @test @inferred(AugmentorBase.operations(p)) === (FlipX(),FlipY(),Rotate90())
end

@testset "ImmutablePipeline with |>" begin
    buf = rand(2,2)

    p = @inferred(FlipX() |> FlipY())
    @test p === AugmentorBase.ImmutablePipeline(FlipX(), FlipY())

    p = @inferred(FlipX() |> buf |> FlipY())
    @test p === AugmentorBase.ImmutablePipeline(FlipX(), CacheImage(buf), FlipY())

    p = @inferred(FlipX() |> CacheImage(buf) |> FlipY())
    @test p === AugmentorBase.ImmutablePipeline(FlipX(), CacheImage(buf), FlipY())

    p = @inferred(FlipX() |> CacheImage() |> FlipY())
    @test p === AugmentorBase.ImmutablePipeline(FlipX(), CacheImage(), FlipY())

    p = @inferred(FlipX() |> FlipY() |> buf)
    @test p === AugmentorBase.ImmutablePipeline(FlipX(), FlipY(), CacheImage(buf))

    p = @inferred(FlipX() |> NoOp() |> FlipY())
    @test p === AugmentorBase.ImmutablePipeline(FlipX(), NoOp(), FlipY())

    p = @inferred((FlipX() |> NoOp()) |> FlipY())
    @test p === AugmentorBase.ImmutablePipeline(FlipX(), NoOp(), FlipY())

    p = @inferred(FlipX() |> (NoOp() |> FlipY()))
    @test p === AugmentorBase.ImmutablePipeline(FlipX(), NoOp(), FlipY())

    p = @inferred(FlipX() |> NoOp() |> FlipY() |> Rotate90())
    @test p === AugmentorBase.ImmutablePipeline(FlipX(), NoOp(), FlipY(), Rotate90())

    p = @inferred((FlipX() |> NoOp()) |> (FlipY() |> Rotate90()))
    @test p === AugmentorBase.ImmutablePipeline(FlipX(), NoOp(), FlipY(), Rotate90())

    p = FlipX() * FlipY() |> Rotate90() * Rotate270()
    @test p == AugmentorBase.ImmutablePipeline(Either(FlipX(),FlipY()), Either(Rotate90(),Rotate270()))

    p = NoOp() * FlipX() * FlipY() |> Rotate90() * Rotate270()
    @test p == AugmentorBase.ImmutablePipeline(Either(NoOp(),FlipX(),FlipY()), Either(Rotate90(),Rotate270()))
end

@testset "Pipeline constructor" begin
    @test_throws MethodError AugmentorBase.Pipeline()
    @test_throws MethodError AugmentorBase.Pipeline(())
    @test_throws MethodError AugmentorBase.Pipeline(1)
    @test_throws MethodError AugmentorBase.Pipeline((1,))

    p = @inferred AugmentorBase.Pipeline(FlipX())
    @test typeof(p) <: AugmentorBase.ImmutablePipeline
    @test p == AugmentorBase.ImmutablePipeline(FlipX())

    p = @inferred AugmentorBase.Pipeline((FlipX(),))
    @test typeof(p) <: AugmentorBase.ImmutablePipeline
    @test p == AugmentorBase.ImmutablePipeline(FlipX())

    p = @inferred AugmentorBase.Pipeline(FlipX(), FlipY())
    @test typeof(p) <: AugmentorBase.ImmutablePipeline
    @test p == AugmentorBase.ImmutablePipeline(FlipX(),FlipY())

    p = @inferred AugmentorBase.Pipeline((FlipX(), FlipY()))
    @test typeof(p) <: AugmentorBase.ImmutablePipeline
    @test p == AugmentorBase.ImmutablePipeline(FlipX(),FlipY())
end


text:
# go back to the nth parent directory
function back2(this_path, n)
return back2_core(this_path, n)
end

# go back to the ith parent directory and precede to directory_ahead
function back2(this_path, n, directory_ahead...)
this_path = back2_core(this_path, n);
return joinpath(this_path,directory_ahead...)
end

function back2_core(this_path, n)
for i = 1:n
	this_path = dirname(this_path);
end
return this_path
end


text:
Mat2x2(T) = T[
  0 0
  0 0
]

Mat3x3(T) = T[
  0 0 0
  0 0 0
  0 0 0
]

Mat4x4(T) = T[
  0 0 0 0
  0 0 0 0
  0 0 0 0
  0 0 0 0
]

const Mat4x4f = SMatrix{4,4,Float32}

const zerosMat4x4f = Mat4x4(Float32)
const eyeMat4x4f = eye(zerosMat4x4f)

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
using IterableTables
using NamedTuples
using Gadfly
using DataTables
using Base.Test

@testset "Gadfly" begin

source_dt = DataTable(a=[1,2,3], b=[4,2,6])

p = Gadfly.plot(source_dt, x=:a, y=:b, Geom.line)

@test isa(p, Gadfly.Plot)

end


text:
# Abstract type: whether the model has random effects
abstract type RandomEffects end
struct hasRandomEffects <: RandomEffects end
struct noRandomEffects <: RandomEffects end

# Model
abstract type LRMoEModel{re <: RandomEffects} end

const LRMoEModelSTD = LRMoEModel{noRandomEffects}
const LRMoEModelRE = LRMoEModel{hasRandomEffects}

# Result
abstract type LRMoEFittingResult end

# Model contains: regression coefficients, component distributions
#   number of iteration, convergence, likelihood, AIC, BIC -> Move to another struct
struct LRMoESTD <: LRMoEModelSTD
    α::Array
    comp_dist::Array
    LRMoESTD(α, comp_dist) = ( size(α)[1] == size(comp_dist)[2] ? new(α, comp_dist) : error("Invalid specification of model.") )
end

struct LRMoESTDFit <: LRMoEFittingResult
    model_fit::LRMoESTD
    converge::Bool
    iter::Integer
    loglik::Real
    loglik_np::Real
    AIC::Real
    BIC::Real
    LRMoESTDFit(model_fit, converge, iter, loglik, loglik_np, AIC, BIC) = new(model_fit, converge, iter, loglik, loglik_np, AIC, BIC)
end

"""
    summary(obj)

Summarizes a fitted LRMoE model.

# Arguments
- `obj`: An object returned by `fit_LRMoE` function.

# Return Values
Prints out a summary of the fitted LRMoE model on screen.
"""
function summary(m::LRMoESTDFit)
    println("Model: LRMoE")
    if m.converge
        println("Fitting converged after $(m.iter) iterations")
    else
        println("Fitting NOT converged after $(m.iter) iterations")
    end
    println("Dimension of response: $(size(m.model_fit.comp_dist)[1])")
    println("Number of components: $(size(m.model_fit.comp_dist)[2])")
    println("Loglik: $(m.loglik)")
    println("Loglik (no penalty): $(m.loglik_np)")
    println("AIC: $(m.AIC)")
    println("BIC: $(m.BIC)")
    println("Fitted α:")
    println("$(m.model_fit.α)")
    println("Fitted component distributions:")
    println("$(m.model_fit.comp_dist)")
end

# tmp = LRMoESTD([1 2 ; 3 4], ["a" "b"])
# tmp1 = LRMoESTDFit(tmp, true, 20, 1.1, 1.2, 1.3)
# summary(tmp1)

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
using Test
using SafeTestsets

@safetestset "Quaternion Tests" begin
    include("quaternion_test.jl")
end

include("diff/diff_test.jl")

@safetestset "Factorization Test" begin
    include("factorization_test.jl")
end

include("initialization/initialization_test.jl")

@safetestset "UI Test" begin
    include("ui_test.jl")
end

@safetestset "URDF Tests" begin
    include("urdf_test.jl")
end

@safetestset "Optionals Tests" begin
    include("optionals_test.jl")
end

@safetestset "Dynamics Tests" begin
    include("dynamics/dynamics_test.jl")
end

@safetestset "Example Tests" begin
    include("example_test.jl")
end


text:
using Base.Test
cd(Pkg.dir("WebSockets", "test"))

@sync yield() # avoid mixing of  output with possible deprecation warnings from .juliarc

@testset "WebSockets" begin
info("\nFragment and unit tests\n")
@testset "Fragment and frame unit tests" begin
    include("frametest.jl");sleep(1)
end

info("\nHttpServer and HTTP handshake\n")
@testset "HttpServer and HTTP handshake" begin
    include("handshaketest.jl");sleep(1)
end

info("\nClient-server test, HTTP client\n")
@testset "Client-server test, HTTP client" begin
    include("client_server_test.jl");sleep(1)
end

info("\nClient test, HTTP client\n")
@testset "Client test, HTTP client" begin
    include("client_test.jl");sleep(1)
end

info("\nAbrupt close & error handling test\n")
@testset "Abrupt close & error handling test" begin
    include("error_test.jl");sleep(1)
end

end
# TODO
# restructure browsertests


text:
function construct_RICE_objective(m::Model,t_choice::Int)

    # Find number of timesteps across model time horizon.
    n_steps = length(dim_keys(m, :time))

    # Pre-allocate matrix to store optimal tax and mitigation rates.
    MIU = ones(n_steps,12)
    MIU[1,:] .= 0

    # Create a function to optimize user-specified model for (i) revenue recycling and (ii) reference case.
    RICE_objective =

    function RICE2010Welfare(MIUvec::Array{Float64,1})
        for j = 1:12
            MIU[2:(t_choice+1),j] = MIUvec[((j-1)*t_choice+1):j*t_choice] #imposes 0 MIU in year 2005.
        end
        set_param!(m, :emissions, :MIU, MIU)
        run(m)
        return(m[:welfare, :welfare])
	end

    # Return the objective function.
    return RICE_objective
end

function vec2array(vec::Array{Float64,1},max_period::Int)
    arr = ones(60,12)
    arr[1,:] .= 0
    for j = 1:12
        arr[2:(max_period+1),j] = vec[((j-1)*max_period+1):j*max_period] #imposes 0 MIU in year 2005.
    end
    return arr
end

function retConstraint(m::Model,max_t)
    run(m)
    backstop = m[:emissions,:pbacktime][2:(max_t+1),:]
    function nlconst(result::Vector,vect::Vector)
        arr = reshape(vect,:,12)
        prices = backstop .* arr .^ 1.8
        result[:] =  std(prices,dims=2) ./ mean(prices,dims=2) .- 0.3
    end
end

#########################################################################################################################################
#all below here is inherited from David/Mimi and used in parameter extraction from Nordhaus' Excel sheet
function getindexfromyear_rice_2010(year)
    baseyear = 2005

    if rem(year - baseyear, 10) != 0
        error("Invalid year")
    end

    return div(year - baseyear, 10) + 1
end

#Function to read a single parameter value from original RICE 2010 model.
function getparam_single(f, range::AbstractString, regions)
    vals= Array{Float64}(undef, length(regions))
    for (i,r) = enumerate(regions)
        data=readxl(f,"$(r)!$(range)")
        vals[i]=data[1]
    end
    return vals
end

#Function to read a time series of parameter values from original RICE 2010 model.
function getparam_timeseries(f, range::AbstractString, regions, T)
    vals= Array{Float64}(undef, T, length(regions))
    for (i,r) = enumerate(regions)
        data=readxl(f,"$(r)!$(range)")
        for n=1:T
            vals[n,i] = data[n]
        end
    end
    return vals
end

max_stars_repo_name: fdennig/RICEupdate

max_issues_repo_name: fdennig/RICEupdate

max_forks_repo_name: fdennig/RICEupdate


text:
export Event, matrix, vector, isfilled

mutable struct Event
    user::Int
    item::Int
    value::Float64 # e.g. rating, 0/1
end

function matrix(m::Int, n::Int, t::Type{<:Number}=Float64)
    Array{Union{Nothing, t}}(nothing, m, n)
end

function vector(m::Int, t::Type{<:Number}=Float64)
    Array{Union{Nothing, t}}(nothing, m)
end

function isfilled(a::AbstractArray)
    nothing ∉ Set(a)
end


text:
const BLOCKSIZE = 1024

Base.mapreduce(f, op, A::AbstractStridedView; dims=:, kw...) =
    Base._mapreduce_dim(f, op, kw.data, A, dims)

Base._mapreduce_dim(f, op, nt::NamedTuple{(:init,)}, A::AbstractStridedView, ::Colon) =
    _mapreduce(f, op, A, nt)
Base._mapreduce_dim(f, op, ::NamedTuple{()}, A::AbstractStridedView, ::Colon) =
    _mapreduce(f, op, A)

Base._mapreduce_dim(f, op, nt::NamedTuple{(:init,)}, A::AbstractStridedView, dims) =
    Base.mapreducedim!(f, op, Base.reducedim_initarray(A, dims, nt.init), A)
Base._mapreduce_dim(f, op, ::NamedTuple{()}, A::AbstractStridedView, dims) =
    Base.mapreducedim!(f, op, Base.reducedim_init(f, op, A, dims), A)

function Base.map(f::F, a1::AbstractStridedView{<:Any,N},
        A::Vararg{AbstractStridedView{<:Any,N}}) where {F,N}
    T = Base.promote_eltype(a1, A...)
    map!(f, similar(a1, T), a1, A...)
end

function Base.map!(f::F, b::AbstractStridedView{<:Any,N}, a1::AbstractStridedView{<:Any,N},
        A::Vararg{AbstractStridedView{<:Any,N}}) where {F,N}
    dims = size(b)

    # Check dimesions
    size(a1) == dims || throw(DimensionMismatch())
    for a in A
        size(a) == dims || throw(DimensionMismatch())
    end

    any(isequal(0), dims) && return b # don't do anything

    _mapreduce_fuse!(f, nothing, nothing, dims, (b, a1, A...))

    return b
end

function _mapreduce(f, op, A::AbstractStridedView{T}, nt = nothing) where {T}
    if length(A) == 0
        b = Base.mapreduce_empty(f, op, T)
        return nt === nothing ? b : op(b, nt.init)
    end

    dims = size(A)
    a = Base.mapreduce_first(f, op, first(A))
    a2 = nt === nothing ? a : op(a, nt.init)
    out = similar(A, typeof(a2), (1,))
    if nt === nothing
        _init_reduction!(out, f, op, a)
    else
        out[ParentIndex(1)] = nt.init
    end
    _mapreducedim!(f, op, nothing, dims, (sreshape(out, one.(dims)), A))
    return out[ParentIndex(1)]
end

@inline function Base.mapreducedim!(f, op, b::AbstractStridedView{<:Any,N},
        a1::AbstractStridedView{<:Any,N}, A::Vararg{AbstractStridedView{<:Any,N}}) where {N}
    outdims = size(b)
    dims = map(max, outdims, map(max, map(size, (a1, A...))...))

    # Check dimesions
    Broadcast.check_broadcast_axes(map(Base.OneTo, dims), b, a1, A...)

    _mapreducedim!(f, op, nothing, dims, (b, a1, A...))
end

function _mapreducedim!(@nospecialize(f), @nospecialize(op), @nospecialize(initop),
        dims::Dims, arrays::Tuple{Vararg{AbstractStridedView}})
    any(isequal(0), dims) && return arrays[1] # don't do anything

    _mapreduce_fuse!(f, op, initop, dims, promoteshape(dims, arrays...))

    return arrays[1]
end

function _mapreduce_fuse!(@nospecialize(f), @nospecialize(op), @nospecialize(initop),
        dims::Dims, arrays::Tuple{Vararg{AbstractStridedView}})
    # Fuse dimensions if possible: assume that at least one array, e.g. the output
    # array in arrays[1], has its strides sorted
    allstrides = map(strides, arrays)
    @inbounds for i = length(dims):-1:2
        merge = true
        for s in allstrides
            if s[i] != dims[i-1]*s[i-1]
                merge = false
                break
            end
        end
        if merge
            dims = TupleTools.setindex(dims, dims[i-1]*dims[i], i-1)
            dims = TupleTools.setindex(dims, 1, i)
        end
    end
    _mapreduce_order!(f, op, initop, dims, allstrides, arrays)
end

function _mapreduce_order!(@nospecialize(f), @nospecialize(op), @nospecialize(initop),
        dims, strides, arrays)
    M = length(arrays)
    N = length(dims)
    # sort order of loops/dimensions by modelling the importance of each dimension
    g = 8*sizeof(Int) - leading_zeros(M+1) # ceil(Int, log2(M+2)) # to account for the fact
    # that there are M arrays, where the first one is counted with a factor 2
    importance = 2 .* ( 1 .<< (g.*(N .- indexorder(strides[1]))))  # first array is output
    # and is more important by a factor 2
    for k = 2:M
        importance = importance .+ ( 1 .<< (g.*(N .- indexorder(strides[k]))))
    end

    importance = importance .* (dims .> 1) # put dims 1 at the back
    p = TupleTools.sortperm(importance, rev = true)
    dims = TupleTools.getindices(dims, p)
    strides = broadcast(TupleTools.getindices, strides, (p,))
    _mapreduce_block!(f, op, initop, dims, strides, arrays)
end

function _mapreduce_block!(@nospecialize(f), @nospecialize(op), @nospecialize(initop),
        dims, strides, arrays)
    offsets = map(offset, arrays)

    if all(l -> l<=BLOCKSIZE, broadcast(_length, (dims,), strides))
        _mapreduce_kernel!(f, op, initop, dims, dims, arrays, strides, offsets)
    else
        minstrides = map(min, strides...)
        mincosts = map(a->ifelse(iszero(a), 1, a << 1), minstrides)
        blocks = _computeblocks(dims, mincosts, strides)

        if Threads.nthreads() == 1 || !use_threads() || prod(dims) < BLOCKSIZE
            _mapreduce_kernel!(f, op, initop, dims, blocks, arrays, strides, offsets)
        elseif _length(dims, strides[1]) == 1 # complete reduction
            T = eltype(arrays[1])
            spacing = isbitstype(T) ? min(1, div(64, sizeof(T))) : 1# to avoid false sharing
            threadblocks, threadoffsets =
                _computethreadblocks(dims, mincosts, strides, offsets)
            threadedout = similar(arrays[1], spacing*length(threadblocks))
            a = arrays[1][ParentIndex(1)]
            if initop !== nothing
                a = initop(a)
            end
            _init_reduction!(threadedout, f, op, a)
            for i = 1:length(threadoffsets)
                threadoffsets[i] = (spacing*(i-1), Base.tail(threadoffsets[i])...)
            end
            _mapreduce_threaded!(threadblocks, threadoffsets, f, op, nothing, blocks,
                (threadedout, Base.tail(arrays)...), strides)
            for i = 1:length(threadblocks)
                a = op(a, threadedout[(i-1)*spacing+1])
            end
            arrays[1][ParentIndex(1)] = a
        else
            mincosts = mincosts .* .!(iszero.(strides[1]))
            # make cost of dimensions with zero stride in output array (reduction
            # dimensions), so that they are not divided in threading (which would lead to
            # race errors)

            threadblocks, threadoffsets =
                _computethreadblocks(dims, mincosts, strides, offsets)
            _mapreduce_threaded!(threadblocks, threadoffsets, f, op, initop, blocks,
                arrays, strides)
        end
    end
    return
end

_init_reduction!(out, f, op::Union{typeof(+),typeof(Base.add_sum)}, a) = fill!(out, zero(a))
_init_reduction!(out, f, op::Union{typeof(*),typeof(Base.mul_prod)}, a) = fill!(out, one(a))
_init_reduction!(out, f, op::typeof(min), a) = fill!(out, a)
_init_reduction!(out, f, op::typeof(max), a) = fill!(out, a)
_init_reduction!(out, f, op::typeof(&), a) = fill!(out, true)
_init_reduction!(out, f, op::typeof(|), a) = fill!(out, false)
_init_reduction!(out, f, op, a) = op(a,a) == a ? fill!(out, a) : error("unknown reduction; incompatible with multithreading")

@noinline function _mapreduce_threaded!(threadblocks, threadoffsets, f, op, initop, blocks, arrays, strides)
    @inbounds Threads.@threads for i = 1:length(threadblocks)
        _mapreduce_kernel!(f, op, initop, threadblocks[i], blocks, arrays, strides,
            threadoffsets[i])
    end
end

@generated function _mapreduce_kernel!(@nospecialize(f), @nospecialize(op),
        @nospecialize(initop), dims::NTuple{N,Int}, blocks::NTuple{N,Int},
        arrays::NTuple{M,AbstractStridedView}, strides::NTuple{M,NTuple{N,Int}},
        offsets::NTuple{M,Int}) where {N,M}

    # many variables
    blockloopvars = Array{Symbol}(undef, N)
    blockdimvars = Array{Symbol}(undef, N)
    innerloopvars = Array{Symbol}(undef, N)
    initblockdimvars = Array{Symbol}(undef, N)
    initvars = Array{Symbol}(undef, N+1)
    for i = 1:N
        blockloopvars[i] = Symbol(:J,i)
        blockdimvars[i] = Symbol(:d,i)
        innerloopvars[i] = Symbol(:j,i)
        initblockdimvars[i] = Symbol(:d′,i)
        initvars[i] = Symbol(:init,i)
    end
    initvars[N+1] = Symbol(:init,N+1)

    Ivars = Array{Symbol}(undef, M)
    Avars = Array{Symbol}(undef, M)
    stridevars = Array{Symbol}(undef, (N,M))
    for j = 1:M
        Ivars[j] = Symbol(:I,j)
        Avars[j] = Symbol(:A,j)
        for i = 1:N
            stridevars[i,j] = Symbol(:stride_,i,:_,j)
        end
    end

    # useful expressions
    pre1ex = Expr(:block)
    pre2ex = Expr(:block)
    pre3ex = Expr(:block)
    for j = 1:M
        push!(pre1ex.args, :($(Avars[j]) = arrays[$j]))
        push!(pre3ex.args, :($(Ivars[j]) = offsets[$j]+1))
        for i = 1:N
            push!(pre2ex.args, :($(stridevars[i,j]) = strides[$j][$i]))
        end
    end

    fcallex = Expr(:call, :f)
    for j = 2:M
        push!(fcallex.args, :($(Avars[j])[ParentIndex($(Ivars[j]))]))
    end
    lhsex = :(A1[ParentIndex($(Ivars[1]))])

    stepstride1ex = Vector{Expr}(undef, N)
    stepstride2ex = Vector{Expr}(undef, N)
    returnstride1ex = Vector{Expr}(undef, N)
    returnstride2ex = Vector{Expr}(undef, N)
    for i = 1:N
        stepstride1ex[i] = :($(Ivars[1]) += $(stridevars[i,1]))
        returnstride1ex[i] = :($(Ivars[1]) -=  $(blockdimvars[i]) * $(stridevars[i,1]))
        stepex = Expr(:block)
        returnex = Expr(:block)
        for j = 2:M
            push!(stepex.args, :($(Ivars[j]) += $(stridevars[i,j])))
            push!(returnex.args, :($(Ivars[j]) -=  $(blockdimvars[i]) * $(stridevars[i,j])))
        end
        stepstride2ex[i] = stepex
        returnstride2ex[i] = returnex
    end

    outerstepstrideex = Vector{Expr}(undef, N)
    outerreturnstrideex = Vector{Expr}(undef, N)
    for i = 1:N
        stepex = Expr(:block)
        returnex = Expr(:block)
        for j = 1:M
            push!(stepex.args, :($(Ivars[j]) +=  $(blockdimvars[i]) * $(stridevars[i,j])))
            push!(returnex.args, :($(Ivars[j]) -=  dims[$i] * $(stridevars[i,j])))
        end
        outerstepstrideex[i] = stepex
        outerreturnstrideex[i] = returnex
    end

    if op == Nothing
        ex = Expr(:(=), lhsex, fcallex)
        exa = Expr(:(=), :a, fcallex)
    else
        ex = Expr(:(=), lhsex, Expr(:call, :op, lhsex, fcallex))
        exa = Expr(:(=), :a, Expr(:call, :op, :a, fcallex))
    end
    i = 1
    if N >= 1
        ex = quote
            if $(stridevars[1,1]) == 0 # explicitly hoist A1[I1] out of loop
                a = $lhsex
                @simd for $(innerloopvars[i]) = Base.OneTo($(blockdimvars[i]))
                    $exa
                    $(stepstride2ex[i])
                end
                $lhsex = a
                $(returnstride2ex[i])
            else
                @simd for $(innerloopvars[i]) = Base.OneTo($(blockdimvars[i]))
                    $ex
                    $(stepstride1ex[i])
                    $(stepstride2ex[i])
                end
                $(returnstride1ex[i])
                $(returnstride2ex[i])
            end
        end
    end
    for outer i = 2:N
        ex = quote
            for $(innerloopvars[i]) = Base.OneTo($(blockdimvars[i]))
                $ex
                $(stepstride1ex[i])
                $(stepstride2ex[i])
            end
            $(returnstride1ex[i])
            $(returnstride2ex[i])
        end
    end

    if initop !== Nothing
        initex = Expr(:(=), lhsex, Expr(:call, :initop, lhsex))
        i = 1
        if N >= 1
            initex = quote
                $(initblockdimvars[i]) = ifelse($(stridevars[i,1]) == 0, 1, $(blockdimvars[i]))
                @simd for $(innerloopvars[i]) in Base.OneTo($(initblockdimvars[i]))
                    $initex
                    $(stepstride1ex[i])
                end
                $(Ivars[1]) -=  $(initblockdimvars[i]) * $(stridevars[i,1])
            end
        end
        for outer i = 2:N
            initex = quote
                $(initblockdimvars[i]) = ifelse($(stridevars[i,1]) == 0, 1, $(blockdimvars[i]))
                for $(innerloopvars[i]) in Base.OneTo($(initblockdimvars[i]))
                    $initex
                    $(stepstride1ex[i])
                end
                $(Ivars[1]) -=  $(initblockdimvars[i]) * $(stridevars[i,1])
            end
        end
        ex = quote
            if $(initvars[1])
                $initex
            end
            $ex
        end
    end

    if initop === Nothing
        for outer i = 1:N
            ex = quote
                for $(blockloopvars[i]) = 1:blocks[$i]:dims[$i]
                    $(blockdimvars[i]) = min(blocks[$i], dims[$i]-$(blockloopvars[i])+1)
                    $ex
                    $(outerstepstrideex[i])
                end
                $(outerreturnstrideex[i])
            end
        end
        ex = quote
            $pre1ex
            $pre2ex
            $pre3ex
            @inbounds $ex
            return A1
        end
    else
        for outer i = 1:N
            ex = quote
                $(initvars[i]) = $(initvars[i+1])
                for $(blockloopvars[i]) = 1:blocks[$i]:dims[$i]
                    $(blockdimvars[i]) = min(blocks[$i], dims[$i]-$(blockloopvars[i])+1)
                    $ex
                    $(initvars[i]) &= $(stridevars[i,1]) > 0
                    $(outerstepstrideex[i])
                end
                $(outerreturnstrideex[i])
            end
        end
        ex = quote
            $pre1ex
            $pre2ex
            $pre3ex
            $(initvars[N+1]) = true
            @inbounds $ex
            return A1
        end
    end
    return ex
end

function indexorder(strides::NTuple{N,Int}) where {N}
    # returns order such that strides[i] is the order[i]th smallest element of strides, not counting zero strides
    # zero strides have order N
    return ntuple(Val(N)) do i
        si = strides[i]
        si == 0 && return 1
        k = 1
        for s in strides
            if s != 0 && s < si
                k += 1
            end
        end
        return k
    end
end

_length(dims::Tuple, strides::Tuple) = ifelse(iszero(strides[1]), 1, dims[1]) * _length(Base.tail(dims), Base.tail(strides))
_length(dims::Tuple{}, strides::Tuple{}) = 1
_maxlength(dims::Tuple, strides::Tuple{Vararg{Tuple}}) = maximum(broadcast(_length, (dims,), strides))

function _lastargmax(t::Tuple)
    i = 1
    for j = 2:length(t)
        @inbounds if t[j] >= t[i]
            i = j
        end
    end
    return i
end

_computeblocks(dims::Tuple{}, costs::Tuple{}, strides::Tuple{Vararg{Tuple{}}}, blocksize::Int = BLOCKSIZE) = ()
function _computeblocks(dims::NTuple{N,Int}, costs::NTuple{N,Int}, strides::Tuple{Vararg{NTuple{N,Int}}}, blocksize::Int = BLOCKSIZE) where {N}
    if _maxlength(dims, strides) <= blocksize
        return dims
    elseif all(isequal(1), first.(indexorder.(strides)))
        return (dims[1], _computeblocks(tail(dims), tail(costs), map(tail, strides), div(blocksize, dims[1]))...)
    elseif blocksize == 0
        return ntuple(n->1, StaticLength(N))
    else
        blocks = dims
        while _maxlength(blocks, strides) >= 2*blocksize
            i = _lastargmax((blocks .- 1) .* costs)
            blocks = TupleTools.setindex(blocks, (blocks[i]+1)>>1, i)
        end
        while _maxlength(blocks, strides) > blocksize
            i = _lastargmax((blocks .- 1) .* costs)
            blocks = TupleTools.setindex(blocks, blocks[i]-1, i)
        end
        return blocks
    end
end

@inbounds function _computethreadblocks(dims::NTuple{N,Int}, costs::NTuple{N,Int}, strides::NTuple{M,NTuple{N,Int}}, offsets::NTuple{M,Int}) where {N,M}
    threadblocks = [dims]
    threadoffsets = [offsets]
    for k in factors
        l = length(threadblocks)
        for j = 1:l
            dims = popfirst!(threadblocks)
            offsets = popfirst!(threadoffsets)
            i = _lastargmax((dims .- (k-1)) .* costs)
            if costs[i] == 0
                push!(threadblocks, dims)
                push!(threadoffsets, offsets)
                return threadblocks, threadoffsets
            end
            ndi = div(dims[i], k)
            newdims = setindex(dims, ndi, i)
            stridesi = getindex.(strides, i)
            for m = 1:k-1
                push!(threadblocks, newdims)
                push!(threadoffsets, offsets)
                offsets = offsets .+ ndi .* stridesi
            end
            ndi = dims[i]-(k-1)*ndi
            newdims = setindex(dims, ndi, i)
            push!(threadblocks, newdims)
            push!(threadoffsets, offsets)
        end
    end
    return threadblocks, threadoffsets
end


text:
@testset "129.sum-root-to-leaf-numbers.jl" begin
    @test sum_numbers(TreeNode{Int}([1, 2 ,3])) == 25
    @test sum_numbers(TreeNode{Int}([4, 9, 0, 5, 1])) == 1026
end


text:
module HodgkinHuxleySingle

using Test
using Conductor, ModelingToolkit, OrdinaryDiffEq, Unitful, IfElse
import ModelingToolkit: isparameter
import Conductor: Na, K
using Unitful: mV, mS, cm, µm, µA, ms, pA

@testset "Hodgkin Huxley Single Compartment" begin

# Symbolic model construction

Vₘ = MembranePotential(-65mV)

nav_kinetics = [
    Gate(AlphaBeta,
         IfElse.ifelse(Vₘ == -40.0, 1.0, (0.1*(Vₘ + 40.0))/(1.0 - exp(-(Vₘ + 40.0)/10.0))),
         4.0*exp(-(Vₘ + 65.0)/18.0), 3, name = :m)
    Gate(AlphaBeta,
         0.07*exp(-(Vₘ+65.0)/20.0),
         1.0/(1.0 + exp(-(Vₘ + 35.0)/10.0)), name = :h)]

kdr_kinetics = [
    Gate(AlphaBeta,
         IfElse.ifelse(Vₘ == -55.0, 0.1, (0.01*(Vₘ + 55.0))/(1.0 - exp(-(Vₘ + 55.0)/10.0))),
         0.125 * exp(-(Vₘ + 65.0)/80.0),
         4, name = :n)]

alphamss = ((4.0exp(-3.6111111111111107 - (0.05555555555555555Vₘ)) +
            IfElse.ifelse(Vₘ == -40.0, 1.0, (4.0 + 0.1Vₘ)*((1.0 - exp(-4.0 - (0.1Vₘ)))^-1)))^-1) *
            IfElse.ifelse(Vₘ == -40.0, 1.0, (4.0 + 0.1Vₘ)*((1.0 - exp(-4.0 - (0.1Vₘ)))^-1))

@test isequal(steadystate(nav_kinetics[1]), alphamss)

@named NaV = IonChannel(Sodium, nav_kinetics, max_g = 120mS/cm^2)
@named Kdr = IonChannel(Potassium, kdr_kinetics, max_g = 36mS/cm^2)
@named leak = IonChannel(Leak, max_g = 0.3mS/cm^2)

@test isparameter(leak.gbar) && !isparameter(get_output(NaV))
channels = [NaV, Kdr, leak]
@test [length(equations(x)) for x in channels] == [3,2,1]
reversals = Equilibria([Na => 50.0mV, K => -77.0mV, Leak => -54.4mV])
@named Iₑ = IonCurrent(NonIonic)
electrode_pulse = Iₑ ~ IfElse.ifelse(t > 100.0, IfElse.ifelse(t < 200.0, ustrip(Float64, µA, 400pA), 0.0), 0.0)

@named neuron = CompartmentSystem(Vₘ, channels, reversals;
                                  geometry = Sphere(radius = 20µm),
                                  stimuli = [electrode_pulse])

@test length.([equations(neuron),
               states(neuron),
               parameters(neuron)]) == [11,13,8]

time = 300.
sim_sys = Simulation(neuron, time = time*ms, return_system = true)

@test length.([equations(sim_sys),
               states(sim_sys),
               parameters(sim_sys)]) == [4,4,8]

expect_syms = [:Vₘ, :NaV₊m, :NaV₊h, :Kdr₊n]
@test all(x -> hasproperty(sim_sys, x), expect_syms)

### Hand-written model setup ###
# Steady-state functions
minf(V0) = inv((4.0 + 0.1V0)*inv(1.0 - exp(-4.0 - 0.1V0)) + 4.0*exp(-3.6111111111111107 - 0.05555555555555555V0)) *
           (4.0 + 0.1V0) * inv(1.0 - exp(-4.0 - 0.1V0))
ninf(V0) = (0.55 + 0.01V0) * inv(1.0 - exp(-5.5 - 0.1V0)) * inv(0.125exp(-0.8125 - 0.0125V0) +
                                                                (0.55 + 0.01V0) * inv(1.0 - exp(-5.5 - 0.1V0)))
hinf(V0) = 0.07*exp(-3.25 - 0.05V0) * inv(0.07*exp(-3.25 - 0.05V0) + inv(1.0 + exp(-3.5 - 0.1V0)) )
hand_pulse(t) = ifelse(100. < t < 200. , 0.0004 , 0.0)

# Initial conditions
V0 = -65.0
u0 = zeros(Float64, 4)
p = zeros(Float64, 8)

# Serialized MTK function initial conditions
u0[1] = V0       # Vₘ
u0[2] = minf(V0) # NaV₊m
u0[3] = hinf(V0) # NaV₊h
u0[4] = ninf(V0) # Kdr₊n

# Parameters
p[1]  = 0.001      # cₘ
p[2]  = 5.02655e-5 # aₘ
p[3]  = 50.0       # ENa
p[4]  = -77.0      # EK
p[5]  = -54.4      # El
p[6]  = 120.0      # NaV₊gbar
p[7]  = 36.0       # Kdr₊gbar
p[8]  = 0.3        # leak₊g

# Cleaned up version of what default MTK output looks like
function hodgkin_huxley!(du, u, p, t)
    @inbounds let Vₘ         = u[1],
                  NaV₊m      = u[2],
                  NaV₊h      = u[3],
                  Kdr₊n      = u[4],

                  cₘ          = p[1],
                  aₘ          = p[2],
                  ENa         = p[3],
                  EK          = p[4],
                  El          = p[5],
                  NaV₊gbar    = p[6],
                  Kdr₊gbar    = p[7],
                  leak₊g      = p[8],
                  Iapp        = hand_pulse(t)
    du[1] = (Iapp - aₘ*leak₊g*(Vₘ - El) - Kdr₊gbar*aₘ*Kdr₊n^4*(Vₘ - EK) -
             NaV₊gbar*aₘ*NaV₊h*NaV₊m^3*(Vₘ - ENa))*inv(aₘ)*inv(cₘ)
    du[2] = ifelse(Vₘ == -40.0 , 1.0 , (4.0 + 0.1Vₘ) * inv(1.0 - exp(-4.0 - 0.1Vₘ))) *
            (1 - NaV₊m) - 4.0*NaV₊m*exp(-3.6111111111111107 + -0.05555555555555555Vₘ)

    du[3] = 0.07*exp(-3.25 - 0.05Vₘ)*(1 - NaV₊h) - NaV₊h*inv(1.0 + exp(-3.5 - 0.1Vₘ))
    du[4] = ifelse(Vₘ == -55.0 , 0.1 , (0.55 + 0.01Vₘ) * inv(1.0 - exp(-5.5 - 0.1Vₘ)))*(1 - Kdr₊n) -
            0.125 * Kdr₊n * exp(-0.8125 - 0.0125Vₘ)
    end
    return nothing
end

byhand_prob = ODEProblem{true}(hodgkin_huxley!, u0, (0.,300.), p)
mtk_prob = ODAEProblem(sim_sys, [], (0., time), [])

byhand_sol = solve(byhand_prob, Rosenbrock23(), reltol=1e-9, abstol=1e-9, saveat=0.025);
current_mtk_sol = solve(mtk_prob, Rosenbrock23(), reltol=1e-9, abstol=1e-9, saveat=0.025);

tsteps = 0.0:0.025:300.0
byhand_out = Array(byhand_sol(tsteps, idxs=1))
current_mtk_out = current_mtk_sol(tsteps)[Vₘ]

@test isapprox(byhand_out, current_mtk_out, rtol=0.0001)

end #testset

end # module


text:
# Linear Regression applied to Housing Prices Data

push!(LOAD_PATH, ".")
push!(LOAD_PATH, "../common/")

using Gadfly
importall linear_regression
importall optimization

# parameters
alpha = 0.2
max_iter = 100
min_err = 0.0

println("Reading data...")
data = float(open(readdlm, "housing.data"))

data = data'
u = mean(data, 2)
r = maximum(data, 2) - minimum(data, 2)
data = (data .- u) ./ r
n, m = size(data)
data = [ones(1,m); data]
n = n + 1

# shuffle
println("Shuffling data...")
data = data[:, shuffle(collect(1:m))]

# seperate train and test data
m_train = 406
m_test = m - m_train
X_train = data[1:n, 1:m_train]
y_train = data[end,1:m_train]
X_test = data[1:n, m_train+1:end]
y_test = data[end, m_train+1:end]

# theta initialization
theta = init_theta(size(X_train, 1))

# training
println("Training data...")
tic()
theta, history, iter = gradient_descent(h, J, grad, X_train, y_train, theta,
                                        alpha, max_iter, min_err, false)
toc()

ms = collect(1:m_test)
y_pred = theta' * X_test
sorted = sortperm(y_test[:])

println("finished")
println("Epochs: ", iter)
println("Cost (trn): ", J(theta, X_train, y_train))
println("Cost (tst): ", J(theta, X_test, y_test))

# plotting
plot(
    layer(x=ms, y=r[end]*y_test[sorted]+u[end], Geom.point, Theme(default_color=colorant"green")),
    layer(x=ms, y=r[end]*y_pred[sorted]+u[end], Geom.point, Theme(default_color=colorant"red")),
    Guide.XLabel("House"),
    Guide.YLabel("Price"),
    Guide.Title("Housing Prices"),
    Guide.manual_color_key("Prices",["Actual", "Predicted"], ["green", "red"])
)

max_stars_repo_name: ilkerkesen/ufldl-tutorial

max_issues_repo_name: ilkerkesen/ufldl-tutorial

max_forks_repo_name: ilkerkesen/ufldl-tutorial


text:
module ActuarialScience

export maleMort, femaleMort,
        MortalityTable,
        qx,px, tpx, tqx,
        tqxy,tpxy, tqx̅y̅, tpx̅y̅,
        ω,w, lx, dx, ex,
        ixVector,
        InterestRate,
        i,vx,tvx,v,
        LifeInsurance,
        Ax, Axn, äx,äxn,
        Dx, Mx

include("interest.jl")
include("decrement.jl")
include("mortality.jl")

# 'actuarial objects' that combine multiple forms of decrements (lapse, interest, death, etc)
abstract type AbstractActuarial end

mutable struct ActuarialHelper
    mort::MortalityTable
    int::InterestRate
end

mutable struct LifeInsurance <: AbstractActuarial
    ah::ActuarialHelper

    function LifeInsurance(mt::MortalityTable,int::InterestRate)
        return new(ActuarialHelper(mt,int))
    end
end

###################
## COMMUTATIONS ###
###################

Dx(ah::ActuarialHelper,x) = tvx(ah.int,x,1) * lx(ah.mort,x)

preMx(ah::ActuarialHelper,x) = tvx(ah.int,x+1,1) * dx(ah.mort,x)

function Nx(ah::ActuarialHelper,x)
range =x:min(ω(ah.mort),ω(ah.int))
    return sum(map((y) -> Dx(ah,y),range))
end

function Mx(ah::ActuarialHelper,x)
    range =  x:(min(ω(ah.mort),ω(ah.int))-1)
    return sum(map((x -> preMx(ah,x)),range))
end
tEx(ah::ActuarialHelper,t,x) = Dx(x+t) / Dx(x)

##################
### Insurances ###
##################

#term insurance on age x for n years
Axn(ins::LifeInsurance,x,n) = (Mx(ins.ah,x) - Mx(ins.ah,x + n) ) / Dx(ins.ah,x)

# whole life insurance
Ax(ins::LifeInsurance,x) = Mx(ins.ah,x) / Dx(ins.ah,x)

# life annuity due
äx(ins::LifeInsurance,x) = Nx(ins.ah,x) / Dx(ins.ah,x)

# finite duration life annuity due
äxn(ins::LifeInsurance,x,n) = (Nx(ins.ah,x) - Nx(ins.ah,x+n) )/ Dx(ins.ah,x)

Dx(ins::LifeInsurance,x) = Dx(ins.ah,x)
Mx(ins::LifeInsurance,x) = Mx(ins.ah,x)

end # module


text:
using Test

function parseInput(input)
    chars = split(input, ",")
    chars = parse.(Int, chars)
    return chars
end

function getVal(val, codes, mode)
    if mode == 0
        return codes[val + 1]
    elseif mode == 1
        return val
    end
end

function part1(codes)
    i = 1
    outputs = []
    while codes[i] != 99
        operation = codes[i]
        numbers = digits(operation, pad=5)
        instruction = numbers[1] + 10*numbers[2]
        mode1 = numbers[3]
        mode2 = numbers[4]
        mode3 = numbers[5]
        if instruction == 1
            val1, val2, dest = codes[i+1:i+3]
            codes[dest+1] = getVal(val1, codes, mode1) + getVal(val2, codes, mode2) #codes[val2]
            i += 4
        elseif instruction == 2
            val1, val2, dest = codes[i+1:i+3]
            codes[dest+1] = getVal(val1, codes, mode1) * getVal(val2, codes, mode2) #codes[val2]
            i += 4
        elseif instruction == 3
            dest = codes[i+1] + 1
            codes[dest] = 1
            i += 2
        elseif instruction == 4
            dest = codes[i+1] + 1
            append!(outputs, codes[dest])
            i += 2
        end
    end
    return codes, outputs
end

function run(s)
    input = parseInput(split(s, "\n")[1])
    # input[2] = 12
    # input[3] = 2
    return part1(input)[2][end]
end

@test part1(parseInput("1002,4,3,4,33"))[1] == [1002, 4, 3, 4, 99]

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time*1000)")
    println(res)
end

main()


text:
export RadialTrajectory, radialNodes, radialDensity

"""
    RadialTrajectory(numProfiles, numSamplingPerProfile
                  ; TE::Float64=0.0
                  , AQ::Float64=1.e-3
                  , angleOffset= :equispaced
                  , kargs...)

returns a 2d radial trajectory.

# Arguments
* `numProfiles::Int64`            - number of profiles
* `numSamplingPerProfile::Int64`  - number of sampling points per profile
* (`TE::Float64=0.0`)             - echo time in s
* (`AQ::Float64=1.e-3`)           - readout duration in s (per profile)
* (`angleOffset= :equispaced`)    - spacing of profile angles (`:equispaced` sampling, `:golden` angle sampling or `:random` sampling)
"""
function RadialTrajectory(::Type{T}, numProfiles, numSamplingPerProfile
                  ; TE::T=0.0
                  , AQ::T=1.e-3
                  , angleOffset= :equispaced
                  , kargs...) where T
  nodes = radialNodes(T, numProfiles, numSamplingPerProfile; angleOffset=angleOffset)
  times = readoutTimes(T, numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ)
  return  Trajectory("Radial", nodes, times, TE, AQ, numProfiles, numSamplingPerProfile, 1, false, true)
end

function radialNodes(::Type{T}, numProfiles::Int64
                          , numSamplingPerProfile::Int64
                          ; angleOffset= :equispaced
                          , kargs...) where T
  nodes = zeros(T, 2, numSamplingPerProfile, numProfiles)
  if angleOffset == :golden
    angles = [i*getGoldenAngleRad()  for i=0:numProfiles-1 ]
  elseif angleOffset == :random
    angles = sort(pi .* rand(numProfiles))
  elseif angleOffset == :equispaced
    angles = collect(pi.*(0:numProfiles-1)/numProfiles)
  end

  pos = collect((0:numSamplingPerProfile-1)/numSamplingPerProfile .- 0.5)
  for l = 1:numProfiles
    for k = 1:numSamplingPerProfile
      nodes[1,k,l] = (-1)^l * pos[k]*cos(angles[l])
      nodes[2,k,l] = (-1)^l * pos[k]*sin(angles[l])
    end
  end
  return reshape(nodes, 2, numSamplingPerProfile*numProfiles)
end

function radialDensity(::Type{T}, numProfiles::Int64, numSamplingPerProfile::Int64) where T
  density = zeros(T, numSamplingPerProfile, numProfiles)
  pos = collect((0:numSamplingPerProfile-1)/numSamplingPerProfile .- 0.5)
  for l = 1:numProfiles
    for k = 1:numSamplingPerProfile
      density[k,l] = abs(pos[k])
    end
  end
  return vec(density)
end


text:
# using Pkg; for p in ("Knet","ArgParse","Images"); haskey(Pkg.installed(),p) || Pkg.add(p); end

"""
Train a Variational Autoencoder with convolutional layers
on the MNIST dataset.
"""
module VAE
using Knet, CUDA, ArgParse, Images, Random, Statistics
import AutoGrad: getval

include(Knet.dir("data","mnist.jl"))
include(Knet.dir("data","imagenet.jl"))

const F = Float32

# global variables reset in main
Atype = gpu() >= 0 ? KnetArray{F} : Array{F}
BINARIZE = false

binarize(x) = convert(Atype, rand(F, size(x))) .< x

function encode(ϕ, x)
    x = reshape(x, (28,28,1,:))

    x = conv4(ϕ[1], x, padding=1)
    x = relu.(x .+ ϕ[2])

    x = conv4(ϕ[3], x, padding=1, stride=2)
    x = relu.(x .+ ϕ[4])
    x = conv4(ϕ[5], x, padding=1)
    x = relu.(x .+ ϕ[6])

    x = conv4(ϕ[7], x, padding=1, stride=2)
    x = relu.(x .+ ϕ[8])

    x = mat(x)
    x = relu.(ϕ[9]*x .+ ϕ[10])

    μ = ϕ[end-3]*x .+ ϕ[end-2]
    logσ² = ϕ[end-1]*x .+ ϕ[end]

    return μ, logσ²
end

function decode(θ, z)
    z = relu.(θ[1]*z .+ θ[2])
    z = relu.(θ[3]*z .+ θ[4])

    filters = size(θ[5], 4)
    width = Int(sqrt(size(z,1) ÷ filters))
    z = reshape(z, (width, width, filters, :))

    z = deconv4(θ[5], z, padding=1, stride=2)
    z = relu.(z .+ θ[6])

    z = deconv4(θ[7], z, padding=1)
    z = relu.(z .+ θ[8])

    z = deconv4(θ[9], z, padding=1, stride=2)
    z = relu.(z .+ θ[10])

    z = deconv4(θ[11], z, padding=1)
    z = sigm.(z .+ θ[12])

    return z
end

function binary_cross_entropy(x, x̂)
    x = reshape(x, size(x̂))
    s = @. x * log(x̂ + F(1e-10)) + (1-x) * log(1 - x̂ + F(1e-10))
    return -sum(s) / length(x)
end

function loss(w, x; samples=1)
    θ, ϕ = w[:decoder], w[:encoder]
    μ, logσ² = encode(ϕ, x)
    σ² = exp.(logσ²)
    σ = sqrt.(σ²)

    KL =  - sum(@. 1 + logσ² - μ*μ - σ²) / 2
    # Normalise by same number of elements as in reconstruction
    KL /= length(x)

    BCE = F(0)
    for s=1:samples
        # ϵ = randn!(similar(μ))
        ϵ = convert(Atype, randn(F, size(μ)))
        z = @. μ + ϵ * σ
        x̂ = decode(θ, z)
        BCE += binary_cross_entropy(x, x̂)
    end
    BCE /= samples

    return BCE + KL
end

lossgradient = grad(loss)

function aveloss(w, xtrn; samples=1, batchsize=100)
    θ, ϕ = w[:decoder], w[:encoder]
    ls = F(0)
    count = 0
    for x in minibatch(xtrn, batchsize; xtype=Atype)
        BINARIZE && (x = binarize(x))
        ls += loss(w, x; samples=samples)
        count += 1
    end
    N = length(xtrn) ÷ size(xtrn, ndims(xtrn))
    return (ls / count) * N
end

function weights(nz, nh)
    θ = [] # z->x

    push!(θ, xavier(nh, nz))
    push!(θ, zeros(nh))

    push!(θ, xavier(64*7^2, nh))
    push!(θ, zeros(64*7^2))

    push!(θ, xavier(4, 4, 32, 64))
    push!(θ, zeros(1,1,32,1))

    push!(θ, xavier(3, 3, 32, 32))
    push!(θ, zeros(1,1,32,1))

    push!(θ, xavier(4, 4, 16, 32))
    push!(θ, zeros(1, 1, 16, 1))

    push!(θ, xavier(3, 3, 1, 16))
    push!(θ, zeros(1,1,1,1))

    θ = map(a->convert(Atype,a), θ)

    ϕ = [] # x->z

    push!(ϕ, xavier(3, 3, 1, 16))
    push!(ϕ, zeros(1, 1, 16, 1))

    push!(ϕ, xavier(4, 4,  16, 32))
    push!(ϕ, zeros(1, 1, 32, 1))

    push!(ϕ, xavier(3, 3, 32, 32))
    push!(ϕ, zeros(1,1, 32, 1))

    push!(ϕ, xavier(4, 4, 32, 64))
    push!(ϕ, zeros(1, 1, 64, 1))

    push!(ϕ, xavier(nh, 64*7^2))
    push!(ϕ, zeros(nh))

    push!(ϕ, xavier(nz, nh)) # μ
    push!(ϕ, zeros(nz))
    push!(ϕ, xavier(nz, nh)) # logσ^2
    push!(ϕ, zeros(nz))

    ϕ = map(a->convert(Atype,a), ϕ)

    return θ, ϕ
end

function plot_dream(w; gridsize=(5,5), scale=1.0, outfile=nothing)
    θ = w[:decoder]
    nh, nz = size(θ[1])
    atype = θ[1] isa KnetArray ? KnetArray : Array
    m, n = gridsize
    nimg = m*n
    z = convert(atype, randn(F, nz, nimg))
    x̂ = Array(decode(θ, z))
    images = map(i->reshape(x̂[:,:,:,i], (28,28,1)), 1:nimg)
    grid = make_image_grid(images; gridsize=gridsize, scale=scale)

    if outfile == nothing
        display(colorview(Gray, grid))
    else
        save(outfile, colorview(Gray, grid))
    end
end

function main(args="")
    s = ArgParseSettings()
    s.description="Variational Auto Encoder on MNIST dataset."
    s.exc_handler=ArgParse.debug_handler
    @add_arg_table s begin
        ("--seed"; arg_type=Int; default=-1; help="random number seed: use a nonnegative int for repeatable results")
        ("--batchsize"; arg_type=Int; default=100; help="minibatch size")
        ("--epochs"; arg_type=Int; default=100; help="number of epochs for training")
        ("--nh"; arg_type=Int; default=400; help="hidden layer dimension")
        ("--nz"; arg_type=Int; default=40; help="encoding dimention")
        ("--atype"; default=(gpu()>=0 ? "KnetArray{F}" : "Array{F}"); help="array type: Array for cpu, KnetArray for gpu")
        ("--infotime"; arg_type=Int; default=1; help="report every infotime epochs")
        ("--binarize"; arg_type=Bool; default=false; help="dinamically binarize during training")
        ("--optim"; default="Adam()"; help="optimizer")
    end
    isa(args, String) && (args=split(args))
    if in("--help", args) || in("-h", args)
        ArgParse.show_help(s; exit_when_done=false)
        return
    end
    o = parse_args(args, s; as_symbols=true)

    global Atype = eval(Meta.parse(o[:atype]))
    global BINARIZE = o[:binarize]
    @info("using $Atype")
    # gc(); knetgc();
    o[:seed] > 0 && Knet.seed!(o[:seed])

    xtrn, ytrn, xtst, ytst = mnist()
    wdec, wenc = weights(o[:nz], o[:nh])
    w = Dict(
        :encoder => wenc,
        :decoder => wdec)
    opt = Dict(
        :encoder => map(wi->eval(Meta.parse(o[:optim])), w[:encoder]),
        :decoder => map(wi->eval(Meta.parse(o[:optim])), w[:decoder]),
    )

    report(epoch) = begin
            println((:epoch, epoch,
                     :trn, aveloss(w, xtrn; batchsize=o[:batchsize]),
                     :tst, aveloss(w, xtst; batchsize=o[:batchsize])))
    end

    report(0)
    @time for epoch=1:o[:epochs]
        @time for x  in minibatch(xtrn, o[:batchsize]; xtype=Atype, shuffle=true)
            BINARIZE && (x = binarize(x))
            dw = lossgradient(w, x)
            update!(w, dw, opt)
        end
        (epoch % o[:infotime] == 0) && report(epoch)
    end

    return w
end

PROGRAM_FILE == "vae_conv_mnist.jl" && main(ARGS)

end # module


text:
#################################################################
#
#    Definition of MCModel types
#
#################################################################

function ispartition(m::Dict, n::Int)
	c = zeros(n)
	for v in values(m)
		c[v[1]:(v[1]+prod(v[2])-1)] = c[v[1]:(v[1]+prod(v[2])-1)] .+ 1
	end
	all(c .== 1.)
end

#### misc functions common to all models  ####
hasgradient{M<:MCModel}(m::M) = m.evalg != nothing
hastensor{M<:MCModel}(m::M) = m.evalt != nothing
hasdtensor{M<:MCModel}(m::M) = m.evaldt != nothing

#### User-facing model creation function  ####

# Currently only a "likelihood" model type makes sense
# Left as is in case other kind of models come up
function model(f::Union(Function, Distribution, Expr); mtype="likelihood", args...)
	if mtype == "likelihood"
		return MCLikelihood(f; args...)
	elseif mtype == "whatever"
	else
	end
end


text:
abstract type SingularityExtractionRule end
regularpart_quadrule(qr::SingularityExtractionRule) = qr.regularpart_quadrule

function momintegrals!(op, g, f, t, s, z, strat::SingularityExtractionRule)

    womps = strat.outer_quad_points

    sop = singularpart(op)
    rop = regularpart(op)

    # compute the regular part
    rstrat = regularpart_quadrule(strat)
    momintegrals!(rop, g, f, t, s, z, rstrat)

    for p in 1 : length(womps)
        x = womps[p].point
        dx = womps[p].weight

        innerintegrals!(sop, x, g, f, t, s, z, strat, dx)
    end # next quadrature point

end


text:
function build_zeroprofit!(m, jm)

    # Add zero profit constraints
    for s in m._productions
        exa = :(
            JuMP.@NLexpression(
                $jm,
                +(
                    $(
                        (:(
                            $(get_jump_variable_for_commodity(jm, input.commodity)) * $(jm[get_comp_demand_name(input)])
                        ) for input in s.inputs)...
                )
                ) -
                +(
                    $(
                        (:(
                            $(get_jump_variable_for_commodity(jm, output.commodity)) * $(jm[get_comp_supply_name(output)])
                        ) for output in s.outputs)...
                    )
                )
            )
        )

        exb = eval(swap_our_param_with_jump_param(jm, exa))

        Complementarity.add_complementarity(jm, get_jump_variable_for_sector(jm, s.sector), exb, string("F_", get_name(s.sector), s.sector.subindex!==nothing ? s.sector.subindex : ""))
        push!(m._nlexpressions, exb)
    end
end