text:
module jlfmt

import Pkg
using ArgParse: @add_arg_table, ArgParseSettings, parse_args, usage_string
using Base: UUID
using JuliaFormatter: JuliaFormatter, format, format_text

function exc_handler(settings, err)
    error(string(err.text, "\n", usage_string(settings)))
end

function parse_settings()
    description = """
    A CLI to run `JuliaFormatter.format`.
    """
    setting = ArgParseSettings(
        prog = "jlfmt",
        description = description,
        exit_after_help = false,
        exc_handler = exc_handler,
    )

    #! format: off
    @add_arg_table setting begin
        "--version"
            dest_name = "version"
            action = :store_true
            help = """
            Print version.
            """
        "--preview"
            dest_name = "overwrite"
            action = :store_false
            help = """
            Set `overwrite=false`.  If it is specified, the formatted
            version of file named `foo.jl` will be written to
            `foo_fmt.jl`.
            """
        "--verbose"
            action = :store_true
            help = """
            Print the name of the file being formatted along with
            relevant details to `stdout`.
            """
        "--always-for-in"
            dest_name = "always_for_in"
            action = :store_true
            help = """
            Always use `in` keyword for `for` loops.
            """
        "--diff"
            action = :store_true
            help = """
            Show diff using `colordiff` or `diff` command instead of writing
            the result to the original files.
            """
        "paths"
            nargs = '*'
            help = """
            Paths to be formatted.  Recursively search Julia files if
            directories are passed.  Read content from stdin if no
            path or a single `-` is specified.
            """
    end
    #! format: on
end

function preprocess_indents(text)
    lines = split(text, "\n")
    lens = [length(match(r" *", l).match) for l in lines if match(r"^ *$", l) === nothing]
    isempty(lens) && return text, 0
    indent = minimum(lens)
    return join((chop(line, head = indent, tail = 0) for line in lines), "\n"), indent
end

function print_diff(a, b, opts = `--unified`; path = nothing)
    prog = something(Sys.which.(("colordiff", "diff"))...)
    if path !== nothing
        opts = `--label a/$path --label b/$path $opts`
    end
    mktempdir() do dir
        pa = joinpath(dir, "a")
        pb = joinpath(dir, "b")
        cmd = `$prog $opts $pa $pb`
        run(`mkfifo $pa $pb`)
        local proc
        @sync begin
            proc = run(pipeline(cmd; stdout = stdout, stderr = stderr); wait = false)
            @async write(pa, a)
            @async write(pb, b)
        end
        wait(proc)
        @assert fetch(proc).exitcode in (0, 1)
        return
    end
end

function print_version(io::IO = stdout)
    deps = Pkg.dependencies()
    info = (
        jlfmt = get(deps, UUID("f15eac7f-9c89-472e-a5bb-5f47caaa2526"), nothing),
        JuliaFormatter = get(deps, UUID("98e50ef6-434e-11e9-1051-2b60c6c9e899"), nothing),
        ArgParse = get(deps, UUID("c7e460c6-2fb9-53a9-8c5b-16f535851c63"), nothing),
    )
    for (k, pkginfo) in pairs(info)
        v = if pkginfo === nothing
            "<unknown version>"
        else
            pkginfo.version
        end
        println(io, k, " ", v)
    end
    println(io, "julia ", VERSION)
    println(io, "in project: ", Base.active_project())
end

function main(args = ARGS)
    arguments = parse_args(args, parse_settings())
    if arguments !== nothing
        if pop!(arguments, "version")
            print_version()
            return
        end

        paths = pop!(arguments, "paths")
        show_diff = pop!(arguments, "diff")
        kwargs = Dict(Symbol(k) => v for (k, v) in arguments)

        if isempty(paths) || paths == ["-"]
            pop!(kwargs, :verbose, nothing)
            pop!(kwargs, :overwrite, nothing)
            text, indent = preprocess_indents(read(stdin, String))
            formatted = format_text(text; margin = 92 - indent, kwargs...)
            show_diff && return print_diff(text, formatted; path = "-")
            lines = (isempty(l) ? l : string(" "^indent, l) for l in split(formatted, "\n"))
            print(join(lines, "\n"))
            return
        end

        if show_diff
            pop!(kwargs, :verbose, nothing)
            pop!(kwargs, :overwrite, nothing)
            for p in paths
                text = read(p, String)
                formatted = format_text(text; kwargs...)
                print_diff(text, formatted; path = p)
            end
            return
        end

        format(paths; kwargs...)
    end
end

end # module

max_stars_repo_name: tkf/jlfmt

max_issues_repo_name: tkf/jlfmt

max_forks_repo_name: tkf/jlfmt


text:
# -*- coding: utf-8 -*-

import FirebirdDB

println("Testing FirebirdDB")
fb = FirebirdDB.connect("testdb1", "sysdba", "password")
println(fb != nothing ? fb.cn :  fb)
@test fb != nothing
@test FirebirdDB.in_transaction(fb, "drop table ttest;")
@test FirebirdDB.in_transaction(fb,
  "create table ttest (id integer, c1 varchar(10), c2 varchar(20));")
@test FirebirdDB.in_transaction(fb,
  "insert into ttest (id, c1, c2) values (1, 'ABC', 'XYZ');")
@test FirebirdDB.in_transaction(fb,
  "insert into ttest (id, c1, c2) values (2, '222', '000');")
@test FirebirdDB.in_transaction(fb,
  "insert into ttest (id, c1, c2) values (4, '333', 'xyz');")
@test FirebirdDB.in_transaction(fb,
  "insert into ttest (id, c1, c2) values (11, 'ABC', 'XYZ');")
@test FirebirdDB.in_transaction(fb,
  "insert into ttest (id, c1, c2) values (22, '222', '000');")
@test FirebirdDB.in_transaction(fb,
  "delete from ttest where c2 in ('abc', 'def', 'xyz');")
@test FirebirdDB.close(fb)
println(fb.cn)
println("ok")

ext: jl


text:
# # Specular Canopy Scattering

# Using packages:
using Plots, Distributions#, PyPlot
using CanopyOptics
pyplot()
theme(:ggplot2)

# Create a specular model with refractive index and "roughness" κ
specularMod = CanopyOptics.SpecularCanopyScattering(nᵣ=1.5, κ=0.2)

# Use standard planophile leaf distribution:
LD = CanopyOptics.planophile_leaves2()

# Create some discretized angles:
# Azimuth (in radians)
ϕ = range(0.0, 2π,  length=200);
# Polar angle as cos(Θ)
μ,w = CanopyOptics.gauleg(180,0.0,1.0);
# Create directional vectors:
dirs = [CanopyOptics.dirVector_μ(a,b) for a in μ, b in ϕ];

# Compute specular reflectance:
R = CanopyOptics.compute_reflection.([specularMod], [dirs[50,1]], dirs, [LD]);

# Polar plot of reflectance
contourf(ϕ, acos.(μ), R,proj=:polar, ylim=(0,π/2), alpha=0.8)

# ### Animation over different Beta leaf distributions
steps = 0.1:0.1:5
α = [collect(steps); 5*ones(length(steps))]
β = [ 5*ones(length(steps)); reverse(collect(steps));]
x = 0:0.01:1
l = @layout [a  b; c d; e f]

anim = @animate for i ∈ eachindex(α)
    LD = CanopyOptics.LeafDistribution(Beta(α[i],β[i]), 2/π)
    R = CanopyOptics.compute_reflection.([specularMod], [dirs[120,1]], dirs, [LD])
    pZ = []
    for m=0:3
        Z⁺⁺, Z⁻⁺ = CanopyOptics.compute_Z_matrices(specularMod, μ, LD, m)
        push!(pZ, contourf(μ, μ, Z⁻⁺, title="Z⁻⁺ Fourier moment $m", xlabel="μꜜ", aspect_ratio = :equal, ylabel="μꜛ"))
    end
    p0 = plot(rad2deg.(π * x/2), pdf.(LD.LD,x), legend=false, ylim=(0,3), title="Leaf angle distribution", xlabel="Θ (degrees)")
    p1 = contourf(ϕ, acos.(μ), R, proj=:polar, ylim=(0,π/2), clims=(0,0.02), alpha=0.8,legend = :none, title="Polar plot R")
    plot(p0, p1, pZ..., layout = l, margin=5Plots.mm)
    plot!(size=(700,700))
end

gif(anim, "anim_fps10.gif", fps = 10)


text:
import Plots
using LoopVectorization

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

@views function SolutionFields_p(mm, mc, rc, er, gr, x, y, in)
    # Outside
    pmat   = (4*mm*rc^2*(mm - mc)*(er*x^2 + gr*x*y - er*y^2))/((x^2 + y^2)^2*(mm + mc))
    # Inside
    pinc   = 0.0
    # Where it's needed
    p      = (in==0)*pmat + (in==1)*pinc
    return p
end

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

@views function SolutionFields_v(mm, mc, rc, er, gr, x, y, in)
    # Outside
    A = mm*(mc-mm)/(mc+mm);
    v1mat = (2*A*er*rc^4*x^3 + 3*A*gr*rc^4*x^2*y - 6*A*er*rc^4*x*y^2 - A*gr*rc^4*y^3 - 4*A*er*rc^2*x^5 -
        4*A*gr*rc^2*x^4*y - 4*A*gr*rc^2*x^2*y^3 + 4*A*er*rc^2*x*y^4 + 2*er*mm*x^7 + 2*gr*mm*x^6*y +
        6*er*mm*x^5*y^2 + 6*gr*mm*x^4*y^3 + 6*er*mm*x^3*y^4 + 6*gr*mm*x^2*y^5 + 2*er*mm*x*y^6 +
        2*gr*mm*y^7)/(2*mm*(x^2 + y^2)^3)
    v2mat = -(A*gr*rc^4*x^3 - 6*A*er*rc^4*x^2*y - 3*A*gr*rc^4*x*y^2 + 2*A*er*rc^4*y^3 + 4*A*er*rc^2*x^4*y +
        4*A*gr*rc^2*x^3*y^2 + 4*A*gr*rc^2*x*y^4 - 4*A*er*rc^2*y^5 + 2*er*mm*x^6*y + 6*er*mm*x^4*y^3 +
        6*er*mm*x^2*y^5 + 2*er*mm*y^7)/(2*mm*(x^2 + y^2)^3)
    # Inside
    v1inc  = (4*er*mm*x + 3*gr*mm*y + gr*mc*y)/(2*(mm + mc));
    v2inc  = -(4*er*mm*y - gr*mm*x + gr*mc*x)/(2*(mm + mc));
    # Where it's needed
    vx      = (in==0)*v1mat + (in==1)*v1inc
    vy      = (in==0)*v2mat + (in==1)*v2inc
    return vx, vy
end

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

@views function SolutionFields_dv(mm, mc, rc, er, gr, x, y, in)
    # Outside
    A = mm*(mc-mm)/(mc+mm);
    dvxdxmat = (- 3*A*er*rc^4*x^4 - 6*A*gr*rc^4*x^3*y + 18*A*er*rc^4*x^2*y^2 + 6*A*gr*rc^4*x*y^3 - 3*A*er*rc^4*y^4 + 2*A*er*rc^2*x^6 + 4*A*gr*rc^2*x^5*y - 10*A*er*rc^2*x^4*y^2 - 10*A*er*rc^2*x^2*y^4 - 4*A*gr*rc^2*x*y^5 + 2*A*er*rc^2*y^6 + er*mm*x^8 + 4*er*mm*x^6*y^2 + 6*er*mm*x^4*y^4 + 4*er*mm*x^2*y^6 + er*mm*y^8)/(mm*(x^2 + y^2)^4);
    dvxdymat = (3*A*gr*rc^4*x^4 - 24*A*er*rc^4*x^3*y - 18*A*gr*rc^4*x^2*y^2 + 24*A*er*rc^4*x*y^3 + 3*A*gr*rc^4*y^4 - 4*A*gr*rc^2*x^6 + 24*A*er*rc^2*x^5*y + 8*A*gr*rc^2*x^4*y^2 + 16*A*er*rc^2*x^3*y^3 + 12*A*gr*rc^2*x^2*y^4 - 8*A*er*rc^2*x*y^5 + 2*gr*mm*x^8 + 8*gr*mm*x^6*y^2 + 12*gr*mm*x^4*y^4 + 8*gr*mm*x^2*y^6 + 2*gr*mm*y^8)/(2*mm*(x^2 + y^2)^4);
    dvydxmat = (A*rc^2*(3*gr*rc^2*x^4 - 24*er*rc^2*x^3*y - 18*gr*rc^2*x^2*y^2 + 24*er*rc^2*x*y^3 + 3*gr*rc^2*y^4 + 8*er*x^5*y + 12*gr*x^4*y^2 - 16*er*x^3*y^3 + 8*gr*x^2*y^4 - 24*er*x*y^5 - 4*gr*y^6))/(2*mm*(x^2 + y^2)^4);
    dvydymat = -(- 3*A*er*rc^4*x^4 - 6*A*gr*rc^4*x^3*y + 18*A*er*rc^4*x^2*y^2 + 6*A*gr*rc^4*x*y^3 - 3*A*er*rc^4*y^4 + 2*A*er*rc^2*x^6 + 4*A*gr*rc^2*x^5*y - 10*A*er*rc^2*x^4*y^2 - 10*A*er*rc^2*x^2*y^4 - 4*A*gr*rc^2*x*y^5 + 2*A*er*rc^2*y^6 + er*mm*x^8 + 4*er*mm*x^6*y^2 + 6*er*mm*x^4*y^4 + 4*er*mm*x^2*y^6 + er*mm*y^8)/(mm*(x^2 + y^2)^4);
    # Inside
    dvxdxinc = (2*er*mm)/(mm + mc);
    dvxdyinc = (gr*(3*mm + mc))/(2*(mm + mc));
    dvydxinc = (gr*(mm - mc))/(2*(mm + mc));
    dvydyinc = -(2*er*mm)/(mm + mc);
    # Where it's needed
    dvxdx      = (in==0)*dvxdxmat + (in==1)*dvxdxinc
    dvxdy      = (in==0)*dvxdymat + (in==1)*dvxdyinc
    dvydx      = (in==0)*dvydxmat + (in==1)*dvydxinc
    dvydy      = (in==0)*dvydymat + (in==1)*dvydyinc
    return dvxdx, dvxdy, dvydx, dvydy
end

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

