text:
cd(ENV["PWD"])

import SeaMass
using Base.Test
using Plots
pyplot()

function test(filename, id, inputSpectrumID, outputSpectrumID)

  # load SMB input spectrum
  smbSpectrumIn = SeaMass.SmbSpectrum(
    "data/out/" * filename * "/1.seamass/" * filename * "." * id * ".smb",
    outputSpectrumID
  )

  # load SMB blurred input spectrum
  smbSpectrumInBlur = false
  try
    smbSpectrumInBlur = SeaMass.SmbSpectrum(
      "data/out/" * filename * "/1.seamass/" * filename * "." * id * ".input.smb",
      outputSpectrumID
    )
  end

  # load SMB reconstructed output spectrum
  smbSpectrumOutReconstruct = SeaMass.SmbSpectrum(
    "data/out/" * filename * "/4.seamass-restore_--reconstruct/" * filename * "." * id * ".smb",
    outputSpectrumID
  )

  # Compare smb
  plot(
    smbSpectrumIn.locations,
    vcat(
      smbSpectrumIn.counts ./ (smbSpectrumIn.locations[2:end] - smbSpectrumIn.locations[1:end-1]),
      0.0,
    ),
    line = :steppost,
    label = "SMB input",
    title = "SMB binned comparison - " * filename,
    xlabel = "m/z (Th)",
    ylabel = "ion count density",
    reuse = false,
  )
  if smbSpectrumInBlur != false
    plot!(
      smbSpectrumInBlur.locations,
      vcat(
        smbSpectrumInBlur.counts ./  (smbSpectrumInBlur.locations[2:end] - smbSpectrumInBlur.locations[1:end-1]),
        0.0,
      ),
      line = :steppost,
      label = "SMB blurred input",
    )
  end
  plot!(
    smbSpectrumOutReconstruct.locations,
    vcat(
      smbSpectrumOutReconstruct.counts ./  (smbSpectrumOutReconstruct.locations[2:end] - smbSpectrumOutReconstruct.locations[1:end-1]),
      0.0,
    ),
    line = :steppost,
    label = "SMB output (seaMass-restore --reconstruct)",
  )
  gui()

  # load SMB output spectrum
  smbSpectrumOut = SeaMass.SmbSpectrum(
    "data/out/" * filename * "/2.seamass-restore/" * filename * "." * id * ".smb",
    outputSpectrumID
  )

  # load SMB deconvolved output spectrum
  smbSpectrumOutDeconvolve = SeaMass.SmbSpectrum(
    "data/out/" * filename * "/3.seamass-restore_--deconvolve/" * filename * "." * id * ".smb",
    outputSpectrumID
  )

  # load SMB centroided output spectrum
  smbSpectrumOutCentroid = SeaMass.SmbSpectrum(
    "data/out/" * filename * "/5.seamass-restore_--centroid/" * filename * "." * id * ".smb",
    outputSpectrumID
  )

  plot(
    smbSpectrumOutReconstruct.locations,
    vcat(
      smbSpectrumOutReconstruct.counts ./  (smbSpectrumOutReconstruct.locations[2:end] - smbSpectrumOutReconstruct.locations[1:end-1]),
      0.0,
    ),
    line = :steppost,
    label = "SMB output (seaMass-restore --reconstruct)",
    title = "SMB proccessed comparison - " * filename,
    xlabel = "m/z (Th)",
    ylabel = "ion count density",
    reuse = false,
  )
  plot!(
    smbSpectrumOutDeconvolve.locations,
    smbSpectrumOutDeconvolve.counts,
    label = "SMB output (seaMass-restore --deconvolve)",
  )
  plot!(
    smbSpectrumOut.locations,
    smbSpectrumOut.counts,
    label = "SMB output (seaMass-restore)",
  )
  sticks!(
    smbSpectrumOutCentroid.locations,
    smbSpectrumOutCentroid.counts,
    label = "SMB output (seaMass-restore --centroid)",
    m = 4,
  )
  gui()

  # load mzMLb reconstructed output spectrum
  mzmlbSpectrumOutReconstruct = SeaMass.MzmlbSpectrum(
    "data/out/" * filename * "/4.seamass-restore_--reconstruct/" * filename * ".mzMLb",
    outputSpectrumID
  )

  # load mzMLb output spectrum
  mzmlbSpectrumOut = SeaMass.MzmlbSpectrum(
    "data/out/" * filename * "/2.seamass-restore/" * filename * ".mzMLb",
    outputSpectrumID
  )

  # load mzMLb centroided output spectrum
  mzmlbSpectrumOutCentroid = SeaMass.MzmlbSpectrum(
    "data/out/" * filename * "/5.seamass-restore_--centroid/" * filename * ".mzMLb",
    outputSpectrumID
  )

  # Compare mzMLb
  plot(
    smbSpectrumIn.locations,
    vcat(
      smbSpectrumIn.counts ./ (smbSpectrumIn.exposure .* (smbSpectrumIn.locations[2:end] - smbSpectrumIn.locations[1:end-1])),
      0.0,
    ),
    line = :steppost,
    label = "SMB input",
    title = "mzMLb comparison - " * filename,
    xlabel = "m/z (Th)",
    ylabel = "intensity density",
    reuse = false,
  )
  plot!(
    mzmlbSpectrumOutReconstruct.mzs,
    mzmlbSpectrumOutReconstruct.intensities,
    label = "mzMLb output (seaMass-restore --reconstruct)",
  )
  plot!(
    mzmlbSpectrumOut.mzs,
    mzmlbSpectrumOut.intensities,
    label = "mzMLb output (seaMass-restore)",
  )
  sticks!(
    mzmlbSpectrumOutCentroid.mzs,
    mzmlbSpectrumOutCentroid.intensities,
    label = "mzMLb output (seaMass-restore --centroid)",
    m = 4,
  )
  gui()

end

filename = "HYE124_TTOF6600_64var_lgillet_I150211_008__index_59994"
id = "p-55-227433333333"
inputSpectrumID = 1
outputSpectrumID = 1
test(filename, id, inputSpectrumID, outputSpectrumID)

filename = "P02U_Swath_1__mzWindow_602_605__scanTime_2300_3500"
id = "p"
inputSpectrumID = (67 - 1) * 51 + 1
outputSpectrumID = 67
test(filename, id, inputSpectrumID, outputSpectrumID)

println("Press <Enter> to finish")
readline()


text:
@testset "Snapshot loading test" begin
    mpSim = MPsim()
    filePath = joinpath( "base", "baseConfig" )
    fileName = string( filePath, ".xlsx" )

    XLSX.openxlsx( fileName ) do xf
        MPIO.readInitPop( mpSim, xf, filePath )
    end  # XLSX.openxlsx( fileName ) do xf

    @test mpSim.orgSize == 50
    @test length( mpSim.attributeList ) == 3
    @test length( mpSim.baseNodeList ) == 5
end  # @testset "Snapshot loading test"

ext: jl


text:
using PyCall
# Routines for interaction with spglib

