text:
"""
Random search (i.e. Monte Carlo) solver parameters.
"""
@with_kw mutable struct RandomSearchSolver
    n_iterations::Int64 = 100
    episode_length::Int64 = 10
    show_progress::Bool = true
end

"""
Random search planner which takes the `solver` parameters and `mdp` problem.
"""
mutable struct RandomSearchPlanner{P<:Union{MDP,POMDP}}
    solver::RandomSearchSolver
    mdp::P
end

# No work performed.
POMDPs.solve(solver::RandomSearchSolver, mdp::Union{POMDP,MDP}) = RandomSearchPlanner(solver, mdp)

# Online work performed here.
function POMDPs.action(planner::RandomSearchPlanner, s)
    # Pull out variables
    mdp::ASTMDP = planner.mdp
    sim::GrayBox.Simulation = mdp.sim
    n_iterations::Int64 = planner.solver.n_iterations
    depth::Int64 = planner.solver.episode_length
    planner.solver.show_progress ? progress = Progress(n_iterations) : nothing

    # Run simulation
    for i in 1:n_iterations
        BlackBox.initialize!(sim)
        AST.go_to_state(mdp, s)
        AST.rollout(mdp, s, depth) # Evaluate system.
        planner.solver.show_progress ? next!(progress) : nothing
    end

    # Pass back action trace if recording is on (i.e. top_k)
    if mdp.params.top_k > 0
        return get_top_path(mdp)
    end
end

"""
    AST.search!(planner::RandomSearchPlanner)

Search using the `RandomSearchPlanner` from an initial AST state. Pass back the best action trace.
"""
function AST.search!(planner::RandomSearchPlanner)
    mdp::ASTMDP = planner.mdp
    Random.seed!(mdp.params.seed) # Determinism
    s = AST.initialstate(mdp)
    return action(planner, s)
end


text:
Flog(x)=(log(x),nothing)

function Flog_inplace(value,auxvalue,x)
    copy!(value,log(x))
end

function Dlog(derivativeIDX,f_c,faux_c,grad_c,grad_n,x)
    axpy!(1.0,grad_c./f_c,grad_n)
end

if PROC=="GPU"

    function Flog(x::CudaArray)
        tmp=CudaArray(Float64,size(x))
        log!(x,tmp)
        return (tmp,nothing) # memory leak here
    end

    function Flog_inplace(value,auxvalue,x::CudaArray)
        log!(x,value)
    end

    function Dlog(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray)
        vdivupdate!(1.0,grad_c,f_c,grad_n)
    end

end

Derivative[Flog]=Dlog
Inplace[Flog]=Flog_inplace

import Base.log

log(n::ADnode)=ADnode(Flog,n)

export log

ext: jl

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
"""
    RecordArray{T,N}(undef, dims; [align])
    RecordVector{T}(undef, length; [align])

Create an array with array-of-structures memory layout.

Optionally accept alignment `align` of each element (i.e.,
`pointer(A:;RecordArray, i)` is divisible by `align` for all `i in
eachindex(A)`).

See also [`RecordArrays.fill`](@ref) and [`RecordArrays.unsafe_zeros`](@ref).

```jldoctest
julia> using RecordArrays

julia> xs = RecordArray{@NamedTuple{a::Float64, b::Char}}(undef, 3);

julia> xs.a .= 1:3;

julia> xs.b .= 'a':'c';

julia> xs
3-element RecordArray{NamedTuple{(:a, :b), Tuple{Float64, Char}},1,…}:
 (a = 1.0, b = 'a')
 (a = 2.0, b = 'b')
 (a = 3.0, b = 'c')
````
"""
(RecordArrays.RecordArray, RecordArrays.RecordVector)

function RecordArrays.RecordArray{T,N}(
    ::UndefInitializer,
    dims::NTuple{N,Integer};
    align = nothing,
) where {T,N}
    check_eltype(T)
    Storage = wrap_union_type(T)
    if align === nothing
    elseif sizeof(Storage) <= Int(align)
    else
        noinline() do
            error("`align = $align` smaller than `sizeof($Storage) = $(sizeof(Storage))`")
        end
    end
    step = Int(something(align, sizeof(T)))
    buffer = Vector{UInt8}(undef, *(step, dims...) + something(align, 0))
    ptr = Ptr{Storage}(alignto(pointer(buffer), align))
    return RecordArray(buffer, ptr, Val(T), dims, Val(step))
end

RecordArrays.RecordArray{T,N}(
    undef::UndefInitializer,
    dims::Vararg{Integer,N};
    kwargs...,
) where {T,N} = RecordArray{T,N}(undef, dims; kwargs...)
RecordArrays.RecordArray{T}(
    undef::UndefInitializer,
    dims::NTuple{N,Integer};
    kwargs...,
) where {T,N} = RecordArray{T,N}(undef, dims; kwargs...)
RecordArrays.RecordArray{T}(
    undef::UndefInitializer,
    dims::Vararg{Integer,N};
    kwargs...,
) where {T,N} = RecordArray{T,N}(undef, dims; kwargs...)

Base.getproperty(A::RecordArray, name::Symbol) = FieldArray{name}(A)
Base.propertynames(A::RecordArray) = fieldnames(eltype(A))

@noinline function Base.setproperty!(A::RecordArray, name::Symbol, _)
    error(
        "$(ndims(A))-dim `RecordArray` does not support `setproperty!`;",
        " use, e.g., `A.$name .= array` instead",
    )
end

function _view_scalar(A, i)
    ptr = pointer(A, i)
    return RecordArray(_buffer(A), ptr, Val(eltype(A)), (), _step(A))
end

Base.view(A::RecordArray, i::Integer) = _view_scalar(A, i)
Base.view(A::RecordArray{T,N}, I::Vararg{Integer,N}) where {T,N} =
    _view_scalar(A, LinearIndices(A)[I...])

"""
    RecordArrays.fill(x::T, dims::NTuple{N,Integer}) -> A::RecordArray{T,N}

Create a `RecordArray` filled with `x`.
"""
RecordArrays.fill
RecordArrays.fill(x, dims::Integer...; kwargs...) = RecordArrays.fill(x, dims; kwargs...)
function RecordArrays.fill(x, dims::Tuple; kwargs...)
    A = RecordArray{typeof(x)}(undef, dims; kwargs...)
    fill!(A, x)
    return A
end

"""
    RecordArrays.unsafe_zeros(T, dims::NTuple{N,Integer}) -> A::RecordArray{T,N}

Create a `RecordArray` filled with zeros for all elements and all fields.
"""
RecordArrays.unsafe_zeros
RecordArrays.unsafe_zeros(T::Type, dims::Integer...; kwargs...) =
    RecordArrays.unsafe_zeros(T, dims; kwargs...)
function RecordArrays.unsafe_zeros(T::Type, dims::Tuple; kwargs...)
    A = RecordArray{T}(undef, dims; kwargs...)
    fill!(_buffer(A), 0)
    return A
end

function Base.summary(io::IO, A::RecordArray{T,N,Vector{UInt8},Storage}) where {T,N,Storage}
    @nospecialize A
    is_standard = T isa Type
    is_standard &= N isa Int
    is_standard &= if T isa Union
        Storage === UnionValue{T}
    else
        T === Storage
    end
    is_standard || return invoke(summary, Tuple{IO,AbstractArray}, io, A)
    print(io, length(A), "-element ")
    print(io, RecordArray)
    print(io, '{')
    show(io, T)
    print(io, ',')
    show(io, N)
    print(io, ",…")
    print(io, '}')
end


text:
# KRYPTO.JL
module Krypto
using Compat
using SHA
using Primes

########  Generic  ########
export encrypt, decrypt, sign, verify

########  Crypto types  ########
export CryptoAlgorithm, SymmetricCrypto, AsymmetricCrypto
export RSA, RLWE, ECC, NTRU
export AES, Blowfish, Twofish, Threefish, Serpent, Salsa20, RABBIT

########  RSA.jl  ########
export RSAKey, RSAPubKey, RSAPrivKey,                         # RSA Key Objects
       RSAKeyGen,                                             # RSA Key Generation Primitive
       RSAEP, RSADP, RSASP1, RSAVP1,                          # RSA Cryptographic Primitives
       RSAES_PKCS1_V1_5_ENCRYPT, RSAES_PKCS1_V1_5_DECRYPT,    # RSA Encryption Scheme PKCS v1.5 functions (encrypt)
       RSASSA_PKCS1_V1_5_SIGN, RSASSA_PKCS1_V1_5_VERIFY,      # RSA Encryption Scheme PKCS v1.5 functions (sign)
       OS2IP, IP2OS, OS2STR,                                  # Octet string convert functions
       RSAEncrypt, RSADecrypt, RSASign, RSAVerify             # RSA crypto function shorthands

########  ECC.jl  --  Elliptic-Curve Cryptography  ########
export ECCKey, ECCPubKey, ECCPrivKey                      # ECC Key objects

########  RLWE.jl  --  Ring Learning with Errors  ########
export RLWEKey, RLWEPubKey, RLWEPrivKey,                 # RLWE Key Objects
       RLWEKeyGen,                                       # RLWE Key Generation Primitive
       RLWEEncrypt, RLWEDecrypt                          # RLWE Cryptographic Primitives

########  CryptoMath.jl  ########
export bitsize, gpow,           # Determine bitsize of N; get first power of E, greater than N
       csrand, csrandbit,       # CSPRNG functions (via RandomDevice)
       e5,                      # Sieve of Atkin prime generation function
       csprime,                 # 'Strong' cryptographical prime generation function
       phi,                     # Euler's totient (phi) function
       gdb, signbin,            # Additional functions, required by ECC
       gprime, gprimemod        # Finds first prime 1) GEQ to N, 2) GEQ to N, while N % M == R

########  ECCMath.jl  ########
export CurveFP,                 # Curve object
       (+), (++), (*), (~),     # Scalar addition (+), fast addition (++), multiplication (*) and negation (~)
       mulf, doublef,           # Optimized (faster) scalar multiplication and doubling
       muladdf, muladdp,        # Hyper-optimized (fastest) scalar multiplication
       curve_q,                 # Find curve's Q % n, having point (x, y) and parameter p
       in,                      # Check if point lies on curve
       norm2proj, proj2norm     # Transfrom xy-coordinates into projective-coordinates, and vice versa

########  LatticeMath.jl  ########
export KRED, KRED2,              # Reduction functions (for NTT)
       (+), (*),
       LatticeModuli, Int16, Int,
       UniformSample,            # Uniformly sample from [-B:B] an polynomial in ring
       NTTK, INTTK,              # Forward Number Theoretic Transform and its inverse
#      GenerateA, GenerateR2,    # Generator functions for testing and crypto-deployment
       genrandpoly

########  Salsa20.jl  ########
export SALSA20,                         # Salsa20 cryptographic primitive
       Salsa20Stream, Salsa20StreamXOR  # Salsa20 stream expander primitives (plain, xor)

########  HMAC.jl  ########
export HMAC               # Hash-based message authentication code

# Load files
include("LUT/rlwe.jl")
include("Util/CryptoTypes.jl")
include("Util/CryptoMath.jl")
include("Util/ECCMath.jl")
include("Util/LatticeMath.jl")
include("Util/HMAC.jl")
include("Asymmetric/RSA.jl")
include("Asymmetric/ECC.jl")
include("Asymmetric/RLWE.jl")
include("Symmetric/Salsa20.jl")

end # module Krypto


text:
module HumanGEM

    import ...Chemostat_Rath2017
    import ...RathData
    import ..Human1
    import MAT
    import Chemostat
    import Chemostat.MetNets
    using ProjAssistant
    @gen_sub_proj

    include("base_intake_info.jl")
    include("const.jl")
    include("exch_met_map.jl")
    include("Hams_medium.jl")
    include("met_readable_ids.jl")
    include("metnet.jl")
    include("mets_map.jl")
    include("Niklas_biomass.jl")

    function __init__()
        @create_proj_dirs
    end
end


text:
@testset "rationalquad" begin
    rng = MersenneTwister(123456)
    x = rand(rng)*2
    v1 = rand(rng, 3)
    v2 = rand(rng, 3)
    @testset "RationalQuadraticKernel" begin
        α = 2.0
        k = RationalQuadraticKernel(α=α)
        @test RationalQuadraticKernel(alpha=α).α == [α]
        @test kappa(k,x) ≈ (1.0+x/2.0)^-2
        @test k(v1,v2) ≈ (1.0+norm(v1-v2)^2/2.0)^-2
        @test kappa(RationalQuadraticKernel(α=α),x) == kappa(k,x)
        @test metric(RationalQuadraticKernel()) == SqEuclidean()
        @test metric(RationalQuadraticKernel(α=2.0)) == SqEuclidean()
        @test repr(k) == "Rational Quadratic Kernel (α = $(α))"
        test_ADs(x->RationalQuadraticKernel(alpha=x[1]),[α])
    end
    @testset "GammaRationalQuadraticKernel" begin
        k = GammaRationalQuadraticKernel()
        @test kappa(k,x) ≈ (1.0+x^2.0/2.0)^-2
        @test k(v1,v2) ≈ (1.0+norm(v1-v2)^4.0/2.0)^-2
        @test kappa(GammaRationalQuadraticKernel(),x) == kappa(k,x)
        a = 1.0 + rand()
        @test GammaRationalQuadraticKernel(alpha=a).α == [a]
        @test repr(k) == "Gamma Rational Quadratic Kernel (α = 2.0, γ = 2.0)"
        #Coherence test
        @test kappa(GammaRationalQuadraticKernel(α=a, γ=1.0), x) ≈ kappa(RationalQuadraticKernel(α=a), x)
        @test metric(GammaRationalQuadraticKernel()) == SqEuclidean()
        @test metric(GammaRationalQuadraticKernel(γ=2.0)) == SqEuclidean()
        @test metric(GammaRationalQuadraticKernel(γ=2.0, α=3.0)) == SqEuclidean()
        # test_ADs(x->GammaRationalQuadraticKernel(α=x[1], γ=x[2]), [a, 2.0])
        @test_broken "All (problem with power operation)"
    end