function Evaluation2D!( ncx, ncy, xc, yc, p, vx, vy, rc, mm, mc )
    gr  = 0;                        # Simple shear: gr=1, er=0
    er  = -1;                       # Strain rate
    @tturbo for i=1:ncx
        for j=1:ncy
            x, y                           = xc[i], yc[j]
            in                             = sqrt(x^2.0 + y^2.0)<=rc
            pa                             = SolutionFields_p( mm, mc, rc, er, gr, x, y, in)
            vxa, vya                       = SolutionFields_v( mm, mc, rc, er, gr, x, y, in)
            dvxdxa, dvxdya, dvydxa, dvydya = SolutionFields_dv(mm, mc, rc, er, gr, x, y, in)
            p[i,j], vx[i,j], vy[i,j]       = pa, vxa, vya
        end
    end
end

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

@views function main()

    # Create domain
    xmin = -3.0
    xmax = 3.0
    ymin = -3.0
    ymax = 3.0
    ncx  = 50
    ncy  = 40
    xc   = LinRange(xmin,xmax,ncx)
    yc   = LinRange(ymin,ymax,ncy)
    xc2d = repeat(xc, 1, length(yc))
    yc2d = repeat(yc, 1, length(xc))'
    p    = zeros(ncx,ncy)
    vx   = zeros(ncx,ncy)
    vy   = zeros(ncx,ncy)
    rc   = 1.0
    mm   = 1.0
    mc   = 100.0

    for itest=1:5
        @time Evaluation2D!( ncx, ncy, xc, yc, p, vx, vy, rc, mm, mc )
    end

    p2 = Plots.heatmap(xc, yc, p', c=:jet1 )
    display(Plots.plot(p2))

end

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

#main()


text:
using OffsetArrays, NCDatasets, Parameters

# Define the contents of the State Struct
# The @with_kw "macro" is used to "unpack" the variables when the struct is passed as an argument to functions
@with_kw struct State{T, S}
    isd        ::Int32
    ied        ::Int32
    jsd        ::Int32
    jed        ::Int32
    is         ::Int32
    ie         ::Int32
    js         ::Int32
    je         ::Int32
    nord       ::Int32
    npx        ::Int32
    npy        ::Int32
    npz        ::Int32
    dt2        ::Float64
    sw_corner  ::Int32
    se_corner  ::Int32
    ne_corner  ::Int32
    nw_corner  ::Int32

    ni         ::Int64
    nj         ::Int64
    nip1       ::Int64
    njp1       ::Int64
    nk         ::Int64
    nk9        ::Int64

    rarea      ::T
    rarea_c    ::T
    sin_sg     ::S
    cos_sg     ::S
    sina_v     ::T
    cosa_v     ::T
    sina_u     ::T
    cosa_u     ::T
    fC         ::T
    rdxc       ::T
    rdyc       ::T
    dx         ::T
    dy         ::T
    dxc        ::T
    dyc        ::T
    cosa_s     ::T
    rsin_u     ::T
    rsin_v     ::T
    rsin2      ::T
    dxa        ::T
    dya        ::T
    delpc      ::S
    delp       ::S
    ptc        ::S
    pt         ::S
    u          ::S
    v          ::S
    w          ::S
    uc         ::S
    vc         ::S
    ua         ::S
    va         ::S
    wc         ::S
    ut         ::S
    vt         ::S
    divg_d     ::S
end

# Assign NetCDF variables,dimensions and attributes to the struct
function State(ds::NCDataset)

    isd =                 ds.attrib["isd"]
    ied =                 ds.attrib["ied"]
    jsd =                 ds.attrib["jsd"]
    jed =                 ds.attrib["jed"]
    is =                  ds.attrib["is"]
    ie =                  ds.attrib["ie"]
    js =                  ds.attrib["js"]
    je =                  ds.attrib["je"]
    nord =                ds.attrib["nord"]
    npx =                 ds.attrib["npx"]
    npy =                 ds.attrib["npy"]
    npz =                 ds.attrib["npz"]
    dt2 =                 ds.attrib["dt2"]
    sw_corner =           ds.attrib["sw_corner"]
    se_corner =           ds.attrib["se_corner"]
    nw_corner =           ds.attrib["nw_corner"]
    ne_corner =           ds.attrib["ne_corner"]
    ni =                  ds.dim["ni"]
    nj =                  ds.dim["nj"]
    nip1 =                ds.dim["nip1"]
    njp1 =                ds.dim["njp1"]
    nk =                  ds.dim["nk"]
    nk9 =                 ds.dim["nk9"]

    # Assign all of the NetCDF variables to a constant 'State' Struct
    return State(
      isd,
      ied,
      jsd,
      jed,
      is,
      ie,
      js,
      je,
      nord,
      npx,
      npy,
      npz,
      dt2,
      sw_corner,
      se_corner,
      nw_corner,
      ne_corner,
      ni,
      nj,
      nip1,
      njp1,
      nk,
      nk9,

      # Assign OffsetArrays
      # The first argument to OffsetArray is the NetCDF variable loaded in memory as a Julia array
      # The second, third (and fourth) args are the non-standardard indexed dimensions
      OffsetArray(Array(ds["rarea"]),        isd:ied,      jsd:jed     ),
      OffsetArray(Array(ds["rarea_c"]),      isd:ied + 1,  jsd:jed + 1 ),
      OffsetArray(Array(ds["sin_sg"]),       isd:ied,      jsd:jed, 1:9),
      OffsetArray(Array(ds["cos_sg"]),       isd:ied,      jsd:jed, 1:9),
      OffsetArray(Array(ds["sina_v"]),       isd:ied,      jsd:jed + 1 ),
      OffsetArray(Array(ds["cosa_v"]),       isd:ied,      jsd:jed + 1 ),
      OffsetArray(Array(ds["sina_u"]),       isd:ied + 1,  jsd:jed     ),
      OffsetArray(Array(ds["cosa_u"]),       isd:ied + 1,  jsd:jed     ),
      OffsetArray(Array(ds["fC"]),           isd:ied + 1,  jsd:jed + 1 ),
      OffsetArray(Array(ds["rdxc"]),         isd:ied + 1,  jsd:jed     ),
      OffsetArray(Array(ds["rdyc"]),         isd:ied,      jsd:jed + 1 ),
      OffsetArray(Array(ds["dx"]),           isd:ied,      jsd:jed + 1 ),
      OffsetArray(Array(ds["dy"]),           isd:ied + 1,  jsd:jed     ),
      OffsetArray(Array(ds["dxc"]),          isd:ied + 1,  jsd:jed     ),
      OffsetArray(Array(ds["dyc"]),          isd:ied,      jsd:jed + 1 ),
      OffsetArray(Array(ds["cosa_s"]),       isd:ied,      jsd:jed     ),
      OffsetArray(Array(ds["rsin_u"]),       isd:ied + 1,  jsd:jed     ),
      OffsetArray(Array(ds["rsin_v"]),       isd:ied,      jsd:jed + 1 ),
      OffsetArray(Array(ds["rsin2"]),        isd:ied,      jsd:jed     ),
      OffsetArray(Array(ds["dxa"]),          isd:ied,      jsd:jed     ),
      OffsetArray(Array(ds["dya"]),          isd:ied,      jsd:jed     ),
      OffsetArray(Array(ds["delpc"]),        isd:ied,      jsd:jed,     1:npz),
      OffsetArray(Array(ds["delp"]),         isd:ied,      jsd:jed,     1:npz),
      OffsetArray(Array(ds["ptc"]),          isd:ied,      jsd:jed,     1:npz),
      OffsetArray(Array(ds["pt"]),           isd:ied,      jsd:jed,     1:npz),
      OffsetArray(Array(ds["u"]),            isd:ied,      jsd:jed + 1, 1:npz),
      OffsetArray(Array(ds["v"]),            isd:ied + 1,  jsd:jed,     1:npz),
      OffsetArray(Array(ds["w"]),            isd:ied,      jsd:jed,     1:npz),
      OffsetArray(Array(ds["uc"]),           isd:ied + 1,  jsd:jed,     1:npz),
      OffsetArray(Array(ds["vc"]),           isd:ied,      jsd:jed + 1, 1:npz),
      OffsetArray(Array(ds["ua"]),           isd:ied,      jsd:jed,     1:npz),
      OffsetArray(Array(ds["va"]),           isd:ied,      jsd:jed,     1:npz),
      OffsetArray(Array(ds["wc"]),           isd:ied,      jsd:jed,     1:npz),
      OffsetArray(Array(ds["ut"]),           isd:ied,      jsd:jed,     1:npz),
      OffsetArray(Array(ds["vt"]),           isd:ied,      jsd:jed,     1:npz),
      OffsetArray(Array(ds["divg_d"]),       isd:ied + 1,  jsd:jed + 1, 1:npz)
    )
end


text:
include(joinpath(pwd(), "src/controller/newton_structure_solver/methods.jl"))

sim = get_simulation("flamingo", "flat_2D_lc", "flat")
model = sim.model
env = sim.env

ref_traj = deepcopy(ContactImplicitMPC.get_trajectory(model, flat_2D_lc,
    joinpath(module_dir(), "src/dynamics/flamingo/gaits/gait_forward_36_4.jld2"),
    load_type = :split_traj_alt))

H_mpc = 15
# obj_mpc = quadratic_objective(model, H_mpc,
#     q = [Diagonal(1e-1 * [3e2, 1e-6, 3e2, 1, 1, 1, 1, 0.1, 0.1]) for t = 1:H_mpc+2],
#     v = [0.0 * Diagonal([1e0,1,1e4,1,1,1,1,1e4,1e4]) for t = 1:H_mpc],
#     u = [Diagonal(3e-1 * [0.1; 0.1; 0.3; 0.3; ones(model.nu-6); 2; 2]) for t = 1:H_mpc-1])

obj_mpc = quadratic_objective(model, H_mpc,
    q = [Diagonal(0.1 * ones(model.nq)) for t = 1:H_mpc+2],
    v = [Diagonal(0.001 * ones(model.nq)) for t = 1:H_mpc],
    u = [Diagonal(0.01 * ones(model.nu)) for t = 1:H_mpc-1])

s = NewtonStructure(sim, H_mpc, ref_traj, obj_mpc, 1.0e-4)
ref_traj.θ
ip_opts = InteriorPointOptions(
			κ_init = 1.0e-4,
			κ_tol = 2.0 * 1.0e-4,
			r_tol = 1.0e-8,
			diff_sol = true,
			solver = :empty_solver)

im_traj = ImplicitTrajectory(ref_traj, sim,
	κ = 1.0e-4,
	mode = :configuration,
	opts=ip_opts)

initialize_trajectories!(s, ref_traj, warm_start_duals = false,
	ref_traj.q[1], ref_traj.q[2])

# compute residual
compute_residual!(s, s.u, s.qa, s.qb, s.ν1, s.ν2, s.H, im_traj, sim.model, sim.env, 1.0e-4)
@show r_norm = residual_norm(s, 1)

# factorize system
factorize!(s)

# solve system
ContactImplicitMPC.solve!(s)

# line search the step direction
α = 1.0
iter = 0
step!(s, α)

# compute candidate residual
compute_residual!(s, s.u_cand, s.qa_cand, s.qb_cand, s.ν1_cand, s.ν2_cand, s.H, im_traj, sim.model, sim.env, 1.0e-4)
residual_norm(s, 1)

newton_solve!(s, sim, im_traj, ref_traj, q0 = ref_traj.q[1], q1 = ref_traj.q[2], warm_start = false, κ = 1.0e-4)

# update!(im_traj, ref_traj_copy, sim, zeros(sim.model.nc), κ = 1.0e-4) #@@@ only keep the rotation stuff not the altitude update.
#
# stride = get_stride(sim.model, ref_traj_copy)
# rot_n_stride!(ref_traj_copy, stride)


text:
#=

@authors: Spencer Lyon, Tom Sargent, John Stachurski

=#

using QuantEcon

# A default Markov chain for the state process
ρ = 0.9
σ = 0.02
n = 25
default_mc = tauchen(n, ρ, σ)

mutable struct AssetPriceModel{TF<:AbstractFloat, TI<:Integer}
    β :: TF            # Discount factor
    γ :: TF            # Coefficient of risk aversion
    mc :: MarkovChain  # State process
    n :: TI            # Number of states
    g :: Function      # Function mapping states into growth rates
end

function AssetPriceModel(;β::AbstractFloat=0.96, γ::AbstractFloat=2.0,
                          mc::MarkovChain=default_mc, g::Function=exp)
    n = size(mc.p)[1]
    return AssetPriceModel(β, γ, mc, n, g)
end

"""
Stability test for a given matrix Q.
"""
function test_stability(ap::AssetPriceModel, Q::Matrix)
    sr = maximum(abs, eigvals(Q))
    if sr >= 1 / ap.β
        msg = "Spectral radius condition failed with radius = $sr"
        throw(ArgumentError(msg))
    end
end

"""
Computes the price-dividend ratio of the Lucas tree.

"""
function tree_price(ap::AssetPriceModel)
    # == Simplify names, set up matrices  == #
    β, γ, P, y = ap.β, ap.γ, ap.mc.p, ap.mc.state_values
    y = reshape(y, 1, ap.n)
    J = P .* ap.g.(y).^(1 - γ)

    # == Make sure that a unique solution exists == #
    test_stability(ap, J)

    # == Compute v == #
    I = eye(ap.n)
    Ones = ones(ap.n)
    v = (I - β * J) \ (β * J * Ones)

    return v
end

"""
Computes price of a consol bond with payoff ζ

"""
function consol_price(ap::AssetPriceModel, ζ::AbstractFloat)
    # == Simplify names, set up matrices  == #
    β, γ, P, y = ap.β, ap.γ, ap.mc.p, ap.mc.state_values
    y = reshape(y, 1, ap.n)
    M = P .* ap.g.(y).^(-γ)

    # == Make sure that a unique solution exists == #
    test_stability(ap, M)

    # == Compute price == #
    I = eye(ap.n)
    Ones = ones(ap.n)
    p = (I - β * M) \ ( β * ζ * M * Ones)

    return p
end

"""
Computes price of a perpetual call option on a consol bond.

"""
function call_option(ap::AssetPriceModel, ζ::AbstractFloat, p_s::AbstractFloat, ϵ=1e-7)

    # == Simplify names, set up matrices  == #
    β, γ, P, y = ap.β, ap.γ, ap.mc.p, ap.mc.state_values
    y = reshape(y, 1, ap.n)
    M = P .* ap.g.(y).^(-γ)

    # == Make sure that a unique console price exists == #
    test_stability(ap, M)

    # == Compute option price == #
    p = consol_price(ap, ζ)
    w = zeros(ap.n, 1)
    error = ϵ + 1
    while (error > ϵ)
        # == Maximize across columns == #
        w_new = max.(β * M * w, p - p_s)
        # == Find maximal difference of each component and update == #
        error = maximum(abs, w - w_new)
        w = w_new
    end

    return w
end


text:
"""
Struct for Asian Fixed Strike Option

		asOption=AsianFixedStrikeOption(T::num1,K::num2,isCall::Bool=true) where {num1 <: Number,num2 <: Number}

Where:\n
		T	=	Time to maturity of the Option.
		K	=	Strike Price of the Option.
		isCall  = true for CALL, false for PUT.
"""
struct AsianFixedStrikeOption{num1 <: Number, num2 <: Number, numtype <: Number} <: AsianPayoff{numtype}
    T::num1
    K::num2
    isCall::Bool
    function AsianFixedStrikeOption(T::num1, K::num2, isCall::Bool = true) where {num1 <: Number, num2 <: Number}
        if T <= 0.0
            error("Time to Maturity must be positive")
        elseif K <= 0.0
            error("Strike Price must be positive")
        else
            zero_typed = zero(num1) + zero(num2)
            return new{num1, num2, typeof(zero_typed)}(T, K, isCall)
        end
    end
end

export AsianFixedStrikeOption;

function payout(S::abstractArray, opt_::AsianFixedStrikeOption) where {abstractArray <: AbstractArray{num_}} where {num_ <: Number}
    iscall = opt_.isCall ? 1 : -1
    zero_typed = zero(eltype(S)) * opt_.K
    return max(iscall * (mean(S) - opt_.K), zero_typed)
end


text:
using ITensors
using Test

@testset "Test debug checks on IndexSet construction" begin
  i = Index(2, "i")

  initially_using_debug_checks = ITensors.using_debug_checks()

  ITensors.disable_debug_checks()
  @test !ITensors.using_debug_checks()
  # Test that no error is thrown in constructor
  @test ITensor(i, i) isa ITensor
  @test ITensor(i, i') isa ITensor

  # Turn on debug checks
  ITensors.enable_debug_checks()
  @test ITensors.using_debug_checks()
  @test_throws ErrorException ITensor(i, i)
  # Test that no error is thrown in constructor
  @test ITensor(i, i') isa ITensor

  # Turn off debug checks
  ITensors.disable_debug_checks()
  @test !ITensors.using_debug_checks()
  # Test that no error is thrown in constructor
  @test ITensor(i, i) isa ITensor
  @test ITensor(i, i') isa ITensor

  # Reset to the initial value
  if !initially_using_debug_checks
    ITensors.disable_debug_checks()
  else
    ITensors.enable_debug_checks()
  end
  @test ITensors.using_debug_checks() == initially_using_debug_checks
end


text:
module DayNNN
export dayNNN

function parseLine(line)
    line
end

function readInput(inputfile="inputs/inputNNN.txt")
    lines = readlines(inputfile)
    lines
    # [parseLine(line) for line in eachline(inputfile)]
end

function dayNNN(inputfile="inputs/inputNNN.txt")
    values = readInput(inputfile)
    (dayNNNPart1(values), dayNNNPart2(values))
end

function dayNNNPart1(values)

end

function dayNNNPart2(values)

end

end #module


text:
using Falcons
using Healpix
using Test

day = 60 * 60 * 24
year = day * 365
prec1 = 60 * 60 * 3

ss = gen_ScanningStrategy(duration=day)

@testset "ScanningStrategy_structure-Test" begin
    @test typeof(ss) <: ScanningStrategy
    @show fieldnames(ScanningStrategy)
    println("Test ScanningStrategy ==> ", ss)
end
println()
@testset "get_pointings_tuple-Test" begin
    pointings = get_pointings_tuple(ss, 0, 30)
    @test typeof(pointings[1]) <: Matrix
    @test typeof(pointings[2]) <: Matrix
    @test typeof(pointings[3]) <: Matrix
    @test typeof(pointings[4]) <: StepRangeLen
    #@show pointings[1][1:10]
    #@show pointings[2][1:10]
    #@show pointings[3][1:10]
    #@show pointings[4][1:10]
end
println()
@testset "get_pointing_pixels-Test" begin
    pointings = get_pointing_pixels(ss, 0, 30)
    @test typeof(pointings[1]) <: Matrix{Int64}
    @test typeof(pointings[2]) <: Matrix{Float64}
    @test typeof(pointings[3]) <: StepRangeLen
    #@show pointings[1][1:10]
    #@show pointings[2][1:10]
    #@show pointings[3][1:10]
end
println()
@testset "get_pointings-Test" begin
    pointings = get_pointings(ss, 0, 30)
    @test typeof(pointings["theta"]) <: Array
    @test typeof(pointings["phi"]) <: Array
    @test typeof(pointings["psi"]) <: Array
    #@show pointings["theta"][1:10]
    #@show pointings["phi"][1:10]
    #@show pointings["psi"][1:10]
end
println()
@testset "ScanningStrategy2map-Test" begin
    n = 1
    out = ScanningStrategy2map(ss, n)
    @test typeof(out) <: Vector
    #@show outmap[2][1:5]
    #@show outmap[3][1:5]
    #@show outmap[4][1:5]
    #@show outmap[5][1:5]
end
println()
@testset "get_psiDataBase-Test" begin
    DB = get_psiDataBase(ss, division=4, idx=1, map_div=4);
    @test typeof(DB) <: Tuple
end


text:
using CLIMA.VariableTemplates

function vars_diagnostic(FT)
    @vars begin
        # state and functions of state
        u::FT
        v::FT
        w::FT
        q_tot::FT
        e_tot::FT
        q_liq::FT
        thd::FT                # θ
        thl::FT                # θ_liq
        thv::FT                # θ_v
        e_int::FT
        h_m::FT
        h_t::FT
        qt_sgs::FT
        ht_sgs::FT

        # vertical fluxes
        vert_eddy_mass_flux::FT # <w′ρ′>
        vert_eddy_u_flux::FT    # <w′u′>
        vert_eddy_v_flux::FT    # <w′v′>
        vert_eddy_qt_flux::FT   # <w'q_tot'>
        vert_qt_flux::FT        # <w q_tot>
        vert_eddy_ql_flux::FT   # <w′q_liq′>
        vert_eddy_qv_flux::FT   # <w′q_vap′>
        vert_eddy_thd_flux::FT  # <w′θ′>
        vert_eddy_thv_flux::FT  # <w′θ_v′>
        vert_eddy_thl_flux::FT  # <w′θ_liq′>

        # variances
        uvariance::FT          # u′u′
        vvariance::FT          # v′v′
        wvariance::FT          # w′w′

        # skewness
        wskew::FT              # w′w′w′

        # turbulent kinetic energy
        TKE::FT
    end
end

num_diagnostic(FT) = varsize(vars_diagnostic(FT))
diagnostic_vars(array) = Vars{vars_diagnostic(eltype(array))}(array)


text:
import Base: isempty, ∈

export AbstractAffineMap,
       matrix, vector, set

"""
    AbstractAffineMap{N, S<:LazySet{N}} <: LazySet{N}

Abstract type for affine maps.

### Notes

See [`AffineMap`](@ref) for a standard implementation of this interface.

Every concrete `AbstractAffineMap` must define the following functions:
- `matrix(::AbstractAffineMap)` -- return the linear map
- `vector(::AbstractAffineMap)` -- return the affine translation vector
- `set(::AbstractAffineMap)` -- return the set that the map is applied to

```jldoctest; setup = :(using LazySets: subtypes)
julia> subtypes(AbstractAffineMap)
6-element Array{Any,1}:
 AffineMap
 ExponentialMap
 ExponentialProjectionMap
 LinearMap
 ResetMap
 Translation
```
"""
abstract type AbstractAffineMap{N, S<:LazySet{N}} <: LazySet{N} end

isoperationtype(::Type{<:AbstractAffineMap}) = true
isconvextype(::Type{<:AbstractAffineMap{N, S}}) where {N, S} = isconvextype(S)

# --- common AbstractAffineMap functions ---

"""
    dim(am::AbstractAffineMap)

Return the dimension of an affine map.

### Input

- `am` -- affine map

### Output

The dimension of an affine map.
"""
function dim(am::AbstractAffineMap)
    return length(vector(am))
end

"""
    σ(d::AbstractVector, am::AbstractAffineMap)

Return the support vector of an affine map.

### Input

- `d`  -- direction
- `am` -- affine map

### Output

The support vector in the given direction.
"""
function σ(d::AbstractVector, am::AbstractAffineMap)
    A = matrix(am)
    return A * σ(_At_mul_B(A, d), set(am)) + vector(am)
end

"""
    ρ(d::AbstractVector, am::AbstractAffineMap)

Return the support function of an affine map.

### Input

- `d`  -- direction
- `am` -- affine map

### Output

The support function in the given direction.
"""
function ρ(d::AbstractVector, am::AbstractAffineMap)
    return ρ(_At_mul_B(matrix(am), d), set(am)) + dot(d, vector(am))
end

"""
    an_element(am::AbstractAffineMap)

Return some element of an affine map.

### Input

- `am` -- affine map

### Output

An element of the affine map. It relies on the `an_element` function of the
wrapped set.
"""
function an_element(am::AbstractAffineMap)
    return matrix(am) * an_element(set(am)) + vector(am)
end

"""
    isempty(am::AbstractAffineMap)

Return whether an affine map is empty or not.

### Input

- `am` -- affine map

### Output

`true` iff the wrapped set is empty.
"""
function isempty(am::AbstractAffineMap)
    return isempty(set(am))
end

"""
    isbounded(am::AbstractAffineMap; cond_tol::Number=DEFAULT_COND_TOL)

Determine whether an affine map is bounded.

### Input

- `am`       -- affine map
- `cond_tol` -- (optional) tolerance of matrix condition (used to check whether
                the matrix is invertible)

### Output

`true` iff the affine map is bounded.

### Algorithm

We first check if the matrix is zero or the wrapped set is bounded.
If not, we perform a sufficient check whether the matrix is invertible.
If the matrix is invertible, then the map being bounded is equivalent to the
wrapped set being bounded, and hence the map is unbounded.
Otherwise, we check boundedness via [`_isbounded_unit_dimensions`](@ref).
"""
function isbounded(am::AbstractAffineMap; cond_tol::Number=DEFAULT_COND_TOL)
    M = matrix(am)
    if iszero(M) || isbounded(set(am))
        return true
    end
    if isinvertible(M; cond_tol=cond_tol)
        return false
    end
    return _isbounded_unit_dimensions(am)
end

"""
    ∈(x::AbstractVector, am::AbstractAffineMap)

Check whether a given point is contained in the affine map of a convex set.

### Input

- `x`  -- point/vector
- `am` -- affine map of a convex set

### Output

`true` iff ``x ∈ am``.

### Algorithm

Note that ``x ∈ M⋅S ⊕ v`` iff ``M^{-1}⋅(x - v) ∈ S``.
This implementation does not explicitly invert the matrix, which is why it also
works for non-square matrices.

### Examples

```jldoctest
julia> am = AffineMap([2.0 0.0; 0.0 1.0], BallInf([1., 1.], 1.), [-1.0, -1.0]);

julia> [5.0, 1.0] ∈ am
false

julia> [3.0, 1.0] ∈ am
true
```

An example with a non-square matrix:

```jldoctest
julia> B = BallInf(zeros(4), 1.);

julia> M = [1. 0 0 0; 0 1 0 0]/2;

julia> [0.5, 0.5] ∈ M*B
true
```
"""
function ∈(x::AbstractVector, am::AbstractAffineMap)
    return matrix(am) \ (x - vector(am)) ∈ set(am)
end

"""
    vertices_list(am::AbstractAffineMap; [apply_convex_hull]::Bool)

Return the list of vertices of a (polyhedral) affine map.

### Input

- `am`                -- affine map
- `apply_convex_hull` -- (optional, default: `true`) if `true`, apply the convex
                         hull operation to the list of vertices transformed by
                         the affine map

### Output

A list of vertices.

### Algorithm

This implementation computes all vertices of `X`, then transforms them through
the affine map, i.e. `x ↦ M*x + v` for each vertex `x` of `X`. By default, the
convex hull operation is taken before returning this list. For dimensions three
or higher, this operation relies on the functionality through the concrete
polyhedra library `Polyhedra.jl`.

If you are not interested in taking the convex hull of the resulting vertices
under the affine map, pass `apply_convex_hull=false` as a keyword argument.

Note that we assume that the underlying set `X` is polyhedral, either concretely
or lazily, i.e. there the function `vertices_list` should be applicable.
"""
function vertices_list(am::AbstractAffineMap; apply_convex_hull::Bool=true)
    # for a zero linear map, the result is just the affine translation
    A = matrix(am)
    b = vector(am)
    if iszero(A)
        return [b]
    end

    # collect vertices list of the wrapped set
    X = set(am)
    @assert applicable(vertices_list, X)  "this function requires that the " *
        "list of vertices is available, but it is not"
    vlist_X = vertices_list(X)

    # create resulting vertices list
    vlist = [A * x + b for x in vlist_X]

    return apply_convex_hull ? convex_hull!(vlist) : vlist
end

"""
    constraints_list(am::AbstractAffineMap)

Return the list of constraints of a (polyhedral) affine map.

### Input

- `am` -- affine map

### Output

The list of constraints of the affine map.

### Notes

We assume that the underlying set `X` is polyhedral, i.e., offers a method
`constraints_list(X)`.

### Algorithm

Falls back to the list of constraints of the translation of a lazy linear map.
"""
function constraints_list(am::AbstractAffineMap)
    return _constraints_list_translation(LinearMap(matrix(am), set(am)),
                                         vector(am))
end

"""
    linear_map(M::AbstractMatrix, am::AbstractAffineMap)

Return the linear map of a lazy affine map.

### Input

- `M`  -- matrix
- `am` -- affine map

### Output

A set corresponding to the linear map of the lazy affine map of a set.
"""
function linear_map(M::AbstractMatrix, am::AbstractAffineMap)
     return translate(linear_map(M * matrix(am), set(am)), M * vector(am))
end


text:
"Map from PCRaster LDD value to a CartesianIndex"
const pcr_dir = [
    CartesianIndex(-1, -1),  # 1
    CartesianIndex(0, -1),  # 2
    CartesianIndex(1, -1),  # 3
    CartesianIndex(-1, 0),  # 4
    CartesianIndex(0, 0),  # 5
    CartesianIndex(1, 0),  # 6
    CartesianIndex(-1, 1),  # 7
    CartesianIndex(0, 1),  # 8
    CartesianIndex(1, 1),  # 9
]

const mv = Float(NaN)

# timestep that the parameter units are defined in
const basetimestep = Second(Day(1))

"Set at indices pit values (default = 5) in a gridded local drainage direction vector"
function set_pit_ldd(pits_2d, ldd, indices; pit = 5)
    pits = pits_2d[indices]
    index = filter(i -> isequal(pits[i], true), 1:length(indices))
    ldd[index] .= pit
    return ldd
end

"Filter upstream neighbors of graph based on logical vector"
function filter_upsteam_nodes(graph, vec_logical::Vector{Bool})
    upstream_nodes = Vector{Int}[]
    for v in topological_sort_by_dfs(graph)
        ups_nodes = inneighbors(graph, v)
        push!(upstream_nodes, filter(i -> !vec_logical[i], ups_nodes))
    end
    return upstream_nodes
end

"""
    active_indices(subcatch_2d, nodata)

Takes a 2D array of the subcatchments. And derive forward and reverse indices.

1: Get a list of `CartesianIndex{2}`` that are active, based on a nodata value.
These map from the 1D internal domain to the 2D external domain.

2: Make a reverse index, a `Matrix{Int}``, which maps from the 2D external domain to
the 1D internal domain, providing an Int which can be used as a linear index. Values of 0
represent inactive cells.
"""
function active_indices(subcatch_2d::AbstractMatrix, nodata)
    A = subcatch_2d
    all_inds = CartesianIndices(size(A))
    indices = filter(i -> !isequal(A[i], nodata), all_inds)

    reverse_indices = zeros(Int, size(A))
    for (i, I) in enumerate(indices)
        reverse_indices[I] = i
    end

    return indices, reverse_indices
end

function active_indices(network, key::Tuple)
    if :reservoir in key
        return network.reservoir.indices_outlet
    elseif :lake in key
        return network.lake.indices_outlet
    elseif :river in key
        return network.river.indices
    elseif :drain in key
        return network.drain.indices
    else
        return network.land.indices
    end
end

function lattometres(lat::Real)
    m1 = 111132.92     # latitude calculation term 1
    m2 = -559.82       # latitude calculation term 2
    m3 = 1.175         # latitude calculation term 3
    m4 = -0.0023       # latitude calculation term 4
    p1 = 111412.84     # longitude calculation term 1
    p2 = -93.5         # longitude calculation term 2
    p3 = 0.118         # longitude calculation term 3

    # Calculate the length of a degree of latitude and longitude in meters
    latlen = m1 + (m2 * cosd(2.0 * lat)) + (m3 * cosd(4.0 * lat)) + (m4 * cosd(6.0 * lat))
    longlen = (p1 * cosd(lat)) + (p2 * cosd(3.0 * lat)) + (p3 * cosd(5.0 * lat))

    return longlen, latlen
end

function cell_lengths(y::AbstractVector, cellength::Real, sizeinmetres::Bool)
    n = length(y)
    xl = fill(mv, n)
    yl = fill(mv, n)
    if sizeinmetres
        xl .= cellength
        yl .= cellength
    else
        for i = 1:n
            longlen, latlen = lattometres(y[i])
            xl[i] = longlen * cellength
            yl[i] = latlen * cellength
        end
    end
    return xl, yl
end

function river_fraction(
    river::AbstractVector,
    riverlength::AbstractVector,
    riverwidth::AbstractVector,
    xl::AbstractVector,
    yl::AbstractVector,
)
    n = length(river)
    riverfrac = fill(mv, n)
    for i = 1:n
        riverfrac[i] = if river[i]
            min((riverlength[i] * riverwidth[i]) / (xl[i] * yl[i]), 1.0)
        else
            0.0
        end
    end
    return riverfrac
end

"""
    set_states(instate_path, model, state_ncnames; <keyword arguments>)

Read states contained in `Dict` `state_ncnames` from NetCDF file located in `instate_path`,
and set states in `model` object. Active cells are selected with the corresponding network's
(`Vector{CartesianIndex}`) from the NetCDF file.

# Arguments
- `type = nothing`: type to convert data to after reading. By default no conversion is done.
"""
function set_states(instate_path, model, state_ncnames; type = nothing)
    @unpack network = model
    # states in NetCDF include dim time (one value) at index 3 or 4, 3 or 4 dims are allowed
    NCDataset(instate_path) do ds
        for (state, ncname) in state_ncnames
            sel = active_indices(network, state)
            dims = length(dimnames(ds[ncname]))
            # 4 dims, for example (x,y,layer,time) where dim layer is an SVector for soil layers
            if dims == 4
                A = read_standardized(ds, ncname, (x = :, y = :, layer = :, time = 1))
                A = permutedims(A[sel, :])
                # note that this array is allowed to have missings, since not every vertical
                # column is `maxlayers` layers deep
                A = replace!(A, missing => NaN)
                # Convert to desired type if needed
                if !isnothing(type)
                    if eltype(A) != type
                        A = map(type, A)
                    end
                end
                # set state in model object
                param(model, state) .= svectorscopy(A, Val{size(A)[1]}())
                # 3 dims (x,y,time)
            elseif dims == 3
                A = read_standardized(ds, ncname, (x = :, y = :, time = 1))
                A = A[sel]
                A = nomissing(A)
                # Convert to desired type if needed
                if !isnothing(type)
                    if eltype(A) != type
                        A = map(type, A)
                    end
                end
                # set state in model object
                param(model, state) .= A
            else
                error(
                    "Number of state dims should be 3 or 4, number of dims = ",
                    string(dims),
                )
            end
        end
    end
end

"""
    ncread(nc, var; <keyword arguments>)

Read `var` from NetCDF file `nc`, type `String` or type `Config` are allowed.  Supports various
keyword arguments to get selections of data in desired types, with or without missing values.

# Arguments
- `sel=nothing`: a selection of indices, such as a `Vector{CartesianIndex}` of active cells,
        to return from the NetCDF. By default all cells are returned.
- `defaults=nothing`: a default value if `var` is not in `nc`. By default it gives an error
    in this case.
- `type=nothing`: type to convert data to after reading. By default no conversion is done.
- `allow_missing=false`: Missing values within `sel` is not allowed by default. Set to
        `true` to allow missing values.
- `fill=nothing`: Missing values are replaced by this fill value if `allow_missing` is `false`.
- `dimname` : name of third dimension of parameter `var`. By default no third dimension is expected.
"""
function ncread(
    nc,
    var;
    sel = nothing,
    defaults = nothing,
    type = nothing,
    allow_missing = false,
    fill = nothing,
    dimname = nothing,
)

    if isnothing(dimname)
        dim_sel = (x = :, y = :)
    else
        @assert dimname == :layer
        dim_sel = (x = :, y = :, layer = :)
    end

    if isnothing(var)
        @assert !isnothing(defaults)
        if !isnothing(type)
            defaults = convert(type, defaults)
        end
        if isnothing(dimname)
            return Base.fill(defaults, length(sel))
        else
            return Base.fill(defaults, (nc.dim[String(dimname)], length(sel)))
        end
    end

    # If var has type Config, input parameters can be changed (through scale, offset and
    # input NetCDF var) or set to a uniform value (providing a value). Otherwise, input
    # NetCDF var is read directly.
    var, mod = ncvar_name_modifier(var)
    if mod.scale != 1.0 || mod.offset != 0.0
        A = read_standardized(nc, var, dim_sel) .* mod.scale .+ mod.offset
    elseif !isnothing(mod.value)
        if isnothing(dimname)
            return Base.fill(mod.value, length(sel))
        else
            return Base.fill(mod.value, (nc.dim[String(dimname)], length(sel)))
        end
    else
        # Read the entire variable into memory, applying scale, offset and
        # set fill_values to missing.
        A = read_standardized(nc, var, dim_sel)
    end

    # Take out only the active cells
    if !isnothing(sel)
        if isnothing(dimname)
            A = A[sel]
        else
            A = permutedims(A[sel, :])
        end
    end

    if !allow_missing
        if isnothing(fill)
            # errors if missings are found
            A = nomissing(A)
        else
            # replaces missings with a fill value
            A = nomissing(A, fill)
        end
    end

    # Convert to desired type if needed
    if !isnothing(type)
        if eltype(A) != type
            A = convert(Array{type}, A)
        end
    end

    return A
end

"""
    set_layerthickness(d::Real, sl::SVector)

Calculate actual soil thickness of layers based on a reference depth (e.g. soil depth or water table depth) `d`,
a SVector `sl` with cumulative soil depth starting at soil surface (0), and a SVector `tl` with actual thickness
per soil layer.
"""
function set_layerthickness(d::Real, sl::SVector, tl::SVector)

    act_d = tl .* mv
    for i = 1:length(act_d)
        if d > sl[i+1]
            act_d = setindex(act_d, tl[i], i)
        elseif d - sl[i] > 0.0
            act_d = setindex(act_d, d - sl[i], i)
        end
    end

    nlayers = length(act_d) - sum(isnan.(act_d))
    return act_d, nlayers
end

"""
    detdrainlength(ldd, xl, yl)

Determines the drainaige length for a non square grid. Input `ldd` (drainage network), `xl` (length of cells in x direction),
`yl` (length of cells in y direction). Output is drainage length.
"""
function detdrainlength(ldd, xl, yl)
    # take into account non-square cells
    # if ldd is 8 or 2 use ylength
    # if ldd is 4 or 6 use xlength
    if ldd == 2 || ldd == 8
        yl
    elseif ldd == 4 || ldd == 6
        xl
    else
        hypot(xl, yl)
    end
end

"""
    detdrainwidth(ldd, xl, yl)

Determines the drainaige width for a non square grid. Input `ldd` (drainage network), `xl` (length of cells in x direction),
`yl` (length of cells in y direction). Output is drainage width.
"""
function detdrainwidth(ldd, xl, yl)
    # take into account non-square cells
    # if ldd is 8 or 2 use xlength
    # if ldd is 4 or 6 use ylength
    slantwidth = (xl + yl) * 0.5
    if ldd == 2 || ldd == 8
        xl
    elseif ldd == 4 || ldd == 6
        yl
    else
        slantwidth
    end
end

# 2.5x faster power method
"Faster method for exponentiation"
pow(x, y) = exp(y * log(x))

"Return the sum of the array `A` at indices `inds`"
function sum_at(A, inds)
    v = zero(eltype(A))
    for i in inds
        v += A[i]
    end
    return v
end

"Return the sum of the function `f` at indices `inds`"
function sum_at(f::Function, inds, T)
    v = zero(T)
    for i in inds
        v += f(i)
    end
    return v
end

# https://juliaarrays.github.io/StaticArrays.jl/latest/pages/api/#Arrays-of-static-arrays-1
function svectorscopy(x::Matrix{T}, ::Val{N}) where {T,N}
    size(x, 1) == N || error("sizes mismatch")
    isbitstype(T) || error("use for bitstypes only")
    copy(reinterpret(SVector{N,T}, vec(x)))
end

"""
    fraction_runoff_toriver(graph, ldd, index_river, slope, n)

Determine ratio `frac` between `slope` river cell `index_river` and `slope` of each upstream neighbor (based on directed acyclic graph
`graph`).
"""
function fraction_runoff_toriver(graph, ldd, index_river, slope, n)
    frac = zeros(n)
    for i in index_river
        nbs = inneighbors(graph, i)
        for j in nbs
            if ldd[j] != ldd[i]
                frac[j] = slope[j] / (slope[i] + slope[j])
            end
        end
    end
    return frac
end

"""
    equal_size_vectors(x)

Used in the structs of arrays to ensure all vectors are of equal length.

`equal_size_vectors(([1,2], [1,2,3]))` would throw an ArgumentError.
`equal_size_vectors(([4,5], [4,5]))` would pass.
`equal_size_vectors((1, [4,5], [4,5]))` would also pass, since `1` is not an AbstractVector.
"""
function equal_size_vectors(x)
    # all vectors in this struct should be the same size
    inds_vec = findall(arg -> isa(arg, AbstractVector), x)
    n = length(x[inds_vec[1]])
    x_vec = x[inds_vec]

    for arr in x_vec
        if length(arr) != n
            throw(ArgumentError("Not all vectors are of equal length"))
        end
    end
    return x
end

"""
    tosecond(x::Period)

Convert a Period into a Float64, which represents the number of seconds. Will fail if this
is not well defined, such as for Month.

# Examples
```julia-repl
julia> tosecond(Day(1))
86400.0
```
"""
tosecond(x::Hour) = Float64(Dates.value(Second(x)))
tosecond(x::Minute) = Float64(Dates.value(Second(x)))
tosecond(x::T) where {T<:DatePeriod} = Float64(Dates.value(Second(x)))
tosecond(x::T) where {T<:TimePeriod} = x / convert(T, Second(1))


text:
using GLPK

"""
    calculate_optimal_flux_distribution(S,[Lv,Uv],[Lx,Ux],c; min_flag=true)

    Computes the optimal metabolic flux distribution given the constraints.

    Inputs:
    `S` - stoichiometric_matrix (M x R)
    `[Lv,Uv]` - (R x 2) array of the flux lower (L) and upper (U) bounds
    `[Lx,Ux]` - (M x 2) array of species lower (L) and upper (U) bounds (if at steady-state, L = U = 0)
    `c` - R x 1 vector holding indexes for objective vector

    Outputs:
    `objective_value` - value of the objective function at the optimum
    `calculated_flux_array` - R x 1 flux array at the optimum
    `dual_value_array` - R x 1 dual values
    `uptake_array` - M x 1 array of S*v
    `exit_flag` = 0 if optimal
    `status_flag` = 5 if optimal
"""
function calculate_optimal_flux_distribution(stoichiometric_matrix::Array{Float64,2}, default_bounds_array::Array{Float64,2}, species_bounds_array::Array{Float64,2}, objective_coefficient_array::Array{Float64,1}; min_flag::Bool = true)

    # set some constants -
    TIME_RESTART_LIM = 60

    # Get the stoichiometric_matrix from data_dictionary -
    (number_of_species,number_of_fluxes) = size(stoichiometric_matrix);

    # # Setup the GLPK problem -
    lp_problem = GLPK.Prob();
    GLPK.set_prob_name(lp_problem, "sample");
    GLPK.set_obj_name(lp_problem, "objective")

    # Set solver parameters
    solver_parameters = GLPK.SimplexParam();
    GLPK.init_smcp(solver_parameters);
    solver_parameters.msg_lev = GLPK.MSG_OFF;

    # Are we doing min -or- max?
    if min_flag == true
    	GLPK.set_obj_dir(lp_problem, GLPK.MIN);
    else
    	GLPK.set_obj_dir(lp_problem, GLPK.MAX);
    end

    # Set the number of constraints and fluxes -
    GLPK.add_rows(lp_problem, number_of_species);
    GLPK.add_cols(lp_problem, number_of_fluxes);

    # Setup flux bounds, and objective function -
    (number_of_fluxes,number_of_bounds) = size(default_bounds_array)
    for flux_index = 1:number_of_fluxes

    	flux_lower_bound = default_bounds_array[flux_index,1]
    	flux_upper_bound = default_bounds_array[flux_index,2]

    	# Check bounds type ... default is DB -
    	if (flux_upper_bound == flux_lower_bound)
    		flux_constraint_type = GLPK.FX
    	else
    		flux_constraint_type = GLPK.DB
    	end

    	# flux symbol? (later use name - for now, fake it)
    	flux_symbol = "R_"*string(flux_index)

    	# Set the bounds in GLPK -
    	GLPK.set_col_name(lp_problem, flux_index, flux_symbol);
    	GLPK.set_col_bnds(lp_problem, flux_index, flux_constraint_type, flux_lower_bound, flux_upper_bound);
    end

    # Setup objective function -
    for (flux_index,obj_coeff) in enumerate(objective_coefficient_array)

    	# Set the objective function value in GLPK -
    	GLPK.set_obj_coef(lp_problem, flux_index, obj_coeff);
    end

    # Setup problem constraints for the metabolites -
    for species_index = 1:number_of_species

    	species_lower_bound = species_bounds_array[species_index,1]
    	species_upper_bound = species_bounds_array[species_index,2]

    	# defualt
    	species_constraint_type = GLPK.FX
    	if (species_lower_bound != species_upper_bound)
    		species_constraint_type = GLPK.DB
    	end

    	# set the symbol -
    	species_symbol = "x_"*string(species_index)

    	# Set the species bounds in GLPK -
    	GLPK.set_row_name(lp_problem, species_index, species_symbol);
    	GLPK.set_row_bnds(lp_problem, species_index, species_constraint_type, species_lower_bound, species_upper_bound);

    end

    # Setup the stoichiometric array -
    counter = 1;
    row_index_array = zeros(Int,number_of_species*number_of_fluxes);
    col_index_array = zeros(Int,number_of_species*number_of_fluxes);
    species_index_vector = collect(1:number_of_species);
    flux_index_vector = collect(1:number_of_fluxes);
    flat_stoichiometric_array = zeros(Float64,number_of_species*number_of_fluxes);
    for species_index in species_index_vector
    	for flux_index in flux_index_vector
    		row_index_array[counter] = species_index;
    		col_index_array[counter] = flux_index;
    		flat_stoichiometric_array[counter] = stoichiometric_matrix[species_index,flux_index];
    		counter+=1;
    	end
    end
    GLPK.load_matrix(lp_problem, number_of_species*number_of_fluxes, row_index_array, col_index_array, flat_stoichiometric_array);

    # Call the solver -
    exit_flag = GLPK.simplex(lp_problem, solver_parameters);

    # Get the objective function value -
    objective_value = GLPK.get_obj_val(lp_problem);

    # Get the calculated flux values from GLPK -
    calculated_flux_array = zeros(Float64,number_of_fluxes);
    for flux_index in flux_index_vector
    	calculated_flux_array[flux_index] = GLPK.get_col_prim(lp_problem, flux_index);
    end

    # Get the dual values -
    dual_value_array = zeros(Float64,number_of_fluxes);
    for flux_index in flux_index_vector
    	dual_value_array[flux_index] = GLPK.get_col_dual(lp_problem, flux_index);
    end

    # is this solution optimal?
    status_flag = GLPK.get_status(lp_problem)

    # Calculate the uptake array -
    uptake_array = stoichiometric_matrix*calculated_flux_array;

    # Formulate the return tuple -
    return (objective_value, calculated_flux_array, dual_value_array, uptake_array, exit_flag, status_flag);
end


text:
"The empty domain with elements of type `T`."
struct EmptySpace{T} <: Domain{T}
end

const AnyEmptySpace = EmptySpace{Any}

EmptySpace() = EmptySpace{Float64}()
EmptySpace(::Type{T}) where {T} = EmptySpace{T}()

emptyspace(d::Domain{T}) where {T} = EmptySpace{T}()

indomain(x::T, d::EmptySpace{T}) where {T} = false

approx_indomain(x, d::EmptySpace, tolerance) = in(x, d)

isempty(d::EmptySpace) = true

# Arithmetic operations

union(d1::EmptySpace, d2::EmptySpace) = d1
union(d1::Domain, d2::EmptySpace) = d1
union(d1::EmptySpace, d2::Domain) = d2

intersect(d1::EmptySpace, d2::EmptySpace) = d1
intersect(d1::Domain, d2::EmptySpace) = d2
intersect(d1::EmptySpace, d2::Domain) = d1

setdiff(d1::EmptySpace, d2::EmptySpace) = d1
setdiff(d1::EmptySpace, d2::Domain) = d1
setdiff(d1::Domain, d2::EmptySpace) = d1

# TODO: verify these - should we restrict x?
(+)(d::EmptySpace, x::Number) = d
(*)(d::EmptySpace, x::Number) = d

==(::EmptySpace, ::EmptySpace) = true

show(io::IO, d::EmptySpace) = print(io, "{} (empty domain)")

"The full space of elements of type `T`."
struct FullSpace{T} <: Domain{T} end

const AnyFullSpace = FullSpace{Any}

FullSpace() = FullSpace{Float64}()
FullSpace(d) = FullSpace{eltype(d)}()

euclideanspace(n::Val{N}) where {N} = euclideanspace(n, Float64)
euclideanspace(::Val{N}, ::Type{T}) where {N,T} = FullSpace{SVector{N,T}}()

indomain(x::T, d::FullSpace{T}) where {T} = true
indomain(x::S, d::FullSpace{T}) where {T,S} = promotes_to(S,T) == Val{true}

approx_indomain(x, d::FullSpace, tolerance) = in(x, d)

# We choose the origin as a point in the full space
point_in_domain(d::FullSpace) = zero(eltype(d))

isempty(::FullSpace) = false # constains zero

# Arithmetic operations

union(d1::FullSpace, d2::FullSpace) = d1
union(d1::Domain, d2::FullSpace) = d2
union(d1::FullSpace, d2::Domain) = d1

intersect(d1::FullSpace, d2::FullSpace) = d1
intersect(d1::Domain, d2::FullSpace) = d1
intersect(d1::FullSpace, d2::Domain) = d2

(+)(d::FullSpace{T}, x::T) where {T} = d

(*)(d::FullSpace, x::Number) = d

show(io::IO, d::FullSpace) = print(io, "{x} (full space)")

convert(::Type{Domain}, ::Type{T}) where T = FullSpace{T}()
convert(::Type{Domain{S}}, ::Type{T}) where {T,S} = convert(Domain{S}, convert(Domain, T))

infimum(d::FullSpace{T}) where {T} = typemin(T)
supremum(d::FullSpace{T}) where {T} = typemax(T)


text:
type Risk
    underlying::FinancialAsset
    hedge::FinancialAsset
    nominalrisk::Float64
end


text:
T=3
tsname="joKron"
@testset "$tsname" begin
for t=1:T # start test loop

    s=rand(3:9,3,2)
    a=rand(ComplexF64,s[1,:]...)
    A=joMatrix(a)
    b=rand(ComplexF64,s[2,:]...)
    B=joMatrix(b;DDT=ComplexF32,RDT=ComplexF64)
    c=rand(ComplexF64,s[3,:]...)
    C=joMatrix(c;DDT=ComplexF64,RDT=ComplexF32)
    k=kron(a,b,c)
    K=joKron(A,B,C)
    Km=K.m
    Kn=K.n

    exp_exp=(1.)/(3.)
    verbose && println("$tsname ($Km,$Kn)")
    @testset "$Km x $Kn" begin
        sn=prod(s[:,2])
        sm=prod(s[:,1])
        vn=rand(ComplexF64,sn);
        vm=rand(ComplexF64,sm);
        mvn=rand(ComplexF64,sn,2);
        mvm=rand(ComplexF64,sm,2);
        @test norm(k*vn-K*vn)<eps(norm(k*vn))^exp_exp
        @test norm(k*mvn-K*mvn)<eps(norm(k*mvn))^exp_exp
        @test norm(transpose(k)*vm-transpose(K)*vm)<eps(norm(transpose(k)*vm))^exp_exp
        @test norm(transpose(k)*mvm-transpose(K)*mvm)<eps(norm(transpose(k)*mvm))^exp_exp
        @test norm(adjoint(k)*vm-adjoint(K)*vm)<eps(norm(adjoint(k)*vm))^exp_exp
        @test norm(adjoint(k)*mvm-adjoint(K)*mvm)<eps(norm(adjoint(k)*mvm))^exp_exp
        @test norm(conj(k)*vn-conj(K)*vn)<eps(norm(conj(k)*vn))^exp_exp
        @test norm(conj(k)*mvn-conj(K)*mvn)<eps(norm(conj(k)*mvn))^exp_exp
        @test norm((-k)*vn-(-K)*vn)<eps(norm((-k)*vn))^exp_exp
        @test norm((-k)*mvn-(-K)*mvn)<eps(norm((-k)*mvn))^exp_exp
    end

end # end test loop
end


text:
using JuMP, GLPK, Cbc

function graph_coloring(; fast = true, benchmark = false)
    m = Model(optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0))

    @variable(m, max_color, Int)
    @variable(m, c[1:49, 1:4], Bin)

    @objective(m, Min, max_color)

    for s in 1:49
        @constraint(m, sum(c[s, 1:4]) == 1)
    end

    if fast
        for s in 1:49
            @constraint(m, sum(c[s, i] * i for i in 1:4) <= max_color)
        end
    else
        for s in 1:49, i in 1:4
            @constraint(m, max_color >= i * c[s, i])
        end
    end

    for i in 1:4
        @constraint(m, c[1, i] + c[11, i] <= 1)
        @constraint(m, c[1, i] + c[8, i] <= 1)
        @constraint(m, c[11, i] + c[8, i] <= 1)
        @constraint(m, c[11, i] + c[22, i] <= 1)
        @constraint(m, c[11, i] + c[24, i] <= 1)
        @constraint(m, c[24, i] + c[22, i] <= 1)
        @constraint(m, c[24, i] + c[36, i] <= 1)
        @constraint(m, c[22, i] + c[8, i] <= 1)
        @constraint(m, c[22, i] + c[23, i] <= 1)
        @constraint(m, c[22, i] + c[36, i] <= 1)
        @constraint(m, c[8, i] + c[2, i] <= 1)
        @constraint(m, c[8, i] + c[6, i] <= 1)
        @constraint(m, c[8, i] + c[23, i] <= 1)
        @constraint(m, c[23, i] + c[6, i] <= 1)
        @constraint(m, c[23, i] + c[31, i] <= 1)
        @constraint(m, c[23, i] + c[41, i] <= 1)
        @constraint(m, c[23, i] + c[36, i] <= 1)
        @constraint(m, c[36, i] + c[31, i] <= 1)
        @constraint(m, c[36, i] + c[41, i] <= 1)
        @constraint(m, c[2, i] + c[4, i] <= 1)
        @constraint(m, c[2, i] + c[5, i] <= 1)
        @constraint(m, c[2, i] + c[6, i] <= 1)
        @constraint(m, c[6, i] + c[5, i] <= 1)
        @constraint(m, c[6, i] + c[15, i] <= 1)
        @constraint(m, c[6, i] + c[31, i] <= 1)
        @constraint(m, c[31, i] + c[15, i] <= 1)
        @constraint(m, c[31, i] + c[33, i] <= 1)
        @constraint(m, c[31, i] + c[37, i] <= 1)
        @constraint(m, c[31, i] + c[41, i] <= 1)
        @constraint(m, c[41, i] + c[37, i] <= 1)
        @constraint(m, c[41, i] + c[40, i] <= 1)
        @constraint(m, c[4, i] + c[10, i] <= 1)
        @constraint(m, c[4, i] + c[5, i] <= 1)
        @constraint(m, c[5, i] + c[10, i] <= 1)
        @constraint(m, c[5, i] + c[13, i] <= 1)
        @constraint(m, c[5, i] + c[15, i] <= 1)
        @constraint(m, c[15, i] + c[13, i] <= 1)
        @constraint(m, c[15, i] + c[35, i] <= 1)
        @constraint(m, c[15, i] + c[33, i] <= 1)
        @constraint(m, c[33, i] + c[35, i] <= 1)
        @constraint(m, c[33, i] + c[37, i] <= 1)
        @constraint(m, c[37, i] + c[46, i] <= 1)
        @constraint(m, c[37, i] + c[40, i] <= 1)
        @constraint(m, c[40, i] + c[46, i] <= 1)
        @constraint(m, c[40, i] + c[47, i] <= 1)
        @constraint(m, c[10, i] + c[7, i] <= 1)
        @constraint(m, c[10, i] + c[13, i] <= 1)
        @constraint(m, c[13, i] + c[7, i] <= 1)
        @constraint(m, c[13, i] + c[26, i] <= 1)
        @constraint(m, c[13, i] + c[35, i] <= 1)
        @constraint(m, c[35, i] + c[26, i] <= 1)
        @constraint(m, c[35, i] + c[32, i] <= 1)
        @constraint(m, c[35, i] + c[39, i] <= 1)
        @constraint(m, c[35, i] + c[46, i] <= 1)
        @constraint(m, c[46, i] + c[39, i] <= 1)
        @constraint(m, c[46, i] + c[43, i] <= 1)
        @constraint(m, c[46, i] + c[47, i] <= 1)
        @constraint(m, c[47, i] + c[43, i] <= 1)
        @constraint(m, c[7, i] + c[26, i] <= 1)
        @constraint(m, c[7, i] + c[49, i] <= 1)
        @constraint(m, c[26, i] + c[21, i] <= 1)
        @constraint(m, c[26, i] + c[32, i] <= 1)
        @constraint(m, c[32, i] + c[21, i] <= 1)
        @constraint(m, c[32, i] + c[25, i] <= 1)
        @constraint(m, c[32, i] + c[29, i] <= 1)
        @constraint(m, c[32, i] + c[34, i] <= 1)
        @constraint(m, c[32, i] + c[39, i] <= 1)
        @constraint(m, c[39, i] + c[34, i] <= 1)
        @constraint(m, c[39, i] + c[38, i] <= 1)
        @constraint(m, c[39, i] + c[44, i] <= 1)
        @constraint(m, c[39, i] + c[42, i] <= 1)
        @constraint(m, c[39, i] + c[43, i] <= 1)
        @constraint(m, c[43, i] + c[42, i] <= 1)
        @constraint(m, c[49, i] + c[21, i] <= 1)
        @constraint(m, c[49, i] + c[25, i] <= 1)
        @constraint(m, c[21, i] + c[25, i] <= 1)
        @constraint(m, c[42, i] + c[44, i] <= 1)
        @constraint(m, c[42, i] + c[48, i] <= 1)
        @constraint(m, c[25, i] + c[17, i] <= 1)
        @constraint(m, c[25, i] + c[29, i] <= 1)
        @constraint(m, c[3, i] + c[12, i] <= 1)
        @constraint(m, c[12, i] + c[9, i] <= 1)
        @constraint(m, c[12, i] + c[14, i] <= 1)
        @constraint(m, c[9, i] + c[14, i] <= 1)
        @constraint(m, c[9, i] + c[16, i] <= 1)
        @constraint(m, c[14, i] + c[16, i] <= 1)
        @constraint(m, c[14, i] + c[19, i] <= 1)
        @constraint(m, c[14, i] + c[18, i] <= 1)
        @constraint(m, c[18, i] + c[19, i] <= 1)
        @constraint(m, c[16, i] + c[18, i] <= 1)
        @constraint(m, c[16, i] + c[17, i] <= 1)
        @constraint(m, c[16, i] + c[20, i] <= 1)
        @constraint(m, c[17, i] + c[20, i] <= 1)
        @constraint(m, c[17, i] + c[28, i] <= 1)
        @constraint(m, c[17, i] + c[30, i] <= 1)
        @constraint(m, c[17, i] + c[29, i] <= 1)
        @constraint(m, c[20, i] + c[28, i] <= 1)
        @constraint(m, c[30, i] + c[27, i] <= 1)
        @constraint(m, c[30, i] + c[29, i] <= 1)
        @constraint(m, c[30, i] + c[34, i] <= 1)
        @constraint(m, c[27, i] + c[34, i] <= 1)
        @constraint(m, c[29, i] + c[34, i] <= 1)
        @constraint(m, c[34, i] + c[38, i] <= 1)
        @constraint(m, c[38, i] + c[45, i] <= 1)
        @constraint(m, c[45, i] + c[44, i] <= 1)
        @constraint(m, c[44, i] + c[48, i] <= 1)
    end

    optimize!(m)
    if !benchmark
        println(JuMP.objective_value(m))
        println(JuMP.termination_status(m))
        _, color_washington = findmax(JuMP.value.(c[1, :]))
        println("color_washington: ", color_washington)
    end
end


text:
@doc raw"""
    BrownianMotionTorus(n::Int)

Returns a hidden state model corresponding to a Brownian motion on an `n`-dimensional torus, with initial condition drawn uniformly at random.
"""
struct BrownianMotionTorus <: HiddenStateModel{Vector{Float64}, ContinuousTime}
    n::Int
end

#####################
### BASIC METHODS ###
#####################

initial_condition(::BrownianMotionTorus) = Uniform(0,2π)

state_dim(model::BrownianMotionTorus) = model.n

"""
    noise_dim(model)

Returns the dimension of the Brownian motion ``W_t`` in the diffusion model ``dX_t = f(X_t)dt + g(X_t)dW_t``.
"""
noise_dim(model::BrownianMotionTorus) = model.n

function initialize(model::BrownianMotionTorus)
    x = rand(initial_condition(model), state_dim(model))
    return x
end

function Base.show(io::IO, ::MIME"text/plain", model::BrownianMotionTorus)
    print(io, "Brownian motion on the ", state_dim(model), "-torus for the hidden state
    type of hidden state:                   ", state_dim(model),"-dimensional vector
    number of independent Brownian motions: ", noise_dim(model),"
    initial condition:                      uniform")
end

function Base.show(io::IO, model::BrownianMotionTorus)
    print(io, "Brownian motion on the ", state_dim(model), "-torus")
end

######################
### TIME EVOLUTION ###
######################

function (model::BrownianMotionTorus)(x::AbstractVector{T}, dt) where T
    return mod2pi.(x .+ sqrt(dt) .* randn(T, noise_dim(model)))
end

function (model::BrownianMotionTorus)(x::AbstractMatrix, dt)
    N   = size(x, 2)
    out = zeros(eltype(x), state_dim(model), N)
    sqr = sqrt(dt)
    @inbounds for b in 1:N
        for a in 1:state_dim(model)
            out[a,b] += mod2pi(x[a,b] + sqr * randn(eltype(x)))
        end
    end
    return out
end

function propagate!(x::AbstractVector, model::BrownianMotionTorus, dt)
    N   = size(x, 2)
    sqr = sqrt(dt)
    @inbounds for a in 1:state_dim(model)
        x[a] +=  sqr * randn(eltype(x))
        x[a]  =  mod2pi(x[a])
    end
    return x
end

function propagate!(x::AbstractMatrix, model::BrownianMotionTorus, dt)
    N   = size(x, 2)
    sqr = sqrt(dt)
    @inbounds for b in 1:N
        for a in 1:state_dim(model)
            x[a,b] += sqr * randn(eltype(x))
            x[a,b]  = mod2pi(x[a,b])
        end
    end
    return x
end

################################
### CONVENIENCE CONSTRUCTORS ###
################################

const BrownianMotionCircle() = BrownianMotionTorus(1)

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
"""
    hamiltonian(L, lH, lHt)

Construct a Hamiltonian from a `Lattice` and operators.
# Arguments
* `L`: Lattice.
* `lH`: local Hamiltonian.
* `lHt`: `Tuple` of the form (t,Op) from which tunnelling terms are built
according to Σ<i,j> (t * Op[i] ⊗ Op[j]' + t' * Op[j] ⊗ Op[i]')
"""
function hamiltonian(L::Lattice, lH::O1, lHt::Tuple{Number,O2}) where {B<:Basis,O1<:AbstractOperator{B,B},O2<:AbstractOperator{B,B}}
    gbasis = CompositeBasis([lH.basis_l for i in 1:nv(L)]...);
    H = begin
            if !isdense(lH) && !isdense(lHt[2])
                SparseOperator(gbasis);
            else
                DenseOperator(gbasis);
            end
        end

    for e in edges(L)
        H.data .+= lHt[1] * embed(gbasis,[e.src, e.dst],[lHt[2]; dagger(lHt[2])]).data;
    end
    H.data .= H.data + H.data';

    for v in vertices(L)
        H.data .+= embed(gbasis, v, lH).data;
    end

    return H;
end
hamiltonian(L::NdLattice{N}, lH::O1, lHt::Tuple{Number,O2}) where {N,B<:Basis,O1<:AbstractOperator{B,B},O2<:AbstractOperator{B,B}} = hamiltonian(L::NdLattice{N}, lH::O1, lHt[1], lHt[2])

#function hamiltonian(L::Lattice, lH::O1, trate::Number, lHt::O2) where {N,Q,B<:Basis,O1<:AbstractOperator{B,B},O2<:AbstractOperator{B,B}}
#    hamiltonian(L, lH, (trate, lHt))
#end

"""
    hamiltonian(L, lH, trate, lHt)

Construct a Hamiltonian from a `Lattice` and operators.
# Arguments
* `L`: `NdLattice{N}`.
* `lH`: local Hamiltonian.
* `trate`: tunnelling rate or N-`Tuple` containing a rate per lattice dimension.
* `lHt`: local operator from which tunnelling terms are built as lHt[i]' ⊗ lHt[j]'.
"""
function hamiltonian(L::NdLattice{N}, lH::O1, trate::Tuple{Vararg{Number,Q}}, lHt::O2) where {N,Q,B<:Basis,O1<:AbstractOperator{B,B},O2<:AbstractOperator{B,B}}
    @assert Q == N "Number of tunnelling rates must match number of dimensions"
    gbasis = CompositeBasis([lH.basis_l for i in 1:nv(L)]...);
    H = begin
            if !isdense(lH) && !isdense(lHt)
                SparseOperator(gbasis);
            else
                DenseOperator(gbasis);
            end
        end
    lids = LinearIndices(L.shape)
    function lattice_slice(d::Int,i::Int)
        symbids = [mod1(i,L.shape[d]);fill(:(:), N-1)]
        Core.eval(CornerSpaceRenorm,Expr(:ref, :($lids), circshift(symbids,d-1)...))[:]
    end
    # E[d]: set of edges along the d-th dimension
    E = [Vector{LightGraphs.SimpleGraphs.SimpleEdge{Int64}}(undef,0) for i in 1:N]
    for d in 1:N
        i_max = L.pbc && (L.shape[d] > 2) ? L.shape[d] : L.shape[d]-1
        E[d] = vcat([[Edge(e) for e in zip(lattice_slice(d,i),lattice_slice(d,i+1))] for i in 1:i_max]...)
        if L.pbc
            if L.shape[d] == 1
                loops = [Edge(v=>v) for v in vertices(L)]
                E[d] = [loops; loops]
            elseif L.shape[d] == 2
                E[d] = [E[d]; E[d]]
            end
        end
    end

    for d in 1:N
        for e in E[d]
            H.data .+= trate[d] * embed_safe(gbasis, e, [lHt; dagger(lHt)]).data;
        end
    end
    H.data .= H.data + H.data';

    for v in vertices(L)
        H.data .+= embed(gbasis, v, lH).data;
    end

    return H;
end

function hamiltonian(L::NdLattice{N}, lH::O1, trate, lHt) where {N,B<:Basis,O1<:AbstractOperator{B,B}}
    gbasis = CompositeBasis([lH.basis_l for i in 1:nv(L)]...);

    _trate = begin
        if typeof(trate) <: Number
            Tuple(fill(ComplexF64(trate), N))
        elseif typeof(trate)<:Tuple && eltype(trate)<:Number
            @assert length(trate) == N "Number of coupling rates must match number of dimensions"
            ComplexF64.(trate)
        else
            error("The coupling rate must either be a Number or a Tuple")
        end
    end

    _lHt = begin
        if typeof(lHt) <: AbstractOperator{B,B}
            Tuple([lHt] for d in 1:N)
        elseif typeof(lHt) <: Vector{O} where {O<:AbstractOperator{B,B}}
            Tuple(deepcopy(lHt) for d in 1:N)
        elseif typeof(lHt)<:Tuple
            @assert length(lHt) == N "Number of coupling operators must match number of dimensions"
            if eltype(lHt)<:AbstractOperator{B,B}
                Tuple([lHt[d]] for d in 1:N)
            elseif eltype(lHt)<:Vector{O} where {O<:AbstractOperator{B,B}}
                Tuple([h for h in lHt[d]] for d in 1:N)
            else
                error("The tuple of coupling operators must contain either local operators or a vector or local operators")
            end
        else
            error("The coupling operator must be a local operator, a tuple of local operators or a tuple of vectors of local operators")
        end
    end

    H = begin
            if !isdense(lH) && !isdense(first(first(_lHt)))
                SparseOperator(gbasis);
            else
                DenseOperator(gbasis);
            end
        end

    lids = LinearIndices(L.shape)
    function lattice_slice(d::Int,i::Int)
        symbids = [mod1(i,L.shape[d]);fill(:(:), N-1)]
        collect(Core.eval(CornerSpaceRenorm,Expr(:ref, :($lids), circshift(symbids,d-1)...)))[:]
    end
    # E[d]: set of edges along the d-th dimension
    E = [Vector{LightGraphs.SimpleGraphs.SimpleEdge{Int64}}(undef,0) for i in 1:N]
    for d in 1:N
        i_max = L.pbc && (L.shape[d] > 2) ? L.shape[d] : L.shape[d]-1
        E[d] = vcat([[Edge(e) for e in zip(lattice_slice(d,i),lattice_slice(d,i+1))] for i in 1:i_max]...)
        if L.pbc
            if L.shape[d] == 1
                loops = [Edge(v=>v) for v in vertices(L)]
                E[d] = [loops; loops]
            elseif L.shape[d] == 2
                E[d] = [E[d]; E[d]]
            end
        end
    end

    for d in 1:N,  e in E[d], h in _lHt[d]
        H.data .+= _trate[d] * embed_safe(gbasis, e, [h, dagger(h)]).data;
    end
    H.data .= H.data + H.data';

    for v in vertices(L)
        H.data .+= embed(gbasis, v, lH).data;
    end

    return H
end

function embed_safe(gb::CompositeBasis, e::AbstractEdge, ops::AbstractVector{O}) where O<:AbstractOperator
    if e.src != e.dst
        embed(gb,[e.src, e.dst],ops)
    else
        embed(gb,e.src,prod(ops))
    end
end

"""
    dissipators(L, J)

Construct the jump operators of a system from a `Lattice` and local jump operators.
# Arguments
* `L`: Lattice.
* `J`: Array of on-site local jump operators.
"""
function dissipators(L::Lattice, J::Vector{O}) where {B<:Basis,O<:AbstractOperator{B,B}}
    gbasis = CompositeBasis([first(J).basis_l for i in 1:nv(L)]...);
    return vcat([[embed(gbasis, v, J[i]) for i in 1:length(J)] for v in vertices(L)]...)
end

isdense(op::DataOperator) = typeof(op.data) <: DenseMatrix
isdense(op::AbstractOperator) = false


text:
import DFT, FFTW
import DFT: fft_kernel_sizes
using Test

# kernel test
@testset "FFT kernel codelet" begin
    for i in fft_kernel_sizes
        data = rand(i)
        @test DFT.fft(data) ≈ FFTW.fft(data)
    end
end


text:
include("../demo_utilities.jl")
include("../../src/cg_modified.jl")
@info "Demo CG Modified Script running at $(@__DIR__)"
println("Chaging pwd to current directory of the script")
cd(@__DIR__)
@info "New pwd: $(pwd())"

function GetObjectiveVals(A, b, xs)
    return xs .|> x -> dot(x, A, x) - 2*dot(b, x)
end

"""
    Bad eigen value distribution, they are distributed from 0.001 to 1 like a
    geometric series.

    * A lot of small eigenvalues are clustered close to zero, few larger ones are
    far away from the other Eigenvalues.

    Parameters:
        rho: Number:
            An number between 0 and 1 that paramaterize the distribution of the
            eigenvalues for the PSD matrix.
        N=20:
            The size of the PSD matrix.

    returns:
        a diagonal matrix.
"""
function GetNastyPSDMatrix(rho::Number, N=20, inverted=false)
    @assert rho <= 1 && rho >= 0
    EigenValues = zeros(N)
    EigenMin, EigenMax = 0.001, 1    # Min Max Eigenvalues.
    EigenValues[1] = EigenMin
    for IdexI in 2:N
        EigenValues[IdexI] = EigenMin +
            ((IdexI - 1)/(N - 1))*(EigenMax - EigenMin)*rho^(N - IdexI)  # formulas
    end
    if inverted
        return Diagonal(1 + EigenMin - EigenValues)
    end
    return Diagonal(EigenValues)
end

N = 512
A = GetNastyPSDMatrix(0.9, N)
A = A^4
b = rand(N)
cg1 = ConjGradModified(A, b)
cg2 = ConjGradModified(A, b)
cg3 = ConjGradModified(A, b)
cg1Soln = Vector{typeof(b)}()
cg2Soln = Vector{typeof(b)}()
cg3Soln = Vector{typeof(b)}()
push!(cg1Soln, cg1.x)
push!(cg2Soln, cg2.x)
push!(cg3Soln, cg3.x)
TurnOffReorthgonalize(cg2)
ChangeStorageLimit(cg3, 32)
for _ in 1:2*N
    cg1(); cg2(); cg3()
    push!(cg1Soln, cg1.x); push!(cg2Soln, cg2.x); push!(cg3Soln, cg3.x)
end

objectiveVals1 = GetObjectiveVals(A, b, cg1Soln)
objectiveVals2 = GetObjectiveVals(A, b, cg2Soln)
objectiveVals3 = GetObjectiveVals(A, b, cg3Soln)
fig = plot(objectiveVals1, label="With Re-orthonalization", title="Convergence of the Objective")
plot!(fig, objectiveVals2, label="without Re-orthogonalization")
plot!(fig, objectiveVals3, label="partial Re-orthogonalization")
savefig(fig, "objectivevals_convergence.png")

cd(".")


text:
using LTC
using Plots
gr()
using BenchmarkTools
using DiffEqSensitivity
using OrdinaryDiffEq
using DiffEqFlux
using GalacticOptim
using Juno
using Cthulhu
using Profile
#using PProf

function generate_data()
    in_features = 2
    out_features = 1
    N = 48
    data_x = [sin.(range(0,stop=3π,length=N)), cos.(range(0,stop=3π,length=N))]
    data_x = [reshape([Float32(data_x[1][i]),Float32(data_x[2][i])],2,1) for i in 1:N]# |> f32
    data_y = [reshape([Float32(y)],1) for y in sin.(range(0,stop=6π,length=N))]# |> f32

    #data_x = [repeat(x,1,20) for x in data_x]
    #data_y = [repeat(x,1,20) for x in data_y]

    data_x, data_y
end

function data(iter; data_x=nothing, data_y=nothing, short=false, noisy=false)
    #noisy_data = Vector{Tuple{Vector{Matrix{Float64}}, Vector{Vector{Float64}}}}([])
    if data_y === nothing
      data_x, data_y = generate_data()
    end
    noisy_data = Vector{Tuple{Vector{Matrix{eltype(data_x[1])}}, Vector{Vector{eltype(data_y[1])}}}}([])
    for i in 1:iter
        x = data_x
        y = data_y
        if short isa Array
          x = x[short[1]:short[2]]
          y = y[short[1]:short[2]]
        end
        push!(noisy_data, (x , noisy ? add_gauss.(y,0.02) : y))
    end
    noisy_data
end

function loss(x,y,m::LTCNet{<:Mapper,<:Mapper,<:LTC.LTCCell,<:AbstractMatrix})
  GalacticOptim.Flux.reset!(m)
  #m = re(θ)
  #ŷ = m.(x)
  #ŷ = map(xi -> m(xi)[end-m.cell.wiring.n_motor+1:end, :], x)
  ŷ = [m(xi)[end-m.cell.wiring.n_motor+1:end, :] for xi in x]
  #ŷ = m.(x)
  #length(ŷ) < length(y) && return Inf
  #sum(Flux.Losses.mse.(ŷ, y; agg=mean))
  sum(sum([(ŷ[i][end,:] .- y[i]) .^ 2 for i in 1:length(y)]))/length(y)#, ŷ
  #sum([sum((ŷ[i] .- y[i]) .^ 2) for i in 1:length(y)])/length(y)
end
function cb(x,y,l,m)
  println(l)
  # pred = m.(x)
  # # isnan(l) && return false
  # fig = plot([ŷ[size(ŷ,1),1] for ŷ in pred])
  # plot!(fig, [yi[size(yi,1),1] for yi in y])
  # display(fig)
  return false
end

function traintest(n, solver=VCABM(), sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)))

  function lg(p,x,y,model)
    # reset_state!(model,p)
    reset!(model)
    # d = train_data[1]
    # x,y = d[1], d[2]
    m = model
    ŷ = [m(xi,p)[end-m.cell.wiring.n_motor+1:end, :] for xi in x]
    #ŷ = m.(x,[p])
    #losses = [Flux.Losses.mse(ŷ[i][end,:], y[i]) for i in 1:length(y)]
    #sum(losses)/length(losses), ŷ
    sum(sum([(ŷ[i][end,:] .- y[i]) .^ 2 for i in 1:length(y)]))/length(y), ŷ
  end
  cbg = function (p,l,pred;doplot=false) #callback function to observe training
    display(l)
    # plot current prediction against data
    if doplot
      fig = plot([ŷ[end,1] for ŷ in pred])
      plot!(fig, [yi[end,1] for yi in y])
      display(fig)
    end
    return false
  end

  x,y = generate_data()
  model = LTC.LTCNet(Wiring(2,1), solver, sensealg)
  #pp, re = Flux.destructure(model)
  # pp = DiffEqFlux.initial_params(model)
  # lower,upper = get_bounds(model)
  lower,upper = [],[]
  #θ = Flux.params(model)
  #θ = Flux.params(pp)

  pp = DiffEqFlux.initial_params(model)
  @show length(pp)
  @show length(pp)

  #@show sum(length.(θ))
  @show length(lower)

  train_data = data(n)

  opt = GalacticOptim.Flux.Optimiser(ClipValue(0.5), ADAM(0.05))

  # Juno.@profiler gs = Flux.Zygote.gradient(θ) do
  #   loss(first(train_data)...,model)
  # end
  # @time  gs = Flux.Zygote.gradient(θ) do
  #   loss(first(train_data)...,model)
  # end
  # @time  gs = Flux.Zygote.gradient(θ) do
  #   loss(first(train_data)...,model)
  # end
  #
  #
  # Juno.@profiler train_loss, back = Flux.Zygote.pullback(() -> loss(x,y,model), θ)
  # Juno.@profiler gs = back(one(train_loss))
  # train_loss, back = Flux.Zygote.pullback(() -> loss(x,y,model), θ)
  # gs = back(one(train_loss))
  # @time train_loss, back = Flux.Zygote.pullback(() -> loss(x,y,model), θ)
  # @time gs = back(one(train_loss))

  # use GalacticOptim.jl to solve the problem
  #adtype = GalacticOptim.AutoZygote()
  #
  #optf = GalacticOptim.OptimizationFunction((x, p) -> lg(x,model), adtype)
  #optfunc = GalacticOptim.instantiate_function(optf, model.cell.p, adtype, nothing)
  #optprob = GalacticOptim.OptimizationProblem(optfunc, model.cell.p, lb=lower, ub=upper)
  #
  #result_neuralode = GalacticOptim.solve(optprob,
  #                                     ParticleSwarm(;lower,upper,n_particles=6),
  #                                     cb = cbg,
  #                                     maxiters = 300)

  ##
  #optfun = OptimizationFunction((x,p,dx,dy)->lg(x,dx,dy,model), GalacticOptim.AutoZygote())
  #optprob = OptimizationProblem(optfun, θ[1], lb=lower, ub=upper)
  ##using IterTools: ncycle
  #res1 = GalacticOptim.solve(optprob, opt, train_data, cb = cbg, maxiters = n)
  #return res1

  #Flux.train!((x,y) -> loss(x,y,model), θ, train_data, opt; cb)

  Profile.clear()
  Profile.clear_malloc_data()

  optfun = OptimizationFunction((θ, p, x, y) -> lg(θ,x,y,model), GalacticOptim.AutoZygote())
  optprob = OptimizationProblem(optfun, pp)
  #using IterTools: ncycle
  #Juno.@profiler GalacticOptim.solve(optprob, opt, train_data, cb = cbg, maxiters = n) C = true
  GalacticOptim.solve(optprob, opt, train_data, cb = cbg, maxiters = 1000)

  # sciml_train(p->lg(p,model), pp, opt, cb = cbg, maxiters=100)

  #Juno.@profiler my_custom_train!(model, (x,y) -> loss(x,y,model), θ, train_data, opt; cb, lower, upper) C = true
end

function cthulu_test()

  m = LTC.LTCNet(Wiring(2,1), VCABM(), InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)))
  x = [rand(Float32,2,1) for _ in 1:10]
  p = initial_params(m)
  m.(x,[p])
  @time m.(x,[p])
  # @profile m.(x)
  #d = Dense(2,2)
  # r = RNN(2,2)
  # @descend_code_warntype r(x)
  #Flux.reset!(m)
  @descend_code_warntype m(x[1],p)