"""
Construct a tuple containing the lattice and the positions of the species
in the convention required to take the place of a `cell` datastructure used in spglib.
"""
function spglib_cell_atommapping_(lattice, atoms)
    lattice = Matrix{Float64}(lattice)  # spglib operates in double precision
    n_attypes = isempty(atoms) ? 0 : sum(length(positions) for (type, positions) in atoms)
    spg_numbers = Vector{Int}(undef, n_attypes)
    spg_positions = Matrix{Float64}(undef, n_attypes, 3)

    offset = 0
    nextnumber = 1
    atommapping = Dict{Int, Any}()
    for (iatom, (type, positions)) in enumerate(atoms)
        atommapping[nextnumber] = type
        for (ipos, pos) in enumerate(positions)
            # assign the same number to all types with this position
            spg_numbers[offset + ipos] = nextnumber
            spg_positions[offset + ipos, :] = pos
        end
        offset += length(positions)
        nextnumber += 1
    end

    # Note: In the python interface of spglib the lattice vectors
    #       are given in rows, but DFTK uses columns
    #       For future reference: The C interface spglib also uses columns.
    (lattice', spg_positions, spg_numbers), atommapping
end
spglib_cell(lattice, atoms) = first(spglib_cell_atommapping_(lattice, atoms))

function spglib_get_symmetry(lattice, atoms; tol_symmetry=1e-5)
    spglib = pyimport("spglib")
    lattice = Matrix{Float64}(lattice)  # spglib operates in double precision

    if isempty(atoms)
        # spglib doesn't like no atoms, so we default to
        # no symmetries (even though there are lots)
        return [Mat3{Int}(I)], [Vec3(zeros(3))]
    end

    # Ask spglib for symmetry operations and for irreducible mesh
    spg_symops = spglib.get_symmetry(spglib_cell(lattice, atoms),
                                     symprec=tol_symmetry)

    # If spglib does not find symmetries give an error
    if spg_symops === nothing
        err_message=spglib.get_error_message()
        error("spglib failed to get the symmetries. Check your lattice, use a " *
              "uniform BZ mesh or disable symmetries. Spglib reported : " * err_message)
    end

    Stildes = [St for St in eachslice(spg_symops["rotations"]; dims=1)]
    τtildes = [rationalize.(τt, tol=tol_symmetry)
               for τt in eachslice(spg_symops["translations"]; dims=1)]
    @assert length(τtildes) == length(Stildes)

    # Checks: (A Stilde A^{-1}) is unitary
    for Stilde in Stildes
        Scart = lattice * Stilde * inv(lattice)  # Form S in cartesian coords
        if maximum(abs, Scart'Scart - I) > tol_symmetry
            error("spglib returned non-unitary rotation matrix")
        end
    end

    # Check (Stilde, τtilde) maps atoms to equivalent atoms in the lattice
    for (Stilde, τtilde) in zip(Stildes, τtildes)
        for (elem, positions) in atoms
            for coord in positions
                diffs = [rationalize.(Stilde * coord + τtilde - pos, tol=tol_symmetry)
                         for pos in positions]

                # If all elements of a difference in diffs is integer, then
                # Stilde * coord + τtilde and pos are equivalent lattice positions
                if !any(all(isinteger, d) for d in diffs)
                    error("Cannot map the atom at position $coord to another atom of the " *
                          "same element under the symmetry operation (Stilde, τtilde):\n" *
                          "($Stilde, $τtilde)")
                end
            end
        end
    end

    Stildes, τtildes
end

function spglib_standardize_cell(lattice::MatT, atoms; correct_symmetry=true,
                                 primitive=false, tol_symmetry=1e-5) where {MatT}
    spglib = pyimport("spglib")
    T = eltype(lattice)

    # Convert lattice and atoms to spglib and keep the mapping between our atoms
    # and spglibs atoms
    cell, atommapping = spglib_cell_atommapping_(lattice, atoms)

    # Ask spglib to standardize the cell (i.e. find a cell, which fits the spglib conventions)
    res = spglib.standardize_cell(spglib_cell(lattice, atoms), to_primitive=primitive,
                                  no_idealize=!correct_symmetry, symprec=tol_symmetry)
    spg_lattice, spg_scaled_positions, spg_numbers = res

    # Note: In the python interface of spglib the lattice vectors
    #       are given in rows, but DFTK uses columns
    #       For future reference: The C interface spglib also uses columns.
    newlattice = MatT(spg_lattice')
    newatoms = [(atommapping[iatom]
                 => T.(spg_scaled_positions[findall(isequal(iatom), spg_numbers), :]))
                for iatom in unique(spg_numbers)]
    newlattice, newatoms
end


text:
import Base: sin, cos, tan, cot, sec, csc, atan

import Base: sinh, cosh, tanh, coth, sech, csch

import Base: abs, exp, log

import Base: +, -, *, /, ^

struct pentacomplex
    i1
    i2
end

function mat(penta::pentacomplex)
    i1, i2 = penta.i1, penta.i2
    return [i1 -i2; i2 i1]
end

function mat2pentacomplex(matrix)
    out1, out2 = matrix[1, 1], matrix[2, 1]
    return pentacomplex(out1, out2)
end

## Basic operators: +, -, *, /, ^

function +(penta1::pentacomplex, penta2::pentacomplex)
    result=mat(penta1)+mat(penta2)
    return mat2pentacomplex(result)
end

+(n, penta::pentacomplex) = mat2pentacomplex(n .+mat(penta))

+(penta::pentacomplex, n) = +(n, penta)

-(penta1::pentacomplex, penta2::pentacomplex) = mat2pentacomplex(mat(penta1)-mat(penta2))
-(n, penta::pentacomplex) = mat2pentacomplex(n .-mat(penta))
-(penta::pentacomplex, n) = -(n, penta)
-(penta::pentacomplex) = mat2pentacomplex(-mat(penta))

*(penta1::pentacomplex, penta2::pentacomplex) = mat2pentacomplex(mat(penta1)*mat(penta2))
*(n, penta::pentacomplex) = mat2pentacomplex(n*mat(penta))
*(penta::pentacomplex, n) = *(n, penta)

/(penta1::pentacomplex, penta2::pentacomplex) = mat2pentacomplex(mat(penta1)/mat(penta2))
/(n, penta::pentacomplex) = mat2pentacomplex(n*inv(mat(penta)))
/(penta::pentacomplex, n) = mat2pentacomplex(mat(penta)/n)

^(penta::pentacomplex, power) = mat2pentacomplex(mat(penta)^power)

## Basic functions: abs, exp, log

abs(penta::pentacomplex) = sqrt(penta.i1^2+penta.i2^2)

function exp(penta::pentacomplex)
    i1, i2 = penta.i1, penta.i2
    out1 = exp(i1)*cos(i2)
    out2 = exp(i1)*sin(i2)
    return pentacomplex(out1, out2)
end

## Basic math functions: sin, cos, tan, cot, sec, csc

function sin(penta::pentacomplex)
    i1, i2 = penta.i1, penta.i2
    out1 = cosh(i2)*sin(i1)
    out2 = sinh(i2)*cos(i1)
    return pentacomplex(out1, out2)
end

function cos(penta::pentacomplex)
    i1, i2 = penta.i1, penta.i2
    out1 = cosh(i2)*cos(i1)
    out2 = -sinh(i2)*sin(i1)
    return pentacomplex(out1, out2)
end

tan(penta::pentacomplex) = sin(penta)/cos(penta)
cot(penta::pentacomplex) = cos(penta)/sin(penta)
sec(penta::pentacomplex) = 1/cos(penta)
csc(penta::pentacomplex) = 1/sin(penta)

## Basic hyperbolic functions: sinh, cosh, tanh, coth, sech, csch

function sinh(penta::pentacomplex)
    i1, i2 = penta.i1, penta.i2
    out1 = cosh(i1)*cos(i2)
    out2 = sinh(i1)*sin(i2)
    return pentacomplex(out1, out2)
end

function cosh(penta::pentacomplex)
    i1, i2 = penta.i1, penta.i2
    out1 = sinh(i1)*cos(i2)
    out2 = cosh(i1)*sin(i2)
    return pentacomplex(out1, out2)
end

tanh(penta::pentacomplex) = sinh(penta)/cosh(penta)
coth(penta::pentacomplex) = cosh(penta)/sinh(penta)
sech(penta::pentacomplex) = 1/cosh(penta)
csch(penta::pentacomplex) = 1/sinh(penta)

function image5(penta::pentacomplex)
    i1, i2 = penta.i1, penta.i2
    return image4(i2)
end

function pentaderivative(f, point, h)
    result = image5(f(pentacomplex(quadcomplex(tricomplex(bicomplex(point+im*h, h), h), h), h)))/h^5
    return result
end


text:
"""
$(TYPEDEF)

Backward Euler's integration scheme with given step-size.

### Fields

- `Δt` -- step-size

"""
struct BackwardEuler{N} <: AbstractSolver
    Δt::N
end

BackwardEuler(; Δt::N) where N = BackwardEuler(Δt)

function _init(alg::BackwardEuler, M, C, K)
    Δt = alg.Δt
    K̂ = K * Δt + C
    return K̂
end

function _solve(alg::BackwardEuler{N},
                ivp::InitialValueProblem{ST, XT},
                NSTEPS::Int; kwargs...) where {N, VT, ST, XT<:Tuple{VT, VT}}

    sys = system(ivp)
    (U₀, _) = initial_state(ivp)

    IMAX = NSTEPS + 1
    M, C, K, R = _unwrap(sys, IMAX)

    iszero(M) || throw(ArgumentError("this method assumes that the system is of first order"))

    K̂ = _init(alg, M, C, K)
    K̂⁻¹ = factorize(K̂)

    # initialize
    U = Vector{VT}(undef, IMAX)
    U[1] = U₀

    Δt = alg.Δt

    @inbounds for i in 1:NSTEPS
        R̂ᵢ₊₁ = R[i+1] * Δt + C*U[i]
        # solve
        U[i+1] = K̂⁻¹ \ R̂ᵢ₊₁
    end

    return _build_solution(alg, U, NSTEPS)
end

function _build_solution(alg::BackwardEuler{N}, U, NSTEPS) where {N}
    t = range(zero(N), step=alg.Δt, length=(NSTEPS+1))
    return Solution(alg, U, t)
end


text:
using ChromNetPaper
using Base.Test

include("evaluate.jl")
include("groupgm.jl")


text:
## Load everything
using Pkg
Pkg.activate(".")
Pkg.instantiate()
include("src/EasyStocks.jl")

## Setup
using CairoMakie # For plotting
mp = ModPar(q=0.005,ψ=0.015,xstar=2.0)
np = NumPar(mp)
MS = ModelSolution(mp,np)
## Run-code

@time SolveModel!(MS)

fig = Figure()
lines(fig[1,1],MS.np.xgrd, MS.α[:,1])
ylims!(-0.05,1.05)
lines(fig[1,2],MS.np.xgrd, MS.V[:,1])
ylims!(-6,1)
lines(fig[2,2],MS.np.xgrd, MS.V[:,2])
ylims!(-1.5,0)
lines(fig[2,1],MS.np.xgrd, MS.b[:,1] + MS.s[:,1])
ylims!(-0.05,5)
linkxaxes!(fig.content...)
fig


text:
using Printf
using TimerOutputs
using FFTW
using Oceananigans
using Oceananigans.Utils

include("benchmark_utils.jl")

#####
##### Benchmark setup and parameters
#####

const timer = TimerOutput()

Nr = 10  # Number of repeats for benchmarking.

# Run benchmark across these parameters.
            Ns = [(256, 256, 256)]
   float_types = [Float64]  # Float types to benchmark.
         archs = [CPU()]             # Architectures to benchmark on.
@hascuda archs = [CPU(), GPU()]      # Benchmark GPU on systems with CUDA-enabled GPUs.

#####
##### Run benchmarks
#####

function dim2xyz(d)
    d == 1 && return :x
    d == 2 && return :y
    d == 3 && return :z
end

for FT in float_types, N in Ns
    Nx, Ny, Nz = N

    #####
    ##### 1D FFT
    #####

    for dim in 1:3
        xyz = dim2xyz(dim)

        bn =  benchmark_name(N, "1D  FFT$xyz", CPU(), FT, npad=3)

        @info @sprintf("Planning transform: %s...\n", bn)
        R = rand(Complex{FT}, Nx, Ny, Nz)
        FFT! = FFTW.plan_fft!(R, dim, flags=FFTW.PATIENT)
        FFT! * R  # warmup

        @info @sprintf("Running benchmark:  %s...\n", bn)
        for _ in 1:Nr
            @timeit timer bn FFT! * R
        end
    end

    #####
    ##### 1D IFFT
    #####

    for dim in 1:3
        xyz = dim2xyz(dim)

        bn =  benchmark_name(N, "1D IFFT$xyz", CPU(), FT, npad=3)

        @info @sprintf("Planning transform: %s...\n", bn)
        R = rand(Complex{FT}, Nx, Ny, Nz)
        IFFT! = FFTW.plan_ifft!(R, 1, flags=FFTW.PATIENT)
        IFFT! * R  # warmup

        @info @sprintf("Running benchmark:  %s...\n", bn)
        for _ in 1:Nr
            @timeit timer bn IFFT! * R
        end
    end

    #####
    ##### 1D DCT
    #####

    bn = benchmark_name(N, "1D  DCTz", CPU(), FT, npad=3)

    @info @sprintf("Planning transform: %s...\n", bn)
    R = rand(FT, Nx, Ny, Nz)
    DCT! = FFTW.plan_r2r!(R, FFTW.REDFT10, 3, flags=FFTW.PATIENT)
    DCT! * R  # warmup

    @info @sprintf("Running benchmark:  %s...\n", bn)
    for _ in 1:Nr
        @timeit timer bn DCT! * R
    end

    #####
    ##### 1D IDCT
    #####

    bn = benchmark_name(N, "1D IDCTz", CPU(), FT, npad=3)

    @info @sprintf("Planning transform: %s...\n", bn)
    R = rand(FT, Nx, Ny, Nz)
    IDCT! = FFTW.plan_r2r!(R, FFTW.REDFT01, 3, flags=FFTW.PATIENT)
    IDCT! * R  # warmup

    @info @sprintf("Running benchmark:  %s...\n", bn)
    for _ in 1:Nr
        @timeit timer bn IDCT! * R
    end

    #####
    ##### 2D FFTxy
    #####

    bn =  benchmark_name(N, "2D  FFTxy ", CPU(), FT, npad=3)

    @info @sprintf("Planning transform: %s...\n", bn)
    R = rand(Complex{FT}, Nx, Ny, Nz)
    FFT! = FFTW.plan_fft!(R, [1, 2], flags=FFTW.PATIENT)
    FFT! * R  # warmup

    @info @sprintf("Running benchmark:  %s...\n", bn)
    for _ in 1:Nr
        @timeit timer bn FFT! * R
    end

    #####
    ##### 2D FFTxy (x then y)
    #####

    bn =  benchmark_name(N, "2D  FFTx,y", CPU(), FT, npad=3)

    @info @sprintf("Planning transform: %s...\n", bn)
    R = rand(Complex{FT}, Nx, Ny, Nz)
    FFTx! = FFTW.plan_fft!(R, 1, flags=FFTW.PATIENT)
    FFTy! = FFTW.plan_fft!(R, 2, flags=FFTW.PATIENT)
    FFTx! * R  # warmup
    FFTy! * R  # warmup

    @info @sprintf("Running benchmark:  %s...\n", bn)
    for _ in 1:Nr
        @timeit timer bn begin
            FFTx! * R
            FFTy! * R
        end
    end
end

#####
##### Print benchmark results
#####

println()
print(versioninfo_with_gpu())
print_timer(timer, title="Transform benchmarks", sortby=:name)
println()


text:
using SafeTestsets

@time begin
    @time @safetestset "IO Tests" begin include("io_test.jl") end
    @time @safetestset "Full Runs Tests" begin include("full_run_test.jl") end
    @time @safetestset "Post-process Tests" begin include("postprocess_test.jl") end
    @time @safetestset "Quadrature Tests" begin include("quadrature_test.jl") end
    @time @safetestset "Struct Tests" begin include("struct_test.jl") end
    @time @safetestset "Tracer Tests" begin include("tracers_test.jl") end
end


text:
_REIF_EQTO_FLOAT_EPSILON = 1.0e-5

"""
Bridges `CP.Reification{CP.DifferentFrom}` to MILP constraints.
"""
struct ReificationDifferentFrom2MILPBridge{T <: Real} <: MOIBC.AbstractBridge
    var_abs::MOI.VariableIndex
    var_abs_int::Union{Nothing, MOI.ConstraintIndex{MOI.VariableIndex, MOI.Integer}}
    con_abs::MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.AbsoluteValue}
    con_bigm::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}
    con_smallm::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}
end

function MOIBC.bridge_constraint(
    ::Type{ReificationDifferentFrom2MILPBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.Reification{CP.DifferentFrom{T}},
) where {T}
    return MOIBC.bridge_constraint(
        ReificationDifferentFrom2MILPBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{ReificationDifferentFrom2MILPBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.Reification{CP.DifferentFrom{T}},
) where {T <: Real}
    f_scalars = MOIU.scalarize(f)

    # Only change with respect to CP.Reification{MOI.EqualTo}: change the
    # sign of the first variable.
    f_scalars[1] = one(T) - f_scalars[1]

    # For this formulation work, both lower and upper bounds are required on
    # the constrained expression (but obviously not on the binary reified
    # variable).
    @assert CP.is_binary(model, f_scalars[1])
    @assert CP.has_lower_bound(model, f_scalars[2])
    @assert CP.has_upper_bound(model, f_scalars[2])

    # If the reified expression is true/false, then the EqualTo constraint
    # must/cannot be satisfied. (If the constraint is satisfied, the reified
    # expression is unconstrained.)
    if T <: Integer
        var_abs, var_abs_int = MOI.add_constrained_variable(model, MOI.Integer())
    else
        var_abs = MOI.add_variable(model)
        var_abs_int = nothing
    end

    con_abs = MOI.add_constraint(
        model,
        MOIU.vectorize(
            MOI.ScalarAffineFunction{T}[
                one(T) * var_abs,
                f_scalars[2] - s.set.value
            ]
        ),
        CP.AbsoluteValue()
    )

    bigm = T(max(
        abs(CP.get_upper_bound(model, f_scalars[2])),
        abs(CP.get_lower_bound(model, f_scalars[2]))
    ))
    con_bigm = MOI.add_constraint(
        model,
        var_abs - bigm * (one(T) - f_scalars[1]),
        MOI.LessThan(zero(T))
    )

    # If the constraint is satisfied, constrain the reified.
    smallm = if T <: Int
        one(T)
    else
        T(_REIF_EQTO_FLOAT_EPSILON)
    end

    con_smallm = MOI.add_constraint(
        model,
        f_scalars[1] - var_abs / smallm,
        MOI.LessThan(zero(T))
    )

    return ReificationDifferentFrom2MILPBridge{T}(var_abs, var_abs_int, con_abs, con_bigm, con_smallm)
end

function MOI.supports_constraint(
    ::Type{ReificationDifferentFrom2MILPBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.Reification{CP.DifferentFrom{T}}},
) where {T <: Real}
    return true
end

function MOIB.added_constrained_variable_types(::Type{ReificationDifferentFrom2MILPBridge{T}}) where {T <: Real}
    return Tuple{Type}[]
end

function MOIB.added_constrained_variable_types(::Type{ReificationDifferentFrom2MILPBridge{T}}) where {T <: Integer}
    return [(MOI.Integer,)]
end

function MOIB.added_constraint_types(::Type{ReificationDifferentFrom2MILPBridge{T}}) where {T <: Real}
    return [
        (MOI.VectorAffineFunction{T}, CP.AbsoluteValue),
        (MOI.VectorAffineFunction{T}, MOI.LessThan{T}),
    ]
end

function MOIB.added_constraint_types(::Type{ReificationDifferentFrom2MILPBridge{T}}) where {T <: Integer}
    return [
        (MOI.VectorAffineFunction{T}, CP.AbsoluteValue),
        (MOI.VectorAffineFunction{T}, MOI.LessThan{T}),
    ]
end

function MOI.get(::ReificationDifferentFrom2MILPBridge{T}, ::MOI.NumberOfVariables) where {T <: Real}
    return 1
end

function MOI.get(
    ::ReificationDifferentFrom2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VariableIndex, MOI.Integer,
    },
) where {T <: Real}
    return 0
end

function MOI.get(
    ::ReificationDifferentFrom2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VariableIndex, MOI.Integer,
    },
) where {T <: Integer}
    return 1