end


text:
module MBTR

using ASE
using JuLIP
using PyCall

function __init__()
    py"""
    import qmmlpack as qmml

    import math
    import numpy as np
    import itertools

    from ase import Atoms
    from ase.data import atomic_numbers, chemical_symbols

    PI = 3.14159265359
    unit_BOHR = 0.52917721067

    # DEFINITONS OF QMML-PACK SUPPORT FUNCTIONS
    #
    def basis_check(cell):
        # Checking lattice vectors whether they all are zeros.
        # If x direction is all zero values,
        # we supply 1,0,0 vector instead.
        # This applies all other directions
        return True if( np.linalg.norm(cell[0]) < 1E-12 and
                        np.linalg.norm(cell[1]) < 1E-12 and
                        np.linalg.norm(cell[2]) < 1E-12 ) else False

    def qmml_mbtr_calc(k, z, r, elms, basis=None, flatten=True):
        repr_ = []
        for t in k:
            accuracy = t[2] # accuracy of parameter expansion
            margins = t[1]  # [1/dmax, 1/dmin] plus margin for t[0][0]==2
                            # [0,Pi] plus margin for t[0][0]==3
            defs = t[0] # Parameters for representation
            (dsigma,) = t[3] # theta parameter for this representation
            definitions = (defs[0],
                           defs[1],
                           defs[2],
                           (defs[3][0], (dsigma,)),
                           defs[4],
                           defs[5],
                           defs[6])
            repr_.append(qmml.many_body_tensor(z, r, margins, definitions,
                         basis=basis, acc=accuracy, elems=elms, flatten=flatten))
        return repr_

    def mbtr_reprf(k, atoms, kernelf=None, kernel_theta=None, flatten=True):
        #MBTR, k=2, inverse distances, quadratic weighting
        #      k=3, angles
        if isinstance(atoms, Atoms):
            ase_atoms = [atoms]
            all_z = atoms.numbers
            all_list = all_z
        else:
            ase_atoms = atoms
            all_z = np.asarray([a.numbers for a in atoms])
            all_list = list(itertools.chain.from_iterable(all_z))
        elms = np.unique(all_list)
        z_basis = []
        r_basis = []
        basis = []
        basis_ids = []
        z_no_basis = []
        r_no_basis = []
        no_basis_ids = []
        # Split basis to None and Atoms.cell
        for ai, a in enumerate(ase_atoms):
            if hasattr(a, 'cell'):
                # check if cell only includes zeros
                if not basis_check(a.cell):
                    z_basis.append(a.numbers)
                    r_basis.append(a.positions)
                    basis.append(a.cell)
                    basis_ids.append(ai)
                else:
                    z_no_basis.append(a.numbers)
                    r_no_basis.append(a.positions)
                    no_basis_ids.append(ai)
            else:
                z_no_basis.append(a.numbers)
                r_no_basis.append(a.positions)
                no_basis_ids.append(ai)

        # representation
        repr_ = []
        all_ids = {}
        if z_basis:
            for i, bi in enumerate(basis_ids):
                all_ids[bi] = i
            repr_.extend(qmml_mbtr_calc(k, z_basis, r_basis,
                                        elms, basis=basis,
                                        flatten=flatten))
        if z_no_basis:
            for i, bi in enumerate(no_basis_ids):
                all_ids[bi] = i + len(basis_ids)
            repr_.extend(qmml_mbtr_calc(k, z_no_basis, r_no_basis,
                                        elms, basis=None,
                                        flatten=flatten))
        repr_ = np.concatenate(repr_, axis=1)
        M = None
        if kernelf is not None:
            # kernel matrix
            M = kernelf(repr_, theta=kernel_theta)
        return repr_, M, all_ids

    def paramf(aa, bb, cc):
        # Parameter splitting function dividing
        # parameters into three blocks of sizes a, b, c.
        #
        # a, b, c are numbers of parameters.
        # a can also be a list of numbers of parameters.
        #
        # paramf([1,2], 0, 1)(['a', 'b', 'c', 'd'])
        #    ->
        # ([('a',), ('b', 'c')], tuple(), ('d',))
        def f(a, b, c, *theta):
            if qmml.is_sequence(a):
                # multiple representation functions
                # as we have only need for cases 2
                # and 3, no need to generalize
                assert len(a) == 2 or len(a) == 3
                resa = [ tuple(theta[:a[0]]),
                         tuple(theta[a[0]:a[0]+a[1]]) ]
                if len(a) == 3: resa.append(
                                tuple(theta[a[0]+a[1]:a[0]+a[1]+a[2]]) )
                a = np.sum(a)
            else:
                # single representation function
                resa = tuple(theta[:a])
            return ( resa, tuple(theta[a:a+b]), tuple(theta[a+b:a+b+c]) )
        return lambda *theta: f(aa, bb, cc, *theta)

    def calc_desc(atoms_list, desc, reprfs, kfunc=None,
                  theta_kernel=None, flatten=True):
        out_desc = []
        # Split dataset into non-PBC and PBC
        # structures and calculate descriptors
        pbc_structs = []
        nonpbc_structs = []
        for aseat in atoms_list:
            if (aseat.pbc[0] == False and
                aseat.pbc[1] == False and
                aseat.pbc[2] == False):
                nonpbc_structs.append(aseat)
            else:
                pbc_structs.append(aseat)
        if pbc_structs:
            # Calculate descriptor for structures with pbc
            rtn_desc, kern, ids = desc(reprfs, pbc_structs, kernelf=kfunc,
                                       kernel_theta=theta_kernel,
                                       flatten=flatten)
            for count, ase_at in enumerate(pbc_structs):
                # Setting database
                out_desc.append({
                    'representation' : rtn_desc[ids[count]],
                    'kernel'         : kern[ids[count]],
                    'symbols'        : ase_at.get_chemical_symbols()
                    })
        if nonpbc_structs:
            # Calculate descriptor for structures with pbc
            rtn_desc, kern, ids = desc(reprfs, nonpbc_structs, kernelf=kfunc,
                                       kernel_theta=theta_kernel,
                                       flatten=flatten)
            for count, ase_at in enumerate(nonpbc_structs):
                # Setting database
                out_desc.append({
                    'representation' : rtn_desc[ids[count]],
                    'kernel'         : kern[ids[count]],
                    'symbols'        : ase_at.get_chemical_symbols()
                    })
        return out_desc

    def MBTR_desc(at, reprs=2, feature='r', kernelf='gaussian',
                  accuracy=0.001, dist_min=0.1, dist_max=1.1,
                  dist_step=100, ang_min=-0.15, ang_max=3.45575,
                  ang_step=100, dist_width=-4.0, ang_width=-3.0,
                  sigma=0.0, stepsize=0.5, reg_strength=-20.0,
                  reg_direction=-1, cutoff=6.0):
        # at            : ASE Atoms
        # reprs         : 0=only distances, 1=only angles, 2=both.
        # feature       : kernel='k' or representation='r'
        # kernelf       : 'gaussian', 'linear', or 'laplacian'.
        # accuracy      : Accuracy for MBTR
        # dist_min      : Minimum distance
        # dist_max      : Maximum distance
        # dist_step     : Number of distance steps
        # ang_min       : Minimum for angle margins
        # ang_max       : Maximum for angle margins
        # ang_step      : Number of angle margin step
        # dist_width    : 2-body MBTR normal distribution width.
        # ang_width     : 3-body MBTR normal distribution width.
        # sigma         : kernel basis sigma.
        # stepsize      : Stepsize for MBTR variables.
        # reg_strength  : Regularization strength.
        # reg_direction : Regularization direction.
        # cutoff        : Reserved if needed. Currently, not used!

        if 'laplacian' in kernelf:
            kfunc = qmml.kernel_laplacian
        elif 'linear' in kernelf:
            kfunc = qmml.kernel_linear
        else:
            kfunc = qmml.kernel_gaussian

        dist_repr = [
            (2, '1/distance', 'identity^2',
                ('normal', 0.), 'identity',
                'noreversals', 'noreversals'
            ),
            (float(dist_min),
             float(dist_max/int(dist_step)),
             int(dist_step)),
            float(accuracy),
            0.
            ]
        ang_repr = [
            (3, 'angle', '1/dotdotdot',
                ('normal', 0.), 'identity',
                'noreversals', 'noreversals'
            ),
            (float(ang_min),
             float(ang_max/int(ang_step)),
             int(ang_step)),
            float(accuracy),
            0.
            ]
        reprfs = []
        if reprs > 1:
            reprfs.append(dist_repr)
        if reprs > 0:
            reprfs.append(ang_repr)
        if reprs < 1:
            reprfs.append(dist_repr)

        variables = (
            { 'value': dist_width, 'priority': 2, 'stepsize': stepsize },
            { 'value': ang_width, 'priority': 1, 'stepsize': stepsize },
            { 'value': sigma, 'priority': 3, 'stepsize': stepsize },
            { 'value': reg_strength, 'priority': 4,
              'direction': reg_direction, 'minimum': reg_strength }
            )
        parmf = paramf([1,1],1,1)
        theta = (dist_width, ang_width, sigma, reg_strength)
        (theta_repr, theta_kernel, theta_method) = parmf(
                                                    *(2.**np.asarray(theta)))
        for ki, k in enumerate(reprfs):
            reprfs[ki][3] = theta_repr[ki]

        fingerprints = calc_desc([at], mbtr_reprf,
                                 reprfs, kfunc,
                                 theta_kernel)
        if feature.startswith('r'):
            fs = np.array(fingerprints[0]['representation'])
        else:
            fs = np.array(fingerprints[0]['kernel'])
        return fs
    """
end

export mbtr

# Please check parameter definitions in python function above.
function mbtr(at; reprs=2, feature="r", kernelf="gaussian",
              accuracy=0.001, dist_min=0.1, dist_max=1.1,
              dist_step=100, ang_min=-0.15, ang_max=3.45575,
              ang_step=100, dist_width=-4.0, ang_width=-3.0,
              sigma=0.0, stepsize=0.5, reg_strength=-20.0,
              reg_direction=-1, cutoff=6.0)
    atom_struct = ASEAtoms(at)
    # Calculate descriptor
    mbtr_rtn = py"MBTR_desc"(atom_struct.po, reprs=reprs,
                             feature=feature, kernelf=kernelf,
                             accuracy=accuracy, dist_min=dist_min,
                             dist_max=dist_max, dist_step=dist_step,
                             ang_min=ang_min, ang_max=ang_max,
                             ang_step=ang_step, dist_width=dist_width,
                             ang_width=ang_width, sigma=sigma,
                             stepsize=stepsize, reg_strength=reg_strength,
                             reg_direction=reg_direction, cutoff=cutoff)
    return mbtr_rtn
end

end # module


text:
#EasyPlot base type, constant & function definitions
#-------------------------------------------------------------------------------

#==Useful constants
===============================================================================#

#=Not used! Informative at the moment.

#const VALID_AXISSCALES = Set([:lin, :log, :dB10, :dB20, :reciprocal])
#const VALID_LINESTYLES = Set([:none, :solid, :dash, :dot, :dashdot])

#In case specified glyph (symbol/marker) not supported...
#implementations should provide default (hopefully not none).
const VALID_GLYPHS = Set([:none,
	:square, :diamond,
	:uarrow, :darrow, :larrow, :rarrow, #usually triangles
	:cross, :+, :diagcross, :x,
	:circle, :o, :star, :*,
])
=#

#==Main types
===============================================================================#
#A plot theme.
mutable struct Theme
	colorscheme::Optional{ColorScheme}

#=Under consideration
	axisline::LineAttributes
	vgridline::LineAttributes
	hgridline::LineAttributes
	dfltline::LineAttributes
	dfltglyph::LineAttributes
	dfltglyphfill
	widthscheme #ex: [1, 3, 5, 7, 9]
=#
end
Theme() = Theme(nothing)

#w & h for plot dimensions (pixels):
abstract type AbstractPlotDimensions; end
struct PlotAutosize <: AbstractPlotDimensions; end
const plotautosize = PlotAutosize()
struct PlotDim <: AbstractPlotDimensions
	w::Int; h::Int
end

abstract type AbstractAxis; end

struct Axis{ID} <: AbstractAxis
	#TODO: Ensure valid id values
end
Axis(id::Symbol) = Axis{id}()

"""`struct FoldedAxis`

Describe how a folded axis is to be displayed.
"""
struct FoldedAxis <: AbstractAxis
	foldinterval::PReal #Interval with which to start overlaying data
	xstart::PReal #Discard data before this point
	xmin::PReal
	xmax::PReal
end
FoldedAxis(foldinterval; xstart=0, xmin=0, xmax=foldinterval) =
	FoldedAxis(foldinterval, xstart, xmin, xmax)

mutable struct Waveform <: AbstractAttributeReceiver
	data::DataMD
	label::String
	line::LineAttributes
	glyph::GlyphAttributes
	strip::Int