end
# cthulu_test()
#Profile.print()

# model = NCP(Wiring(2,1), VCABM(), InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)))
# ppp = Flux.params(model)
# sum(length.(ppp))
# Flux.trainable(model)

#@time traintest(10, AutoTsit5(Rosenbrock23()))
@time traintest(40)
#@time traintest(10)
#00:53 - 01:02 = 9 min compilation time

#@time traintest(300, VCABM(), ForwardDiffSensitivity())
#@time traintest(300, VCABM(), ReverseDiffAdjoint())
# @time traintest(300, Euler(), ForwardDiffSensitivity())
#@time traintest(300, VCABM(), InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)))
#@time traintest(300, VCABM(), InterpolatingAdjoint(autojacvec=false))
#@time traintest(300, VCABM(), InterpolatingAdjoint(autojacvec=ZygoteVJP()))
#@time traintest(300, VCABM(), InterpolatingAdjoint(checkpointing=true))
#traintest(300, AutoTsit5(Rosenbrock23()), InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)))
#@btime traintest(3, CVODE_BDF(linear_solver=:GMRES), InterpolatingAdjoint(autojacvec=ZygoteVJP()))
#@time traintest(300, CVODE_BDF(linear_solver=:GMRES), InterpolatingAdjoint(autojacvec=ZygoteVJP()))
#@time traintest(300, CVODE_BDF(linear_solver=:GMRES), InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)))