end

function MOI.get(
    ::ReificationDifferentFrom2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, CP.AbsoluteValue,
    },
) where {T <: Real}
    return 1
end

function MOI.get(
    ::ReificationDifferentFrom2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.LessThan{T},
    },
) where {T <: Real}
    return 2
end

function MOI.get(
    b::ReificationDifferentFrom2MILPBridge{T},
    ::MOI.ListOfVariableIndices,
) where {T <: Real}
    return [b.var_abs]
end

function MOI.get(
    b::ReificationDifferentFrom2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VariableIndex, MOI.Integer,
    },
) where {T <: Integer}
    return [b.var_abs_int]
end

function MOI.get(
    b::ReificationDifferentFrom2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, CP.AbsoluteValue,
    },
) where {T <: Real}
    return [b.con_abs]
end

function MOI.get(
    b::ReificationDifferentFrom2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.LessThan{T},
    },
) where {T <: Real}
    return [b.con_bigm, b.con_smallm]
end


text:
module TowerField12Test

using Test

using TowerFields.TowerField12: create_frobmap, mulroot
using TowerFields.BigIntField: powmod, frobmod, mulmod

const Q = 13
const frobmap = create_frobmap(Q)

@testset "frobmod: FQ" begin
    @test powmod(Q, mulroot, 7, Q) == frobmod(Q, frobmap, mulroot, 7)
end

@testset "frobmod: FQ2" begin
    @test powmod(Q, mulroot, (7, 11), Q) == frobmod(Q, frobmap, mulroot, (7, 11))
end

@testset "frobmod: FQ3" begin
    e = ((7, 11), (6, 10), (5, 9))
    @test powmod(Q, mulroot, e, Q) == frobmod(Q, frobmap, mulroot, e)
end

@testset "mulroot: FQ" begin
    @test mulroot(123) == -123
end

@testset "mulroot: FQ2" begin
    @test mulroot((3, 4)) == (-1, 7)
end

@testset "mulroot: FQ6" begin
    x = ((7, 11), (6, 10), (5, 9))
    @test mulroot(x) == ((-4, 14), (7, 11), (6, 10))
end

end


text:
using MathProgBase
function MathProgBase.initialize(d::SSMOptim, requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in [:Grad, :Jac])
            error("Unsupported feature $feat")
        end
    end
end

MathProgBase.features_available(d::SSMOptim) = [:Grad, :Jac]

function MathProgBase.eval_f(d::SSMOptim, x)
    ssm_obj_fun(d, x)
end

function MathProgBase.eval_g(d::SSMOptim, g, x)
end

function MathProgBase.eval_grad_f(d::SSMOptim, grad_f, x)
    ForwardDiff.gradient!(grad_f, y->ssm_obj_fun(d, y), x, d.d.jcfg, Val{false}())
end

MathProgBase.jac_structure(d::SSMOptim) = Int[],Int[]
# lower triangle only
MathProgBase.hesslag_structure(d::SSMOptim) = Int[]

function MathProgBase.eval_jac_g(d::SSMOptim, J, x)
end

function MathProgBase.eval_hesslag(d::SSMOptim, H, x, σ, μ)
end


text:
"""
The exported names in this module are:

$(EXPORTS)
"""
module Io

using DocStringExtensions
using JLD2
using FileIO

using Gridap.Helpers
import JSON
using BSON

export to_dict
export from_dict

export to_json
export to_json_file
export to_bson_file
export to_jld2_file
export from_json
export from_json_file
export from_bson_file
export from_jld2_file

include("IoInterfaces.jl")

include("Json.jl")

include("Bson.jl")

include("JLD2.jl")
end


text:
## INSTRUCTIONS:
## RUN IN LINE WITH SIMULATION MAIN FLOW CODE
## REQUIRES "data_file"

### UPDATE FUNCTION FOR DATA FILE ARRAYS

function UpdateArray(data_dictionary,string,alterations)
  # access array
  #println(data_dictionary)
  #println(alterations)
  data_array = data_dictionary[string];
  # set modifications for all conditions here
  for i in 1:Int(size(alterations)[1])
    updates = alterations[Int(i)];
    index = updates[1];
    value = updates[2];
    data_array[Int(index)]=value;
  data_dictionary[string] = data_array;
  end
  return data_array
end

##

#tgfb forward rate

forward_master = 0.5 #0.5

tgfb_add_forward_rate = forward_master #4.0
vegf_add_forward_rate = forward_master*10
vivit_add_forward_rate = forward_master*10

##
#### set values for multiple variables

# # # TGFB SIGNALING RATES
tgfb12_to_ap1_sp1_gain = 1.0 #1.5
tgfb3_to_smad_gain_rate = 2.0 *1.0 # TGFB3_R_phosphorylation_Smad: Smad =(TGFB3_R)=> pSmad # 2.0
tgfb12_to_smad_rate = 0.25
MAPK_activation_smad_rate = 0.2
# # # BETA-CATEININ RATES
bcat_deg_rate = 0.3 # default 0.1
# bcat_deg_gsk3_rate = 0.05 # default 0.1

bcat_trans_rate = 0.5-0.1
# # # BETA-CATEININ COMPLEXING
ecad_complex_gain = 12.0 *2 # 6 or 8
ecad_complex_reverse = 1.0 *2
ecad_bcat_degradation_rate = 0.3 + 0.3
lef1_complex_gain = 4.0
tcf4_complex_gain = 4.0 + 4
apc_complex_forward_rate=11.0 # 11.0
apc_bcat_deg_rate = 0.1
tgfb3_pik3_rate = 0.5  #5.0
gsk3_p_rate = 0.5 # 3.0
#ecad_bais_factor = 1.0 # multiply ecadherin-bcat only
ecad_trans_rate = 1.0+0.7+0.5# 1.5
ecad_deg_rate= 1.3 -0.9 -0.2 # 0.1

ecad_transc_rate = 1.0+0.4 # 1.5
ecad_mrna_deg_rate= 0.1 + 0.1 # 0.1

pi3k_translation_rate = 0.5 # 0.5
gsk3_translation_rate = 0.5 # 0.5

apc_complex_trans_rate = 0.6
lef1_trans_rate = 1.0 # 1.0
tcf4_trans_rate = 1.0 # 0.5

snail_slug_trans_rate = 1.0 # increases lef1 mrna too much

yreg_deg_rate = 1.5 # 1.5
# for ox/ko
yreg_ox_deg_rate = 0.15
yreg_ko_deg_rate_protein = 15

yreg_translation_rate = 1

### SELECTED DATA FILE INDEXES

# INDEX VALUES CHANGE WITH UPDATES TO MODEL.NET

# initial conditions
tgfb12_initial_cond_index = 88 # TGFB12
R_initial_cond_index = 87
vegf_initial_cond_index = 93 # VEGF

yreg_translation_rate_index = 152 # translation_mRNA_YREG1: mRNA_YREG1 = YREG1
yreg_protein_deg_rate_index = 48 # YREG1_degradation: YREG1 =([])=> []
# TGFB signaling
tgfb3_binding_rate_index = 3 # TGFB3_binding: TGFB3+R =([])=> TGFB3_R
R_deg_rate_index = 57 # degradation_R: R =([])=> []
TGFB12_R_activation_RAS_rate_index = 16 # TGFB12_R_activation_RAS: RAS =(TGFB12_R)=> RAS_GTP
RAS_activation_RAF_rate_index  = 17 # RAS_GTP_phosphorylation_RAF: RAF =(RAS_GTP)=> RAF_P
RAF_activation_MAPK_rate_index = 18 # RAF_P_phosphorylation_MAPK: MAPK =(RAF_P)=> pMAPK
MAPK_activation_AP1_SP1_rate_index = 20 # MAPK_phosphorylation_AP1_SP1: AP1_SP1 =(pMAPK)=> AP1_SP1_P
MAPK_activation_smad_rate_index = 19 # MAPK_phosphorylation_Smad: Smad =(pMAPK)=> pSmad
smad_deg_rate_index = 45 # Smad_degradation: Smad =([])=> []
psmad_deg_rate_index = 66 # pSmad_degradation: pSmad =([])=> []

# Bcatenin complex formation
bcat_ecad_forward_index = 21 # Binding_Ecadherin_Bcatenin: Ecadherin+Bcatenin =([])=> Ecadherin_Bcatenin
bcat_ecad_reverse_index = 22 # Binding_Ecadherin_Bcatenin: Ecadherin_Bcatenin =([])=> Ecadherin+Bcatenin (reverse)
#bcat_apc_reverse_index = 23
bcat_tcf4_forward_index = 25 # TCF4_Bcatenin_complex_formation: TCF4+Bcatenin =([])=> Active_TCF4
bcat_tcf4_reverse_index = 26 # TCF4_Bcatenin_complex_formation: Active_TCF4 =([])=> TCF4+Bcatenin (reverse)
bcat_lef1_forward_index = 27 # LEF1_Bcatenin_complex_formation: LEF1+Bcatenin =([])=> Active_LEF1
bcat_lef1_reverse_index = 28 # LEF1_Bcatenin_complex_formation: Active_LEF1 =([])=> LEF1+Bcatenin (reverse)
bcat_deg_rate_index = 43 # Bcatenin_degradation: Bcatenin =([])=> []
bcat_trans_rate_index = 149 # translation_mRNA_Bcatenin: mRNA_Bcatenin = Bcatenin

ecad_bcat_degradation_rate_index = 58 # Ecadherin_Bcatenin_degradation: Ecadherin_Bcatenin =([])=> []

# VEGF signaing
PKC_activation_RAF_rate_index = 11 # PKC_activation_RAF: RAF =(Active_PKC)=> RAF_P

snail_slug_trans_rate_index = 143 # induction_gene_SNAIL_SLUG: gene_SNAIL_SLUG = mRNA_SNAIL_SLUG
lef1_trans_rate_index = 153 # translation_mRNA_LEF1: mRNA_LEF1 = LEF1

apc_complex_forward_rate_index = 23 # Binding_APC_AXIN_forw: APC_AXIN+Bcatenin =([])=> APC_AXIN_Bcatenin
tgfb3_pik3_rate_index = 14 # TGFB3_R_phosphorylation_PI3K: PI3K =(TGFB3_R)=> PI3K_P
apc_bcat_deg_rate_index = 77 # degradation_APC_AXIN_Bcatenin: APC_AXIN_Bcatenin =([])=> []
gsk3_p_rate_index = 15 # PI3K_phosphorylation_GSK3: GSK3 =(PI3K_P)=> GSK3_P
pi3k_translation_rate_index = 169 # translation_mRNA_PI3K: mRNA_PI3K = PI3K
gsk3_translation_rate_index = 168	# translation_mRNA_GSK3: mRNA_GSK3 = GSK3
apc_complex_trans_rate_index = 167 # translation_mRNA_APC_AXIN: mRNA_APC_AXIN = APC_AXIN

#ecad bcat balance
ecad_trans_rate_index = 148 # translation_mRNA_Ecadherin: mRNA_Ecadherin = Ecadherin
ecad_deg_rate_index = 44 # Ecadherin_degrdation: Ecadherin =([])=> []
ecad_mrna_deg_rate_index = 102 # degradation_mRNA_Ecadherin: mRNA_Ecadherin = []
ecad_transc_rate_index = 125 # induction_gene_Ecadherin: gene_Ecadherin = mRNA_Ecadherin
# TGFB dosing
tgfb12_forward_rate_index = 91 # addition_TGFb12: [] =([])=> TGFB12
tgfb12_reverse_rate_index = 92 # addition_TGFb12: TGFB12 =([])=> [] (reverse)
tgfb12_binding_rate_index = 1 # TGFB12_binding: TGFB12+R =([])=> TGFB12_R