end
_Waveform(data::DataMD, label::String, strip::Int) =
	Waveform(data, label, LineAttributes(), GlyphAttributes(), strip)
Waveform(data::DataMD, args...; label::String="", strip=1, kwargs...) =
	_apply(_Waveform(data, label, strip), args, kwargs)

mutable struct YStrip
	scale::Symbol
	striplabel::String
	axislabel::String
	ext::Extents1D
	grid::AbstractGrid
end
_YStrip() = YStrip(:lin, "", "", _Extents1D(), GridCartesian())

mutable struct Plot <: AbstractAttributeReceiver
	title::String
	xlabel::String
	legend::Bool
	xaxis::AbstractAxis
	xext::Extents1D
	ystriplist::Vector{YStrip}
	wfrmlist::Vector{Waveform}
	annot::Vector{AbstractPlotAnnotation}
end
_Plot(title::String, legend::Bool) =
	Plot(title, "", legend, Axis(:lin), _Extents1D(), YStrip[_YStrip()], Waveform[], AbstractPlotAnnotation[])
Plot(args...; title::String="", legend=true, kwargs...) =
	_apply(_Plot(title, legend), args, kwargs)

mutable struct PlotCollection <: AbstractAttributeReceiver
	title::String
	ncolumns::Int #TODO: Create a more flexible system
	plotlist::Vector{Plot}
	bblist::NullOr{Vector{BoundingBox}} #User can specify where to draw plots
	theme::Theme
	opt::Any #Optional data/instructions to relay to builder
end
_PlotCollection(title::String, ncolumns::Int) =
	PlotCollection(title, ncolumns, Plot[], nothing, Theme(), nothing)
PlotCollection(args...; title::String="", ncolumns::Int=1, kwargs...) =
	_apply(_PlotCollection(title, ncolumns), args, kwargs)

#==Register constructors with cons() interface
===============================================================================#
cons(::DS{:plot_collection}, args...; kwargs...) = PlotCollection(args...; kwargs...)
cons(::DS{:plotcoll}, args...; kwargs...) = PlotCollection(args...; kwargs...)
cons(::DS{:plot}, args...; kwargs...) = Plot(args...; kwargs...)
cons(::DS{:wfrm}, args...; kwargs...) = Waveform(args...; kwargs...)
cons(::DS{:fldaxis}, args...; kwargs...) = FoldedAxis(args...; kwargs...)

#==Accessors
===============================================================================#
Base.Symbol(a::Axis{T}) where T = Symbol(T)
Base.Symbol(a::FoldedAxis) = :lin
isauto(d::AbstractPlotDimensions) = false
isauto(d::PlotAutosize) = true

#==Copy constructors
===============================================================================#
Base.copy(src::Waveform) =
	Waveform(src.data, src.label, deepcopy(src.line), deepcopy(src.glyph), src.strip)

#==Helper functions
===============================================================================#
#Create a new Extents1D with conditionally overwritten min/max values.
function _new_overwrite(srcext::Extents1D, _min, _max)
	newmin = srcext.min; newmax = srcext.max
	NoOverwrite(_min) || (newmin = _min)
	NoOverwrite(_max) || (newmax = _max)
	return Extents1D(newmin, newmax)
end

function _get_strip(p::Plot, istrip::Int)
	nstrips = length(p.ystriplist)
	if istrip < 1 || istrip > nstrips
		throw(ArgumentError("set(`::Plot`, ystrip=...): strip index=$istrip, but strip count=$nstrips"))
	end
	return p.ystriplist[istrip]
end

_setgrid(::DS, p::Plot, istrip::Int; kwargs...) =
	throw(ArgumentError("set(`::Plot`, grid=set(fmt=...), ...): fmt must be one of: {:cartesian}"))

function _setgrid(::DS{:cartesian}, p::Plot, istrip::Int; kwargs...)
	strip = _get_strip(p, istrip)
	strip.grid = GridCartesian(; kwargs...)
	return p
end

#==Methods to resolve user-supplied values
===============================================================================#
_resolve_xaxis(xscale::Symbol) = Axis(xscale)
_resolve_xaxis(a::FoldedAxis) = a

#_resolve_colorant(c::Symbol) = getcolor(c)
#_resolve_colorant(c::Colorant) = c

#Resolve to appropriate ColorRef value:
_resolve_ColorRef(c::Nothing) = c
_resolve_ColorRef(c::Symbol) = getcolor(c)
_resolve_ColorRef(c::Int) = c
_resolve_ColorRef(c::Colorant) = c

#==Implement _apply() methods
===============================================================================#
function _apply(w::Waveform, ::DS{:line}; style=nooverwrite, width=nooverwrite, color=nooverwrite)
	_style = w.line.style; _width = w.line.width; _color = w.line.color
	NoOverwrite(style) || (_style = style)
	NoOverwrite(width) || (_width = width)
	NoOverwrite(color) || (_color = _resolve_ColorRef(color))
	w.line = LineAttributes(style=_style, width=_width, color=_color)
	return w
end

function _apply(w::Waveform, ::DS{:glyph}; shape=nooverwrite, size=nooverwrite,
	color=nooverwrite, fillcolor=nooverwrite)
	g=w.glyph; _shape=g.shape; _size=g.size; _color=g.color; _fillcolor=g.fillcolor
	NoOverwrite(shape) || (_shape = shape)
	NoOverwrite(size) || (_size = size)
	NoOverwrite(color) || (_color = _resolve_ColorRef(color))
	NoOverwrite(fillcolor) || (_fillcolor = _resolve_ColorRef(fillcolor))
	w.glyph = GlyphAttributes(shape=_shape, size=_size, color=_color, fillcolor=_fillcolor)
	return w
end

#Configure Plot with x/y axes (single y-strip) and change scales:
function _apply(p::Plot, ::DS{:xyaxes}; xscale=nooverwrite, yscale=nooverwrite,
	xmin=nooverwrite, xmax=nooverwrite, ymin=nooverwrite, ymax=nooverwrite)
	if length(p.ystriplist) > 0
		resize!(p.ystriplist, 1)
	else
		push!(p.ystriplist, _YStrip())
	end
	strip = p.ystriplist[1]
	NoOverwrite(xscale) || (p.xaxis = _resolve_xaxis(xscale))
	NoOverwrite(yscale) || (strip.scale = yscale)
	p.xext = _new_overwrite(p.xext, xmin, xmax)
	strip.ext = _new_overwrite(strip.ext, ymin, ymax)
	return p
end

#Change plot labels:
function _apply(p::Plot, ::DS{:labels}; title=nooverwrite,
	xaxis=nooverwrite, yaxis=nooverwrite)
	if !NoOverwrite(yaxis) && length(p.ystriplist) < 1
		push!(p.ystriplist, _YStrip())
	end
	NoOverwrite(title) || (p.title = title)
	NoOverwrite(xaxis) || (p.xlabel = xaxis)
	NoOverwrite(yaxis) || (p.ystriplist[1].axislabel = yaxis)
	return p
end

#Change x-axis parameters of a plot:
function _apply(p::Plot, ::DS{:xaxis}; scale=nooverwrite,
	min=nooverwrite, max=nooverwrite, label=nooverwrite)
	NoOverwrite(scale) || (p.xaxis = _resolve_xaxis(scale))
	NoOverwrite(label) || (p.xlabel = label)
	p.xext = _new_overwrite(p.xext, min, max)
	return p
end

#Change number of y-strips in a plot:
function _apply(p::Plot, ::DS{:nstrips}, nstrips::Int)
	if length(p.ystriplist) > nstrips
		resize!(p.ystriplist, 1)
	else
		newstrips = nstrips - length(p.ystriplist)
		for i in 1:newstrips
			push!(p.ystriplist, _YStrip())
		end
	end
	return p
end

#Change misc. strip properties:
function _apply(p::Plot, ::DS{:ystrip}, istrip::Int;
	min=nooverwrite, max=nooverwrite, scale=nooverwrite,
	axislabel=nooverwrite, striplabel=nooverwrite
)
	strip = _get_strip(p, istrip)
	strip.ext = _new_overwrite(strip.ext, min, max)
	NoOverwrite(scale) || (strip.scale = scale)
	NoOverwrite(axislabel) || (strip.axislabel = axislabel)
	NoOverwrite(striplabel) || (strip.striplabel = striplabel)
	return p
end
#Strip aliases to avoid parameter name collisions (can only set ystrip=set(...) once):
_apply(p::Plot, ::DS{:ystrip1}; kwargs...) = _apply(p, DS(:ystrip), 1; kwargs...)
_apply(p::Plot, ::DS{:ystrip2}; kwargs...) = _apply(p, DS(:ystrip), 2; kwargs...)
_apply(p::Plot, ::DS{:ystrip3}; kwargs...) = _apply(p, DS(:ystrip), 3; kwargs...)
_apply(p::Plot, ::DS{:ystrip4}; kwargs...) = _apply(p, DS(:ystrip), 4; kwargs...)
_apply(p::Plot, ::DS{:ystrip5}; kwargs...) = _apply(p, DS(:ystrip), 5; kwargs...)
_apply(p::Plot, ::DS{:ystrip6}; kwargs...) = _apply(p, DS(:ystrip), 6; kwargs...)
_apply(p::Plot, ::DS{:ystrip7}; kwargs...) = _apply(p, DS(:ystrip), 7; kwargs...)
_apply(p::Plot, ::DS{:ystrip8}; kwargs...) = _apply(p, DS(:ystrip), 8; kwargs...)
_apply(p::Plot, ::DS{:ystrip9}; kwargs...) = _apply(p, DS(:ystrip), 9; kwargs...)

function _apply(p::Plot, ::DS{:xfolded}, foldinterval::PReal;
	xstart=0.0, xmin=0.0, xmax=foldinterval
)
	p.xaxis = FoldedAxis(foldinterval, xstart=xstart, xmin=xmin, xmax=xmax)
	return p
end

function _apply(g::GridCartesian; vmajor=nooverwrite, vminor=nooverwrite,
	hmajor=nooverwrite, hminor=nooverwrite
)
	NoOverwrite(vmajor) || (g.vmajor = vmajor)
	NoOverwrite(vminor) || (g.vminor = vminor)
	NoOverwrite(hmajor) || (g.hmajor = hmajor)
	NoOverwrite(hminor) || (g.hminor = hminor)
	return g
end

function _apply(p::Plot, ::DS{:grid}, istrip::Int=1; fmt=nooverwrite, kwargs...)
	if NoOverwrite(fmt)
		strip = _get_strip(p, istrip)
		_apply(strip.grid; kwargs...)
	elseif typeof(fmt) != Symbol
		_setgrid(DS(:throwerror), p, istrip)
	else
		_setgrid(DS(:cartesian), p, istrip; kwargs...)
	end
	return p
end

#Grid aliases to avoid parameter name collisions (can only set ystrip=set(...) once):
_apply(p::Plot, ::DS{:grid1}; kwargs...) = _apply(p, DS(:grid), 1; kwargs...)
_apply(p::Plot, ::DS{:grid2}; kwargs...) = _apply(p, DS(:grid), 2; kwargs...)
_apply(p::Plot, ::DS{:grid3}; kwargs...) = _apply(p, DS(:grid), 3; kwargs...)
_apply(p::Plot, ::DS{:grid4}; kwargs...) = _apply(p, DS(:grid), 4; kwargs...)
_apply(p::Plot, ::DS{:grid5}; kwargs...) = _apply(p, DS(:grid), 5; kwargs...)
_apply(p::Plot, ::DS{:grid6}; kwargs...) = _apply(p, DS(:grid), 6; kwargs...)
_apply(p::Plot, ::DS{:grid7}; kwargs...) = _apply(p, DS(:grid), 7; kwargs...)
_apply(p::Plot, ::DS{:grid8}; kwargs...) = _apply(p, DS(:grid), 8; kwargs...)
_apply(p::Plot, ::DS{:grid9}; kwargs...) = _apply(p, DS(:grid), 9; kwargs...)

#==push! interface
===============================================================================#
function Base.push!(c::PlotCollection, p::Plot, args...)
	push!(c.plotlist, p)
	if length(args)>0; push!(c, args...); end
	return c
end
function Base.push!(p::Plot, w::Waveform, args...)
	push!(p.wfrmlist, w)
	if length(args)>0; push!(p, args...); end
	return p
end
function Base.push!(p::Plot, a::AbstractPlotAnnotation, args...)
	push!(p.annot, a)
	if length(args)>0; push!(p, args...); end
	return p
end

#Last line


text:
function jsonify_inputs(inputs, is_value)
    out = Array{String, 1}()
    for input in inputs
        if !ismissing(is_value) && input == is_value
            push!(out, "'$input': value")
        else
            push!(out, "'$input' : this.$(input)")
        end
    end
    join(out, ", ")
end

vue_models(m::UIModel) = "$(m.id): $(m.default)"

function matte_builtin_models()
"""
            matte_notconnected_overlay: true,
            error_snackbar: false,
            matte_error_msg: "",
            session_id: uuidv1()"""
end

function vue_models(content::NTuple{N, UIModel}) where N
    join([matte_builtin_models(), [vue_models(x) for x in content]...], ",\n")
end