#  601.084 ms (1153411 allocations: 98.26 MiB)
#  8.510164 seconds (35.02 M allocations: 3.450 GiB, 5.96% gc time)

#@time traintest(300, VCABM(), InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)))

#@time traintest(3, AutoTsit5(Rosenbrock23()), InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)))
#@time traintest(300, AutoTsit5(Rosenbrock23()), InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)))

#  2.118 s (6965793 allocations: 461.57 MiB)

#ltc = Flux.Chain(Dense(2,5),Flux.LSTM(5,5),Flux.Dense(5,1))

# Flux.reset!(ltc)
# opt = GalacticOptim.Flux.Optimiser(ClipValue(0.1), ADAM(0.001))
#
# my_custom_train!(ltc, (x,y) -> lossf(x,y)[1], Flux.params(ltc), data(3), opt; cb=()->cbf(data_x,data_y,ltc),lower,upper)
# my_custom_train!(ltc, (x,y) -> lossf(x,y)[1], Flux.params(ltc), data(1000), opt; cb=()->cbf(data_x,data_y,ltc),lower,upper)

# my_custom_train!(ltc, (x,y) -> lossf(x,y)[1], Flux.params(ltc), data(100), opt; data_range=[1,8], cb=()->cbf(data_x,data_y,ltc),lower,upper)
# my_custom_train!(ltc, (x,y) -> lossf(x,y)[1], Flux.params(ltc), data(100), opt; data_range=[15,20], cb=()->cbf(data_x,data_y,ltc),lower,upper)
# my_custom_train!(ltc, (x,y) -> lossf(x,y)[1], Flux.params(ltc), data(100), opt; data_range=[10,30], cb=()->cbf(data_x,data_y,ltc),lower,upper)
# my_custom_train!(ltc, (x,y) -> lossf(x,y)[1], Flux.params(ltc), data(100), opt; data_range=[1,35], cb=()->cbf(data_x,data_y,ltc),lower,upper)
# my_custom_train!(ltc, (x,y) -> lossf(x,y)[1], Flux.params(ltc), data(100, short=[1,20]), ADAM(0.001); cb=()->cbf(data_x,data_y,ltc),lower,upper)
# my_custom_train!(ltc, (x,y) -> lossf(x,y)[1], Flux.params(ltc), data(100, short=[1,30]), ADAM(0.001); cb=()->cbf(data_x,data_y,ltc),lower,upper)
# my_custom_train!(ltc, (x,y) -> lossf(x,y)[1], Flux.params(ltc), data(3000), ADAM(0.001); cb=()->cbf(data_x,data_y,ltc),lower,upper)
#
# Flux.train!((x,y)->lossf(x,y)[1],Flux.params(ltc),data(200),ADAM(0.02); cb = ()->cbf(data_x,data_y,ltc))