# VEGF dosing
vegf_forward_rate_index = 93 # addition_VEGF: [] =([])=> VEGF
vegf_reverse_rate_index = 94 # addition_VEGF: VEGF =([])=> [] (reverse)
vegf_binding_rate_index = 5 # VEGF_binding: VEGF+VEGFR =([])=> VEGF_VEGFR

# TGFB3_Ab dosing
TGFB3_Ab_foward_rate_index = 86 # addition_TGFB3_Ab: [] =([])=> TGFB3_Ab (reverse)
TGFB3_Ab_reverse_rate_index = 85 # addition_TGFB3_Ab: TGFB3_Ab =([])=> []

# U0126 dosing
U0126_foward_rate_index = 98 # addition_U0126: [] =([])=> U0126 (reverse)
U0126_reverse_rate_index = 97 # addition_U0126: U0126 =([])=> []

# TGFB3 dosing
TGFB3_foward_rate_index = 95 # addition_TGFB3: [] =([])=> TGFB3
TGFB3_reverse_rate_index = 96 # addition_TGFB3: TGFB3 =([])=> [] (reverse)
tgfb12_to_smad_rate_index = 12 # TGFB12_R_phosphorylation_Smad: Smad =(TGFB12_R)=> pSmad
tgfb3_to_smad_gain_rate_index = 13 # TGFB3_R_phosphorylation_Smad: Smad =(TGFB3_R)=> pSmad

# vivit

vivit_forward_rate_index  = 100 # addition_vivit: [] =([])=> vivit (reverse)
vivit_reverse_rate_index = 99 # addition_vivit: vivit =([])=> []

# Snail slug transfeciton
ss_virus_forward_rate_index = 88 # addition_ss_virus: [] =([])=> ss_virus (reverse)
ss_virus_reverse_rate_index = 87 # addition_ss_virus: ss_virus =([])=> []

# Ecadherin transfection
ecad_virus_forward_rate_index = 90 # addition_ecad_virus: [] =([])=> ecad_virus (reverse)
ecad_virus_reverse_rate_index = 89 # addition_ecad_virus: ecad_virus =([])=> []
# DN smad
# DN SMAD temoprary
smad_induction_rate_index = 124 # induction_gene_Smad: gene_Smad = mRNA_Smad
# redundant psmad_deg_rate_index = 66 # pSmad_degradation: pSmad =([])=> []
# reverse convention based on product term
#smad_virus_forward_rate_index = 146 # addition_smad_virus: [] =([])=> smad_virus (reverse)
#smad_virus_reverse_rate_index = 122 # addition_smad_virus: smad_virus =([])=> []

# DN 1
lef1_induction_rate_index = 130 # induction_gene_LEF1: gene_LEF1 = mRNA_LEF1
# reverse convention based on product term
#LEF1_virus_forward_rate_index = 18 # addition_lef1_virus: [] =([])=> lef1_virus (reverse)
#LEF1_virus_reverse_rate_index = 111 # addition_lef1_virus: lef1_virus =([])=> []

lef1_deg_rate_index = 50 # LEF1_degradation: LEF1 =([])=> []
Active_LEF1_deg_rate_index = 51 # Active_LEF1_degradation: Active_LEF1 =([])=> []

#yreg_deg_rate_index = 122 # degradation_mRNA_YREG1: mRNA_YREG1 = []
yreg_deg_rate_index = 106 # degradation_mRNA_YREG1: mRNA_YREG1 = []

tcf4_trans_rate_index = 150 # translation_mRNA_TCF4: mRNA_TCF4 = TCF4

# MAPK inhibitor
# ?
raf_mapk_phos_rate_index = 18 # RAF_P_phosphorylation_MAPK: MAPK =(RAF_P)=> pMAPK

#### THESE SETTINGS ARE FOR STEADY STATE AND ARE DEFAULT ####

# initialize arrays to hold updates for each array in data_file
initial_cond_update_array = [];
rate_constant_update_array = [];
control_parameter_update_array = [];

# initial conditions
push!(initial_cond_update_array,[tgfb12_initial_cond_index,0.0]) # TGFB12 initial conditions
push!(initial_cond_update_array,[R_initial_cond_index,6.0]) # R initial condition
push!(initial_cond_update_array,[vegf_initial_cond_index,0.0]) # VEGF

# rates
# TGFB dosing
push!(rate_constant_update_array,[tgfb12_forward_rate_index,0.0]) # addition_TGFb12: [] =([])=> TGFB12
push!(rate_constant_update_array,[tgfb12_reverse_rate_index,0.0]) # addition_TGFb12: TGFB12 =([])=> [] (reverse)
push!(rate_constant_update_array,[tgfb12_binding_rate_index,3.0]) # TGFB12_binding: TGFB12+R =([])=> TGFB12_R

# VEGF dosing
push!(rate_constant_update_array,[vegf_forward_rate_index,0.0]) # addition_VEGF: [] =([])=> VEGF
push!(rate_constant_update_array,[vegf_reverse_rate_index,0.0]) # addition_VEGF: VEGF =([])=> [] (reverse)
push!(rate_constant_update_array,[vegf_binding_rate_index,3.0]) # VEGF_binding: VEGF+VEGFR =([])=> VEGF_VEGFR

# ecad deg and trans
push!(rate_constant_update_array,[ecad_deg_rate_index,ecad_deg_rate]) #default 0.1
push!(rate_constant_update_array,[ecad_trans_rate_index ,ecad_trans_rate]) #default 1.0

## Inhibitor dosing

# SNAIL_SLUG dosing
# push!(rate_constant_update_array,[snail_slug_forward_rate_index,0.0]) # addition_VEGF: [] =([])=> VEGF
# push!(rate_constant_update_array,[snail_slug_reverse_rate_index,0.0]) # addition_VEGF: VEGF =([])=> [] (reverse)

# TGFB3_Ab dosing
push!(rate_constant_update_array,[TGFB3_Ab_foward_rate_index,0.0])
push!(rate_constant_update_array,[TGFB3_Ab_reverse_rate_index,0.0])

# U0126 dosing
push!(rate_constant_update_array,[U0126_foward_rate_index ,0.0])
push!(rate_constant_update_array,[U0126_reverse_rate_index,0.0])

# vivit dosing
push!(rate_constant_update_array,[vivit_forward_rate_index,0.0])
push!(rate_constant_update_array,[vivit_reverse_rate_index,0.0])

## Virus dosing

push!(rate_constant_update_array,[ss_virus_forward_rate_index,0.0])
push!(rate_constant_update_array,[ss_virus_reverse_rate_index,0.0])

push!(rate_constant_update_array,[ecad_virus_forward_rate_index,0.0])
push!(rate_constant_update_array,[ecad_virus_reverse_rate_index,0.0])

#push!(rate_constant_update_array,[smad_virus_forward_rate_index,0.0])
#push!(rate_constant_update_array,[smad_virus_reverse_rate_index,0.0])

#push!(rate_constant_update_array,[LEF1_virus_forward_rate_index,0.0])
#push!(rate_constant_update_array,[LEF1_virus_reverse_rate_index,0.0])

# TGFB3 dosing

push!(rate_constant_update_array,[TGFB3_foward_rate_index,0.0])
push!(rate_constant_update_array,[TGFB3_reverse_rate_index,0.0])

##
# TGFB signaling
push!(rate_constant_update_array,[tgfb3_binding_rate_index,3.5]) # TGFB3_binding: TGFB3+R =([])=> TGFB3_R
push!(rate_constant_update_array,[tgfb3_to_smad_gain_rate_index,tgfb3_to_smad_gain_rate]) # TGFB3_binding: TGFB3+R =([])=> TGFB3_R

push!(rate_constant_update_array,[tgfb12_to_smad_rate_index,tgfb12_to_smad_rate]) # TGFB3_binding: TGFB3+R =([])=> TGFB3_R

push!(rate_constant_update_array,[R_deg_rate_index,0.01]) # R_degradation: R =([])=> []
push!(rate_constant_update_array,[TGFB12_R_activation_RAS_rate_index,1.0]) #TGFB12_R_activation_RAS: RAS =(TGFB12_R)=> RAS_GTP
push!(rate_constant_update_array,[RAS_activation_RAF_rate_index,tgfb12_to_ap1_sp1_gain]) # RAS_GTP_phosphorylation_RAF: RAF =(RAS_GTP)=> RAF_P
push!(rate_constant_update_array,[RAF_activation_MAPK_rate_index,tgfb12_to_ap1_sp1_gain]) # RAF_P_phosphorylation_MAPK: MAPK =(RAF_P)=> pMAPK
push!(rate_constant_update_array,[MAPK_activation_AP1_SP1_rate_index,tgfb12_to_ap1_sp1_gain]) # MAPK_phosphorylation_AP1_SP1: AP1_SP1 =(pMAPK)=> AP1_SP1_P
push!(rate_constant_update_array,[MAPK_activation_smad_rate_index,MAPK_activation_smad_rate]) # MAPK_phosphorylation_AP1_SP1: AP1_SP1 =(pMAPK)=> AP1_SP1_P

# VEGF signaing
push!(rate_constant_update_array,[PKC_activation_RAF_rate_index,0.0]) # PKC_activation_RAF: RAF =(Active_PKC)=> RAF_P

# Bcatenin complex formation
push!(rate_constant_update_array,[bcat_ecad_forward_index,ecad_complex_gain]) # 19	Binding_Ecadherin_Bcatenin: Ecadherin+Bcatenin =([])=> Ecadherin_Bcatenin
push!(rate_constant_update_array,[bcat_ecad_reverse_index,ecad_complex_reverse]) # 20	Binding_Ecadherin_Bcatenin: Ecadherin_Bcatenin =([])=> Ecadherin+Bcatenin (reverse)
push!(rate_constant_update_array,[bcat_tcf4_forward_index,tcf4_complex_gain]) # 21	TCF4_Bcatenin_complex_formation: TCF4+Bcatenin =([])=> Active_TCF4
push!(rate_constant_update_array,[bcat_tcf4_reverse_index,1.0]) # 22	TCF4_Bcatenin_complex_formation: Active_TCF4 =([])=> TCF4+Bcatenin (reverse)
push!(rate_constant_update_array,[bcat_lef1_forward_index,lef1_complex_gain]) # 23	LEF1_Bcatenin_complex_formation: LEF1+Bcatenin =([])=> Active_LEF1
push!(rate_constant_update_array,[bcat_lef1_reverse_index,1.0]) # 24	LEF1_Bcatenin_complex_formation: Active_LEF1 =([])=> LEF1+Bcatenin (reverse)

push!(rate_constant_update_array,[bcat_deg_rate_index,bcat_deg_rate]) # Bcatenin_degradation: Bcatenin =([])=> []
#push!(rate_constant_update_array,[bcat_deg_gsk3_mediated_index,bcat_deg_gsk3_rate])
#push!(rate_constant_update_array,[pmapk_inact_gsk3_index,pmapk_inact_gsk3_rate])

push!(rate_constant_update_array,[yreg_deg_rate_index,yreg_deg_rate])
push!(rate_constant_update_array,[bcat_trans_rate_index,bcat_trans_rate])
push!(rate_constant_update_array,[ecad_bcat_degradation_rate_index,ecad_bcat_degradation_rate])

# LEF1
push!(rate_constant_update_array,[lef1_trans_rate_index,lef1_trans_rate])

#
push!(rate_constant_update_array,[apc_complex_forward_rate_index,apc_complex_forward_rate])
push!(rate_constant_update_array,[tgfb3_pik3_rate_index,tgfb3_pik3_rate])
push!(rate_constant_update_array,[gsk3_p_rate_index,gsk3_p_rate])
push!(rate_constant_update_array,[apc_bcat_deg_rate_index,apc_bcat_deg_rate])
push!(rate_constant_update_array,[pi3k_translation_rate_index,pi3k_translation_rate])
push!(rate_constant_update_array,[gsk3_translation_rate_index,gsk3_translation_rate])

push!(rate_constant_update_array,[apc_complex_trans_rate_index,apc_complex_trans_rate])

push!(rate_constant_update_array,[tcf4_trans_rate_index ,tcf4_trans_rate])

push!(rate_constant_update_array,[yreg_translation_rate_index ,yreg_translation_rate])

push!(rate_constant_update_array,[ecad_transc_rate_index,ecad_transc_rate]) # addition_TGFb12: [] =([])=> TGFB12
push!(rate_constant_update_array,[ecad_mrna_deg_rate_index,ecad_mrna_deg_rate]) # addition_TGFb12: [] =([])=> TGFB12

# control parameter array
#num_regs = 15; # number of inducing or repressing intereactions (could have multiple for species)
# do forced update for now
#for i in 1:num_regs
  #push!(control_parameter_update_array[i,1],GAIN_MASTER) # Gain for transcriptional activity