function fetch_update_method(input_id, reverse_dependencies, dependency_tree)
    out = "
    fetch_update_$input_id: function(value) {\n"
    for output in reverse_dependencies
        out *= "        Matte.request_update(\"$output\", {$(jsonify_inputs(dependency_tree[output], input_id))})\n"
    end
    out *= "    }"

    out
end

function watch_function(input_id, content)
    "
    $input_id: function(value) {
        $content
    }"
end

function watch_functions(rev_dep, watch_extras)
    # Define the set of watch functions for inputs that affected outputs
    watch_dict = Dict{String, String}(input_id => "this.fetch_update_$input_id(value)" for input_id in keys(rev_dep))

    # define manually added watch statements
    for extra in watch_extras
        if haskey(watch_dict, extra.id)
            watch_dict[extra.id] *= "\n" * extra.code
        else
            watch_dict[extra.id] = extra.code
        end
    end

    join([watch_function(id, code) for (id, code) in watch_dict], ",\n")
end

function on_connected(dep_tree)
    mount_calls = join(["Matte.request_update(\"$output_id\", {$(jsonify_inputs(dependencies, missing))})" for (output_id, dependencies) in dep_tree], ";\n")
    """
    on_connected: function() {
        this.matte_notconnected_overlay = false;
        $mount_calls
    }"""
end

function update_model()
    """
    update_model: function(response) {
        if (Object.prototype.hasOwnProperty.call(response, "matte_error_msg") && !(response["matte_error_msg"] === null)) {
            this.matte_error_msg = response["matte_error_msg"]
            this.error_snackbar = true
        } else if (!(response['value'] === null)) {
            this[response['id']] = response['value']
        }
    }"""
end

function vue_methods(models, rev_dep, dep_tree)
    for model in models
        if !haskey(rev_dep, model.id)
            rev_dep[model.id] = []
        end
    end
    methods = [fetch_update_method(input_id, outputs_affected, dep_tree) for (input_id, outputs_affected) in rev_dep]
    pushfirst!(methods, on_connected(dep_tree))
    pushfirst!(methods, update_model())
    join(methods, ",\n")
end

function vue_watch(rev_dep, watch)
    watch_functions(rev_dep, watch)
end

function vue_js(rev_dep, dep_tree, models, watch)
    """
    new Vue({
      el: '#app',
      vuetify: new Vuetify(),
      data: {
        $(vue_models(models))
      },
      methods: {
        $(vue_methods(models, rev_dep, dep_tree))
      },
      watch: {
        $(vue_watch(rev_dep, watch))
      },
      mounted: function() {
        Matte.load_channels(this);
      }
    })
    """
end


text:
module MacroDatabase
using ..Fxcrd.ShellCompiler
using ..Fxcrd.ShellInterop
export Macro, Macros, get_macro, macrocall

struct Macro
  name::String
  func::Function
  eval::Bool
  meta::Bool
end
const Macros = Dict{String,Macro}()

function get_macro(name)
  @assert haskey(Macros, name) "macro $name not defined"
  return Macros[name]
end

function macrocall(m::Macro, args::String...; input::String="")::String
  if m.eval
    shell() do sh
      args = map(args) do arg
        read_sh(sh, "echo $arg")
      end
    end
  end

  res = m.func(args..., input=input)
  if m.meta
    shell() do sh
      res = read_sh(sh, res)
    end
  end

  return res
end

end

max_stars_repo_name: fcard/Fxcrd

max_issues_repo_name: fcard/Fxcrd

max_forks_repo_name: fcard/Fxcrd


text:
using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--folder_name"
            help = "The name of the folder containing the images and other data"
            arg_type = String
            default = "labrun5"
        "--previous"
            help = "Compute within a previous results folder."
            arg_type = String
            default = ""
        "--jldfile"
            help = "Which JLD file to load from the previous results folder"
            arg_type = String
            default = ""
        "--failsafe"
            help = "Single process and recursive Bayes tree solving, slower by but more robust with better error stack traces."
            action = :store_true
        "--cxoffset"
            help = "Center x offset"
            arg_type = Float64
            default = 0.0
        "--cyoffset"
            help = "Center y offset"
            arg_type = Float64
            default = 0.0
        "--focalscale"
            help = "Focal distance scaling"
            arg_type = Float64
            default = 1.0
        "--tagscale"
            help = "Tag sized scaling"
            arg_type = Float64
            default = 1.0
        "--pose_trigger_distance"
            help = "Required distance travelled to trigger new pose."
            arg_type = Float64
            default = 0.4
        "--pose_trigger_rotate"
            help = "Required rotation to trigger new pose."
            arg_type = Float64
            default = pi/4
        "--naive_frac"
            help = "Fraction of particles to use naive odo model."
            arg_type = Float64
            default = 0.6
        "--ADAM_step"
            help = "ADAM training step size (might be annealed locally during training)."
            arg_type = Float64
            default = 0.1
        "--batch_size"
            help = "Batch solve every n poses"
            arg_type = Int64
            default = 20
        "--epochsFlux"
            help = "How many epochs to induce during Flux training."
            arg_type = Int64
            default = 15
        "--fluxGenerations"
            help = "How many times to repeat the Flux training cycles."
            arg_type = Int64
            default = 5
        "--trainingNumOffset"
            help = "Offset training folder number."
            arg_type = Int64
            default = 0
        "--rndSkip"
            help = "Flux training number of interpose starting points to skip; can/should be used in combination with rndChords.  -1 is randomized"
            arg_type = Int64
            default = -1
        "--rndChord"
            help = "Flux training number of interpose chords.  Int[] is randomized"
            arg_type = String
            default = "Int[]"
        "--loadInitModels"
            help = "Load initial model weights from factor graph object"
            arg_type = String
            default = ""
        "--fgpathsflux"
            help = "Point to a text file containing which results (savedfg.tar.gz destinations) to use for training."
            arg_type = String
            default = ""
        "--numFGDatasets"
            help = "How many fgpathsflux datasets to include in FG for Flux training; default -1 is use all available."
            arg_type = Int64
            default = -1
        "--localprocs"
            help = "How many local Julia processes to use"
            arg_type = Int64
            default = 4
        "--remoteprocs"
            help = "How many remote Julia processes to use"
            arg_type = Int64
            default = 0
        "--remoteserver"
            help = "Which remote server to use for processing"
            arg_type = String
            default = "JL_CLUSTER_HY"
        "--iterposes"
            help = "maximum number of poses to process"
            arg_type = Int64
            default = Int64(2)^63-Int64(1)
        "--msgloops"
            help = "Maximum number or messages to consume in the main loop (ROS/LCM/etc)"
            arg_type = Int64
            default = Int64(2)^63-Int64(1)
        "--show"
            help = "Show pdfs during computation using `evince`"
            action = :store_true
        "--vis2d"
            help = "Draw 2D visualizations"
            action = :store_true
        "--vis3d"
            help = "Draw 3D visualizations"
            action = :store_true
        "--imshow"
            help = "Show camera visualization"
            action = :store_true
        "--batch_resolve"
            help = "Resolve the entire factor graph at the end"
            action = :store_true
        "--report_factors"
            help = "Resolve the entire factor graph at the end"
            action = :store_true
        "--savedfg"
            help = "Store DFG objects to file"
            action = :store_true
        "--dbg"
            help = "Enable debugging flags"
            action = :store_true
        "--usesimmodels"
            help = "Use sim trained models"
            action = :store_true
    end

    return parse_args(s)
end

parsed_args = parse_commandline()

# convert from strings to types
parsed_args["rndChord"] = eval(Meta.parse(parsed_args["rndChord"]))


text:
module Reference

include("../compat.jl")

using OpenGene

import DataStructures:OrderedDict

include("Genome/Genome.jl")
include("Gencode/Gencode.jl")

end


text: error("oops")


text:
include("../Resistances/heat.jl")

using Plots, Measurements

# Boundary model(s)
T_a = 95.6
T_b = 273.15

@named boundaryT_a =Temperature(T=T_a)
@named boundaryT_b =Temperature(T=T_b)

@variables Q_total(t)
D = Differential(t)

# Sphere model
r_inner = 0.1
th_1 = 0.0025
th_ins = 0.01
th_2 = 0.0025

r_ins_inner = r_inner + th_1
r_ins_outer = r_ins_inner + th_ins
r_outer = r_ins_outer + th_2

@named convectionInner = Convection(surfaceArea=4*3.14159*r_inner^2,alpha=150)
@named linerInner = Sphere(r_inner=r_inner,r_outer=r_ins_inner,lambda=15)
@named contact_1 = Contact(surfaceArea=4*3.14159*r_ins_inner^2,Rc_pp=0.003)
@named insulation = Sphere(r_inner=r_ins_inner,r_outer=r_ins_outer,lambda=0.033)
@named contact_2 = Contact(surfaceArea=4*3.14159*r_ins_outer^2,Rc_pp=0.003)
@named linerOuter = Sphere(r_inner=r_ins_outer,r_outer=r_outer,lambda=15)
@named radiationOuter = Radiation(surfaceArea=4*3.14159*r_outer^2,epsilon=0.7)
@named convectionOuter = Convection(surfaceArea=4*3.14159*r_outer^2,alpha=6.0)

eqs =
[
    connect(convectionInner.port_a,boundaryT_a.port)
    connect(convectionInner.port_b,linerInner.port_a)
    connect(linerInner.port_b,contact_1.port_a)
    connect(contact_1.port_b,insulation.port_a)
    connect(insulation.port_b,contact_2.port_a)
    connect(contact_2.port_b,linerOuter.port_a)
    connect(linerOuter.port_b,radiationOuter.port_a,convectionOuter.port_a)
    connect(convectionOuter.port_b, radiationOuter.port_b,boundaryT_b.port)
    D(Q_total) ~ convectionInner.port_a.Q_flow
]

@named model = ODESystem(eqs, t, [Q_total],[], systems=[boundaryT_a,boundaryT_b,convectionInner,linerInner,contact_1,insulation,contact_2,linerOuter,radiationOuter,convectionOuter])

sys = structural_simplify(model)
u0 = [Q_total=>0.0]
prob = ODAEProblem(sys, u0, (0, 1.0), [T_a, T_b, r_inner, th_1, r_ins_inner, th_ins, r_ins_outer, th_2, r_outer])
sol = solve(prob, Tsit5())
y_min, t_min, y_max, t_max = reportMinMax("Q_total","Thermal Resistances Model",sol)
errorCheck(-69.43455, y_min)
plot(sol)


text:
import HoTRG_lemon.SpinModule: isClassical
#=
LatticeInfo(latticecode::AbstractString)
getDimension(latticeinfo::LatticeInfo)
getGeometry(latticeinfo::LatticeInfo)
getCoarserate(latticeinfo::LatticeInfo)
getQuantumOrClassical(latticeinfo::LatticeInfo)
isClassical(latticeinfo::LatticeInfo)
getDimension(latticeinfo::LatticeInfo, dimension::Int)
isTwoDimension(latticeinfo::LatticeInfo)
isSquareLattice(latticeinfo::LatticeInfo)
=#

type LatticeInfo
	quantumOrClassical::AbstractString
	dimension::Int
	geometry::AbstractString
	coarserate::Int
	legextension::Int

	function LatticeInfo(latticecode::AbstractString)
		this = new()
		initiateLattice!(latticecode, this)
		return this
	end
end

#---
# auxiliary function for constructor

function initiateLattice!(latticecode::AbstractString, latticeinfo::LatticeInfo)
	lattice = split(latticecode,"_")
	latticeinfo.quantumOrClassical = lowercase(lattice[1])
	latticeinfo.dimension = parse(Int,match(r"(?<dimension>\d+)d",lattice[2])[:dimension])
	latticeinfo.geometry = lowercase(lattice[3])
	if getGeometry(latticeinfo) == "square"
		latticeinfo.coarserate = 2
		latticeinfo.legextension = 0
	elseif getGeometry(latticeinfo) == "fractal"
		latticeinfo.legextension =
			try
				parse(Int,lattice[4])
			catch
				1
			end
	end
	return latticeinfo
end

#---
# functions for latticeinfo
function getDimension(latticeinfo::LatticeInfo)
	return latticeinfo.dimension
end

function getGeometry(latticeinfo::LatticeInfo)
	return latticeinfo.geometry
end

function getLegextension(latticeinfo::LatticeInfo)
	return latticeinfo.legextension
end

function getCoarserate(latticeinfo::LatticeInfo)
	return latticeinfo.coarserate
end

function getQuantumOrClassical(latticeinfo::LatticeInfo)
	return latticeinfo.quantumOrClassical
end

function isClassical(latticeinfo::LatticeInfo)
	if getQuantumOrClassical(latticeinfo) == "classical"
		return true
	else
		return false
	end
end

function getDimension(latticeinfo::LatticeInfo, dimension::Int)
	if getDimension(latticeinfo) == dimension
		return true
	else
		return false
	end
end

function isTwoDimension(latticeinfo::LatticeInfo)
	return getDimension(latticeinfo, 2)
end

function isSquareLattice(latticeinfo::LatticeInfo)
	if getGeometry(latticeinfo) == "square"
		return true
	else
		return false
	end
end

function isFractalLattice(latticeinfo::LatticeInfo)
	if getGeometry(latticeinfo) == "fractal"
		return true
	else
		return false
	end
end
#---
# test function

function testlatticeinfo()
	latticecode = "quantum_2d_square"
	latticeinfo = LatticeInfo(latticecode)
	println(isTwoDimension(latticeinfo))
	println(isSquareLattice(latticeinfo))
	println(getDimension(latticeinfo,4))
end


text:
"""
`GasliftValves`: a type to define a string of gas lift valves for valve & pressure calculations.

Constructor: `GasliftValves(md::Array, PTRO::Array, R::Array, port::Array)`

Port sizes must be in integer increments of 64ths inches.

Indicate orifice valves with an R-value and PTRO of 0.
"""
struct GasliftValves

    md::Array{Float64,1}
    PTRO::Array{Float64,1}
    R::Array{Float64,1}
    port::Array{Int64,1}

    function GasliftValves(md::Array{T} where T <: Real, PTRO::Array{T} where T <: Real, R::Array{T} where T <: AbstractFloat, port::Array{T} where T <: Union{Real, Int})

        ports = try
            convert(Array{Int64,1}, port)
        catch
            throw(ArgumentError("Specify port sizes in integer 64ths inches, e.g. 16 for a quarter-inch port."))
        end

        if any(R .> 1) || any(R .< 0)
            throw(ArgumentError("R-values are the area ratio of the port to the bellows and must be in [0, 1]."))
        elseif any(R .> 0.2)
            @info "Large R-value(s) entered--validate valve entry data."
        end

        new(convert(Array{Float64,1}, md), convert(Array{Float64,1}, PTRO), convert(Array{Float64,1}, R), ports)
    end
end

#printing for gas lift valves
Base.show(io::IO, valves::GasliftValves) = print(io, "Valve design with $(length(valves.md)) valves and bottom valve at $(valves.md[end])' MD.")

"""
`Wellbore`: type to define a flow path as an input for pressure drop calculations

See `read_survey` for helper method to create a Wellbore object from deviation survey files.

# Fields
- `md::Array{Float64, 1}`: measured depth for each segment in feet
- `inc::Array{Float64, 1}`: inclination from vertical for each segment in degrees, e.g. true vertical = 0°
- `tvd::Array{Float64, 1}`: true vertical depth for each segment in feet
- `id::Array{Float64, 1}`: inner diameter for each pip segment in inches

# Constructors
By default, negative depths are disallowed, and a 0 MD / 0 TVD point is added if not present, to allow graceful handling of outlet pressure definitions.
To bypass both the error checking and convenience feature, pass `true` as the final argument to the constructor.

`Wellbore(md, inc, tvd, id::Array{Float64, 1}, allow_negatives = false)`: defines a new Wellbore object from a survey with inner diameter defined for each segment. Lengths of each input array must be equal.

`Wellbore(md, inc, tvd, id::Float64, allow_negatives = false)`: defines a new Wellbore object with a uniform ID along the entire flow path.

`Wellbore(md, inc, tvd, id, valves::GasliftValves, allow_negatives = false)`: defines a new Wellbore object and adds interpolated survey points for each gas lift valve.
"""
struct Wellbore

    md::Array{Float64, 1}
    inc::Array{Float64, 1}
    tvd::Array{Float64, 1}
    id::Array{Float64, 1}

    function Wellbore(md, inc, tvd, id::Array{Float64, 1}, allow_negatives::Bool = false)

        lens = length.([md, inc, tvd, id])

        if !( count(x -> x == lens[1], lens) == length(lens) )
            throw(DimensionMismatch("Mismatched number of wellbore elements used in wellbore constructor."))
        end

        if !allow_negatives
            if minimum(md) < 0 || minimum(tvd) < 0
                throw(ArgumentError("Survey contains negative measured or true vertical depths. Pass the `allow_negatives` constructor flag if this is intentional."))
            end

            #add the origin/outlet reference point if missing
            if !(md[1] == tvd[1] <= 0)
                md = vcat(0, md)
                inc = vcat(0, inc)
                tvd = vcat(0, tvd)
                id = vcat(id[1], id)
            end
        end

        new(md, inc, tvd, id)
    end
end #struct Wellbore

#convenience constructor for uniform tubulars
Wellbore(md, inc, tvd, id::Float64, allow_negatives::Bool = false) = Wellbore(md, inc, tvd, repeat([id], inner = length(md)), allow_negatives)

#convenience constructors to add reference depths for valves so that they can be used as injection points
function Wellbore(md, inc, tvd, id, valves::GasliftValves, allow_negatives::Bool = false)

    well = Wellbore(md, inc, tvd, id, allow_negatives)

    for v in 1:length(valves.md)

        upper_index = searchsortedlast(well.md, valves.md[v])

        if well.md[upper_index] != valves.md[v]
            lower_index = upper_index + 1 #also the target insertion position

            x1, x2 = well.md[upper_index], well.md[lower_index]
            for property in [well.inc, well.tvd, well.id]
                y1, y2 = property[upper_index], property[lower_index]
                interpolated_value = y1 + (y2 - y1)/(x2 - x1) * (valves.md[v] - x1)
                insert!(property, lower_index, interpolated_value)
            end

            insert!(well.md, lower_index, valves.md[v])
        end
    end

    return well
end

#handle argument defaults in read_survey
function Wellbore(md, inc, tvd, id, valves::Nothing, allow_negatives::Bool = false)
    Wellbore(md, inc, tvd, id, allow_negatives)
end

#Printing for Wellbore structs
Base.show(io::IO, well::Wellbore) = print(io,
    "Wellbore with $(length(well.md)) points.\n",
    "Ends at $(well.md[end])' MD / $(well.tvd[end])' TVD.\n",
    "Max inclination $(maximum(well.inc))°. Average ID $(round(sum(well.id)/length(well.id), digits = 3)) in.")

#model struct. ONLY applies to wrapper functions.
"""
`WellModel`: Makes it easier to iterate well models

`pressure_and_temp(;model::WellModel)`

Develop pressure traverse in psia and temperature profile in °F from wellhead down to datum for a WellModel object. Requires the following fields to be defined in the model:

Returns a pressure profile as an Array{Float64,1} and a temperature profile as an Array{Float64,1}, referenced to the measured depths in the original Wellbore object.

Pressure correlation functions available:
- `BeggsAndBrill` with Payne correction factors
- `HagedornAndBrown` with Griffith and Wallis bubble flow correction

## Required
- `well::Wellbore`: Wellbore object that defines segmentation/mesh, with md, tvd, inclination, and hydraulic diameter
- `roughness`: pipe wall roughness in inches
- `temperature_method = "linear"`: temperature method to use; "Shiu" for Ramey method with Shiu relaxation factor, "linear" for linear interpolation
- `WHT = missing`: wellhead temperature in °F; required for `temperature_method = "linear"`
- `geothermal_gradient = missing`: geothermal gradient in °F per 100 ft; required for `temperature_method = "Shiu"`
- `BHT` = bottomhole temperature in °F
- `WHP`: absolute outlet pressure (wellhead pressure) in **psig**
- `dp_est`: estimated starting pressure differential (in psi) to use for all segments--impacts convergence time
- `q_o`: oil rate in stocktank barrels/day
- `q_w`: water rate in stb/d
- `GLR`: **total** wellhead gas:liquid ratio, inclusive of injection gas, in scf/bbl
- `APIoil`: API gravity of the produced oil
- `sg_water`: specific gravity of produced water
- `sg_gas`: specific gravity of produced gas

## Optional
- `injection_point = missing`: injection point in MD for gas lift, above which total GLR is used, and below which natural GLR is used
- `naturalGLR = missing`: GLR to use below point of injection, in scf/bbl
- `pressurecorrelation::Function = BeggsAndBrill: pressure correlation to use
- `error_tolerance = 0.1`: error tolerance for each segment in psi
- `molFracCO2 = 0.0`, `molFracH2S = 0.0`: produced gas fractions of hydrogen sulfide and CO2, [0,1]
- `pseudocrit_pressure_correlation::Function = HankinsonWithWichertPseudoCriticalPressure`: psuedocritical pressure function to use
- `pseudocrit_temp_correlation::Function = HankinsonWithWichertPseudoCriticalTemp`: pseudocritical temperature function to use
- `Z_correlation::Function = KareemEtAlZFactor`: natural gas compressibility/Z-factor correlation to use
- `gas_viscosity_correlation::Function = LeeGasViscosity`: gas viscosity correlation to use
- `solutionGORcorrelation::Function = StandingSolutionGOR`: solution GOR correlation to use
- `bubblepoint::Union{Function, Real} = StandingBubblePoint`: either bubble point correlation or bubble point in **psia**
- `oilVolumeFactor_correlation::Function = StandingOilVolumeFactor`: oil volume factor correlation to use
- `waterVolumeFactor_correlation::Function = GouldWaterVolumeFactor`: water volume factor correlation to use
- `dead_oil_viscosity_correlation::Function = GlasoDeadOilViscosity`: dead oil viscosity correlation to use
- `live_oil_viscosity_correlation::Function = ChewAndConnallySaturatedOilViscosity`: saturated oil viscosity correction function to use
- `frictionfactor::Function = SerghideFrictionFactor`: correlation function for Darcy-Weisbach friction factor
- `outlet_referenced = true`: whether to use outlet pressure (WHP) or inlet pressure (BHP) for
"""
mutable struct WellModel

    wellbore::Wellbore; roughness
    valves::Union{GasliftValves, Missing}
    temperatureprofile::Union{Array{T,1}, Missing} where T <: Real
    temperature_method; WHT; geothermal_gradient; BHT; casing_temp_factor
    pressurecorrelation::Function; outlet_referenced::Bool
    WHP; CHP; dp_est; dp_est_inj; error_tolerance; error_tolerance_inj
    q_o; q_w; GLR
    injection_point; naturalGLR
    APIoil; sg_water; sg_gas; sg_gas_inj
    molFracCO2; molFracH2S; molFracCO2_inj; molFracH2S_inj
    pseudocrit_pressure_correlation::Function
    pseudocrit_temp_correlation::Function
    Z_correlation::Function
    gas_viscosity_correlation::Function
    solutionGORcorrelation::Function
    bubblepoint::Union{Function, Real}
    oilVolumeFactor_correlation::Function
    waterVolumeFactor_correlation::Function
    dead_oil_viscosity_correlation::Function
    live_oil_viscosity_correlation::Function
    frictionfactor::Function

    function WellModel(;wellbore, roughness, valves = missing, temperatureprofile = missing,
                        temperature_method = "linear", WHT = missing, geothermal_gradient = missing, BHT = missing, casing_temp_factor = 0.85,
                        pressurecorrelation = BeggsAndBrill, outlet_referenced = true,
                        WHP, CHP = missing, dp_est, dp_est_inj = 0.1 * dp_est, error_tolerance = 0.1, error_tolerance_inj = 0.05,
                        q_o, q_w, GLR, injection_point = missing, naturalGLR = missing,
                        APIoil, sg_water, sg_gas, sg_gas_inj = sg_gas,
                        molFracCO2 = 0.0, molFracH2S = 0.0, molFracCO2_inj = molFracCO2, molFracH2S_inj = molFracH2S,
                        pseudocrit_pressure_correlation = HankinsonWithWichertPseudoCriticalPressure,
                        pseudocrit_temp_correlation = HankinsonWithWichertPseudoCriticalTemp,
                        Z_correlation = KareemEtAlZFactor, gas_viscosity_correlation = LeeGasViscosity,
                        solutionGORcorrelation = StandingSolutionGOR, bubblepoint = StandingBubblePoint, oilVolumeFactor_correlation = StandingOilVolumeFactor,
                        waterVolumeFactor_correlation = GouldWaterVolumeFactor,
                        dead_oil_viscosity_correlation = GlasoDeadOilViscosity, live_oil_viscosity_correlation = ChewAndConnallySaturatedOilViscosity,
                        frictionfactor = SerghideFrictionFactor)

        new(wellbore, roughness, valves, temperatureprofile, temperature_method, WHT, geothermal_gradient, BHT, casing_temp_factor,
            pressurecorrelation, outlet_referenced, WHP, CHP, dp_est, dp_est_inj, error_tolerance, error_tolerance_inj,
            q_o, q_w, GLR, injection_point, naturalGLR, APIoil, sg_water, sg_gas, sg_gas_inj, molFracCO2, molFracH2S, molFracCO2_inj, molFracH2S_inj,
            pseudocrit_pressure_correlation, pseudocrit_temp_correlation, Z_correlation, gas_viscosity_correlation, solutionGORcorrelation, bubblepoint,
            oilVolumeFactor_correlation, waterVolumeFactor_correlation, dead_oil_viscosity_correlation, live_oil_viscosity_correlation, frictionfactor)
    end

end

#Printing for model structs
function Base.show(io::IO, model::WellModel)

    fields = fieldnames(WellModel)
    values = map(f -> getfield(model, f), fields) |> list -> map(x -> !(x isa Array) ? string(x) : "$(length(x)) points from $(maximum(x)) to $(minimum(x)).", list)
    msg = string.(fields) .* " : " .* values

    println(io, "Well model: ")
    for item in msg
        println(io, item)
    end
end


text:
moons = html_div() do
    html_div(id="fields") do
        html_div(className="container-slider") do
            html_label("Samples: "),
            html_div(
                dcc_slider(
                        id = "n_samples",
                        min = 1,
                        max = 10000,
                        marks = Dict(i => "$i" for i = 0:1000:10000),
                        value = 100,
                ),
            style = (width = "70%", display = "inline-block"),
            className = "slider")
        end,
        html_div(className="container-slider") do
            html_label("noise: "),
            html_div(
                dcc_slider(
                        id = "noise",
                        min = 0.00,
                        max = 1.0,
                        step = 0.01,
                        marks = Dict(i => "$i" for i = 0.0:0.1:1.0),
                        value = 0.1,
                ),
            style = (width = "70%", display = "inline-block"),
            className = "slider")
        end
    end,
    html_div(className="container-graph") do
        dcc_graph(
            id = "plot_moons",
            figure = gen_moons(1000, 0.25)
        )
    end
end
callback!(app, Output("plot_moons", "figure"),[Input("n_samples", "value"), Input("noise", "value")]) do n_samples, noise
    gen_moons(n_samples,noise)
end


text:
"""
A generic interface for single-player games and two-player zero-sum games.

Stochastic games and intermediate rewards are supported. By convention,
rewards are expressed from the point of view of the player called _white_.
In two-player zero-sum games, we call `black` the player trying to minimize the reward.
"""
module GameInterface

export AbstractGameSpec, AbstractGameEnv

using ..AlphaZero: Util

#####
##### Game environments and game specifications
#####

# NAMING CONVENTION:
# - A game specification is usually written: `game_spec` or `gspec`.
# - A game environment is usually written: `game` or `g`.
#   We usually avoid `env` as it is too generic and can clash with the name for an
#   MCTS environment or training environment.
# - A state is usually written: `state` or `s`.
# - An action is usually written: `action` or `a`.

"""
    AbstractGameSpec

Abstract type for a game specification.

The specification holds all _static_ information about a game, which does not
depend on the current state.
"""
abstract type AbstractGameSpec end

"""
    AbstractGameEnv

Abstract base type for a game environment.

Intuitively, a game environment holds a game specification and a current state.
"""
abstract type AbstractGameEnv end

"""
    init(::AbstractGameSpec) :: AbstractGameEnv

Create a new game environment in a (possibly random) initial state.
"""
function init end

"""
    spec(game::AbstractGameEnv) :: AbstractGameSpec

Return the game specification of an environment.
"""
function spec end

#####
##### Queries on specs
#####

"""
    two_players(::AbstractGameSpec) :: Bool

Return whether or not a game is a two-players game.
"""
function two_players end

"""
    actions(::AbstractGameSpec)

Return the vector of all game actions.
"""
function actions end

"""
    vectorize_state(::AbstractGameSpec, state) :: Array{Float32}

Return a vectorized representation of a given state.
"""
function vectorize_state end

#####
##### Operations on envs
#####

"""
    set_state!(game::AbstractGameEnv, state)

Modify the state of a game environment in place.
"""
function set_state! end

"""
    current_state(game::AbstractGameEnv)

Return the game state.

!!! warn

    The state returned by this function may be stored (e.g. in the MCTS tree) and must
    therefore either be fresh or persistent. If in doubt, you should make a copy.

"""
function current_state end

# TODO: maybe MCTS should make the copy itself. The performance cost should not be great
# and it would probably avoid people a lot of pain.

"""
    game_terminated(::AbstractGameEnv)

Return a boolean indicating whether or not the game is in a terminal state.
"""
function game_terminated end

"""
    white_playing(::AbstractGameEnv) :: Bool

Return `true` if white is to play and `false` otherwise.

For a one-player game, this function must always return `true`.
"""
function white_playing end

"""
    actions_mask(::AbstractGameEnv)

Return a boolean mask indicating what actions are available.

The following identities must hold:

  - `game_terminated(game) || any(actions_mask(game))`
  - `length(actions_mask(game)) == length(actions(spec(game)))`
"""
function actions_mask end

"""
    play!(game::AbstractGameEnv, action)

Update the game environment by making the current player perform `action`.
Note that this function does not have to be deterministic.
"""
function play! end

"""
    white_reward(game::AbstractGameEnv)

Return the intermediate reward obtained by the white player after the last
transition step. The result is undetermined when called at an initial state.
"""
function white_reward end

"""
    heuristic_value(game::AbstractGameEnv)

Return a heuristic estimate of the state value for the current player.

The given state must be nonfinal and returned values must belong to the
``(-∞, ∞)`` interval.

This function is not needed by AlphaZero but it is useful for building
baselines such as minmax players.
"""
function heuristic_value end

#####
##### Symmetries
#####

"""
    symmetries(::AbstractGameSpec, state)

Return the vector of all pairs `(s, σ)` where:
  - `s` is the image of `state` by a nonidentical symmetry
  - `σ` is the associated actions permutation, as an integer vector of
     size `num_actions(game)`.

A default implementation is provided that returns an empty vector.

Note that the current state of the passed environment is ignored by this function.

# Example

In the game of tic-tac-toe, there are eight symmetries that can be
obtained by composing reflexions and rotations of the board (including the
identity symmetry).

# Property

If `(s2, σ)` is a symmetry for state `s1`, then `mask2 == mask1[σ]` must hold where
`mask1` and `mask2` are the available action masks for `s1` and `s2` respectively.
"""
function symmetries(::AbstractGameSpec, state)
  return Tuple{typeof(state), Vector{Int}}[]
end

#####
##### Interface for interactive exploratory tools
#####

"""
    render(game::AbstractGameEnv)

Print the game state on the standard output.
"""
function render end

"""
    action_string(::AbstractGameSpec, action) :: String

Return a human-readable string representing the provided action.
"""
function action_string end

"""
    parse_action(::AbstractGameSpec, str::String)

Return the action described by string `str` or `nothing` if `str` does not
denote a valid action.
"""
function parse_action end

"""
    read_state(game_spec::AbstractGameSpec)

Read a state from the standard input.
Return the corresponding state (with type `state_type(game_spec)`)
or `nothing` in case of an invalid input.
"""
function read_state end

#####
##### Derived spec functions
#####

"""
    state_type(::AbstractGameSpec)

Return the state type associated to a game.

State objects must be persistent or appear as such as they are stored into
the MCTS tree without copying. They also have to be comparable and hashable.
"""
function state_type(game_spec::AbstractGameSpec)
  return typeof(current_state(init(game_spec)))
end

"""
    state_dim(::AbstractGameSpec)

Return a tuple that indicates the shape of a vectorized state representation.
"""
function state_dim(game_spec::AbstractGameSpec)
  state = current_state(init(game_spec))
  return size(vectorize_state(game_spec, state))
end

"""
    state_memsize(::AbstractGameSpec)

Return the memory footprint occupied by a state of the given game.

The computation is based on a random initial state, assuming that all states have an
identical footprint.
"""
function state_memsize(game_spec::AbstractGameSpec)
  state = current_state(init(game_spec))
  return Base.summarysize(state)
end

"""
    action_type(::AbstractGameSpec)

Return the action type associated to a game.
"""
function action_type(game_spec::AbstractGameSpec)
  return eltype(actions(game_spec))
end

"""
    num_actions(::AbstractGameSpec)

Return the total number of actions associated with a game.
"""
num_actions(game_spec::AbstractGameSpec) = length(actions(game_spec))

"""
    init(::AbstractGameSpec, state) :: AbstractGameEnv

Create a new game environment, initialized in a given state.
"""
function init(gspec::AbstractGameSpec, state)
  env = init(gspec)
  set_state!(env, state)
  return env
end

#####
##### Derived env functions
#####

"""
    clone(::AbstractGameEnv)

Return an independent copy of the given environment.
"""
function clone(game::AbstractGameEnv)
  return init(spec(game), current_state(game))
end

"""
    available_actions(::AbstractGameEnv)

Return the vector of all available actions.
"""
function available_actions(game::AbstractGameEnv)
  mask = actions_mask(game)
  return actions(spec(game))[mask]
end

"""
    apply_random_symmetry!(::AbstractGameEnv)

Update a game environment by applying a random symmetry
to the current state (see [`symmetries`](@ref)).
"""
function apply_random_symmetry!(game::AbstractGameEnv)
  gspec = spec(game)
  syms = symmetries(gspec, current_state(game))
  @assert !isempty(syms) "no symmetries were declared for this game"
  symstate, _ = rand(syms)
  set_state!(game, symstate)
  return
end

end


text:
# Stolen from examples/polymake.jl
# Needs polymake (obviously)

# Returns all vectors v such that Av == b and Cv >= 0.
function solve_mixed(A::fmpz_mat, b::fmpz_mat, C::fmpz_mat)
  bA = Array{BigInt, 2}(hcat(-b, A))
  z = findall(i->!iszero_row(C, i), 1:nrows(C))
  zI = Array{BigInt, 2}(hcat(zero_matrix(FlintZZ, nrows(C), 1), C))[z, :]
  p = Polymake.@pm Polytope.Polytope(:EQUATIONS => bA, :INEQUALITIES => zI)

  res = Vector{fmpz_mat}(undef, p.N_LATTICE_POINTS)
  if p.N_LATTICE_POINTS == 0
    return res
  end

  inner_points = p.INTERIOR_LATTICE_POINTS
  boundary_points = p.BOUNDARY_LATTICE_POINTS

  for i = 1:Int(size(boundary_points)[1]) # "==" nrows(boundary_points)
    @assert boundary_points[i, 1] == 1
    c = zero_matrix(FlintZZ, 1, ncols(A))
    for j = 1:ncols(A)
      c[1, j] = boundary_points[i, j + 1]
    end
    res[i] = c
  end

  for i = 1:Int(size(inner_points)[1])
    @assert inner_points[i, 1] == 1
    c = zero_matrix(FlintZZ, 1, ncols(A))
    for j = 1:ncols(A)
      c[1, j] = inner_points[i, j + 1]
    end
    res[i + Int(size(boundary_points)[1])] = c
  end
  return res
end


text:
using Test
using FixedPoint

const R = Rational{BigInt}
const R1 = Rational{Int64}

# Exhaustive tests for Int8, just because we can

function nullary(::Type{F8}) where F8
    f = precision(F8)
    typemax8 = typemax(Int8) >> f
    inrange8(x) = abs(x) <= typemax8
    eps8 = R(1, 2^f)
    round8(x) = round(x * 2^f) / 2^f

    @testset "Nullary functions, $f fraction bits" begin
        @test isequal(R(zero(F8)), zero(R))
        @test isequal(R(one(F8)), one(R))
        @test isequal(R(eps(F8)), inv(R(2^f)))
        @test isequal(R(typemin(F8)), -R(typemax(F8)))
        @test isequal(R(typemax(F8)), R(typemax(Int8)) / R(2^f))
    end
end

function basic(::Type{F8}) where {F8}
    f = precision(F8)
    typemax8 = typemax(Int8) >> f
    inrange8(x) = abs(x) <= typemax8
    eps8 = R(1, 2^f)
    round8(x) = round(x * 2^f) / 2^f

    @testset "Basic properties, $f fraction bits" begin
        for i in -typemax(Int8) : typemax(Int8)
            x = F8(FixedPoint.Internal(), i)
            r = R(i, 2^f)

            @test isequal(eltype(x), Int8)
            @test isequal(precision(x), f)
            @test isequal(numerator(x), i)
            @test isequal(denominator(x), Int8(2)^f)
        end
    end
end

function conversions(::Type{F8}) where {F8}
    f = precision(F8)
    typemax8 = typemax(Int8) >> f
    inrange8(x) = abs(x) <= typemax8
    eps8 = R(1, 2^f)
    round8(x) = round(x * 2^f) / 2^f

    @testset "Type conversions, $f fraction bits" begin
        for i in -typemax(Int8) : typemax(Int8)
            x = F8(FixedPoint.Internal(), i)
            r = R(i, 2^f)

            @test isequal(F8(x), x)
            @test isequal(R(x), r)
            @test isequal(F8(r), x)
            @test isequal(R(Float32(x)), r)
            @test isequal(F8(Float32(r)), x)
            @test isequal(R(float(x)), float(r))
            @test isequal(F8(float(r)), x)
        end
    end
end

function unary(::Type{F8}) where {F8}
    f = precision(F8)
    typemax8 = typemax(Int8) >> f
    inrange8(x) = abs(x) <= typemax8
    eps8 = R(1, 2^f)
    round8(x) = round(x * 2^f) / 2^f

    @testset "Unary functions, $f fraction bits" begin
        for i in -typemax(Int8) : typemax(Int8)
            x = F8(FixedPoint.Internal(), i)
            r = R(i, 2^f)

            @test isequal(R(+x), +r)
            @test isequal(R(-x), -r)
            @test isequal(R(*(x)), *(r))
            @test isequal(R(abs(x)), abs(r))
            @test hash(x) === hash(x)
            @test hash(x) !== hash(i)
            @test hash(x) !== hash(r)
            if r != 0 && inrange8(inv(r))
                # y = inv(r)
                # dy = max(abs(r), abs(y)) * eps8
                # # TODO: Instead of allowing for an error here, fix the
                # # implementation and require things to be rounded
                # # correctly
                # if inrange8(y) && abs(y) > dy
                #     @test abs(R(inv(x)) - round8(inv(r))) <= dy
                # end
                if !isequal(R(inv(x)), round8(inv(r)))
                    @show f i x r R(inv(x)) round8(inv(r))
                end
                @test isequal(R(inv(x)), round8(inv(r)))
            end
            @test isequal(R(sign(x)), sign(r))
            for a in -10:10
                if r != 0
                    y = r^a
                    dy = 4 * max(1, abs(r), abs(y)) * eps8
                    if inrange8(y) && (a >= 0 || inrange8(inv(y))) &&
                            abs(y) > dy && (a >= 0 || abs(inv(y)) > dy)
                        @test abs(R(x^a) - y) <= dy
                    end
                end
            end
        end
    end
end

function binary(::Type{F8}) where {F8}
    f = precision(F8)
    typemax8 = typemax(Int8) >> f
    inrange8(x) = abs(x) <= typemax8
    eps8 = R(1, 2^f)
    round8(x) = round(x * 2^f) / 2^f

    @testset "Binary functions, $f fraction bits" begin
        for i in -typemax(Int8) : typemax(Int8),
                j in -typemax(Int8) : typemax(Int8)
            x = F8(FixedPoint.Internal(), i)
            y = F8(FixedPoint.Internal(), j)
            r = R(i, 2^f)
            s = R(j, 2^f)

            @test (r == s) === (hash(x) == hash(y))

            if inrange8(r + s)
                @test isequal(R(x + y), r + s)
            end
            if inrange8(r - s)
                @test isequal(R(x - y), r - s)
            end
            if inrange8(r * s)
                @test isequal(R(x * y), round8(r * s))
            end
            if s != 0 && inrange8(r / s)
                @test isequal(R(x / y), round8(r / s))
            end
            if r != 0 && inrange8(r \ s)
                @test isequal(R(x \ y), round8(r \ s))
            end
            @test isequal(R(copysign(x, y)), copysign(r, s))
            @test isequal(R(flipsign(x, y)), flipsign(r, s))
            @test isequal(R(max(x, y)), max(r, s))
            @test isequal(R(min(x, y)), min(r, s))
        end
    end
end

function comparisons(::Type{F8}) where {F8}
    f = precision(F8)
    typemax8 = typemax(Int8) >> f
    inrange8(x) = abs(x) <= typemax8
    eps8 = R(1, 2^f)
    round8(x) = round(x * 2^f) / 2^f

    @testset "Comparison operators, $f fraction bits" begin
        for i in -typemax(Int8) : typemax(Int8),
                j in -typemax(Int8) : typemax(Int8)
            x = F8(FixedPoint.Internal(), i)
            y = F8(FixedPoint.Internal(), j)
            r = R(i, 2^f)
            s = R(j, 2^f)
            @test isequal(x == y, r == s)
            @test isequal(x != y, r != s)
            @test isequal(x < y, r < s)
            @test isequal(x <= y, r <= s)
            @test isequal(x > y, r > s)
            @test isequal(x >= y, r >= s)
        end
    end
end

function ternary(::Type{F8}) where {F8}
    f = precision(F8)
    typemax8 = typemax(Int8) >> f
    inrange8(x) = abs(x) <= typemax8
    eps8 = R(1, 2^f)
    round8(x) = round(x * 2^f) / 2^f

    @testset "Ternary functions, $f fraction bits" begin
        for i in -typemax(Int8) : typemax(Int8),
                j in -typemax(Int8) : typemax(Int8),
                k in -typemax(Int8) : typemax(Int8)
            x = F8(FixedPoint.Internal(), i)
            y = F8(FixedPoint.Internal(), j)
            z = F8(FixedPoint.Internal(), k)
            r = R1(i, 2^f)
            s = R1(j, 2^f)
            t = R1(k, 2^f)

            if inrange8(r * s) && inrange8(muladd(r, s, t))
                q = muladd(r, s, t)
                dq = max(1, abs(r), abs(s), abs(t), abs(q)) * R1(1, 2^f)
                @test abs(R1(muladd(x, y, z)) - q) <= dq
            end
            if inrange8(r * s) && inrange8(fma(r, s, t))
                @test isequal(R1(fma(x, y, z)), round8(fma(r, s, t)))
            end
        end
    end
end

function elementary(::Type{F8}) where {F8}
    f = precision(F8)
    typemax8 = typemax(Int8) >> f
    inrange8(x) = abs(x) <= typemax8
    eps8 = R(1, 2^f)
    round8(x) = round(x * 2^f) / 2^f

    @testset "Elementary functions, $f fraction bits" begin
        for i in -typemax(Int8) : typemax(Int8)
            x = F8(FixedPoint.Internal(), i)
            r = R(i, 2^f)

            if r >= 0
                @test isequal(R(sqrt(x)), round8(sqrt(r)))
            end
        end
    end
end

for f in 0:6
    F8 = Fixed{f, Int8}
    nullary(F8)
    basic(F8)
    conversions(F8)
    unary(F8)
    binary(F8)
    comparisons(F8)
    ternary(F8)
    elementary(F8)
end


text:
"""
    finalize_jsi()

Finalize the package JustSayIt.

See also: [`init_jsi`](@ref)
"""
function finalize_jsi()

end


text:
using NeighbourLists
using JuLIP
using BenchmarkTools
using PyCall
using ASE
# using ProfileView

include("test_aux.jl")

# @pyimport matscipy.neighbours as matscipy_neighbours
matscipy_neighbours = pyimport("matscipy.neighbours")

function matscipy_nlist(at, cutoff)
   pycall(matscipy_neighbours["neighbour_list"],
          NTuple{5, PyArray}, "ijdDS", at.po, cutoff)
end

println("# Threads = ", Base.Threads.nthreads())

println("Neighbourlist assembly benchmark")
for L in [2, 4, 10, 30]
   print("L = $L")
   # si, non-cubic cell, mixed bc
   at = bulk(:Si, cubic=true) * L
   println(", N = $(length(at))")
   set_pbc!(at, (true, true, true))
   C = JMat(cell(at))
   X = positions(at)
   perbc = JVec(pbc(at))
   cutoff = 3.5 * rnn(:Si)

   println("Julia Nlist")
   @btime PairList($X, $cutoff, $C, $perbc, sorted = false)
   println("Matscipy Nlist")
   @btime matscipy_nlist($(ASEAtoms(at)), $cutoff)
   println("------------------------------------------")
end

print("N-body energy and forces benchmark: ")
at = bulk(:Cu, cubic=true) * 10
println("bulk :Cu with nat = $(length(at))")
r0 = rnn(:Cu)
rcut = 2.1 * r0
X = positions(at)
C = cell(at)
f, f_d = gen_fnbody(rcut, r0)
nlist = PairList(X, rcut, C, (false, false, false), sorted = true)

for M in [2, 3] #, 4, 5]
   println("M = $M")
   print("  energy: ")
   @btime n_body($f, $M, $nlist)
   print("  forces: ")
   @btime grad_n_body($f_d, $M, $nlist)
end


text:
module Rsvg

depsjl = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
isfile(depsjl) ? include(depsjl) : error("Rsvg not properly ",
    "installed. Please run\nPkg.build(\"Rsvg\")")

using Cairo
using Gtk

include("types.jl")
include("calls.jl")

export RsvgDimensionData, RsvgHandle

Rsvg.set_default_dpi(72.0)

end # module


text:
export AbstractDevice, Device
export setup_ϵᵣ!, setup_src!
export Mode, add_mode!

struct Mode
    pol::Polarization
    dir::Direction
    neff::Number
    pt::Point
    width::Number
end

abstract type AbstractDevice{D} end

Base.size(d::AbstractDevice) = size(d.grid);
Base.size(d::AbstractDevice, i::Int) = size(d.grid, i);
Base.length(d::AbstractDevice) = length(d.grid);

mutable struct Device{D} <: AbstractDevice{D}
	grid::Grid{D}
	ϵᵣ::Array{Complex}
    src::Array{Complex}
    ω::AbstractVector{Float}
    modes::Array{Mode}
end

"    Device(grid::Grid, ω::AbstractVector{Float})"
function Device(grid::Grid, ω::AbstractVector{Float})
	ϵᵣ = ones(Complex, size(grid));
    src = zeros(Complex, size(grid));
	return Device(grid, ϵᵣ, src, ω, Array{Mode}(undef, 0))
end

"    Device(grid::Grid, ω::Float)"
Device(grid::Grid, ω::Float) = Device(grid, [ω]);

# ============================================================================ #

"    normalize_parameters(g::Grid)"
normalize_parameters(g::Grid) = (ϵ₀*g.L₀, μ₀*g.L₀, c₀/g.L₀);

"    normalize_parameters(d::AbstractDevice)"
normalize_parameters(d::AbstractDevice) = normalize_parameters(d.grid);

# ============================================================================ #

function _compose_shapes!(pixels::AbstractArray, grid::Grid, shapes::AbstractVector{<:Shape})
    kd = KDTree(shapes);
    for i in eachindex(pixels)
        x = xc(grid, i);
        y = yc(grid, i);
        shape = findfirst([x, y, 0.0], kd);
        if shape != nothing
            if iscallable(shape.data)
                pixels[i] = shape.data(x, y);
            else
                pixels[i] = shape.data;
            end
        end
    end
end

function _mask_values!(pixels::AbstractArray, grid::Grid, region, value)
    if ndims(grid) == 2
        mask = [region(x, y) for x in xc(grid), y in yc(grid)];
        if iscallable(value)
            value_assigned = [value(x, y) for x in xc(grid), y in yc(grid)];
            pixels[mask]   = value_assigned[mask];
        else
            pixels[mask]  .= value;
        end
    elseif ndims(grid) == 1
        mask = [region(x) for x in xc(grid)];
        if iscallable(value)
            value_assigned = [value(x) for x in xc(grid)];
            pixels[mask]   = value_assigned[mask];
        else
            pixels[mask]  .= value;
        end
    else
        error("Unkown device dimension!")
    end;
end

"    setup_ϵᵣ!(d::AbstractDevice, shapes::AbstractVector{<:Shape})"
setup_ϵᵣ!(d::AbstractDevice, shapes::AbstractVector{<:Shape}) = _compose_shapes!(d.ϵᵣ, d.grid, shapes)

"    setup_ϵᵣ!(d::AbstractDevice, region, value)"
setup_ϵᵣ!(d::AbstractDevice, region, value) = _mask_values!(d.ϵᵣ, d.grid, region, value)

"    setup_src!(d::AbstractDevice, region, value)"
setup_src!(d::AbstractDevice, region, value) = _mask_values!(d.src, d.grid, region, value)

"    setup_src!(d::AbstractDevice, pt::Point)"
function setup_src!(d::AbstractDevice, pt::Point)
    (indx, indy) = coord2ind(d.grid, pt);
    d.src[indx, indy] = 1im;
end

"    setup_src!(d::AbstractDevice, pt::Point, srcnormal::Direction)"
function setup_src!(d::AbstractDevice, pt::Point, srcnormal::Direction)
    (indx, indy) = coord2ind(d.grid, pt);
    if srcnormal == x̂
        d.src[indx, :] = 1im;
    elseif srcnormal == ŷ
        d.src[:, indy] = 1im;
    end
end

# ============================================================================ #

"    add_mode!(d::AbstractDevice, mode::Mode)"
function add_mode!(d::AbstractDevice, mode::Mode)
    append!(d.modes, [mode]);
end

"    setup_mode!(d::AbstractDevice, pol::Polarization, ω::Float, neff::Number, pt::Point, srcnormal::Direction, srcwidth::Number)"
function setup_mode!(d::AbstractDevice, pol::Polarization, ω::Float, neff::Number, pt::Point, srcnormal::Direction, srcwidth::Number)
    (_, vector, indx, indy) = get_modes(d, pol, ω, neff, 1, pt, srcnormal, srcwidth)
    d.src[indx, indy] += normalize(abs.(vector[:]));
end

"    get_modes(d::AbstractDevice, pol::Polarization, ω::Float, neff::Number, nmodes::Int, pt::Point, slicenormal::Direction, slicewidth::Number)"
function get_modes(d::AbstractDevice, pol::Polarization, ω::Float, neff::Number, nmodes::Int, pt::Point, slicenormal::Direction, slicewidth::Number)
    (indx, indy) = coord2ind(d.grid, pt);

    slicenormal == x̂ && (srcpoints = Int(round(slicewidth/dy(d.grid))));
    slicenormal == ŷ && (srcpoints = Int(round(slicewidth/dx(d.grid))));
    srcpoints % 2 == 0 && (srcpoints += 1); # source points is odd

    M = Int((srcpoints-1)/2);
    srcpoints = 2*M+1;

    if slicenormal == x̂
        indx = indx;
        indy = indy .+ (-M:M);
        dh = dy(d.grid);
    elseif slicenormal == ŷ
        indx = indx+(-M:M);
        indy = indy;
        dh = dx(d.grid);
    end

    g1D   = Grid(srcpoints, [0 srcpoints*dh], L₀=d.grid.L₀);
    dev1D = Device(g1D, ω);
    dev1D.ϵᵣ = d.ϵᵣ[indx, indy];
    (β, vectors) = eigenmode(dev1D, pol, neff, nmodes);
    return (β, vectors, indx, indy)
end


text:
abstract type AbstractGainPhaseMismCrosstalk{T} end

"""
(a * t^e + b) / (t^e + 1)
"""
struct AsymptoticGainPhaseMismCrosstalk{
        N,
        T <: AbstractFloat,
        E <: Union{Float64, Int64},
        CT <: SMatrix{N, N, Complex{T}}
    } <: AbstractGainPhaseMismCrosstalk{T}
    a::SVector{N, T}
    b::SVector{N, T}
    e::E
    C::CT
    t::typeof(1.0s)
end

function AsymptoticGainPhaseMismCrosstalk(
    a::SVector{N, T},
    b::SVector{N, T};
    e = 2,
    C = SMatrix{N,N,Complex{T}}(I),
    t = 0.0s
) where {N, T <: AbstractFloat}
    AsymptoticGainPhaseMismCrosstalk(a, b, e, C, t)
end

function propagate(gpmc::AsymptoticGainPhaseMismCrosstalk, Δt, rng)
    t = gpmc.t + Δt
    AsymptoticGainPhaseMismCrosstalk(gpmc.a, gpmc.b, gpmc.e, gpmc.C, t)
end

function get_gain_phase_mism_crosstalk(gpmc::AsymptoticGainPhaseMismCrosstalk{N,T}) where {N,T}
    phase_misms = cis.(
        (gpmc.a * (gpmc.t / s) ^ gpmc.e + gpmc.b) ./ ((gpmc.t / s) ^ gpmc.e + 1)
    )
    SMatrix{N,N,Complex{T}}(phase_misms .* gpmc.C)
end

@inline function propagate(gain_phase_mism_crosstalk, Δt, rng)
    gain_phase_mism_crosstalk
end

@inline function get_gain_phase_mism_crosstalk(gain_phase_mism_crosstalk)
    gain_phase_mism_crosstalk
end

"""
$(SIGNATURES)

Normalizes the gain and phase mismatch and crosstalk function so that the norm over columns is 1.

"""
function normalize_gain_phase_mism_and_crosstalk(gain_and_phase_mism_and_crosstalk)
    gain_and_phase_mism_and_crosstalk_norm = map(
        norm,
        eachcol(gain_and_phase_mism_and_crosstalk)
    )'
    gain_and_phase_mism_and_crosstalk ./ gain_and_phase_mism_and_crosstalk_norm
end

struct RandomWalkGainPhaseMismCrosstalk{
    N,
    T <: AbstractFloat,
    P <: SMatrix,
    G <: SVector,
    CP <: SMatrix,
    CA <: SVector
} <: AbstractGainPhaseMismCrosstalk{T}
    phase_mism::P
    gain_mism::G
    crosstalk_phase::CP
    crosstalk_ampl::CA
    phase_mism_vel_std::Float64
    crosstalk_phase_vel_std::Float64
end

function RandomWalkGainPhaseMismCrosstalk(
    N,
    T = Float64;
    gain_mism_std = 0.1,
    relative_crosstalk_ampl_std = 0.1,
    phase_mism_vel_std,
    crosstalk_phase_vel_std,
    crosstalk_ampl,
)
    gain_mism = ones(SVector{N, T}) + randn(SVector{N, T}) * gain_mism_std
    phase_mism = zeros(SMatrix{2, N, T})
    crosstalk_phase = zeros(SMatrix{2, N^2 - N, T})
    crosstalk_ampl = (ones(SVector{N^2 - N, T}) + randn(SVector{N^2 - N, T}) * relative_crosstalk_ampl_std) * crosstalk_ampl
    RandomWalkGainPhaseMismCrosstalk{N, T, typeof(phase_mism), typeof(gain_mism), typeof(crosstalk_phase), typeof(crosstalk_ampl)}(
        phase_mism,
        gain_mism,
        crosstalk_phase,
        crosstalk_ampl,
        phase_mism_vel_std,
        crosstalk_phase_vel_std,
    )
end

@inline function propagate(gpmc::RandomWalkGainPhaseMismCrosstalk{N, T, P, G, CP, CA}, Δt, rng) where {N, T, P, G, CP, CA}
    process = get_process(Order(2), upreferred(Δt / s))
    C = cholesky(get_process_covariance(Order(2), upreferred(Δt / s))).L
    next_phase_mism = process * gpmc.phase_mism + C * randn(rng, SMatrix{2, N, T}) * gpmc.phase_mism_vel_std
    next_crosstalk_phase = process * gpmc.crosstalk_phase + C * randn(rng, SMatrix{2, N^2 - N, T}) * gpmc.crosstalk_phase_vel_std
    RandomWalkGainPhaseMismCrosstalk{N, T, P, G, CP, CA}(next_phase_mism, gpmc.gain_mism, next_crosstalk_phase, gpmc.crosstalk_ampl, gpmc.phase_mism_vel_std, gpmc.crosstalk_phase_vel_std)
end

@inline function get_gain_phase_mism_crosstalk(gpmc::RandomWalkGainPhaseMismCrosstalk{N, T}) where {N, T}
    C = MMatrix{N,N,Complex{T}}(undef)
    crosstalk_idx = 1
    for i = 1:N, j = 1:N
        if i == j
            C[i,i] = cis(gpmc.phase_mism[1, i]) * gpmc.gain_mism[i]
        else
            C[i,j] = cis(gpmc.crosstalk_phase[1, crosstalk_idx]) * gpmc.crosstalk_ampl[crosstalk_idx]
            crosstalk_idx += 1
        end
    end
    SMatrix(C)
end


text:
type ActivationDesc
  ptr::Ptr{Void}
end

function ActivationDesc(mode, relu_nanopt=CUDNN_NOT_PROPAGATE_NAN, relu_ceiling=1.0)
  p = Ptr{Void}[0]
  cudnnCreateActivationDescriptor(p)
  cudnnSetActivationDescriptor(p[1], mode, relu_nanopt, relu_ceiling)
  desc = ActivationDesc(p[1])
  finalizer(desc, cudnnDestroyActivationDescriptor)
  desc
end

Base.unsafe_convert(::Type{Ptr{Void}}, desc::ActivationDesc) = desc.ptr


text:
using JuliaFormatter

function main()
    perfect = format(joinpath(@__DIR__, ".."); style=YASStyle())
    if perfect
        @info "Linting complete - no files altered"
    else
        @info "Linting complete - files altered"
        run(`git status`)
    end
    return nothing
end

main()


text:
abstract type Input end

"""
    get_value(input::Constant, time)

Get input value at time t.
"""
function get_value end

"""
    ConstantInput(val)

Constant value input.
"""
struct ConstantInput <: Input
    val::Float64
end

"""
    StepInput(val, tini, tfin)

Step input. Takes value from tini to tfin and 0 otherwise.
"""
struct StepInput <: Input
    val::Float64
    tini::Float64
    tfin::Float64
end

"""
    DoubletInput(val, tini, tfin)

Symmetric doublet input with val peak to peak amplitude. First positive.
"""
DoubletInput(val, tini, tfin) = vcat(
    StepInput(0.5*val, tini, 0.5*(tini+tfin)),
    StepInput(-0.5*val, 0.5*(tini+tfin), tfin)
)

"""
    InverseDoubletInput(val, tini, tfin)

Symmetric doublet input with val peak to peak amplitude. First negative.
"""
InverseDoubletInput(val, tini, tfin) = vcat(
    StepInput(-0.5*val, tini, 0.5*(tini+tfin)),
    StepInput(0.5*val, 0.5*(tini+tfin), tfin)
)

"""
    RampInput(val, tini, tend)

Ramp input. Takes value at tend increasing linearly from tini to tend.
"""
struct RampInput <: Input
    val::Float64
    tini::Float64
    tfin::Float64
end

"""
    SinusoidalInput(A, f, ϕ, tini, tfin)

Sinusoidal input. Takes a A/2 * sin(2πf (t-tini) + ϕ) value between tini and tfin.
Value is 0 otherwise.
"""
struct SinusoidalInput <: Input
    "Peak to peak amplitude"
    A::Float64
    "Frequency (1/s)"
    f::Float64
    "Phase (rad)"
    ϕ::Float64
    tini::Float64
    tfin::Float64
end

# --------- get_value methods ---------
get_value(input::ConstantInput, time) = input.val
get_value(input::StepInput, time) = input.tini <= time < input.tfin ? input.val : 0.

function get_value(input::RampInput, time)
    if input.tini <= time < input.tfin
        return input.val * (time - input.tini) / (input.tfin - input.tini)
    else
        return 0.
    end
end

function get_value(input::SinusoidalInput, time)
    if input.tini <= time < input.tfin
        return 0.5 * input.A * sin(2π * input.f * (time - input.tini) + input.ϕ)
    else
        return 0.
    end
end

function get_value(input_arr::Array{T,1}, t) where T<:Input
    value = 0.0
    for input in input_arr
        value = value + get_value(input, t)
    end
    return value
end


text:
# Custom pretty-printing

function Base.show(io::IO, ::MIME"text/plain", a::Union{VectorDifference,MatrixDifference})
    println(io, nameof(typeof(a)), " with indices:")
    print(io, " modified: ", a.modinds, "\n added: ", a.addinds, "\n removed: ", a.reminds)
    println(io, "\nand values:")
    println(io, " common: ", a.modvals)
    println(io, " added: ", a.addvals)
    print(io, " removed: ", a.remvals)
end

function Base.show(
    io::IO, ::MIME"text/plain", a::Union{DictDifference,NamedTupleDifference}
)
    println(io, nameof(typeof(a)), " with values:")
    println(io, " modified: ", a.modvals)
    println(io, " added: ", a.addvals)
    print(io, " removed: ", a.remvals)
end

function Base.show(io::IO, ::MIME"text/plain", a::SetDifference)
    println(io, nameof(typeof(a)), " with values:")
    println(io, " common: ", a.comvals)
    println(io, " added: ", a.addvals)
    print(io, " removed: ", a.remvals)
end


text:
using Plots
unicodeplots()

xs = -10:10
p1 = plot(xs,xs)
p2 = plot(xs,xs.^2)
p3 = plot(xs,xs.^3)
p4 = plot(xs,xs.^4)

plot(p1,p2,p3,p4,layout=(2,2),legend=false)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
###############
# data
################

type class_variables
	# private variables
  _v::Array{Float64,1}
	_x_ind::UnitRange{Int64}
	_y_ind::UnitRange{Int64}
	_s_ind::UnitRange{Int64}
	_tau_ind::Int64
	_kappa_ind::Int64

	function class_variables(n::Int64, m::Int64)
		this = new();
		# construct indicies
		this._x_ind = 1:n;
		this._y_ind = (n+1):(n+m);
		this._s_ind = (n+m+1):(2*n + m);
		this._tau_ind = 2*n + m + 1;
		this._kappa_ind = 2*n + m + 2;

		# intialize
		intial_val = 1.0;
		this._v = zeros(2*(n+1)+m);
		this._v[this._x_ind] = intial_val;
		this._v[this._s_ind] = intial_val;
		this._v[this._tau_ind] = intial_val;
		this._v[this._kappa_ind] = intial_val;

		return this;
	end
end

##################################
# methods
##################################

function x(vars::class_variables)
    return vars._v[vars._x_ind]
end

function x(vars::class_variables, new_vals::Array{Float64})
    return vars._v[vars._x_ind] = new_vals;
end

function y(vars::class_variables)
    return vars._v[vars._y_ind]
end

function y(vars::class_variables, new_vals::Array{Float64})
    return vars._v[vars._y_ind] = new_vals;
end

function s(vars::class_variables)
    return vars._v[vars._s_ind]
end

function s(vars::class_variables, new_vals::Array{Float64})
    return vars._v[vars._s_ind] = new_vals;
end

function tau(vars::class_variables)
    return vars._v[vars._tau_ind]
end

function tau(vars::class_variables, new_vals::Float64)
    return vars._v[vars._tau_ind] = new_vals;
end

function kappa(vars::class_variables)
    return vars._v[vars._kappa_ind]
end

function kappa(vars::class_variables, new_vals::Float64)
    return vars._v[vars._kappa_ind] = new_vals;
end

function x_scaled(vars::class_variables)
     return x(vars)/tau(vars)
end

function y_scaled(vars::class_variables)
     return y(vars)/tau(vars)
end

function n(vars::class_variables)
    return length(vars._x_ind)
end

function m(vars::class_variables)
    return length(vars._y_ind)
end

function mu(vars::class_variables)
    return mean([x(vars) .* s(vars); tau(vars) * kappa(vars)])
end

# point = point + alpha * direction
function move!(vars::class_variables, alpha::Float64, direction::class_variables)
    try
      vars._v = vars._v + alpha * direction._v;
      check_positive(vars);
    catch e
        println("ERROR class_variables/move!")
        throw(e)
    end
end

#function ipopt_start()
#	 this.v[this.tau_ind] = 1.0;
#	 this.v[this.kappa_ind] = 1e-20;
#end

# check for errors

# check x > 0
function check_positive(vars::class_variables)
      try
        for x_i in x(vars)
           @assert(x_i > 0)
        end

        for s_i in s(vars)
           @assert(s_i > 0)
        end

        @assert(tau(vars) > 0)
        @assert(kappa(vars) > 0)

    catch e
        println("ERROR in check_positive")
        throw(e)
    end

     check_for_wrong_vals(vars)
end

function check_for_wrong_vals(vars::class_variables)
      try
          # check for infinity
          for x_i in x(vars)
            @assert(x_i != -Inf)
            @assert(x_i != Inf)
          end

          for s_i in s(vars)
            @assert(s_i != -Inf)
            @assert(s_i != Inf)
          end

          @assert(tau(vars) != -Inf)
          @assert(tau(vars) != Inf)
          @assert(kappa(vars) != -Inf)
          @assert(kappa(vars) != Inf)

          # check for NaNs
          for x_i in x(vars)
            @assert(!isnan(x_i))
          end

          for s_i in s(vars)
            @assert(!isnan(s_i))
          end

          @assert(!isnan(tau(vars)))
          @assert(!isnan(kappa(vars)))
      catch e
          println("ERROR in check_for_wrong_vals")
          throw(e)
      end
end

max_forks_repo_name: ohinder/Everest


text:
jl_name2id(m::jlModel, objtype, name::String) = Int(MJCore.mj_name2id(m, objtype, name) + 1)
jl_id2name(m::jlModel, objtype, id::Integer) = MJCore.mj_id2name(m, objtype, Int32(id - 1))

jl_enabled(m::jlModel, flag::MJCore.mjtEnableBit) = _flagtest(m.opt.enableflags, flag)
jl_disabled(m::jlModel, flag::Union{Integer, MJCore.mjtDisableBit}) = _flagtest(m.opt.disableflags, flag)
_flagtest(flags, flag) = Integer(flags) & Integer(flag) == Integer(flag)