text:
using Logging, Test, HSM

struct HsmStateInfoState <: HSM.AbstractHsmState
    # Actual contents are not necessary for tests, so they're left out for
    # convenience.
end

@testset "HsmStateInfo()" begin
    # Test default constructor.
    info = HSM.HsmStateInfo(nothing)
    @test !isnothing(info)
    @test isnothing(info.parent_state)
    @test isnothing(info.active_substate)

    something = HsmStateInfoState()
    info = HSM.HsmStateInfo(something)
    @test !isnothing(info)
    @test info.parent_state == something
    @test isnothing(info.active_substate)

    # Test multiple-argument constructor does not exist.
    @test_throws MethodError HSM.HsmStateInfo(nothing, nothing)
end


text:
# Comment: all the function have been adapted from a different repository.
# A precise reference will be inserted once that repository will become public

struct ODE{P}
    poly_ring::MPolyRing
    x_vars::Array{P, 1}
    y_vars::Array{P, 1}
    u_vars::Array{P, 1}
    parameters::Array{P, 1}
    x_equations::Dict{P, <: Union{P, Generic.Frac{P}}}
    y_equations::Dict{P, <: Union{P, Generic.Frac{P}}}
    function ODE{P}(
            x_eqs::Dict{P, <: Union{P, Generic.Frac{P}}},
            y_eqs::Dict{P, <: Union{P, Generic.Frac{P}}},
            inputs::Array{P, 1}
        ) where {P <: MPolyElem{<: FieldElem}}
        # Initialize ODE
        # x_eqs is a dictionary x_i => f_i(x, u, params)
        # y_eqs is a dictionary y_i => g_i(x, u, params)

        num, den = unpack_fraction(collect(values(x_eqs))[1])
        poly_ring = parent(num)
        if !all(isascii.(string.(gens(poly_ring))))
            nonascii_chars = filter(g->!isascii(g), string.(gens(poly_ring)))
            st = join(nonascii_chars, ", ")
            @warn "Non-ascii characters are not supported by Singular: " * st
        end
        x_vars = collect(keys(x_eqs))
        y_vars = collect(keys(y_eqs))
        u_vars = inputs
        parameters = filter(v -> (!(v in x_vars) && !(v in u_vars) && !(v in y_vars)), gens(poly_ring))
        new{P}(poly_ring, x_vars, y_vars, u_vars, parameters, x_eqs, y_eqs)
    end