#end
#push!(control_parameter_update_array,[7,GAIN_TCF4_TGFB3])
#push!(control_parameter_update_array,[6,GAIN_BCAT_TGFB3])
#push!(control_parameter_update_array,[5,GAIN_ACTIVE_LEF1_ECAD])

# push modifications to RAM
#println(data_dictionary["INITIAL_CONDITION_ARRAY"])
#println(data_dictionary)
#UpdateArray(data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
#UpdateArray(data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)
#UpdateArray(data_dictionary,"CONTROL_PARAMETER_ARRAY",control_parameter_update_array)
###  I can't figure out how to get data_dictionary to be passed into this streamed file,
### It would seem that the stream opens a new variable space and then merges it with the main variable space
### Hence, these will need to be called post-stream

#


text:
println("Hello world!");

@everywhere using Printf;

@show workdir = if length(ARGS) > 0
  ARGS[1];
else
  ".";
end

@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x));
@everywhere fmt_int(x) = @sprintf("%03d", x);

@everywhere function prefix(case)
  "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))";
end

cases = Any[];
E0s = [1e-1; 1.0; 10.0];
K1s = [0.0; 1e-1; 0.5; 1.0; 2.0];
K2s = [0.0; 1e-1; 0.5; 1.0; 2.0];
kTs = [1.0];
Fzs = [0.0; 0.5; 1.0; 2.0];
Fxs = [0.0; 0.5; 1.0; 2.0];
ns = Int[100; 200];
bs = [0.5; 1.0; 2.0];
for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, E0 in E0s, K1 in K1s, K2 in K2s
  if K1==K2; continue; end
  push!(cases, Dict(:E0 => E0, :K1 => K1, :K2 => K2,
                    :kT => kT, :Fz => Fz,
                    :Fx => Fx, :n => n, :b => b));
end

@info "total number of cases to run: $(length(cases))";