end

# ------------------------------------------------------------------------------
"""
    func set_parameter_values(ode::ODE{P}, param_values::OrderedDict{P,T}) where {T <: FieldElem,P <: MPolyElem{T}}

Substitute parameters with numerical values.

## Input:
- ode, an ODE as above
- param_values, values for (some of) the parameters as dictionary parameter => value

## Output:
- new ode with the parameters in param_values plugged with the given numbers
"""
function set_parameter_values(ode::ODE{P}, param_values::OrderedDict{P,T}) where {T <: FieldElem,P <: MPolyElem{T}}
    new_vars = map(var_to_str, [v for v in gens(ode.poly_ring) if !(v in keys(param_values))])
    small_ring, small_vars = Nemo.PolynomialRing(base_ring(ode.poly_ring), new_vars)
    eval_dict = OrderedDict(str_to_var(v, ode.poly_ring) => str_to_var(v, small_ring) for v in new_vars)
    merge!(eval_dict, OrderedDict(p => small_ring(val) for (p, val) in param_values))

    return ODE{P}(
        OrderedDict{P,Union{P,Generic.Frac{P}}}(eval_at_dict(v, eval_dict) => eval_at_dict(f, eval_dict) for (v, f) in ode.x_equations),
         OrderedDict{P,Union{P,Generic.Frac{P}}}(eval_at_dict(v, eval_dict) => eval_at_dict(f, eval_dict) for (v, f) in ode.y_equations),
        [eval_at_dict(u, eval_dict) for u in ode.u_vars]
    )
end

# ------------------------------------------------------------------------------

"""
    func print_for_SIAN(ode::ODE{P}, outputs::Array{P,1}) where P <: MPolyElem{<: FieldElem}

Prints the ODE in the format accepted by SIAN (https://github.com/pogudingleb/SIAN)
"""
function print_for_SIAN(ode::ODE{P}, outputs::Array{P,1}) where P <: MPolyElem{<: FieldElem}
    vars_str = OrderedDict(x => var_to_str(x) * "(t)" for x in vcat(ode.x_vars, ode.u_vars))
    merge!(vars_str, OrderedDict(p => var_to_str(p) for p in ode.parameters))
    R_print, vars_print = Nemo.PolynomialRing(base_ring(ode.poly_ring), [vars_str[v] for v in gens(ode.poly_ring)])
    result = ""

    function _lhs_to_str(lhs)
        num, den = unpack_fraction(lhs)
        result = string(evaluate(num, vars_print))
        if den != 1
            result = "($result) / ($(evaluate(den, vars_print)))"
        end
        return result
    end

    for (x, f) in ode.equations
        result = result * "diff(" * var_to_str(x) * "(t), t) = $(_lhs_to_str(f)), \n"
    end
    for (y_ind, g) in enumerate(outputs)
        result = result * "y_var_$y_ind(t) = $(_lhs_to_str(g)), \n"
    end
    return result