pmap(case -> begin;
  outfile = joinpath(workdir, "$(prefix(case)).out");
  if !isfile(outfile)
    println("Running case: $case.");
    command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type interacting -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 250000 -v 2 --prefix $(joinpath(workdir, prefix(case)))`;
    output = read(command, String);
    write(outfile, output);
  else
    println("Case: $case has already been run.");
  end
end, cases);

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats


text:
immutable Token
    kind::Int
    contents::UTF8String
end

immutable Variable
	varname::UTF8String
	vartype::UTF8String
	description::UTF8String
end

type DocumentationChunk
    name::UTF8String
    description::UTF8String
    args::Vector{Variable}
    fields::Vector{Variable}
    returns::Vector{Variable}
    examples::UTF8String
    exported::Bool
end

function DocumentationChunk()
    DocumentationChunk("",
    	               "",
					   Variable[],
					   Variable[],
					   Variable[],
					   "",
					   false)
end


text:
 # Predicts whether a student will pass a course using the course's passrate
# as a probability.

module PassRate
    # Train the model
    function train(curriculum)
        for course in curriculum.courses
            model = Dict()
            model[:passrate] = course.passrate
            course.model = model
        end

        curriculum.stopoutModel[:rates] = [0.0838, 0.1334, 0.0465, 0.0631, 0.0368, 0.0189, 0.0165] * 100
    end

    # Predict grade
    function predict_grade(course, student)
        roll = rand()

        if roll <= course.model[:passrate]
            return 4.0
        else
            return 0.0
        end
    end

    # Predict stopout
    function predict_stopout(student, currentTerm, model)
        if currentTerm > 7
            return false
        else
            roll = rand(1:100)
            return roll <= model[:rates][currentTerm]
        end
    end
end

# Helper functions

# Sets all course model passrates to given value
function setPassrates(courses, passrate)
    for course in courses
        course.passrate = passrate
    end
end

max_stars_repo_name: M-S-H/CASL

max_issues_repo_name: M-S-H/CASL

max_forks_repo_name: M-S-H/CASL


text:
module ObjectPools

export new!, recycle!, ArrayPool

using ArgCheck: @argcheck
using Parameters: @unpack
using DocStringExtensions: SIGNATURES, TYPEDEF

####
#### general interface
####

"""
    new!(pool, specifications...)

Return a yet unused object from `pool` according to `specifications`. The interpretation of
specifications depends on the pool.
"""
function new! end

####
#### specific pools
####

"""
$(TYPEDEF)

This type is used for internal implementation only.
"""
mutable struct HomogeneousArrayPool{T,N,D}
    arrays::Vector{Array{T,N}}
    used::Int
    """
    $(SIGNATURES)

    An object pool for homogeneous arrays with element type `T` and dimensions `D`.
    """
    function HomogeneousArrayPool{T}(D::NTuple{N,Int}) where {T,N}
        @argcheck all(d -> d isa Int && d ≥ 0, D)
        new{T,N,D}(Vector{Array{T,N}}(), 0)
    end
end

"""
$(SIGNATURES)

A new array of the type and dimensions determined by the pool.
"""
function _new!(pool::HomogeneousArrayPool{T,N,D}) where {T,N,D}
    @unpack arrays = pool
    pool.used += 1
    if length(arrays) ≥ pool.used
        arrays[pool.used]
    else
        A = Array{T}(undef, D)
        push!(arrays, A)
        A
    end
end

function _recycle!(pool::HomogeneousArrayPool)
    pool.used = 0
    nothing
end

struct ArrayPool
    arrays::Dict{Any,Any}
    @doc """
    $(SIGNATURES)

    An object pool of `Array`s with heterogeneous element types and dimensions.
    """
    ArrayPool() = new(Dict{Any,Any}())
end

"""
    new!(pool, Array{T}, dims)
    new!(pool, Array{T}, dims...)

A new `Array{T}` with dimensions `dims` (can also be specified in a splat syntax).
"""
function new!(pool::ArrayPool, ::Type{<:Array{T,N}}, dims::NTuple{N,Int}) where {T,N}
    _pool = get!(() -> HomogeneousArrayPool{T}(dims), pool.arrays, (T, dims...))
    _new!(_pool)::Array{T,N}
end

new!(pool::ArrayPool, T, dims::Int...) = new!(pool, T, dims)

function new!(pool::ArrayPool, ::Type{Array{T}}, dims::NTuple{N,Int}) where {T,N}
    new!(pool, Array{T,N}, dims)
end

function recycle!(pool::ArrayPool)
    foreach(_recycle!, values(pool.arrays))
    nothing
end

end # module

@doc read(joinpath(dirname(@__DIR__), "README.md"), String) ObjectPools


text:
export initialize, terminate

function initialize()
    err = ccall((:Pa_Initialize, portaudio), Int32, ())
    check_pa_error(err)
end

function terminate()
    err = ccall((:Pa_Terminate, portaudio), Int32, ())
    check_pa_error(err)
end

function initialize(f::Function)
    initialize()
    f()
    terminate()
end


text:
# Bogumił Kamiński, 2019-03-25

# Start Julia run with the following command:
#     julia -p auto
#
# .int file format specification:
#     * Each image is encoded by 1+28^2 bytes.
#     * Images are stored consecutively.
#     * First byte for each image is its class from 0 to 9
#     * Then each of remaining 28^2 has value from 0 to 255
#       and encodes one entry in a 28x28 image

using Statistics
using Distributed

try
    using PyPlot
    @info "PyPlot loaded"
catch
    @warn "PyPlot missing"
end

@everywhere const TRAIN = (class=UInt8[], image=Matrix{Int32}[])
@everywhere const TEST = (class=UInt8[], image=Matrix{Int32}[])

@everywhere for (filename, data) in [("mnist_train.int", TRAIN),
                                     ("mnist_test.int", TEST)]
    open(filename) do f
        while !eof(f)
            c = read(f, UInt8)
            v = read(f, 28^2)
            push!(data.class, c)
            push!(data.image, reshape(v, 28, 28))
        end
    end
end

@everywhere function distance(a, b)
    d = 0
    @simd for i in 1:length(a)
        @inbounds d += (a[i] - b[i]) ^ 2
    end
    d
end

@everywhere function knnacc(i)
    dist = distance.(Ref(TEST.image[i]), TRAIN.image)
    knn_loc = partialsortperm(dist, 1:20)
    TRAIN.class[knn_loc]
end

@time neivec = pmap(knnacc, axes(TEST.class, 1))
const NEIGHBOURS = reduce(hcat, neivec)

function acc(classes, testclass)
    v = zeros(Int, 10)
    for c in classes
        v[c + 1] += 1
    end
    m = maximum(v)
    (v[testclass + 1] == m) / count(==(m), v)
end

function evalk(k)
    mean(axes(NEIGHBOURS, 2)) do i
        classes = view(NEIGHBOURS, 1:k, i)
        testclass = TEST.class[i]
        acc(classes, testclass)
    end
end

evalk.(axes(NEIGHBOURS, 1))

max_stars_repo_path: 6. MNIST/mnist_1_knn.jl

max_issues_repo_path: 6. MNIST/mnist_1_knn.jl

max_forks_repo_path: 6. MNIST/mnist_1_knn.jl


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -6.148474362391325 <= q <= 10.677081718106185)

                     add_NL_constraint(m, :(sigmoid(-0.2518902526786948 + 0.9847866884384731*sigmoid(0.2752793536861313 + -0.1568397657479923*sigmoid(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + -0.1062303426643516*sigmoid(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + 0.5642457669318222*sigmoid(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))) + -0.7649756572019379*sigmoid(0.9918790053549298 + 0.85981905523253*sigmoid(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + 0.9639515148457134*sigmoid(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + -0.08816238542180388*sigmoid(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))) + 0.10087195663512549*sigmoid(-0.8202416682813327 + -0.17766510212211006*sigmoid(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + -0.1793020696087071*sigmoid(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + -0.41892665263312834*sigmoid(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5])))) + sigmoid(0.638435553115452 + 0.09453389081965424*sigmoid(0.2752793536861313 + -0.1568397657479923*sigmoid(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + -0.1062303426643516*sigmoid(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + 0.5642457669318222*sigmoid(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))) + -0.07927014075141203*sigmoid(0.9918790053549298 + 0.85981905523253*sigmoid(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + 0.9639515148457134*sigmoid(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + -0.08816238542180388*sigmoid(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))) + 0.2624391124261729*sigmoid(-0.8202416682813327 + -0.17766510212211006*sigmoid(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + -0.1793020696087071*sigmoid(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + -0.41892665263312834*sigmoid(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5])))) + sigmoid(0.5273095272255199 + 0.4025366346817978*sigmoid(0.2752793536861313 + -0.1568397657479923*sigmoid(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + -0.1062303426643516*sigmoid(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + 0.5642457669318222*sigmoid(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))) + -0.585290488444234*sigmoid(0.9918790053549298 + 0.85981905523253*sigmoid(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + 0.9639515148457134*sigmoid(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + -0.08816238542180388*sigmoid(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))) + 0.1647823489958915*sigmoid(-0.8202416682813327 + -0.17766510212211006*sigmoid(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + -0.1793020696087071*sigmoid(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + -0.41892665263312834*sigmoid(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#RSDeltaSigmaPort: A port of Richard Schreier's Delta Sigma Toolbox
#-------------------------------------------------------------------------------
#__precompile__(false) #Bit faster to develop when using ControlSystems.jl
#=Conventions
 - Function name `_stairs__()`
   - Leading `_`: Private function (not meant to be called externally from this
     module).
   - Trailing `__`: Implementation for a specific function (do not call unless
     from within that parent function).
=#

"""`RSDeltaSigmaPort`: A port of Richard Schreier's Delta Sigma Toolbox

# Sample usage
```julia-repl
julia> using RSDeltaSigmaPort #Will take a while to load, compile, etc...
julia> import RSDeltaSigmaPort: @runsample

julia> @runsample("dsdemo1.jl")
julia> @runsample("dsdemo2.jl")
julia> @runsample("dsdemo3.jl")
julia> @runsample("dsdemo4_audio.jl")
julia> @runsample("dsdemo5.jl")
julia> @runsample("dsdemo6.jl")
julia> @runsample("dsexample1.jl")
julia> @runsample("dsexample2.jl")
julia> @runsample("demoLPandBP.jl")
```

# See also:
[`simulateDSM`](@ref), [`simulateMS`](@ref), [`simulateSNR`](@ref), [`simulateHBF`](@ref) |
[`synthesizeNTF`](@ref), [`realizeNTF`](@ref), [`realizeNTF_ct`](@ref) |
[`calculateSNR`](@ref), [`peakSNR`](@ref), [`predictSNR`](@ref) |
[`calculateTF`](@ref), [`evalTF`](@ref), [`evalTFP`](@ref) |
[`stuffABCD`](@ref), [`scaleABCD`](@ref), [`mapABCD`](@ref), [`partitionABCD`](@ref) |
[`mapCtoD`](@ref), [`mapQtoR`](@ref) |
[`exampleHBF`](@ref), [`designHBF`](@ref), [`simulateHBF`](@ref) |
[`pulse`](@ref), [`impL1`](@ref) |
[`lollipop`](@ref), [`logsmooth`](@ref) |
[`documentNTF`](@ref), [`plotExampleSpectrum`](@ref)
"""
module RSDeltaSigmaPort

const rootpath = realpath(joinpath(@__DIR__, ".."))

#=Flags:
	-TODO
	-VERIFYME: is this the correct behaviour???
   -VERSIONSWITCH: Algorithm originally depends on software version (Needs to be checked)
   -NEEDSTOOLKIT
	-REQUESTHELP
=#

import Random
import Statistics: mean
import LinearAlgebra
import LinearAlgebra: norm, diagm, eigen, cond
import SpecialFunctions: erfinv
import Interpolations
import FFTW: fft, fftshift
import DSP: conv, filt, remez
import Optim
import Optim: optimize, GoldenSection, IPNewton, GradientDescent, LBFGS
import Optim: Fminbox, TwiceDifferentiableConstraints
import Polynomials
import Polynomials: Polynomial
import Printf: @sprintf
using CMDimData
using CMDimData.MDDatasets
using CMDimData.EasyPlot
using CMDimData.EasyPlot.Colors #Should be ok with whatever version it needs
import CMDimData.EasyPlot: BoundingBox #Use this one to avoid version conflicts with Graphics.
using InspectDR
CMDimData.@includepkg EasyPlotInspect

import ControlSystems
import ControlSystems: isdiscrete
#using ControlSystems: zpk, zpkdata

function throw_unreachable()
	msg = "Code expected to be unreachable with exception handling"
	error(msg)
end

function throw_notimplemented()
	msg = "Feature not implemented"
	error(msg)
end

const j = im
include("base.jl")
include("arrays.jl")
include("mlfun.jl")
include("mlfun_control.jl")
include("text.jl")
include("datasets.jl")
include("signals_time.jl")
include("power.jl")
include("timedomain.jl")
include("windowing.jl")
include("snr.jl")
include("optimize.jl")
include("transferfunctions.jl")
include("statespace.jl")
include("mapCtoD.jl")
include("synthesizeNTF.jl")
include("quantizer.jl")
include("simulate_base.jl")
include("simulateDSM.jl")
include("simulateMS.jl")
include("realizeNTF.jl")
include("realizeNTF_ct.jl")
include("filter_base.jl")
include("exampleHBF.jl")
include("designHBF.jl")
include("designLCBP.jl")
include("simulateHBF.jl")
include("calc_spectrum.jl")
include("plot_base.jl")
include("plot_transient.jl")
include("plot_spectrum.jl")
include("plot_NTF.jl")
include("plot_zplane.jl")
include("plot_SNR.jl")
include("plot_state.jl")
include("display.jl")

#==Convenience macros
===============================================================================#
macro runsample(filename::String)
	fnstr = "$filename"
	fpath = realpath(joinpath(rootpath, "sample", fnstr))
	@info("Running RSDeltaSigmaPort sample:\n$fpath...")
	m = quote
		include($fpath)
	end
	return esc(m) #esc: Evaluate in calling module
end

#==Exported interface
===============================================================================#

#Simple calculations:
export dbm, dbp, dbv, rms
export undbm, undbp, undbv

#Compatibility:
#WARN: Might eventually collide with names from other packages:
#NOTE: Some use "_" prefix to avoid collisions with similar functions in other packages
export cplxpair
export _ss, _zpk, _zpkdata
export _zp2ss, _zp2tf, _freqz, _tf, _minreal, _impulse
export interp1_lin, interp1_cubic #Specialized names
export interp
export squeeze
export eye, orth, eig

#Data/arrays:
export AbstractDSM, RealDSM, QuadratureDSM
export padl, padr, padt, padb
export array_round #Converts a floating-point value to an (Int) index usable by Julia arrays.
export mapQtoR

#Accessors
export isquadrature

#Signal generators:
export ds_f1f2, default_ftest, default_fband
export genTestTone, genTestTone_quad, genTestTone_sin
export pulse, impL1

#Windowing/time domain:
export ds_hann, applywnd
export ds_therm
export sinc_decimate

#Misc. calculations:
export calcSpecInfo
export calcSNRInfo, calculateSNR, peakSNR, predictSNR

#Transfer functions:
export synthesizeNTF, realizeNTF, realizeNTF_ct
export calculateTF, evalTF, evalTFP
export rmsGain, cancelPZ

#State space:
export stuffABCD, scaleABCD, mapABCD, partitionABCD
export mapCtoD

#Filters:
export exampleHBF, designHBF, designLCBP

#Simulations:
export simulateDSM, simulateMS
export simulateSNR, simulateHBF

#Plotting:
export waveform, wfrm_stairs, lollipop, logsmooth
export plotPZ, plotPZ!, plotNTF, plotNTF!
export documentNTF
export plotSNR, plotSNR!, plotStateMaxima, plotUsage
export plotModTransient, plotLollipop
export plotSpec, plotSpec!, plotModSpectrum, plotModSpectrum!
export plotSpectrum, plotSpectrum!, plotExampleSpectrum

#Displaying plots:
export inlinedisp, saveimage, displaygui

#Display/text:
export ds_orderString, str_modulatortype

end # module


text:
using KernelFunctions
using Test
using Flux

@testset "Params" begin
    ν = 2.0; c = 3.0; d = 2.0; γ = 2.0; α = 2.5
    kc = ConstantKernel(c=c)
    @test all(params(kc) .== params([c]))
    km = MaternKernel(ν=ν)
    @test all(params(km) .== params([ν]))
    kl = LinearKernel(c=c)
    @test all(params(kl) .== params([c]))
    kge = GammaExponentialKernel(γ=γ)
    @test all(params(kge) .== params([γ]))
    kgr = GammaRationalQuadraticKernel(γ=γ, α=α)
    @test all(params(kgr) .== params([α], [γ]))
    kp = PolynomialKernel(c=c, d=d)
    @test all(params(kp) .== params([d], [c]))
    kr = RationalQuadraticKernel(α=α)
    @test all(params(kr) .== params([α]))

    k = km + kc
    @test all(params(k) .== params([k.weights], km, kc))

    k = km * kc
    @test all(params(k) .== params(km, kc))

    s = 2.0
    k = transform(km, s)
    @test all(params(k) .== params([s], km))

    v = [2.0]
    k = transform(kc, v)
    @test all(params(k) .== params(v, kc))

    P = rand(3, 2)
    k = transform(km,LowRankTransform(P))
    @test all(params(k) .== params(P, km))

    k = transform(km, LowRankTransform(P) ∘ ScaleTransform(s))
    @test all(params(k) .== params([s], P, km))

    c = Chain(Dense(3, 2))
    k = transform(km, FunctionTransform(c))
    @test all(params(k) .== params(c, km))

end


text:
using TravelingSalesmanHeuristics
using Base.Test
using Distances

###
# helpers
###

# generate a Euclidean distance matrix for n points in the unit square
function generate_planar_distmat(n)
	pts = rand(2, n)
	return pairwise(Euclidean(), pts, pts)
end

# test that a path is acceptable:
# - at most one city appears twice, and if so must be first and last
# - all values 1:n are present where n is the maximum city index, and no other values
function testpathvalidity(path, iscycle)
	if iscycle
		@test path[1] == path[end]
		pop!(path)
	end
	path = sort(path)
	for i in 1:length(path)
		@test path[i] == i
	end
end

###
# main tests
###

function test_nearest_neighbor()
	distmats = cell(2)
	distmats[1] = generate_planar_distmat(10)
	distmats[2] = generate_planar_distmat(2)

	for dm in distmats
		n = size(dm, 1)
		randstartcity = rand(1:n)
		# standard
		path, cost = nearest_neighbor(dm)
		@test cost > 0
		testpathvalidity(path, true)
		# repetitive
		path, cost = nearest_neighbor(dm, repetitive = true)
		@test cost > 0
		testpathvalidity(path, true)
		# no loop, 2 opt
		path, cost = nearest_neighbor(dm, closepath = false)
		@test cost > 0
		testpathvalidity(path, false)
		# no loop, no 2 opt
		path, cost = nearest_neighbor(dm, closepath = false, do2opt = false)
		@test cost > 0
		testpathvalidity(path, false)
		# fixed start, no loop
		path, cost = nearest_neighbor(dm, closepath = false, firstcity = Nullable(randstartcity))
		@test cost > 0
		testpathvalidity(path, false)
	end

	dm_bad = rand(3,2)
	@test_throws ErrorException nearest_neighbor(dm_bad)
end

function test_cheapest_insertion()
	# default random start
	dm = generate_planar_distmat(8)
	path, cost = cheapest_insertion(dm)
	@test cost > 0
	testpathvalidity(path, true) # should be a closed path
	# repetitive start
	path, cost = cheapest_insertion(dm, repetitive = true)
	@test cost > 0
	testpathvalidity(path, true)

	# bad input
	dmbad = rand(2,3)
	@test_throws ErrorException cheapest_insertion(dmbad)
end

function test_simulated_annealing()
	dm = generate_planar_distmat(14)
	# single start
	path, cost = simulated_annealing(dm)
	@test cost > 0
	testpathvalidity(path, true) # closed path
	# multi-start
	dm = generate_planar_distmat(8)
	path, cost = simulated_annealing(dm, num_starts = 3)
	@test cost > 0
	testpathvalidity(path, true) # also closed
	# given init path
	init_path = collect(1:8)
	push!(init_path, 1)
	reverse!(init_path, 2, 6)
	path, cost = simulated_annealing(dm; init_path = Nullable(init_path))
	@test cost > lowerbound(dm)
	testpathvalidity(path, true) # still closed
end

function test_bounds()
	dm = generate_planar_distmat(2)
	path, cost = solve_tsp(dm)
	lb = lowerbound(dm)
	@test lb <= cost
	dm = generate_planar_distmat(30)
	path, cost = solve_tsp(dm)
	lb = lowerbound(dm)
	@test lb <= cost
end

function test_path_costs()
	dm = [1 2 3; 4 5 6; 7 8 9]
	p1 = [1, 2, 3, 1]
	@test TravelingSalesmanHeuristics.pathcost(dm, p1) == 2 + 6 + 7
	# various reverse indices to test
	revs = [(1,4), (2,3), (2,4), (3,3)]
	for rev in revs
		path_rev = reverse(p1, rev[1], rev[2])
		@test TravelingSalesmanHeuristics.pathcost(dm, path_rev) ==
			  TravelingSalesmanHeuristics.pathcost_rev(dm, p1, rev[1], rev[2])
	end
end

###
# run
###
srand(47)
test_nearest_neighbor()
test_cheapest_insertion()
test_simulated_annealing()
test_bounds()
test_path_costs()
println("Done testing.")


text:
module MultiDictsTest

using MultiDicts
using Base.Test

function test{K,V}(::Type{K}, ::Type{V})
    @assert Char <: K
    @assert Int <: V
    d = MultiDict{K,V}()
    @test isempty(d)
    d['a'] = 1
    @test d['a'] == 1
    d['b'] = 2
    @test d['b'] == 2
    d['a'] = 1
    @test d['a'] == 1
    addindex!(d, 'b')
    @test d['b'] == 2
    i = get!(d, 'a', 2)
    @test i == 1
    j = get!(d, 'c', 3)
    @test j == 3
    delete!(d, 'a')
    @test d['a'] == 1
    delete!(d, 'a')
    @test !haskey(d, 'a')
    @test haskey(d, 'c')
    delete!(d, 'c')
    @test !haskey(d, 'c')
    @test !isempty(d)
    empty!(d)
    @test isempty(d)
end

function main()
    test(Char, Int)
    test(Char, Any)
    test(Any, Int)
    test(Any, Any)
end

end


text:
# License for this file: MIT
# Copyright 2018-2019, DLR Institute of System Dynamics and Control
#
"""
   module Fluid4_Components

Date: Feb. 25, 2019

Fluid component models with ModiaMedia based on the new approach described in

- Zimmer, Bender, Pollok (2018): Robust Modeling of Directed Thermofluid Flows in Complex Networks.
  2nd Japanese Modelica Conference, Tokyo, Proceedings, pp 39-48. Download:
  https://www.modelica.org/events/modelica2018japan/conference-proceedings/modelica-final-proceedings-2018-Japan.pdf

Difference to Fluid1_Components: Connector with (Medium,m_flow,r,state) and Medium is also in state

# Main Authors of this module
```

      /|      Dirk Zimmer, Martin Otter
  ___/_|___   German Aerospace Center (DLR)
 /  /  /  /   Institute of System Dynamics and Control
/__/__/__/    https://www.dlr.de/sr/en/
   | /
   |/

```
"""
module Fluid4_Components

using  Modia
using  ModiaMedia
import ModiaMath

using  ModiaMedia.StaticArrays # included via ModiaMedia, to avoid requirement to add it in the standard environment
using  ModiaMedia.Unitful      # included via ModiaMedia, to avoid requirement to add it in the standard environment

export add!
export BaseProperties
export StaticPressure, InertialPressure, SpecificEnthalpy, Temperature, MassFlowRate
export InPlug, OutPlug, FixedSource_pT, FixedSink, ShortPipe, Splitter, Junction, Compressor, ClosedVolume

# Variable types
StaticPressure(  ; args...) = Variable(; start=1.0e5, nominal=1.0e5, size = (), T = u"Pa"  , info = "Pressure for static mass-flow", args...)
InertialPressure(; args...) = Variable(; start=0.0  , nominal=1.0e5, size = (), T = u"Pa"  , info = "Pressure for accelerating the fluid", args...)
SpecificEnthalpy(; args...) = Variable(; start=1.0e4, nominal=1.0e4, size = (), T = u"J/kg", info = "Specific enthalpy", args...)
Temperature(     ; args...) = Variable(; start=300.0, nominal=300.0, size = (), T = u"K"   , info = "Temperature", args...)
MassFlowRate(    ; args...) = Variable(; start=0.0  ,                size = (), T = u"kg/s", info = "Mass flow rate", args...)

const MediumVariable() = Var(size=())
const MediumState()    = Var(size=())

###### Define state-dependent variables for plotting, depending on medium
const dependentVariables = Dict{AbstractString,Function}("p"   => ModiaMedia.pressure,
                                                         "T"   => ModiaMedia.temperature,
                                                         "d"   => ModiaMedia.density,
                                                         "h"   => ModiaMedia.specificEnthalpy,
                                                         "u"   => ModiaMedia.specificInternalEnergy,
                                                         "Cp"  => ModiaMedia.specificHeatCapacityCp,
                                                         "eta" => ModiaMedia.dynamicViscosity)

ModiaMath.variablesDependingOnStruct(state::ModiaMedia.SimpleMediumState)         = dependentVariables
ModiaMath.variablesDependingOnStruct(state::ModiaMedia.SimpleIdealGasMediumState) = dependentVariables

###### Define BaseProperties models supported here (user can provide further definitions)

@model SimpleMedium_BaseProperties begin
    p_start = 1.0e5                # parameter
    T_start = 300.0                # parameter

    Medium = MediumVariable()      # input
    p      = StaticPressure(  start=p_start)  # input
    h      = SpecificEnthalpy(start=1.0)      # input
    T      = Temperature(     start=T_start)  # output
    d      = Float(start=1.0)      # output
    u      = Float(start=1.0)      # output
    der_d  = Float(start=1.0)      # output
    der_u  = Float(start=1.0)      # output
    state  = MediumState()         # output

@equations begin
    d     = density(Medium)
    u     = specificInternalEnergy_T(Medium,T)
    state = setState_pT(Medium,p,T)
    h     = specificEnthalpy_T(Medium,T)
    der_d = 0.0
    der_u = specificInternalEnergy_T_der_2(Medium,T)*der(T)
    end
end

@model SimpleIdealGasMedium_BaseProperties begin
    p_start = 1.0e5                # parameter
    T_start = 300.0                # parameter

    Medium = MediumVariable()      # input
    p      = StaticPressure(  start=p_start)  # input
    h      = SpecificEnthalpy(start=1.0)      # input
    T      = Temperature(     start=T_start)  # output
    d      = Float(start=1.0)      # output
    u      = Float(start=1.0)      # output
    der_d  = Float(start=1.0)      # output
    der_u  = Float(start=1.0)      # output
    state  = MediumState()         # output

@equations begin
    d     = density_pT(Medium,p,T)
    u     = specificInternalEnergy_T(Medium,T)
    h     = specificEnthalpy_T(Medium,T)
    der_d = density_pT_der_2(Medium,p,T)*der(p) + density_pT_der_3(Medium,p,T)*der(T)
    der_u = specificInternalEnergy_T_der_2(Medium,T)*der(T)
    state = setState_pT(Medium,p,T)
    end
end

BaseProperties(Medium::ModiaMedia.SimpleMedium        ; p_start=1e5, T_start=300.0, args...) = SimpleMedium_BaseProperties(        Medium=Medium, p_start=p_start, T_start=T_start, args...)
BaseProperties(Medium::ModiaMedia.SimpleIdealGasMedium; p_start=1e5, T_start=300.0, args...) = SimpleIdealGasMedium_BaseProperties(Medium=Medium, p_start=p_start, T_start=T_start, args...)

####### Define Fluid components

"""
    connector InPlug - Uni-directional flow of fluid (mass flow is entering the component)
"""
@model InPlug(:connector) begin
    m_flow = MassFlowRate(flow = true)
    r      = InertialPressure()
    state  = MediumState()
end

"""
    connector OutPlug - Uni-directional flow of fluid (mass flow is leaving the component)
"""
@model OutPlug(:connector) begin
    m_flow = MassFlowRate(flow = true)
    r      = InertialPressure()
    state  = MediumState()
end

"""
    model FixedSource_pT - Fixed pressure and fixed specific enthalpy at source
"""
@model FixedSource_pT begin
    Medium  = MediumVariable()    # Medium MUST be redefined when instanting Source, otherwise error
    outPlug = OutPlug()
    p0 = StaticPressure(variability=parameter, info = "Fixed static pressure at source")
    T0 = Temperature(   variability=parameter, info = "Fixed temperature at source")
    @equations begin
        outPlug.r      = 0.0
        outPlug.state  = setState_pT(Medium, p0, T0)
    end
end

"""
    model FixedSink - Fixed pressure at sink
"""
@model FixedSink begin
    inPlug = InPlug()
    p0 = StaticPressure(variability=parameter, info = "Fixed static pressure at sink")
    @equations begin
        pressure(inPlug.state) + inPlug.r = p0
    end
end

"""
    model ShortPipe - Simple pipe (without volume)
"""
@model ShortPipe begin
    inPlug  = InPlug()
    outPlug = OutPlug()

    k = Parameter(1.0e5, info="Quadratic pressure drop coefficient")
    l = Parameter(1.0  , info="Length of pipe")
    A = Parameter(1.0  , info="Area of pipe")    # should be taken from a global world model or propagated through the connection structure

    dp     = StaticPressure(  info = "delta pressure in flow direction")
    dr     = InertialPressure(info = "delta inertial pressure in flow direction")
    m_flow = MassFlowRate()

    eta   = Float(size=())
    @equations begin
        # mass flow balance
        m_flow = inPlug.m_flow
        inPlug.m_flow + outPlug.m_flow = 0

        # pressure drop
        dp = -m_flow*abs(m_flow) * k

        # Propagation of specific quantities
        outPlug.r     = inPlug.r + dr;
        outPlug.state = isenthalpicState(inPlug.state, dp)

        # inertial component of pressure
        der(m_flow)/A*l = -dr

        # compute dynamic viscosity to test dispatch of function
        eta = dynamicViscosity(inPlug.state)
    end
end

"""
    model Splitter - one inflow, two outflows
"""
@model Splitter begin
    inPlugA  = InPlug()
    outPlugB = OutPlug()
    outPlugC = OutPlug()

    eps = 0.000001

    @equations begin
        # mass flow balance
        inPlugA.m_flow + outPlugB.m_flow + outPlugC.m_flow = 0.0

        # propagation of specific quantities
        outPlugB.r     = inPlugA.r
        outPlugB.state = inPlugA.state

        outPlugC.r     = inPlugA.r
        outPlugC.state = inPlugA.state
    end
end

"""
    model Junction - two inflows, one outflow
"""
@model Junction begin
    inPlugA  = InPlug()
    inPlugB  = InPlug()
    outPlugC = OutPlug()

    eps = 0.000001

    # local variables
    pA = StaticPressure()
    pB = StaticPressure()
    pC = StaticPressure()
    hA = SpecificEnthalpy()
    hB = SpecificEnthalpy()
    hC = SpecificEnthalpy()

    @equations begin
        # mass flow balance
        inPlugA.m_flow + inPlugB.m_flow + outPlugC.m_flow = 0.0

        # Get p and h from inflowing fluids
        pA = pressure(        inPlugA.state)
        pB = pressure(        inPlugB.state)
        hA = specificEnthalpy(inPlugA.state)
        hB = specificEnthalpy(inPlugB.state)

        # total pressure balance
        pA + inPlugA.r = pC + outPlugC.r
        pB + inPlugB.r = pC + outPlugC.r

        # mixing of fluid
            # formulation safeguarded against zero-massflow and backward flow
            (abs(outPlugC.m_flow) + 2*eps)*hC = (abs(inPlugA.m_flow) + eps)*hA + (abs(inPlugB.m_flow) + eps)*hB

            # for simplicity m_flow is taken for the weighted mean of p but it should be V_flow
            (abs(outPlugC.m_flow) + 2*eps)*pC = (abs(inPlugA.m_flow) + eps)*pA + (abs(inPlugB.m_flow) + eps)*pB

        # construct outPlug state
        outPlugC.state = setState_ph(getMedium(inPlugA.state), pC, hC)
    end
end

"""
    model Compressor - Simple compressor
"""
@model Compressor begin
    inPlug  = InPlug()
    outPlug = OutPlug()

    pressureRatio = Parameter(2.0, info="Fixed pressure ratio")
    efficiency    = Parameter(0.9)
    kappa         = Parameter(1.4)
    R             = Parameter(8.314, info="Gas constant")

    k = Parameter(1.0e5, info="Quadratic pressure drop coefficient")
    l = Parameter(1.0  , info="Length of compressor")
    A = Parameter(1.0  , info="Area of compressor")    # should be taken from a global world model or propagated through the connection structure

    dr     = InertialPressure(info = "delta inertial pressure in flow direction")
    dh     = SpecificEnthalpy(info = "delta in spec. enthalpy in flow direction")
    m_flow = MassFlowRate()
    y      = Float(info="Specific energy inserted into the mass-flow")
    k_frac = Float()

    pin = StaticPressure()
    hin = SpecificEnthalpy()
    @equations begin
        k_frac = kappa/(kappa-1.0)

        # mass flow balance
        m_flow = inPlug.m_flow
        inPlug.m_flow + outPlug.m_flow = 0

        # Change of h due to compressor
        y  = k_frac*R*(pressureRatio^(1/k_frac)-1)
        dh = y/efficiency;
        pin = pressure(        inPlug.state)
        hin = specificEnthalpy(inPlug.state)

        # Propagation of specific quantities
        outPlug.r     = inPlug.r + dr
        outPlug.state = setState_ph(getMedium(inPlug.state), pressureRatio*pin, hin+dh)

        # inertial component of pressure
        der(m_flow)/A*l = -dr
    end
end

"""
    model ClosedVolume - Simple volume
"""
@model ClosedVolume begin
    Medium  = AbstractMedium    # Medium MUST be redefined when instanting ClosedVolume, otherwise error

    inPlug  = InPlug()
    outPlug = OutPlug()

    V  = 1.0       # Fixed volume
    p0 = Medium.infos.p_default     # initial pressure
    T0 = Medium.infos.T_default     # initial temperature

    M     = Float(info="Mass of volume"      , start=1.0)
    der_M = Float(info="Total mass flow rate", start=0.0)

    #medium = getBaseProperties(Medium)(p_start=p0, T_start=T0)
    medium = BaseProperties(Medium; p_start=p0, T_start=T0)

    pin  = StaticPressure()
    hin  = SpecificEnthalpy()
    hout = SpecificEnthalpy()
    h    = SpecificEnthalpy()
@equations begin
    # Input quantities
    pin  = pressure(        inPlug.state)
    hin  = specificEnthalpy(inPlug.state)
    hout = specificEnthalpy(outPlug.state)

    # Mass balance
    der_M          = inPlug.m_flow + outPlug.m_flow
    medium.der_d*V = der_M

    # Energy balance
    M = medium.d*V
    der_M*medium.u + M*medium.der_u = inPlug.m_flow*(hin + inPlug.r/medium.d) + outPlug.m_flow*hout

    # Pressure at input
    pin + inPlug.r = medium.p

    # Properties at outPlug
    outPlug.r     = 0.0
    outPlug.state = medium.state

    # For testing
    h = specificEnthalpy(medium.state)
    end
end

end


text:
export DispersiveMode
export DispersiveField
export calc_disp

# Calculate the unnormalized Ω from normalized ω.
function ω2Ω(ω::AbstractRange{<:Real}, ω₀::Real, unit::MaxwellUnit)
    ωmin, ωmax = ω[1], ω[end]
    ωmin ≤ ω₀ ≤ ωmax || @error "ω₀ = $ω₀ should be within range of ω."

    # Convert the units of dm's quantities to the SI units.
    ω₀ *= unit.ω
    ωmin *= unit.ω
    ωmax *= unit.ω
    Nω = length(ω)

    # Construct Ω without catastrophic cancellation.
    Ωmin, Ωmax = ωmin-ω₀, ωmax-ω₀
    Ω = range(Ωmin, Ωmax, length=Nω)

    return Ω
end

# The type of ω is limited to AbstractRange in order to avoid catastrophic cancellation in
# Ω = ω .- ω₀.  Ω is constructed by range(ω[1]-ω₀, ω[end]-ω₀, length(ω)).
const DispersiveMode{K,Kₑ,Kₘ,K₊₁,VR<:AbstractRange{<:Real}} = ParametrizedMode{K,Kₑ,Kₘ,1,K₊₁,NamedTuple{(:ω,),Tuple{VR}}}

DispersiveMode(ω::AbstractRange{<:Real}, mdl::Model) = (θ = (ω=ω,); ParametrizedMode(θ,mdl))

struct DispersiveField{K}
    Fᵣ_fun::Array{Spline1D,K}
    Fᵢ_fun::Array{Spline1D,K}
end

function DispersiveField(Ω::AbstractRange{<:Real},
                         F::AbsArrNumber{K₊₁};  # F[Ω,x,y,z]
                         α::Number=1.0  # scale factor
                         ) where {K₊₁}
    size(F,1)==length(Ω) || @error "size(F,1) = $(size(F,1)) and length(Ω) = $(length(Ω)) should be same."

    K = K₊₁ - 1
    dim_shp = SVec(ntuple(k->k+1, Val(K)))  # 2:K+1
    sz_shp = size(F)[dim_shp]
    CI = CartesianIndices(sz_shp)

    bc = "extrapolate"
    Fᵣ_fun = Array{Spline1D}(undef, sz_shp)
    Fᵢ_fun = Array{Spline1D}(undef, sz_shp)
    for ci = CI
        Fᵣ_fun[ci] = Spline1D(Ω, real.(α .* @view(F[:,ci])); bc)
        Fᵢ_fun[ci] = Spline1D(Ω, imag.(α .* @view(F[:,ci])); bc)
    end

    return DispersiveField(Fᵣ_fun, Fᵢ_fun)
end

Base.size(df::DispersiveField) = size(df.Fᵣ_fun)

# Interpolate a given field at Ω.
function (df::DispersiveField)(Ω::Real)
    sz_shp = size(df)
    F = ArrComplexF(undef, sz_shp)

    CI = CartesianIndices(F)
    for ci = CI
        F[ci] = df.Fᵣ_fun[ci](Ω) + im * df.Fᵢ_fun[ci](Ω)
    end

    return F
end

# Interpolate a given tuple of fields at Ω.
(df::NTuple{N,DispersiveField})(Ω::Real) where {N} = ntuple(k->df[k](Ω), Val(N))

# Create interpolated quantities related to dispersion.  These quantities are usually
# calculated around a specific frequency ω₀, which is provided to the function.
#
# This often generates the final results for publication, so the outputs are described in
# the SI units.
#
# The function is used after constructing dm by scan_param().
function calc_disp(dm::DispersiveMode, ω₀::Real, unit::MaxwellUnit)
    # Convert the units of dm's quantities to the SI units.
    ω = dm.θ.ω
    Ω = ω2Ω(ω, ω₀, unit)
    ω *= unit.ω  # ω is newly allocated because *= instead of .*= is used

    β = real.(dm.β) ./ unit.L  # β is newly allocated
    neff = c₀ .* β ./ ω  # phase velocity is c = c₀ / neff = ω / β

    # Construct the interpolators for neff and β.
    bc = "extrapolate"
    neff_fun = Spline1D(Ω, neff; bc)
    β_fun = Spline1D(Ω, β; bc)

    # Construct the interpolators for phase mismatch ∆β.
    β₀ = β_fun(0.0)
    β₁₀ = derivative(β_fun, 0.0)
    ∆β = β .- β₀ .- β₁₀ .* Ω
    ∆β_fun = Spline1D(Ω, ∆β; bc)

    # Construct the interpolator for group velocity vg.
    # Another option is to define a function that evaluates the derivatives on the fly as
    # vg_fun(Ω) = 1 / derivative(β_fun, Ω), but this turns out to be very slow.  So we
    # evaluate the derivatives at the sampling points and construct a new interpolator for
    # the derivatives.
    β₁ = derivative(β_fun, Ω)
    β₁_fun = Spline1D(Ω, β₁; bc)

    vg = 1 ./ β₁
    vg_fun = Spline1D(Ω, vg; bc)

    # Construct the interpolator for GVD parameter.
    β₂ = derivative(β_fun, Ω, nu=2)  # nu: order of derivative
    β₂_fun = Spline1D(Ω, β₂; bc)

    D = -ω.^2 .* β₂ ./ 2πc₀  # D = -(ω²/2πc₀) β₂
    D_fun = Spline1D(Ω, D; bc)

    # dm's E and H are normalized to have unit longitudinal energy flux, but the energy flux
    # is in the units of unit.P, which is not 1.  Once E, H, and length in transverse
    # dimension are multiplied by unit.E, unit.H, and unit.L, energy flux is multiplied by
    # unit.P.  In order to normalize E and H in the SI units whose longitudinal energy flux
    # is unit.P, we need to normalize both fields by the same factor, which is √unit.P.
    E_fun = DispersiveField.(Ref(Ω), dm.E, α=unit.E/√unit.P)
    H_fun = DispersiveField.(Ref(Ω), dm.H, α=unit.H/√unit.P)

    # Do not return ω and ω₀ to leave them normalized outside this function.
    return (Ω=Ω, neff=neff_fun, β=β_fun, β₁=β₁_fun, β₂=β₂_fun, ∆β=∆β_fun, vg=vg_fun, D=D_fun, E=E_fun, H=H_fun)
end


text:
if environ("DASH_APP_LOCATION") != "ABSOLUTE"
    import .server.app
else
    import server.app
end

@doc """
This is a helper function which returns an empty string if
an environment variable corresponding to `var` cannot be found,
or the value of this variable when it exists.
"""
function environ(var)
  if var in keys(ENV)
    return ENV[var]
  else
    return ""
  end
end

@doc """
This function simplifies prepending relative paths
with (optional) paths derived from the `DASH_DOCS_URL_PREFIX`
environment variable. In Dash Enterprise, the help content is
typically hosted within `/Docs`, so `relpath("/somepath")`
would return "/Docs/somepath".
"""ex
function dash_relpath(path)
    if startswith(path, "/") && "DASH_DOCS_URL_PREFIX" in keys(ENV) && !startswith(path, string("/", strip(environ("DASH_DOCS_URL_PREFIX"), ['/'])))
        # In enterprise docs, all assets are under `/Docs`
        return string(rstrip(environ("DASH_DOCS_URL_PREFIX"), ['/']), path)
    end
    return path
end

@doc """
This is a helper function which catches errors, and prints a
useful console message describing the path of the file containing
the code which threw the error. The original error is then
rethrown. Julia offers macros, which are more powerful than
Python's decorators, but here we can just wrap a function by
prepending `exception_handler`, and appending the arguments, e.g.
`exception_handler(load_example)("/somepath", true)`
"""
function exception_handler(func)
    function wrapper(path, kwargs...)
        try
            return(func(path, kwargs...))
        catch e
            println(string("\nError running ", path, "\n", '='^76))
            rethrow()
        end
    end
    return wrapper
end

@doc """
`load_examples` is a helper function for assembling an
array of examples which are referenced in the page
rendered from `index_filename`. It works in concert with
`load_example` below.
"""
function load_examples(index_filename, omit = [])
    dir = dirname(relpath(index_filename))
    example_dir = joinpath(dir, "examples")
    try
        example_filenames = readdir(example_dir)
    catch
        return []
    end

    examples = []
    for filename in example_filenames
        full_filename = joinpath(example_dir, filename)
        if !(filename in omit) && isfile(full_filename) && endswith(filename, ".py")
            push!(examples, load_example(full_filename))
        end
    end
    return examples
end

function load_example(path, relative_path=false)
    if relative_path
        thispath = @__DIR__
        # in Julia, if path is absolute, it will omit
        # thispath in the statement below, so should
        # probably use lstrip to avoid this behaviour
        path = joinpath(thispath, lstrip(path, ['/']))
    end
    open(path) do file
        _source = read(file, string)
        _example = _source

        if !occursin("app = dash()", _example) && !occursin("app = CustomDash()", _example)
            throw("Didn't declare app")
        end

        _example = replace(_example, r"app = dash()" => "# app = dash()")

        commented_configs = [
            "app.scripts.config.serve_locally"
            "app.css.config.serve_locally"
        ]

        for config in commented_configs
            _example = replace(_example, Regex("# " * config))
        end

        if !occursin("using Dash\n", _example)
            throw("Didn't import dash")
        end

        if !occursin("app.layout = ", _example)
            throw("app.layout not assigned")
        end

        _example = replace(_example, r"app.layout =" => "layout =")

        if !occursin(r"run_server\(.*?\)", _example)
            throw("run_server(app) missing")
        end

        if occursin(r"$tools", _example)
            _example = replace(_example, r"$tools" => dirname(@__FILE__))
        end

        # replace remote datasets with local ones
        # so that the app can run in internet-less environments
        find_and_replace = Dict(
            "https://raw.githubusercontent.com/plotly/datasets/master/gapminderDataFiveYear.csv"=>
            "datasets/gapminderDataFiveYear.csv",

            "https://raw.githubusercontent.com/plotly/datasets/master/gapminder2007.csv"=>
            "datasets/gapminder2007.csv",

            "https://raw.githubusercontent.com/plotly/datasets/master/solar.csv"=>
            "datasets/solar.csv",

            "https://raw.githubusercontent.com/plotly/datasets/master/Emissions%20Data.csv"=>
            "datasets/Emissions%20Data.csv",

            "https://raw.githubusercontent.com/plotly/datasets/master/1962_2006_walmart_store_openings.csv"=>
            "datasets/1962_2006_walmart_store_openings.csv",

            "https://upload.wikimedia.org/wikipedia/commons/e/e4/Mitochondria%2C_mammalian_lung_-_TEM_%282%29.jpg"=>
            "datasets/mitochondria.jpg",

            "https://gist.githubusercontent.com/chriddyp/cb5392c35661370d95f300086accea51/raw/8e0768211f6b747c0db42a9ce9a0937dafcbd8b2/indicators.csv"=>
            "datasets/indicators.csv",

            "https://gist.githubusercontent.com/chriddyp/c78bf172206ce24f77d6363a2d754b59/raw/c353e8ef842413cae56ae3920b8fd78468aa4cb2/usa-agricultural-exports-2011.csv"=>
            "datasets/usa-agricultural-exports-2011.csv",

            "https://gist.githubusercontent.com/chriddyp/5d1ea79569ed194d432e56108a04d188/raw/a9f9e8076b837d541398e999dcbac2b2826a81f8/gdp-life-exp-2007.csv"=>
            "datasets/gdp-life-exp-2007.csv",

            "https://plotly.github.io/datasets/country_indicators.csv"=>
            "datasets/country_indicators.csv",

            "https://github.com/plotly/datasets/raw/master/26k-consumer-complaints.csv"=>
            "datasets/26k-consumer-complaints.csv",

            "https://js.cytoscape.org/demos/colajs-graph/data.json"=>
            "datasets/colajs-graph-data.json",

            "https://js.cytoscape.org/demos/colajs-graph/cy-style.json"=>
            "datasets/colajs-graph-cy-style.json",

            "https://www.publicdomainpictures.net/pictures/60000/nahled/flower-outline-coloring-page.jpg"=>
            dash_relpath("/assets/images/gallery/flower-outline-coloring-page.jpg"),

            "https://raw.githubusercontent.com/plotly/datasets/master/mitochondria.jpg"=>
            dash_relpath("/assets/images/gallery/mitochondria.jpg")
        )

        for key in find_and_replace
            if occursin(key, _example)
                _example = replace(_example, Regex("# " * config) => find_and_replace[key])
            end
        end

        try
            include_string(app, _example)
        catch e
            print(_example)
            rethrow()
        end

        return(
            string("```python \n", _source, "```")
        )
    end
end

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
module Justitia

import GZip
import MathOptInterface
const MOI = MathOptInterface

include("interface.jl")

include("implementations/csv_record.jl")
include("implementations/milp_result.jl")
include("implementations/moi_model.jl")
include("implementations/moi_backed_instance.jl")
include("implementations/moi_based_approach.jl")
include("implementations/miplib_instance.jl")

include("driver.jl")

end


text:
using Serialization
using Test

using ComparativeAutograder
using ComparativeAutograder: TestSuite, TestCase, launch_runner_proc, readpipe!

const ADD_INT_FUNCTION_SRC = """
function add(x::Int, y::Int)
    x + y
end
"""

@testset "Test subprocess runner" begin
    mktemp() do sub_path, sub_io
        write(sub_io, ADD_INT_FUNCTION_SRC)
        close(sub_io)

        suite = TestSuite([
            TestCase(1, 2),
            TestCase(3, 4),
            TestCase(1.23, 1.45),
        ])

        proc = launch_runner_proc(suite, sub_path, "add"; debug="all")
        if proc.exitcode != 0
            @info "Printing runner stdout."
            println(proc.stdout)
            @info "Printing runner stderr."
            println(proc.stderr)
        end
        @test proc.exitcode == 0
        @test proc.stderr != ""

        result::TestSuiteResult = proc.result.result
        @test result.results[1].result == 3
        @test result.results[2].result == 7
        @test result.results[3].result === nothing
        @test occursin("MethodError", result.results[3].exception)
    end
end