end

# ------------------------------------------------------------------------------
"""
    func macrohelper_extract_vars(equations::Array{Expr,1})

A helper-function for a macro used in extracting variables from equations.
"""
function macrohelper_extract_vars(equations::Array{Expr,1})
    funcs, x_vars, all_symb = Array{Any}(undef, 0), Array{Any}(undef, 0), Array{Any}(undef, 0)
    aux_symb = Set([:(+), :(-), :(=), :(*), :(^), :t, :(/), :(//)])
    for eq in equations
        MacroTools.postwalk(
            x -> begin
            if @capture(x, f_'(t))
                push!(x_vars, f)
                push!(all_symb, f)
            elseif @capture(x, f_(t))
                push!(funcs, f)
            elseif (x isa Symbol) && !(x in aux_symb)
                push!(all_symb, x)
            end
            return x
        end,
            eq
        )
    end
    io_vars = setdiff(funcs, x_vars)
    all_symb = vcat(x_vars, io_vars, setdiff(all_symb, funcs))
    return x_vars, io_vars, all_symb
end

# ------------------------------------------------------------------------------
"""
    func macrohelper_clean(ex::Expr)

A cleanup helper for the macro.
"""
function macrohelper_clean(ex::Expr)
    ex = MacroTools.postwalk(x -> @capture(x, f_'(t)) ? f : x, ex)
    ex = MacroTools.postwalk(x -> @capture(x, f_(t)) ? f : x, ex)
    ex = MacroTools.postwalk(x -> x == :(/) ? :(//) : x, ex)
    return ex
end

# ------------------------------------------------------------------------------
"""
    macro ODEmodel(ex::Expr...)

Macro for creating an ODE from a list of equations and injecting all variables into the global scope.

Example:

```
ode = @ODEmodel(
    x1'(t) = - a * x1(t),
    y1(t) = x1(t),
)
```
"""
macro ODEmodel(ex::Expr...)
    equations = [ex...]
    x_vars, io_vars, all_symb = macrohelper_extract_vars(equations)

    # creating the polynomial ring
    vars_list = :([$(all_symb...)])
    R = gensym()
    vars_aux = gensym()
    exp_ring = :(($R, $vars_aux) = SIAN.Nemo.PolynomialRing(SIAN.Nemo.QQ, map(string, $all_symb)))
    assignments = [:($(all_symb[i]) = $vars_aux[$i]) for i in 1:length(all_symb)]

    # preparing equations
    equations = map(macrohelper_clean, equations)
    x_dict = gensym()
    y_dict = gensym()
    y_vars = Array{Any}(undef, 0)
    x_dict_create_expr = :($x_dict = Dict{SIAN.Nemo.fmpq_mpoly,Union{SIAN.Nemo.fmpq_mpoly,SIAN.Nemo.Generic.Frac{SIAN.Nemo.fmpq_mpoly}}}())
    y_dict_create_expr = :($y_dict = Dict{SIAN.Nemo.fmpq_mpoly,Union{SIAN.Nemo.fmpq_mpoly,SIAN.Nemo.Generic.Frac{SIAN.Nemo.fmpq_mpoly}}}())
    eqs_expr = []
    for eq in equations
        if eq.head != :(=)
            throw("Problem with parsing at $eq")
        end
        lhs, rhs = eq.args[1:2]
        loc_all_symb = macrohelper_extract_vars([rhs])[3]
        to_insert = undef
        if lhs in x_vars
            to_insert = x_dict
        elseif lhs in io_vars
            to_insert = y_dict
            push!(y_vars, lhs)
        else
            throw("Unknown left-hand side $lhs")
        end
        if isempty(loc_all_symb)
            push!(eqs_expr, :($to_insert[$lhs] = $R($rhs)))
        else
            push!(eqs_expr, :($to_insert[$lhs] = ($rhs)))
        end
    end

    u_vars = setdiff(io_vars, y_vars)
    params = setdiff(all_symb, union(x_vars, y_vars, u_vars))
    print("Summary of the model:\n")
    print("State variables: [", join(map(string, x_vars), ", "), "]\n")
    print("Parameter: [", join(map(string, params), ", "), "]\n")
    print("Inputs: [", join(map(string, u_vars), ", "), "]\n")
    print("Outputs: [", join(map(string, y_vars), ", "), "]\n")

    # creating the ode object
    ode_expr = :(ODE{SIAN.Nemo.fmpq_mpoly}($x_dict, $y_dict, Array{SIAN.Nemo.fmpq_mpoly}([$(u_vars...)])))

    result = Expr(
        :block,
        exp_ring, assignments...,
        x_dict_create_expr, y_dict_create_expr, eqs_expr...,
        ode_expr
    )
    return esc(result)
end

# ------------------------------------------------------------------------------
"""
    func generate_replica(ode::ODE{P}, r::Int) where P <: MPolyElem

Generate a replica of the original input system as per <Theorem here>.
Returns `ode_r`, and r-fold replica of the original ode.
States, outputs, and inputs are replicated, parameters are not.
"""
function generate_replica(ode::ODE{P}, r::Int) where P <: MPolyElem
    """
    Returns ode_r, and r-fold replica of the original ode.
    States, outputs, and inputs are replicated, parameters are not
    """
    new_varnames = Array{String}(undef, 0)
    for v in vcat(ode.x_vars, ode.y_vars, ode.u_vars)
        append!(new_varnames, [var_to_str(v) * "_r$i" for i in 1:r])
    end
    append!(new_varnames, map(string, ode.parameters))
    new_ring, new_vars = Nemo.PolynomialRing(base_ring(ode.poly_ring), new_varnames)
    new_x_eqs = OrderedDict{P,Union{P,Generic.Frac{P}}}()
    new_y_eqs = OrderedDict{P,Union{P,Generic.Frac{P}}}()
    new_us = Array{P,1}()
    for i in 1:r
        eval = merge(
            OrderedDict(v => str_to_var(var_to_str(v) * "_r$i", new_ring) for v in vcat(ode.x_vars, ode.y_vars, ode.u_vars)),
            OrderedDict(p => switch_ring(p, new_ring) for p in ode.parameters)
        )
        eval_vec = [eval[v] for v in gens(ode.poly_ring)]
        new_x_eqs = merge(
            new_x_eqs,
            OrderedDict{P,Union{P,Generic.Frac{P}}}(evaluate(x, eval_vec) => evaluate(f, eval_vec) for (x, f) in ode.x_equations)
        )
        new_y_eqs = merge(
            new_y_eqs,
            OrderedDict{P,Union{P,Generic.Frac{P}}}(evaluate(x, eval_vec) => evaluate(f, eval_vec) for (x, f) in ode.y_equations)
        )
        append!(new_us, [str_to_var(var_to_str(u) * "_r$i", new_ring) for u in ode.u_vars])
    end
    return ODE{P}(new_x_eqs, new_y_eqs, new_us)
end

# ------------------------------------------------------------------------------
"""
    func _reduce_poly_mod_p(poly::MPolyElem{Nemo.fmpq}, p::Int)

Reduces a polynomial modulo p.
"""
function _reduce_poly_mod_p(poly::MPolyElem{Nemo.fmpq}, p::Int)
    """
    Reduces a polynomial over Q modulo p
    """
    den = denominator(poly)
    num = change_base_ring(Nemo.ZZ, den * poly)
    if Nemo.GF(p)(den) == 0
        throw(Base.ArgumentError("Prime $p divides the denominator of $poly"))
    end
    return change_base_ring(Nemo.GF(p), num) * (1 // Nemo.GF(p)(den))
end

# ------------------------------------------------------------------------------

max_stars_repo_name: alexeyovchinnikov/SIAN-Julia

max_issues_repo_name: alexeyovchinnikov/SIAN-Julia

max_forks_repo_name: alexeyovchinnikov/SIAN-Julia


text:
s = open("input.txt") do file
    readlines(file);
end

function RunProgramUntilRepeat(program)
    seen = Set()

    i = 1;
    acc = 0;
    while !in(i, seen) && i <= length(program)
        push!(seen, i);
        parts = split(program[i], " ");

        cmd = parts[1];
        num = parse(Int, parts[2]);

        if cmd == "nop"
            i += 1;
        elseif cmd == "acc"
            acc += num;
            i += 1;
        elseif cmd == "jmp"
            i += num
        end
    end

    if i == length(program) + 1
        return acc
    end

    return -999999;
end

function FixProgram(program)
    for i = 1:length(program)
        if occursin("nop", program[i])
            program[i] = replace(program[i], "nop"=>"jmp");
            val = RunProgramUntilRepeat(program);

            if val != -999999
                return val
            end
            program[i] = replace(program[i], "jmp"=>"nop");

        elseif occursin("jmp", program[i])
            program[i] = replace(program[i], "jmp"=>"nop");
            val = RunProgramUntilRepeat(program);

            if val != -999999
                return val
            end
            program[i] = replace(program[i], "nop"=>"jmp");
        end
    end
end

FixProgram(s)


text:
module TestDoctest

using Documenter
using Test
using Transducers

@testset "/docs" begin
    include("../docs/utils.jl")
    transducers_rm_examples()
    doctest(Transducers; manual=true)
end

@testset "/test/doctests" begin
    doctest(joinpath((@__DIR__), "doctests"), Module[])
end

end  # module