text:
struct DisjointSet
    p::Vector{Int}
    rank::Vector{Int}
    function DisjointSet(n::Int)
        p = [i for i in 1:n]
        rank = zeros(Int, n)
        new(p, rank)
    end
end

function find_set(i::Int, ds::DisjointSet)::Int
    while i != ds.p[i]
        i = ds.p[i]
    end
    return ds.p[i]
end

function union_set!(i::Int, j::Int, ds::DisjointSet)
    i = find_set(i, ds)
    j = find_set(j, ds)
    if ds.rank[i] > ds.rank[j]
        ds.p[j] = i
    else
        if ds.rank[i] == ds.rank[j]
            ds.rank[j] += 1
        end
        ds.p[i] = j
    end
end


text:
using Aligner

rawsent1 = "i played in the jazz combo all four years"
rawsent2 = "i played in the jazz calm bowl all four years"

sent1 = split(rawsent1, " ")
sent2 = split(rawsent2, " ")

mat = generic_distance_matrix(sent1, sent2, 1, 1, normalized_word_distance)
printmat(mat)

alignment = align(sent1, sent2, 1, 1, normalized_word_distance)
println(alignment.res1)
println(alignment.res2)
printmat(alignment.mat)
println(alignment.idx1)

ext: jl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "eccodes"
version = v"2.18.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://confluence.ecmwf.int/download/attachments/45757960/eccodes-2.18.0-Source.tar.gz", "d88943df0f246843a1a062796edbf709ef911de7269648eef864be259e9704e3"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd eccodes-2.18.0-Source
if [ ${target} = "x86_64-w64-mingw32" ] || [ ${target} = "i686-w64-mingw32" ] ; then
    chmod +x cmake/ecbuild_windows_replace_symlinks.sh
    atomic_patch -p1 /workspace/srcdir/patches/windows.patch
else
    atomic_patch -p1 /workspace/srcdir/patches/unix.patch
fi
cd ..
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DENABLE_NETCDF=OFF -DENABLE_PNG=ON -DENABLE_PYTHON=OFF -DENABLE_FORTRAN=OFF ../eccodes-2.18.0-Source/
make -j${nproc}
make install
install_license ../eccodes-2.18.0-Source/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="libpng_jll", uuid="b53b4c65-9356-5827-b1ea-8c7a1a84506f"))
    Dependency(PackageSpec(name="OpenJpeg_jll", uuid="643b3616-a352-519d-856d-80112ee9badc"))
]

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


text:
function Kprec( pw_grid::PWGrid, psi::Array{ComplexF64,2} )
    Nstates = size(psi)[2]
    G2 = pw_grid.G2
    Npoints = pw_grid.Npoints
    Kpsi = zeros( ComplexF64, size(psi) )
    for ist = 1:Nstates
        for ip = 1:Npoints
            Kpsi[ip,ist] = psi[ip,ist] / ( 1.0 + G2[ip] )
        end
    end
    return Kpsi
end

# no preconditioning
function Kprec(psi)
  return psi
end


text:
#
# Conjugate gradient
#
# This is an independent implementation of:
#   W. W. Hager and H. Zhang (2006) Algorithm 851: CG_DESCENT, a
#     conjugate gradient method with guaranteed descent. ACM
#     Transactions on Mathematical Software 32: 113–137.
#
# Code comments such as "HZ, stage X" or "HZ, eqs Y" are with
# reference to a particular point in this paper.
#
# Several aspects of the following have also been incorporated:
#   W. W. Hager and H. Zhang (2012) The limited memory conjugate
#     gradient method.
#
# This paper will be denoted HZ2012 below.
#
# There are some modifications and/or extensions from what's in the
# paper (these may or may not be extensions of the cg_descent code
# that can be downloaded from Hager's site; his code has undergone
# numerous revisions since publication of the paper):
#
# cgdescent: the termination condition employs a "unit-correct"
#   expression rather than a condition on gradient
#   components---whether this is a good or bad idea will require
#   additional experience, but preliminary evidence seems to suggest
#   that it makes "reasonable" choices over a wider range of problem
#   types.
#
# linesearch: the Wolfe conditions are checked only after alpha is
#   generated either by quadratic interpolation or secant
#   interpolation, not when alpha is generated by bisection or
#   expansion. This increases the likelihood that alpha will be a
#   good approximation of the minimum.
#
# linesearch: In step I2, we multiply by psi2 only if the convexity
#   test failed, not if the function-value test failed. This
#   prevents one from going uphill further when you already know
#   you're already higher than the point at alpha=0.
#
# both: checks for Inf/NaN function values
#
# both: support maximum value of alpha (equivalently, c). This
#   facilitates using these routines for constrained minimization
#   when you can calculate the distance along the path to the
#   disallowed region. (When you can't easily calculate that
#   distance, it can still be handled by returning Inf/NaN for
#   exterior points. It's just more efficient if you know the
#   maximum, because you don't have to test values that won't
#   work.) The maximum should be specified as the largest value for
#   which a finite value will be returned.  See, e.g., limits_box
#   below.  The default value for alphamax is Inf. See alphamaxfunc
#   for cgdescent and alphamax for linesearch_hz.

macro cgtrace()
    quote
        if tracing
            dt = Dict()
            if o.extended_trace
                dt["x"] = copy(x)
                dt["g(x)"] = copy(g)
                dt["Current step size"] = alpha
            end
            g_norm = vecnorm(g, Inf)
            update!(tr,
                    iteration,
                    f_x,
                    g_norm,
                    dt,
                    o.store_trace,
                    o.show_trace,
                    o.show_every,
                    o.callback)
        end
    end
end

immutable ConjugateGradient{T} <: Optimizer
    eta::Float64
    P::T
    precondprep!::Function
    linesearch!::Function
end

function ConjugateGradient(;
                           linesearch!::Function = hz_linesearch!,
                           eta::Real = 0.4,
                           P::Any = nothing,
                           precondprep! = (P, x) -> nothing)
    ConjugateGradient{typeof(P)}(Float64(eta),
                                 P, precondprep!,
                                 linesearch!)
end

function optimize{T}(df::DifferentiableFunction,
                     initial_x::Array{T},
                     mo::ConjugateGradient,
                     o::OptimizationOptions)
    # Print header if show_trace is set
    print_header(o)

    # Maintain current state in x and previous state in x_previous
    x, x_previous = copy(initial_x), copy(initial_x)

    # Count the total number of iterations
    iteration = 0

    # Track calls to function and gradient
    f_calls, g_calls = 0, 0

    # Count number of parameters
    n = length(x)

    # Maintain current gradient in g and previous gradient in g_previous
    g, g_previous = similar(x), similar(x)

    # Maintain the preconditioned gradient in pg
    pg = similar(x)

    # The current search direction
    s = similar(x)

    # Buffers for use in line search
    x_ls, g_ls = similar(x), similar(x)

    # Intermediate value in CG calculation
    y = similar(x)
    py = similar(x)

    # Store f(x) in f_x
    f_x = df.fg!(x, g)
    @assert typeof(f_x) == T
    f_x_previous = convert(T, NaN)
    f_calls, g_calls = f_calls + 1, g_calls + 1
    copy!(g_previous, g)

    # Keep track of step-sizes
    alpha = alphainit(one(T), x, g, f_x)

    # TODO: How should this flag be set?
    mayterminate = false

    # Maintain a cache for line search results
    lsr = LineSearchResults(T)

    # Trace the history of states visited
    tr = OptimizationTrace{typeof(mo)}()
    tracing = o.store_trace || o.show_trace || o.extended_trace || o.callback != nothing
    @cgtrace

    # Output messages
    if !isfinite(f_x)
        error("Must have finite starting value")
    end
    if !all(isfinite(g))
        @show g
        @show find(!isfinite(g))
        error("Gradient must have all finite values at starting point")
    end

    # Determine the intial search direction
    #    if we don't precondition, then this is an extra superfluous copy
    #    TODO: consider allowing a reference for pg instead of a copy
    mo.precondprep!(mo.P, x)
    A_ldiv_B!(pg, mo.P, g)
    scale!(copy!(s, pg), -1)

    # Assess multiple types of convergence
    x_converged, f_converged = false, false
    g_converged = vecnorm(g, Inf) < o.g_tol

    # Iterate until convergence
    converged = g_converged
    while !converged && iteration < o.iterations
        # Increment the number of steps we've had to perform
        iteration += 1

        # Reset the search direction if it becomes corrupted
        dphi0 = vecdot(g, s)
        if dphi0 >= 0
            @simd for i in 1:n
                @inbounds s[i] = -pg[i]
            end
            dphi0 = vecdot(g, s)
            if dphi0 >= 0
                break
            end
        end

        # Refresh the line search cache
        clear!(lsr)
        @assert typeof(f_x) == T
        @assert typeof(dphi0) == T
        push!(lsr, zero(T), f_x, dphi0)

        # Pick the initial step size (HZ #I1-I2)
        alpha, mayterminate, f_update, g_update =
          alphatry(alpha, df, x, s, x_ls, g_ls, lsr)
        f_calls, g_calls = f_calls + f_update, g_calls + g_update

        # Determine the distance of movement along the search line
        alpha, f_update, g_update =
          mo.linesearch!(df, x, s, x_ls, g_ls, lsr, alpha, mayterminate)
        f_calls, g_calls = f_calls + f_update, g_calls + g_update

        # Maintain a record of previous position
        copy!(x_previous, x)

        # Update current position # x = x + alpha * s
        LinAlg.axpy!(alpha, s, x)

        # Maintain a record of the previous gradient
        copy!(g_previous, g)

        # Update the function value and gradient
        f_x_previous, f_x = f_x, df.fg!(x, g)
        f_calls, g_calls = f_calls + 1, g_calls + 1

        x_converged,
        f_converged,
        g_converged,
        converged = assess_convergence(x,
                                       x_previous,
                                       f_x,
                                       f_x_previous,
                                       g,
                                       o.x_tol,
                                       o.f_tol,
                                       o.g_tol)

        # Check sanity of function and gradient
        if !isfinite(f_x)
            error("Function must finite function values")
        end

        # Determine the next search direction using HZ's CG rule
        #  Calculate the beta factor (HZ2012)
        # -----------------
        # Comment on py: one could replace the computation of py with
        #    ydotpgprev = vecdot(y, pg)
        #    vecdot(y, py)  >>>  vecdot(y, pg) - ydotpgprev
        # but I am worried about round-off here, so instead we make an
        # extra copy, which is probably minimal overhead.
        # -----------------
        mo.precondprep!(mo.P, x)
        dPd = dot(s, mo.P, s)
        etak::T = mo.eta * vecdot(s, g_previous) / dPd
        @simd for i in 1:n
            @inbounds y[i] = g[i] - g_previous[i]
        end
        ydots = vecdot(y, s)
        copy!(py, pg)        # below, store pg - pg_previous in py
        A_ldiv_B!(pg, mo.P, g)
        @simd for i in 1:n     # py = pg - py
           @inbounds py[i] = pg[i] - py[i]
        end
        betak = (vecdot(y, pg) - vecdot(y, py) * vecdot(g, s) / ydots) / ydots
        beta = max(betak, etak)
        @simd for i in 1:n
            @inbounds s[i] = beta * s[i] - pg[i]
        end

        @cgtrace
    end

    return MultivariateOptimizationResults("Conjugate Gradient",
                                           initial_x,
                                           x,
                                           Float64(f_x),
                                           iteration,
                                           iteration == o.iterations,
                                           x_converged,
                                           o.x_tol,
                                           f_converged,
                                           o.f_tol,
                                           g_converged,
                                           o.g_tol,
                                           tr,
                                           f_calls,
                                           g_calls)
end


text:
using NCDatasets
using DataStructures
using ArgParse, JSON
using DataStructures

println("""
This program is run after a CESM case has been setup before execute.
It reads env_build.xml and env_run.xml to fill in paths and files.
""")

function runOneCmd(cmd)
    println(">> ", string(cmd))
    run(cmd)
end

function pleaseRun(cmd)
    if isa(cmd, Array)
        for i = 1:length(cmd)
            runOneCmd(cmd[i])
        end
    else
        runOneCmd(cmd)
    end
end

function getCESMConfig(path, ids)

    println("CESM config path: $(path)")

    wdir = pwd()

    cd(path)

    d = OrderedDict()
    for id in ids
        d[id] = readchomp(`./xmlquery $(id) -silent -valonly`)
    end

    cd(wdir)

    return d
end

#=
function getXML(filename)
    println("Reading XML file $(filename)")
    xml_doc = parse_file(filename)
    entries = get_elements_by_tagname(root(xml_doc), "entry")

    d = OrderedDict()
    for entry in entries
        d[attribute(entry, "id")] = attribute(entry, "value")
    end

    return d
end
=#

function parse_commandline()

    s = ArgParseSettings()

    @add_arg_table s begin

        "--domain-file"
            help = "Horizontal domain file containing z_w_top, z_w_bot in meters with z=0 is the surface."
            arg_type = String
            required = true

        "--zdomain-file"
            help = "Vertical domain file containing z_w_top, z_w_bot in meters with z=0 is the surface."
            arg_type = String
            required = true

        "--topo-file"
            help = "File containing Nz_bot."
            arg_type = String
            required = true

        "--forcing-file-HMXL"
            help = "The forcing file"
            arg_type = String
            default = ""

        "--forcing-file-TEMP"
            help = "The forcing file"
            arg_type = String
            default = ""

        "--forcing-file-SALT"
            help = "The forcing file"
            arg_type = String
            default = ""

        "--forcing-file-QFLXT"
            help = "The forcing file"
            arg_type = String
            default = ""

        "--forcing-file-QFLXS"
            help = "The forcing file"
            arg_type = String
            default = ""

        "--forcing-time"
            help = "Beg, end, align"
            arg_type = String
            nargs = 3
            default = ["", "", ""]

        "--init-file"
            help = "The init ocean file."
            arg_type = String
            required = true

        "--ocn-model"
            help = "The name of ocean model. Can be: `SOM`, `MLM`, `EOM` and `EMOM`"
            arg_type = String
            required = true

        "--output-file"
            help = "The output file."
            arg_type = String
            default = "config.toml"

        "--finding-QFLX-timescale"
            help = "If this is set to a real number, finding QFLX mode is on. And timescale is in days."
            arg_type = Float64
            default = NaN

        "--casename"
            help = "Casename. Default: UNSET"
            arg_type = String
            default = "UNSET"

        "--caseroot"
            help = "Caseroot. Default: UNSET"
            arg_type = String
            default = "UNSET"

        "--caserun"
            help = "Caserun. Default: UNSET"
            arg_type = String
            default = "UNSET"

        "--archive-root"
            help = "Archive root. Default: UNSET"
            arg_type = String
            default = "UNSET"

        "--calendar"
            help = "Calendar type. Two possible choices: NO_LEAP, GREGORIAN. Default: NO_LEAP"
            arg_type = String
            default = "NO_LEAP"

    end

    return parse_args(s)
end

parsed = DataStructures.OrderedDict(parse_commandline())
JSON.print(parsed,4)

if parsed["ocn-model"] == "SOM"
    advection_scheme = "static"
elseif parsed["ocn-model"] == "MLM"
    advection_scheme = "static"
elseif parsed["ocn-model"] in ["EOM", "EMOM"]
    advection_scheme = "ekman_CO2012"
else
    throw(ErrorException("Error: Unknown ocean model `$(parsed["ocn-model"])`."))
end

if parsed["finding-QFLX-timescale"] > 0.0
    println("QFLX finding mode on. Timescale = $(parsed["finding-QFLX-timescale"]) days.")
    Qflx = "off"
    τwk_SALT = 86400.0 * parsed["finding-QFLX-timescale"]
    τwk_TEMP = 86400.0 * parsed["finding-QFLX-timescale"]
elseif isnan(parsed["finding-QFLX-timescale"])
    println("QFLX finding mode off. Weak restoring timescale = 100 years for TEMP and Salt.")
    Qflx = "on"
    τwk_SALT = 86400.0 * 365 * 100
    τwk_TEMP = 86400.0 * 365 * 100
else
    throw(ErrorException("Unknown scenario: got `finding-QFLX-timescale` = $(parsed["finding-QFLX-timescale"])"))
end

calendar = Dict(
    "NO_LEAP"   => "DateTimeNoLeap",
    "GREGORIAN" => "DateTimeProlepticGregorian",
)[parsed["calendar"]]

var_file_map = Dict()

for var in ["HMXL", "TEMP", "SALT", "QFLXT", "QFLXS"]
    file = parsed["forcing-file-$(var)"]
    if file != ""
        var_file_map[var] = file
    end
end

config = DataStructures.OrderedDict{Any, Any}(

    "DRIVER" => DataStructures.OrderedDict(
        "casename"           => parsed["casename"],
        "caseroot"           => parsed["caseroot"],
        "caserun"            => parsed["caserun"],
        "archive_root"       => parsed["archive-root"],
    ),

    "MODEL_MISC" => DataStructures.OrderedDict(
        "timetype"               => calendar,
        "init_file"              => parsed["init-file"],
        "rpointer_file"          => "rpointer.iom",
        "daily_record"           => [],
        "monthly_record"         => ["{ESSENTIAL}", "QFLXT", "QFLXS"],
        "enable_archive"         => true,
    ),

    "MODEL_CORE" => DataStructures.OrderedDict(

        "domain_file"                  => parsed["domain-file"],
        "topo_file"                    => parsed["topo-file"],
        "cdata_var_file_map"           => var_file_map,

        "cdata_beg_time"               => parsed["forcing-time"][1],
        "cdata_end_time"               => parsed["forcing-time"][2],
        "cdata_align_time"             => parsed["forcing-time"][3],

        "z_w"                          => nothing,

        "substeps"                     => 8,
        "MLD_scheme"                   => "datastream",
        "Qflx"                         => Qflx,
        "convective_adjustment"        => "on",
        "advection_scheme"             => advection_scheme,

        "weak_restoring"               => "on",
        "τwk_TEMP"                     => τwk_TEMP,
        "τwk_SALT"                     => τwk_SALT,

        "Ekman_layers"                 => 5,
        "Returnflow_layers"            => 28,

        "transform_vector_field"       => true,
    ),

)

Dataset(parsed["zdomain-file"], "r") do ds

    z_w_top = nomissing(ds["z_w_top"][:], NaN)
    z_w_bot = nomissing(ds["z_w_bot"][:], NaN)

    z_w      = zeros(Float64, length(z_w_top)+1)
    z_w[1:end-1]   = z_w_top
    z_w[end] = z_w_bot[end]

    global config["MODEL_CORE"]["z_w"]  = z_w
end

using TOML

output_config = parsed["output-file"]
println("Output file: $(output_config)")
open(output_config, "w") do io
    TOML.print(io, config; sorted=true)
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
@testset "Correlator" begin

    @testset "Early prompt late correlator" begin

        correlator = @inferred EarlyPromptLateCorrelator()

        @test @inferred(get_early(correlator)) == 0.0
        @test @inferred(get_prompt(correlator)) == 0.0
        @test @inferred(get_late(correlator)) == 0.0

        correlator = @inferred EarlyPromptLateCorrelator(NumAnts(1))

        @test @inferred(get_early(correlator)) == 0.0
        @test @inferred(get_prompt(correlator)) == 0.0
        @test @inferred(get_late(correlator)) == 0.0

        correlator = @inferred EarlyPromptLateCorrelator(NumAnts(2))

        @test @inferred(get_early(correlator)) == SVector(0.0 + 0.0im, 0.0 + 0.0im)
        @test @inferred(get_prompt(correlator)) == SVector(0.0 + 0.0im, 0.0 + 0.0im)
        @test @inferred(get_late(correlator)) == SVector(0.0 + 0.0im, 0.0 + 0.0im)

        correlator = @inferred EarlyPromptLateCorrelator(
            SVector(1.0 + 0.0im, 1.0 + 0.0im),
            SVector(1.0 + 0.0im, 1.0 + 0.0im),
            SVector(1.0 + 0.0im, 1.0 + 0.0im)
        )

        @test @inferred(zero(correlator)) == EarlyPromptLateCorrelator(
            SVector(0.0 + 0.0im, 0.0 + 0.0im),
            SVector(0.0 + 0.0im, 0.0 + 0.0im),
            SVector(0.0 + 0.0im, 0.0 + 0.0im)
        )

        correlator = @inferred EarlyPromptLateCorrelator(
            1.0 + 0.0im,
            1.0 + 0.0im,
            1.0 + 0.0im
        )

        @test @inferred(zero(correlator)) == EarlyPromptLateCorrelator(
            0.0 + 0.0im,
            0.0 + 0.0im,
            0.0 + 0.0im
        )

        correlator = @inferred EarlyPromptLateCorrelator(
            SVector(1.0 + 0.0im, 1.0 + 0.0im),
            SVector(1.0 + 0.0im, 1.0 + 0.0im),
            SVector(1.0 + 0.0im, 1.0 + 0.0im)
        )
        filtered_correlator = @inferred Tracking.filter(x -> x[1], correlator)
        @test filtered_correlator == EarlyPromptLateCorrelator(
            1.0 + 0.0im,
            1.0 + 0.0im,
            1.0 + 0.0im
        )

        correlator = @inferred EarlyPromptLateCorrelator(
            1.0 + 0.0im,
            1.0 + 0.0im,
            1.0 + 0.0im
        )
        filtered_correlator = @inferred Tracking.filter(x -> x, correlator)
        @test filtered_correlator == EarlyPromptLateCorrelator(
            1.0 + 0.0im,
            1.0 + 0.0im,
            1.0 + 0.0im
        )

        early_late_sample_shift = @inferred get_early_late_sample_shift(
            GPSL1,
            correlator,
            4e6Hz,
            0.5
        )
        @test early_late_sample_shift ≈ round(0.5 * 4e6 / 1023e3)
        @test typeof(early_late_sample_shift) <: Integer

        normalized_correlator = @inferred Tracking.normalize(correlator, 10)
        @test normalized_correlator == EarlyPromptLateCorrelator(
            0.1 + 0.0im,
            0.1 + 0.0im,
            0.1 + 0.0im
        )

        correlator = EarlyPromptLateCorrelator(
            SVector(1.0 + 0.0im, 1.0 + 0.0im),
            SVector(1.0 + 0.0im, 1.0 + 0.0im),
            SVector(1.0 + 0.0im, 1.0 + 0.0im)
        )
        normalized_correlator = @inferred Tracking.normalize(correlator, 10)
        @test normalized_correlator == EarlyPromptLateCorrelator(
            SVector(0.1 + 0.0im, 0.1 + 0.0im),
            SVector(0.1 + 0.0im, 0.1 + 0.0im),
            SVector(0.1 + 0.0im, 0.1 + 0.0im)
        )

        signal = 1.0 + 0.0im
        carrier = cis(0.0)
        early_late_sample_shift = 2
        # Late Prompt Early
        early_code = 1
        code_register = 1 << 4 + 0 << 2 + 1
        correlator = EarlyPromptLateCorrelator(1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im)

        next_correlator = @inferred Tracking.correlate_iteration(
            GPSL1,
            correlator,
            signal,
            carrier,
            code_register,
            early_late_sample_shift,
            early_code
        )
        @test get_early(next_correlator) ≈ 1.0 + signal * conj(carrier) * 1
        @test get_prompt(next_correlator) ≈ 1.0 + signal * conj(carrier) * -1
        @test get_late(next_correlator) ≈ 1.0 + signal * conj(carrier) * 1

        early_code = -1
        code_register = 0 << 4 + 1 << 2 + 0

        next_correlator = @inferred Tracking.correlate_iteration(
            GPSL1,
            correlator,
            signal,
            carrier,
            code_register,
            early_late_sample_shift,
            early_code
        )
        @test get_early(next_correlator) ≈ 1.0 + signal * conj(carrier) * -1
        @test get_prompt(next_correlator) ≈ 1.0 + signal * conj(carrier) * 1
        @test get_late(next_correlator) ≈ 1.0 + signal * conj(carrier) * -1

    end

end


text:
#-------------------------------------------------------------------------------
# File: move_insertion.jl
# Description: This file contains all functions relatives to the insertion move.
# Date: November 4, 2019
# Author: Jonathan Fontaine, Killian Fretaud, Rémi Garcia,
#         Boualem Lamraoui, Benoît Le Badezet, Benoit Loger
#-------------------------------------------------------------------------------

"""
    move_insertion!(solution::Solution, old_index::Int,
                    new_index::Int, instance::Instance)
Replaces car at `old_index` at `new_index` and updates the solution's matrices.
"""
function move_insertion!(solution::Solution, old_index::Int,
                         new_index::Int, instance::Instance)

    # #DEBUG and TODO: Comment / uncomment for debbuging.
    # #                Ensure the quality of the solution.
    # update_matrices!(solution, instance)
    # initialize_batches!(solution, instance)

    # #DEBUG and TODO: Comment / uncomment for debbuging.
    # #                Test the quality of the given solution.
    # s1 = deepcopy(solution)
    #
    # update_matrices!(solution, instance)
    # if !(solution.colors === nothing)
    #     initialize_batches!(solution, instance)
    # end
    #
    # error1 = false
    # if !(solution.colors === nothing)
    #     for i in 1:solution.length
    #         bool_print = false
    #         if s1.colors[i].start != solution.colors[i].start
    #             #println("\tBad start batch (",i,") - ", s1.colors[i].start," and ",solution.colors[i].start)
    #             error1 = true
    #             bool_print= true
    #         end
    #         if s1.colors[i].width != solution.colors[i].width
    #             #println("\tBad width batch (",i,") - ", s1.colors[i].width," and ",solution.colors[i].width)
    #             error1 = true
    #             bool_print= true
    #         end
    #         if bool_print
    #             #println("\t", instance.color_code[i] ,"-", instance.color_code[car_pos_a]," and ",instance.color_code[car_pos_b])
    #         end
    #     end
    # end
    #
    # error2 = false
    # for o in 1:instance.nb_HPRC+instance.nb_LPRC
    #     for i in 1:solution.length
    #         if s1.M1[o,i] != solution.M1[o,i]
    #             error2 = true
    #         end
    #         if s1.M2[o,i] != solution.M2[o,i]
    #             error2 = true
    #         end
    #         if s1.M3[o,i] != solution.M3[o,i]
    #             error2 = true
    #         end
    #     end
    # end
    # if error1 || error2
    #     print("\t\t\tIn move_insertion for given sol ", old_index, " at ", new_index)
    #     if error1
    #         print(" --  Error batch")
    #     end
    #     if error2
    #         print(" -- Error matrices")
    #     end
    #     println()
    # end

    #Update sequence (better complexity ?)
    car_inserted = solution.sequence[old_index]
    if old_index < new_index
        for car_moved_pos in old_index:(new_index-1)
            solution.sequence[car_moved_pos] = solution.sequence[car_moved_pos+1]
        end
        solution.sequence[new_index] = car_inserted
    end
    if old_index > new_index
        for car_moved_pos in old_index:-1:(new_index+1)
            solution.sequence[car_moved_pos] = solution.sequence[car_moved_pos-1]
        end
        solution.sequence[new_index] = car_inserted
    end

    # TODO: COMPLEXITY NON-OPTIMAL
    update_matrices!(solution, instance)
    if !(solution.colors === nothing)
        initialize_batches!(solution, instance)
    end

    # #DEBUG and TODO: Comment / uncomment for debbuging.
    # #                Test the quality of the move.
    # s1 = deepcopy(solution)
    #
    # update_matrices!(solution, instance)
    # if !(solution.colors === nothing)
    #     initialize_batches!(solution, instance)
    # end
    #
    # error1 = false
    # for i in 1:solution.length
    #     bool_print = false
    #     if s1.colors[i].start != solution.colors[i].start
    #         println("\tBad start batch (",i,") - ", s1.colors[i].start," and ",solution.colors[i].start)
    #         error1 = true
    #         bool_print= true
    #     end
    #     if s1.colors[i].width != solution.colors[i].width
    #         println("\tBad width batch (",i,") - ", s1.colors[i].width," and ",solution.colors[i].width)
    #         error1 = true
    #         bool_print= true
    #     end
    #     if bool_print
    #         println("\t", instance.color_code[i] ,"-", instance.color_code[old_index]," to ",instance.color_code[new_index])
    #     end
    # end
    #
    # error2 = false
    # for o in 1:instance.nb_HPRC+instance.nb_LPRC
    #     for i in 1:solution.length
    #         if s1.M1[o,i] != solution.M1[o,i]
    #             error2 = true
    #         end
    #         if s1.M2[o,i] != solution.M2[o,i]
    #             error2 = true
    #         end
    #         if s1.M3[o,i] != solution.M3[o,i]
    #             error2 = true
    #         end
    #     end
    # end
    # if error1 || error2
    #     println("\nIn move_insertion for ", old_index, " to ", new_index)
    #     if error1
    #         println("\nError batch")
    #     end
    #     if error2
    #         println("\nError matrices")
    #     end
    # end

    return solution
end

"""
    cost_move_insertion(solution::Solution, car_pos_a::Int,
                        instance::Instance, objectives::BitArray{1})
Returns the cost of the insertion of the car `car_pos_a` for all valid positions
with respect to objectives `objectives`. A negative cost means that
the move is interesting with respect to treated objectives.
"""
function cost_move_insertion(solution::Solution, position::Int,
                             instance::Instance, objectives::BitArray{1})
    @assert length(objectives) == 3

    # # DEBUG and TODO: Comment / uncomment for debbuging.
    # #                 Ensure to test every variation.
    # objectives = trues(3)

    cost_on_objective = zeros(Int, solution.length, 3)
    b0 = instance.nb_late_prec_day+1

    delta1_for_first = zeros(Int, solution.length)
    delta1_for_second = zeros(Int, solution.length)
    delta2_for_first = zeros(Int, solution.length)
    delta2_for_second = zeros(Int, solution.length)

        #---------------------------------------------------------- Initialization of M1

    sequence = copy(solution.sequence)
    deleteat!(sequence, position)
    M1 = copy(solution.M1)

    #violations_caused_on_X removing the car will cause a variation on number of violations
    if objectives[1]
        M1, violations_caused_on_first = update_lines_remove!(
            solution, instance, # instance
            M1, position, # calcul
            1, instance.nb_HPRC # lines modified
        )
    end

    if objectives[2]
        M1, violations_caused_on_second = update_lines_remove!(
            solution, instance, # instance
            M1, position, # calcul
            instance.nb_HPRC+1, instance.nb_HPRC+instance.nb_LPRC # lines modified
        )
    end

        #---------------------------------------------------------- Cost of insertion

    # Cost on HPRC
    if objectives[1]
        #What about insertion at b0 (first place available)
        delta2_for_first[b0] = compute_delta2_for_b0(
            solution, instance, # instance
            M1, sequence, position, # calcul
            1, instance.nb_HPRC # for options
        )
        delta1_for_first[b0] = compute_delta1(
            solution, instance, # instance
            M1, sequence, b0, position, # calcul
            1, instance.nb_HPRC # for options
        )

        # for all other valid positions
        for index in (b0+1):solution.length
            # Delta 1 is compute as for b0
            delta1_for_first[index] = compute_delta1(solution, instance,
                                                     M1, sequence, index, position,
                                                     1, instance.nb_HPRC)
        end
        for index in (b0+1):solution.length
            delta2_for_first = compute_delta2(solution, instance,
                                              M1, sequence, index, position,
                                              1, instance.nb_HPRC,
                                              delta2_for_first)
        end
        # The cost is variation of deletion + delta1 (new sequence) + delta2 (modified sequences)
        for i in b0:solution.length
            cost_on_objective[i, 1] = delta1_for_first[i] + delta2_for_first[i] + violations_caused_on_first
        end
    end

    # Cost on LPRC
    if objectives[2]
        delta2_for_second[b0] = compute_delta2_for_b0(
            solution, instance,
            M1, sequence, position,
            instance.nb_HPRC+1,
            instance.nb_HPRC+instance.nb_LPRC)
        delta1_for_second[b0] = compute_delta1(
            solution, instance,
            M1, sequence, b0, position,
            instance.nb_HPRC+1, instance.nb_HPRC+instance.nb_LPRC)
        for index in (b0+1):solution.length
            # Delta 1 is compute as for b0
            delta1_for_second[index] = compute_delta1(
                solution, instance,
                M1, sequence, index, position,
                instance.nb_HPRC+1, instance.nb_HPRC+instance.nb_LPRC)
        end
        for index in (b0+1):solution.length
            delta2_for_second = compute_delta2(
                solution, instance,
                M1, sequence, index, position,
                (instance.nb_HPRC+1), (instance.nb_HPRC+instance.nb_LPRC),
                delta2_for_second)
        end

        # The cost is variation of deletion + delta1 (new sequence) + delta2 (modified sequences)
        for i in b0:solution.length
            cost_on_objective[i, 2] = delta1_for_second[i] + delta2_for_second[i] + violations_caused_on_second
        end
    end

    # Cost on PCC
    if objectives[3]
        for index in b0:solution.length
            car = solution.sequence[position]

            # Due to deletion
            if position > b0
                if instance.color_code[car] != instance.color_code[solution.sequence[position-1]]
                    cost_on_objective[index, 3] -= 1
                end

                if position < solution.length
                    if instance.color_code[solution.sequence[position-1]] != instance.color_code[solution.sequence[position+1]]
                        cost_on_objective[index, 3] += 1
                    end
                end
            end
            if position < solution.length
                if instance.color_code[car] != instance.color_code[solution.sequence[position+1]]
                    cost_on_objective[index, 3] -= 1
                end
            end

            # Due to insertion
            if index > b0 && index < solution.length
                if instance.color_code[sequence[index-1]] != instance.color_code[sequence[index]]
                    cost_on_objective[index, 3] -= 1
                end
            end
            if index > b0
                if instance.color_code[sequence[index-1]] != instance.color_code[car]
                    cost_on_objective[index, 3] += 1
                end
            end
            if index < solution.length
                if instance.color_code[sequence[index]] != instance.color_code[car]
                    cost_on_objective[index, 3] += 1
                end
            end
        end
    end

    # #DEBUG and TODO: Comment / uncomment for debbuging.
    # #                Test the quality of the variation computed.
    # s1 = deepcopy(solution)
    # s1.sequence[car_pos_a], s1.sequence[car_pos_b] = s1.sequence[car_pos_b], s1.sequence[car_pos_a]
    # update_matrices!(s1, instance)
    #
    # real_cost = cost(s1, instance, objectives) - cost(solution, instance, objectives)
    # if real_cost != cost_on_objective
    #     println("\nIn cost_move_exchange for ", car_pos_a, " and ", car_pos_b)
    #     println("\tError (real vs computed) : ", real_cost, " - ", cost_on_objective)
    # end

    return cost_on_objective
end

"""
    cost_move_insertion(solution::Solution, car_pos_a::Int,
                        instance::Instance, objective::Int)
Returns the cost of the insertion of the car `car_pos_a` for all valid positions
with respect to objectives 1 to `objective`. A negative cost means that the move
is interesting with respect to treated objective.
"""
cost_move_insertion(solution::Solution, position::Int, instance::Instance, objective::Int) =
    cost_move_insertion(solution, position, instance, [trues(objective) ; falses(3-objective)])

        #-------------------------------------------------------#
        #                                                       #
        #                   Factorization                       #
        #         #TODO: Need refactoring and checking          #
        #-------------------------------------------------------#

"""
    update_lines_remove!(solution::Solution, instance::Instance,
                         M1::Array{Int,2}, car_pos_a::Int,
                         first_line::Int, last_line::Int)
For factorization
"""
function update_lines_remove!(solution::Solution, instance::Instance,
                              M1::Array{Int,2}, car_pos_a::Int,
                              first_line::Int, last_line::Int)

    violations_caused = 0 # Removing car_pos_a will change the number of violations.
    due_to_sequence_removed = 0

    b0 = instance.nb_late_prec_day+1
    # forall option there is 1 to last_sequence_intact sequences unchanged
    for option in first_line:last_line
        last_sequence_intact = car_pos_a - instance.RC_q[option]
        if last_sequence_intact < 0
            last_sequence_intact = 0
        end
        # for car in 1:last_sequence_intact : do nothing

        # forall sequences from last_sequence_intact+1 to pos_a-1,
        for index in (last_sequence_intact+1):(car_pos_a-1)
            # one less violation if car has the option
            if instance.RC_flag[solution.sequence[car_pos_a], option]
                if M1[option, index] > instance.RC_p[option] #
                    violations_caused -= 1
                end
                M1[option, index] -= 1
            end
            # one more violation if new car reached has the option
            new_index_reached = index + instance.RC_q[option]
            if new_index_reached <= solution.length #TODO: is it not < ? Don't think so
                if instance.RC_flag[solution.sequence[new_index_reached], option]
                    if M1[option, index] >= instance.RC_p[option]
                        violations_caused += 1
                    end
                    M1[option, index] += 1
                end
            end
        end

        # sequence pos_a is deleted, every violation of this sequence disappear
        violations_caused -= max(0, solution.M1[option, car_pos_a] - instance.RC_p[option])

        if (car_pos_a == b0)
            due_to_sequence_removed -= (max(0, solution.M1[option, car_pos_a] - instance.RC_p[option]))
        end

        # forall sequences from pos_a+1 to n are intact but shifted,
        for index in (car_pos_a):(instance.nb_cars-1)
            M1[option, index] = M1[option, index+1]
        end

        M1[option, instance.nb_cars] = 0 # There is no car here
    end

    return M1, violations_caused
end

"""
    compute_delta2_for_b0(solution::Solution, instance::Instance,
                          M1::Array{Int,2}, sequence::Array{Int,1},
                          car_pos_a::Int, first_line::Int, last_line::Int)
For factorization
"""
function compute_delta2_for_b0(solution::Solution, instance::Instance,
                               M1::Array{Int,2}, sequence::Array{Int,1},
                               car_pos_a::Int, first_line::Int, last_line::Int)

    b0 = instance.nb_late_prec_day+1
    delta2 = 0

    C_in = solution.sequence[car_pos_a]

    for option in first_line:last_line

        #there is q(oj)-1 modified sequences...
        first_modified_sequence = (b0 - instance.RC_q[option]) + 1
        if first_modified_sequence < 1 # must be a valid index
            first_modified_sequence = 1
        end

        for modified_sequence in first_modified_sequence:(b0-1)
            new_unreached_index = modified_sequence + instance.RC_q[option] - 1

            # car force another one to quit the sequence
            if new_unreached_index <= solution.length-1 # must be a valid index
                C_out = sequence[new_unreached_index]
                # C_in has different option than C_out -> variation in violations number
                if xor(instance.RC_flag[C_in, option], instance.RC_flag[C_out, option])
                    if instance.RC_flag[C_out, option] && M1[option, modified_sequence] > instance.RC_p[option]
                        delta2 -= 1
                    end
                    if instance.RC_flag[C_in, option] && M1[option, modified_sequence] >= instance.RC_p[option]
                        delta2 += 1
                    end
                end
            else # it is an ending sequence, there is no new unreached car
                if instance.RC_flag[C_in, option] && M1[option, modified_sequence] >= instance.RC_p[option]
                    delta2 += 1
                end
            end
        end
    end

    return delta2
end

"""
    compute_delta1(solution::Solution, instance::Instance,
                   M1::Array{Int,2}, sequence::Array{Int,1}, cursor::Int,
                   car_pos_a::Int, first_line::Int, last_line::Int)
For factorization
"""
function compute_delta1(solution::Solution, instance::Instance,
                        M1::Array{Int,2}, sequence::Array{Int,1}, cursor::Int,
                        car_pos_a::Int, first_line::Int, last_line::Int)

    delta1 = 0
    for option in first_line:last_line
        # ... and one new sequence at cursor
        # This is basically sequence at cursor without the (cursor + instance.RC_q[option] - 1)-th index (ejected by car)
        first_unreached_index = cursor + instance.RC_q[option] - 1
        # has different option than unreached -> variation in violations number
        # car is in ...
        variation = 0
        if instance.RC_flag[solution.sequence[car_pos_a], option]
            variation += 1
        end
        # ... and first_unreached_index (if it exists) is not
        if first_unreached_index <= (solution.length-1)
            if instance.RC_flag[sequence[first_unreached_index], option]
                variation -= 1
            end
        end
        delta1 += max(0, (M1[option, cursor] + variation) - instance.RC_p[option])
    end

    return delta1
end

"""
    compute_delta2(solution::Solution, instance::Instance,
                   M1::Array{Int,2}, sequence::Array{Int,1}, index::Int,
                   car_pos_a::Int, first_line::Int, last_line::Int,
                   delta2_for_objective)
For factorization
"""
function compute_delta2(solution::Solution, instance::Instance,
                        M1::Array{Int,2}, sequence::Array{Int,1}, index::Int,
                        car_pos_a::Int, first_line::Int, last_line::Int,
                        delta2_for_objective)
    #TODO: How to handle this case?
    if index == solution.length
        C_insert = solution.sequence[car_pos_a]
        for option in first_line:last_line
            sequence_unreaching_it = index - instance.RC_q[option]
            if sequence_unreaching_it > 0
                for i in (sequence_unreaching_it+1):solution.length
                    if (M1[option, i] >= instance.RC_p[option] && instance.RC_flag[C_insert, option])
                        delta2_for_objective[index] += 1
                    end
                end
            end
        end
    else
        delta2_for_objective[index] = delta2_for_objective[index-1]

        C_insert = solution.sequence[car_pos_a] # the car we want to insert

        for option in first_line:last_line
            # What is called in the article: delta_{b-q(o_j)}
            # there is one sequence too far from us now
            sequence_unreaching_it = index - instance.RC_q[option]
            if sequence_unreaching_it > 0 && index < solution.length
                C_in = sequence[index-1] # It was previously push out (shifted to index in the last calcul of d2)
                if xor(instance.RC_flag[C_insert, option], instance.RC_flag[C_in, option])
                    # The last sequence is now more violated than before ?
                    if instance.RC_flag[C_in, option]
                        if M1[option, sequence_unreaching_it] > instance.RC_p[option]
                            delta2_for_objective[index] += 1
                        end
                    end
                    # Car was counted as a violation
                    if instance.RC_flag[C_insert, option]
                        if M1[option, sequence_unreaching_it] >= instance.RC_p[option]
                            delta2_for_objective[index] -= 1
                        end
                    end
                end
            end

            # What is called in the article: delta_{b-1}
            new_modified_sequence = index - 1
            first_unreached_index = index + instance.RC_q[option] - 2
            if new_modified_sequence > 0 && first_unreached_index <= (solution.length-1)
                C_out = sequence[first_unreached_index]
                if xor(instance.RC_flag[C_out, option], instance.RC_flag[C_insert, option])
                    # The new sequence is now more violated than before
                    if instance.RC_flag[C_out, option]
                        if M1[option, new_modified_sequence] > instance.RC_p[option]
                            delta2_for_objective[index] -= 1
                        end
                    end
                    # Car dropped one violation
                    if instance.RC_flag[C_insert, option]
                        if M1[option, new_modified_sequence] >= instance.RC_p[option]
                            delta2_for_objective[index] += 1
                        end
                    end
                end
            else # sequence was at the end, no-one is getting out
                # Car may cause one violation
                if instance.RC_flag[C_insert, option]
                    if M1[option, new_modified_sequence] >= instance.RC_p[option]
                        delta2_for_objective[index] += 1
                    end
                end
            end
        end
    end

    return delta2_for_objective
end


text:
# Empty matrix and vector for function evaluation
const emptyfvec = Array{Float64}(undef, 0)
const emptyfmat = Array{Float64}(undef, 0, 0)

# Define Abstract NLP Evaluator for STO problem
abstract type STOev <: MathOptInterface.AbstractNLPEvaluator end

# Linear Case
mutable struct linSTOev <: STOev
  # Parameters
  x0::Array{Float64,1}                    # Initial State x0
  nx::Int                                 # State dimension
  A::Array{Float64,3}                     # Dynamics Matrices
  N::Int                                  # Number of Switching Times
  t0::Float64                             # Initial Time
  tf::Float64                             # Final Time (from setup)
  tfdelta::Float64                        # Final Time (from current delta)
  Q::Array{Float64,2}                     # State Cost
  E::Array{Float64,2}                    # Final State Cost

  # Grid Variables
  ngrid::Int64                             # Number of grid points
  tgrid::Array{Float64,1}                  # Array of the grid
  tvec::Array{Float64,1}                   # Array of the complete grid with sw times
  tauIdx::Array{Int,1}                     # Array of the position of switching times in the complete grid vector tvec
  tgridIdx::Array{Int,1}                   # Array of the position of grid points in the complete grid vector tvec
  deltacomplete::Array{Float64,1}          # Array of the complete sw intervals

  # Precomputed Values
  V::Array{Complex{Float64}, 3}            # Dynamics Matrices decomp V
  invV::Array{Complex{Float64}, 3}         # Dynamics Matrices decomp V^-1
  D::Array{Complex{Float64}, 2}            # Dynamics Matrices decomp D
  isDiag::Array{Bool,1}                    # Vector of boolean variables stating if dynamics matrices are diagonalizable

  IndTril::Array{Int, 1}                   # Single Element of Lower Triangular Matrices (Hessian)
  Itril::Array{Int, 1}                     # Double Element Indeces of Lower Triangular Matrices (Hessian)
  Jtril::Array{Int, 1}                     # Double Element Indeces of Lower Triangular Matrices (Hessian)
  Ag::Array{Float64, 2}                    # Linear Constraints Matrix
  Ig::Array{Int, 1}                        # Index Linear Constraints
  Jg::Array{Int, 1}                        # Index Linear Constraints
  Vg::Array{Float64, 1}                    # Value Linear Constraints
  bg::Array{Float64, 1}                    # Constant Term Linear Constraints

  # Shared Data Between Functions
  deltafun_prev::Array{Float64, 1}            # Previous switching intervals for function evaluation
  deltagrad_prev::Array{Float64, 1}             # Previous switching intervals for gradient evaluation
  deltahess_prev::Array{Float64, 1}             # Previous switching intervals for hessian evaluation
  xpts::Array{Float64, 2}                  # States at Switching Times
  expMat::Array{Float64, 3}                # Matrix Exponentials
  Phi::Array{Float64, 4}                   # Matrix of State Transition Matrices
  M::Array{Float64, 3}                     # Integrals over Switching Intervals
  S::Array{Float64, 3}                     # S Matrices for each interval
  C::Array{Float64, 3}                     # C Matrices for each interval

  # Store Evaluations
  obj::Array{Float64, 1}                   # Store Objective Value
  deltaval::Array{Float64, 2}              # Store Iterates
  nobjeval::Int                            # Number of objective function evaluations
  ngradeval::Int                           # Number of gradient evaluations
  nhesseval::Int                           # Number of hessian evaluations
end

# Nonlinear Case
mutable struct nlinSTOev <: STOev
  # Parameters
  x0::Array{Float64,1}                     # Initial State x0
  nx::Int                                  # State dimension
  A::Array{Float64,3}                      # Linearized dynamics matrices
  N::Int                                   # Number of switching times
  t0::Float64                              # Initial Time
  tf::Float64                              # Final Time (from setup)
  tfdelta::Float64                         # Final Time (from current delta)
  Q::Array{Float64,2}                      # State Cost
  E::Array{Float64,2}                     # Final State Cost
  uvec::Array{Float64,2}                   # Array of switching inputs
  # nartsw::Int64                          # Number of artificial switches

  # Grid Variables
  ngrid::Int64                             # Number of grid points
  tgrid::Array{Float64,1}                  # Array of the grid
  tvec::Array{Float64,1}                   # Array of the complete grid with sw times
  tauIdx::Array{Int,1}                     # Array of the position of switching times in the complete grid vector tvec
  tgridIdx::Array{Int,1}                   # Array of the position of grid points in the complete grid vector tvec
  deltacomplete::Array{Float64,1}          # Array of the complete sw intervals including the grid ones

  # Nonlinear Dynamics and Derivatives Functions
  nonlin_dyn::Function
  nonlin_dyn_deriv::Function

  # Precomputed Values
  IndTril::Array{Int, 1}                  # Single Element Indeces of Lower Triangular Matrices (Hessian)
  Itril::Array{Int, 1}                    # Double Element Indeces of Lower Triangular Matrices (Hessian)
  Jtril::Array{Int, 1}                    # Double Element Indeces of Lower Triangular Matrices (Hessian)
  Ag::Array{Float64, 2}                   # Linear Constraints (Vector)
  Ig::Array{Int, 1}                       # Index Linear Constraint
  Jg::Array{Int, 1}                       # Index Linear Constraint
  Vg::Array{Float64, 1}                   # Value Linear Constraints
  bg::Array{Float64, 1}                   # Constant Term Linear Constraints

  # Shared Data Between Functions
  deltafun_prev::Array{Float64, 1}            # Previous switching intervals for function evaluation
  deltagrad_prev::Array{Float64, 1}             # Previous switching intervals for gradient evaluation
  deltahess_prev::Array{Float64, 1}             # Previous switching intervals for hessian evaluation
  xpts::Array{Float64, 2}                 # States at Switching Times
  expMat::Array{Float64, 3}               # Matrix Exponentials
  Phi::Array{Float64, 4}                  # Matrix of State Transition Matrices
  M::Array{Float64, 3}                    # Integrals over Switching Intervals
  S::Array{Float64, 3}                    # S Matrices for each interval
  C::Array{Float64, 3}                     # C Matrices for each interval

  # Store Evaluations
  obj::Array{Float64, 1}                   # Store Objective Value
  deltaval::Array{Float64, 2}              # Store Iterates
  nobjeval::Int                            # Number of objective function evaluations
  ngradeval::Int                           # Number of gradient evaluations
  nhesseval::Int                           # Number of hessian evaluations
end

# Create switching time optimization (STO) abstract type
abstract type STO end

mutable struct linSTO <: STO  # Linear STO type
  model::MathOptInterface.AbstractOptimizer   # Nonlinear Program Model
  STOev::linSTOev                             # NLP Evaluator for linear STO

  # Data Obtained after Optimization
  delta::Array{Float64,1}                     # Optimal Switching Intervals
  delta_MOI::Array{MathOptInterface.VariableIndex,1}
  tau::Array{Float64,1}                       # Optimal Switching Times
  objval::Float64                             # Optimal Value of Cost Function
  stat::MathOptInterface.TerminationStatusCode# Status of Opt Problem
  soltime::Float64                            # Time Required to solve Opt

  # Inner Contructor for Incomplete Initialization
  linSTO(model, STOev, delta, delta_MOI) = new(model, STOev, delta, delta_MOI)

end

mutable struct nlinSTO <: STO  # Nonlinear STO type
  model::MathOptInterface.AbstractOptimizer  # Nonlinear Program Model
  STOev::nlinSTOev                            # NLP Evaluator for nonlinear STO

  # Data Obtained After Optimization
  delta::Array{Float64,1}                     # Optimal Switching Intervals
  delta_MOI::Array{MathOptInterface.VariableIndex,1}
  tau::Array{Float64,1}                       # Optimal Switching Times
  objval::Float64                             # Optimal Value of Cost Function
  stat::MathOptInterface.TerminationStatusCode# Status of Opt Problem
  soltime::Float64                            # Time Required to solve Opt

  # Inner Contructor for Incomplete Initialization
  nlinSTO(model, STOev, delta, delta_MOI) = new(model, STOev, delta, delta_MOI)
end


text:
# Test relief.jl

# Tests could be more comprehensive (of course).  In some cases all
# that is done is to execute the function in some way just to make
# sure no exceptions are thrown as a result of something breaking with
# a new version of Julia.

println("testing relief")

img = sineramp()

rimg = relief(img)
az = 45
el = 45
gradscale = 2
rimg = relief(img, az, el, gradscale, cmap("L10")[1])
rimg = relief(img, az, el, gradscale)


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

import Base: launch, manage, connect, exit

mutable struct UnixDomainCM <: ClusterManager
    np::Integer
end

function launch(manager::UnixDomainCM, params::Dict, launched::Array, c::Condition)
#    println("launch $(manager.np)")
    cookie = Base.cluster_cookie()
    for i in 1:manager.np
        sockname = tempname()
        try
            cmd = `$(params[:exename]) --startup-file=no $(@__FILE__) udwrkr $sockname $cookie`
            io, pobj = open(cmd, "r")

            wconfig = WorkerConfig()
            wconfig.userdata = Dict(:sockname=>sockname, :io=>io, :process=>pobj)
            push!(launched, wconfig)
            notify(c)
        catch e
            println(e)
        end
    end
end

function connect(manager::UnixDomainCM, pid::Int, config::WorkerConfig)
    if myid() == 1
#        println("connect_m2w")
        config.connect_at = get(config.userdata)[:sockname] # This will be useful in the worker-to-worker connection setup.

        print_worker_stdout(get(config.userdata)[:io], pid)
    else
#        println("connect_w2w")
        sockname = get(config.connect_at)
        config.userdata = Dict{Symbol, Any}(:sockname=>sockname)
    end

    t = time()
    while true
        try
            address = get(config.userdata)[:sockname]
            if isa(address, Tuple)
                sock = connect(address...)
            else
                sock = connect(address)
            end
            return (sock, sock)
        catch e
            if (time() - t) > 30.0
                rethrow(e)
            else
                sleep(0.1)
            end
        end
    end

end

# WORKER
function start_worker(sockname, cookie)
    Base.init_worker(cookie, UnixDomainCM(0))

    srvr = listen(sockname)
    while true
        sock = accept(srvr)
        Base.process_messages(sock, sock)
    end
end

function manage(manager::UnixDomainCM, id::Int, config::WorkerConfig, op)
    # Does not seem to be required, filesystem entry cleanup is happening automatically on process exit
#     if op == :deregister
#         try
#             rm(get(config.userdata)[:sockname])
#         end
#     end
    nothing
end

function print_worker_stdout(io, pid)
    @schedule while !eof(io)
        line = readline(io)
        println("\tFrom worker $(pid):\t$line")
    end
end

if (length(ARGS) > 0) && (ARGS[1] == "udwrkr")
    # script has been launched as a worker
    start_worker(ARGS[2], ARGS[3])
end


text:
# StrICU.jl - Wrapper for ICU (International Components for Unicode) library

# Some content of the documentation strings was derived from the ICU header files.
# (Those portions copyright (C) 1996-2015, International Business Machines Corporation and others)

"""
    StrICU (International Components for Unicode) Wrapper
"""
module StrICU

using ModuleInterfaceTools

using ICU_jll

# for now, just leave this hardcoded, until I can see how to get the version from ICU_jll
const suffix = "_68"
const iculib = ICU_jll.libicuuc
const iculibi18n = ICU_jll.libicui18n

@api extend! StrBase

finalizer(o, f::Function) = Base.finalizer(f, o)

import Base: parse, get, close

export ICU
const ICU = StrICU

const cvt_utf8 = utf8
const cvt_utf16 = utf16
export cvt_utf8, cvt_utf16

const ByteStr   = Union{ASCIIStr, UTF8Str, String}

const WordStringCSE = Union{UCS2CSE, _UCS2CSE, UTF16CSE}
const WordStrings = Str{<:WordStringCSE}

export set_locale!

function __init__()
    set_locale!("")
end

_libicu(s, lib, p) = ( Symbol(string(p, s, suffix)), lib )

const UBool      = Int8
const UChar      = UInt16
const UErrorCode = Int32
const U_PARSE_CONTEXT_LEN = 16
const U_PARSE_TUPLE = ntuple((i)->0x20, U_PARSE_CONTEXT_LEN)

struct UParseError
    line::UErrorCode                                 # The line on which the error occured
    offset::UErrorCode                               # The character offset to the error
    preContext::NTuple{U_PARSE_CONTEXT_LEN, UInt8}   # Textual context before the error
    postContext::NTuple{U_PARSE_CONTEXT_LEN, UInt8}  # The error itself and/or textual context after the error
end
UParseError() = UParseError(0, 0, U_PARSE_TUPLE, U_PARSE_TUPLE)

FAILURE(x::Integer) = x > 0
SUCCESS(x::Integer) = x <= 0
U_BUFFER_OVERFLOW_ERROR = 15

const locale   = ASCIIStr[""]

include("utext.jl")
include("ustring.jl")
include("ubrk.jl")
include("ucnv.jl")
include("ucol.jl")
include("ucsdet.jl")
include("udat.jl")
include("ucal.jl")
include("ucasemap.jl")

@api freeze

end # module StrICU


text:
function templateTestTmapfoldr1(arr::AbstractArray)::Bool
	init = randInit();
	return mapfoldr(x->2x, +, arr; init = init) == tmapfoldr(x->2x, +, arr; init = init);
end

function testTmapfoldr1Stable()::Bool
	return templateTestTmapfoldr1(randArray());
end

function testTmapfoldr1Unstable()::Bool
	arr = randStrArray();
	str = lipsum();
	init = randstring(rand(1:1000));
	return mapfoldr(x->str[1:x], *, arr; init = init) == tmapfoldr(x->str[1:x], *, arr; init = init);
end

function testTmapfoldr1Offset()::Bool
	return templateTestTmapfoldr1(offsetArray());
end

function testTmapfoldr1Empty()::Bool
	return templateTestTmapfoldr1([]);
end

#####

function templateTestTmapfoldr2(arr::AbstractArray)::Bool
	return mapfoldr(x->2x, +, arr) == tmapfoldr(x->2x, +, arr);
end

function testTmapfoldr2Stable()::Bool
	return templateTestTmapfoldr2(randArray());
end

function testTmapfoldr2Unstable()::Bool
	arr = randStrArray();
	str = lipsum();
	return mapfoldr(x->str[1:x], *, arr) == tmapfoldr(x->str[1:x], *, arr);
end

function testTmapfoldr2Offset()::Bool
	return templateTestTmapfoldr2(offsetArray());
end

function testTmapfoldr2Empty()
	tmapfoldr(x->2x, (x, y)->x, []);
end


text:
## test data available for use in PopGen
export @nancycats, @gulfsharks, dataset

"""
    dataset(::String)
Load an example dataset from either `"gulfsharks"` (SNP) or `"nancycats"` (microsatellite).
Can also use `"sharks"` and `"cats"` as shorthands. Use `?@nancycats` and
`?@gulfsharks` to learn more about these datasets. Use `PopGenCore.dataset` if using
`PopGenCore.jl` directly.

### Example
```
ncats = PopGen.dataset("nancycats")
gsharks = PopGen.dataset("sharks")
```
"""
function dataset(name::String)
    if lowercase(name) in ["nancycats", "cats"]
        filename = normpath(joinpath(@__DIR__,"..","data", "nancycats.gen"))
    elseif lowercase(name) in ["gulfsharks", "sharks"]
        filename = normpath(joinpath(@__DIR__,"..","data", "gulfsharks.csv"))
    else
        throw(ArgumentError("Please choose either the \"nancycats\" or \"gulfsharks\" datasets"))
    end
    read(filename, silent = true, allow_monomorphic = true)
end

"""
    @nancycats
Returns `PopData` of corresponding "nancycats" dataset as featured in
the R package `adegenet`. This dataset is composed of 9 microsatellite
loci in 237 individuals across 17 populations.

Example:
```
ncats = @nancycats
```
"""
macro nancycats()
    return esc(quote
        filename = normpath(joinpath(@__DIR__,"..","data", "nancycats.gen"))
        genepop(filename, silent = true, allow_monomorphic = true)
    end)
end

"""
    @gulfsharks
Returns `PopData` corresponding the Blacknose shark dataset as used in
Dimens et al. 2019. This is a mid-sized SNP dataset of 2209 haplotypes
in 212 individuals, across 7 populations.

Example:

```
sharks = @gulfsharks
```
"""
macro gulfsharks()
    return esc(quote
        filename = normpath(joinpath(@__DIR__,"..","data", "gulfsharks.csv"))
        delimited(filename, silent = true, allow_monomorphic = true)
    end)
end


text:
@testset "active pattern" begin
    @testset "regular case" begin
        @active LessThan0(x) begin
            if x > 0
                nothing
            else
                x
            end
        end

        @test (@match 15 begin
            LessThan0(_) => :a
            _ => :b
        end) === :b

        @test (@match -15 begin
            LessThan0(a) => a
            _ => 0
        end) == -15
    end
    @testset "parametric case" begin

        @active Re{r :: Regex}(x) begin
             match(r, x)
        end

        @test (@match "123" begin
            Re{r"\d+"}(x) => x.match
            _ => @error ""
        end) == "123"

        @test_skip @match "abc" begin
            Re{r"\d+"}(x) => x
            _ => @error ""
        end
    end

    @testset "custom pattern for given structs" begin
        @eval struct Interval end

        @active internal Interval{a, b}(arg) begin
            a <= arg <= b
        end

        @use Enum

        @active visible in (@__MODULE__) IsEven(a) begin
            a % 2 === 0
        end

        function parity(x)
            @match x begin
                IsEven => :even
                _ => :odd
            end
        end
        @test :even === parity(4)
        @test :odd === parity(3)

        @test 2 == @match 3 begin
            Interval{1, 2} => 1
            Interval{3, 4} => 2
            Interval{5, 6} => 3
            Interval{7, 8} => 4
        end

    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CUDNN_CUDA111_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("CUDNN_CUDA111")
JLLWrappers.@generate_main_file("CUDNN_CUDA111", UUID("2839a42b-739d-52bf-a6f0-df8e4ebc59da"))
end  # module CUDNN_CUDA111_jll


text:
using Documenter, DocumenterLaTeX

makedocs(
	sitename = "MimiFUND.jl",
	pages = [
		"Home" => "index.md",
		"Science" => "science.md",
		"Tables" => "tables.md"]
)

deploydocs(
    repo = "github.com/fund-model/MimiFUND.jl.git"
)


text:
using Plots
push!(LOAD_PATH, "./src")
using Basis2DTri
using Basis1D
using LaTeXStrings
theme(:wong)
gr(aspect_ratio=1,legend=nothing,axis=nothing,border=nothing,ticks=nothing,
   markerstrokewidth=0,markersize=1)

r1D, w1D = gauss_lobatto_quad(0,0,4)

plot([-1;-1],[-1;1],seriescolor=:black,lw=3,label = nothing)
plot!([-1;1],[1;1],seriescolor=:black,lw=3,label = nothing)
plot!([1;1],[1;-1],seriescolor=:black,lw=3,label = nothing)
plot!([1;-1],[-1;-1],seriescolor=:black,lw=3,label = nothing)
scatter!(r1D,ones(size(r1D)),markersize=7,markercolor=:royalblue1,markerstrokecolor=:black,markerstrokewidth=2)
scatter!(r1D,r1D[2]*ones(size(r1D)),markersize=7,markercolor=:royalblue1,markerstrokecolor=:black,markerstrokewidth=2)
scatter!(r1D,zeros(size(r1D)),markersize=7,markercolor=:royalblue1,markerstrokecolor=:black,markerstrokewidth=2)
scatter!(r1D,r1D[4]*ones(size(r1D)),markersize=7,markercolor=:royalblue1,markerstrokecolor=:black,markerstrokewidth=2)
scatter!(r1D,-1*ones(size(r1D)),markersize=7,markercolor=:royalblue1,markerstrokecolor=:black,markerstrokewidth=2)

plot!([-1;-1].+2.5,[-1;1],seriescolor=:black,lw=3,label = nothing)
plot!([-1;1].+2.5,[1;1],seriescolor=:black,lw=3,label = nothing)
plot!([1;1].+2.5,[1;-1],seriescolor=:black,lw=3,label = nothing)
plot!([1;-1].+2.5,[-1;-1],seriescolor=:black,lw=3,label = nothing)
scatter!(r1D.+2.5,ones(size(r1D)),markersize=7,markercolor=:darkorange1,markerstrokecolor=:black,markerstrokewidth=2)
scatter!(r1D.+2.5,r1D[2]*ones(size(r1D)),markersize=7,markercolor=:darkorange1,markerstrokecolor=:black,markerstrokewidth=2)
scatter!(r1D.+2.5,zeros(size(r1D)),markersize=7,markercolor=:darkorange1,markerstrokecolor=:black,markerstrokewidth=2)
scatter!(r1D.+2.5,r1D[4]*ones(size(r1D)),markersize=7,markercolor=:darkorange1,markerstrokecolor=:black,markerstrokewidth=2)
scatter!(r1D.+2.5,-1*ones(size(r1D)),markersize=7,markercolor=:darkorange1,markerstrokecolor=:black,markerstrokewidth=2)

savefig("LGL_spectral_elem_2.png")

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
export simulate_bb

"""
Makes a random move. It does NOT add a new tile
"""
function randmove(bitboard::Bitboard)::Bitboard
    dirs = collect(DIRS)

    for i in 4:-1:1
        j = rand(1:i)
        dir = dirs[j]
        new_bitboard = move(bitboard, dir)
        if new_bitboard != bitboard
            return new_bitboard
        end
        dirs[j] = dirs[i]
    end
    return bitboard
end

"""
Simulate a game given the bitboard state till the end
"""
function simulate_bb(bitboard::Bitboard=initbboard())
    while true
        new_bitboard = randmove(bitboard)
        if new_bitboard != bitboard
            bitboard = add_tile(new_bitboard)
            @assert bitboard != new_bitboard
        else
            # lost if no move can be made
            return bitboard
        end
    end
end


text:
using StaticArrays: SMatrix

using ..Continuations: AbstractContinuationProblem, AbstractContinuationSolver,
    AbstractContinuationCache,
    as, SweepSetup, ContinuationSweep, ContinuationSolution,
    ContinuationCache, ContinuationOptions, ContinuationSolver,
    residual_jacobian!
import ..Continuations: step!, new_sweep!

using ..BifurcationsBase: timekind, Continuous, Discrete
import ..BifurcationsBase: TimeKind

using ..BifurcationsBase
using ..BifurcationsBase: SpecialPoint, SpecialPointInterval,
    BifurcationSweep, BifurcationSolution, BifurcationSolver,
    BifurcationProblem, BifurcationCache,
    eigvals_prototype, allocate_sweep!, check_sweep_length, record!
import ..BifurcationsBase: analyze!, re_analyze!

"""
Codimension 2 bifurcation points and special points.

These points are supported:

- [Cusp bifurcation]
- [Bogdanov-Takens bifurcation]
- [Bautin Bifurcation]
- [Fold-Hopf Bifurcation]
- [Hopf-Hopf Bifurcation]

[Cusp bifurcation]:
http://www.scholarpedia.org/article/Cusp_bifurcation
[Bogdanov-Takens bifurcation]:
http://www.scholarpedia.org/article/Bogdanov-Takens_bifurcation
[Bautin Bifurcation]:
http://www.scholarpedia.org/article/Bautin_bifurcation
[Fold-Hopf Bifurcation]:
http://www.scholarpedia.org/article/Fold-Hopf_bifurcation
[Hopf-Hopf Bifurcation]:
http://www.scholarpedia.org/article/Hopf-Hopf_bifurcation

* http://www.scholarpedia.org/article/Bifurcation
* https://www.encyclopediaofmath.org/index.php/Codimension-two_bifurcations
"""
module PointTypes
@enum(
    PointType,
    none,
    # initial,
    cusp,
    bautin,
    bogdanov_takens,
    fold_hopf,
    hopf_hopf,
)
end  # module
using .PointTypes: PointType

const AbstractCodim2SpecialPoint{tkind <: TimeKind} =
    AbstractSpecialPoint{tkind, PointType}

const Codim2Sweep{tkind <: TimeKind, ptType <: PointType} =
    BifurcationSweep{tkind, ptType}
const Codim2Solution{S <: ContinuationSolution, W <: Codim2Sweep} =
    BifurcationSolution{S, W}
# The constraints `<: TimeKind` and `<: ContinuationSolution` are
# important.
# See: [[../codim1/solver.jl::Codim1Sweep]]

BifurcationsBase.point_type_type(::Codim2Problem) = PointType

BifurcationsBase.regular_point(::Type{PointType}) = PointTypes.none

function BifurcationsBase.eigvals_prototype(prob::Codim2Problem,
                                            cache::ContinuationCache)
    N = length(cache.u) ÷ 2 - 1
    return cache.u[1:N] .+ 1im
end
# TODO: improve it for SVector

mutable struct Codim2Cache{P, C <: ContinuationCache{P},
                           JType, eType,
                           } <: BifurcationCache{P}
    super::C
    J::JType
    det::Float64
    prev_det::Float64
    eigvals::eType
    prev_eigvals::eType
    point_type::PointType

    """
    The quadratic coefficient ``a(0)``.
    See: http://www.scholarpedia.org/article/Saddle-node_bifurcation
    """
    quadratic_coefficient::Float64
    prev_quadratic_coefficient::Float64

    """
    The first Lyapunov coefficient ``l_1(0)``.
    """
    lyapunov_coefficient::Float64
    prev_lyapunov_coefficient::Float64
end
# TODO: Maybe rename Codim2Cache to SNContCache for something.
# Or put bifurcation type-specific cache in a sub-cache?

function Codim2Cache(super::C,
                     J::JType,
                     eigvals::eType,
                     point_type::PointType = PointTypes.none,
                     ) where {P, C <: ContinuationCache{P},
                              JType, eType}
    return Codim2Cache{P, C, JType, eType}(super, J,
                                           NaN, NaN,
                                           eigvals, copy(eigvals),
                                           point_type,
                                           NaN, NaN,
                                           NaN, NaN)
end
# TODO: Remove this constructor after removing the type parameter `P`.

Codim2Cache(prob::Codim2Problem, super::ContinuationCache) =
    Codim2Cache(
        super,
        ds_jacobian(super),
        copy(eigvals_prototype(prob, super)),
    )

BifurcationsBase.BifurcationCache(prob::Codim2Problem,
                                  super::ContinuationCache) =
    Codim2Cache(prob, super)

const Codim2Solver{
        R <: ContinuationSolver,
        P <: Codim2Problem,
        C <: Codim2Cache,
        S <: Codim2Solution,
        } =
    BifurcationSolver{R, P, C, S}

# TODO: trait
function BifurcationsBase.contkind(prob::DiffEqCodim2Problem)
    if eltype(prob.v0) <: Complex
        return HopfCont()
    else
        return SaddleNodeCont()
    end
end

function re_analyze!(solver::Codim2Solver, u::AbstractVector)
    residual_jacobian!(as(solver.cache, ContinuationCache), u)
    analyze!(solver.cache, solver.opts)

    # Suppress special point recording:
    # It's a bit ugly hack... (communicate by sharing!)
    solver.cache.point_type = PointTypes.none  # TODO: FIX!
    solver.cache.lyapunov_coefficient = NaN    # TODO: FIX!

    record!(solver.sol, solver.cache)
end

function analyze!(cache::Codim2Cache, opts)
    cnt = contkind(cache)
    cache.J = J = ds_jacobian(cache)
    cache.prev_det = cache.det
    cache.det = det(J)
    cache.prev_eigvals = cache.eigvals
    cache.eigvals = ds_eigvals(timekind(cache), J)
    set_quadratic_coefficient!(cnt, cache)
    set_lyapunov_coefficient!(cnt, cache)
    cache.point_type = guess_point_type(cnt, timekind(cache), cache, opts)
    set_augsys_cache!(cache)
end

function BifurcationsBase.post_step!(solver::Codim2Solver)
    cache = solver.cache
    if cache.point_type == PointTypes.bogdanov_takens &&
            as(cache, ContinuationCache).simple_bifurcation
        # Switching at Bogdanov-Takens point requires to switch
        # problem type and so cannot be handled by the continuation
        # solver.  As a quick solution, I'm removing the last
        # SimpleBifurcationInterval (though communication by sharing
        # is ugly...)
        pop!(as(solver, ContinuationSolver).sol.sweeps[end].simple_bifurcation)
        # TODO: Find a better solution to avoid registering
        # Bogdanov-Takens point as a simple bifurcation point.
    end
end

struct FakeCache{P, TU, TJ}
    prob_cache::P
    u::TU
    J::TJ
end

const Cachish = Union{ContinuationCache, FakeCache}

ds_jacobian(solver) = ds_jacobian(as(solver, ContinuationSolver).cache)
ds_jacobian(cache::Codim2Cache) = ds_jacobian(as(cache, ContinuationCache))

function ds_jacobian(cache::Cachish)
    prob = cache.prob_cache.prob  # TODO: interface
    return ds_jacobian(prob, cache.J)
end

ds_jacobian(prob::DiffEqCodim2Problem, J) = ds_jacobian(contkind(prob), J)

function ds_jacobian(ckind::ContinuationKind, J::AbstractMatrix)
    d = dims_from_augsys(J, ckind)
    return @view J[1:d.ds_dim, 1:d.ds_dim]
end

@generated function ds_jacobian(ckind::ContinuationKind,
                                 J::SMatrix{S1, S2, T}
                                 ) where {S1, S2, T}
    d = dims_from_augsys(S2, ckind())
    N = d.ds_dim
    values = [:(J[$i, $j]) for j in 1:N for i in 1:N]
    quote
        return SMatrix{$N, $N, $T}($(values...))
    end
end

function ds_state(cache::Cachish)
    prob = cache.prob_cache.prob  # TODO: interface
    return ds_state(prob, cache.u)
end

function ds_eigvec(cache::Cachish)
    prob = cache.prob_cache.prob  # TODO: interface
    return ds_eigvec(prob, cache.u)
end

# TODO: define ds_f!
ds_f(x, prob, cache::Codim2Cache) = ds_f(x, as(cache, ContinuationCache))

function ds_f(x, cache::Cachish)
    prob = cache.prob_cache.prob  # TODO: interface
    f = prob.de_prob.f  # TODO: interface
    p = modified_param!(prob, cache.u)
    return f(x, p, 0)
end

set_quadratic_coefficient!(::HopfCont, ::Codim2Cache) = nothing

function set_quadratic_coefficient!(::SaddleNodeCont, cache::Codim2Cache)
    cache.prev_quadratic_coefficient = cache.quadratic_coefficient
    cache.quadratic_coefficient = sn_quadratic_coefficient(cache)
end

set_lyapunov_coefficient!(::SaddleNodeCont, ::Codim2Cache) = nothing

function set_lyapunov_coefficient!(::HopfCont, cache::Codim2Cache)
    cache.prev_lyapunov_coefficient = cache.lyapunov_coefficient
    cache.lyapunov_coefficient = first_lyapunov_coefficient(cache)
end


text:
#=

Author: Shunsuke Hori

=#

using QuantEcon
using NLsolve
using NLopt

mutable struct Para{TF <: AbstractFloat,
                    TM <: AbstractMatrix{TF},
                    TV <: AbstractVector{TF}}
    beta::TF
    Pi::TM
    G::TV
    Theta::TV
    transfers::Bool
    U::Function
    Uc::Function
    Ucc::Function
    Un::Function
    Unn::Function
    n_less_than_one::Bool
end

"""
Class returns planner's allocation as a function of the multiplier
on the implementability constraint mu
"""
struct Planners_Allocation_Sequential{TP <: Para,
                                      TI <: Integer,
                                      TV <: AbstractVector}
    para::TP
    mc::MarkovChain
    S::TI
    cFB::TV
    nFB::TV
    XiFB::TV
    zFB::TV
end

"""
Initializes the class from the calibration Para
"""
function Planners_Allocation_Sequential(para::Para)
    beta, Pi, G, Theta =
        para.beta, para.Pi, para.G, para.Theta
    mc = MarkovChain(Pi)
    S = size(Pi, 1) # number of states
    #now find the first best allocation
    cFB, nFB, XiFB, zFB = find_first_best(para, S, 1)

    return Planners_Allocation_Sequential(para, mc, S, cFB, nFB, XiFB, zFB)
end

"""
Find the first best allocation
"""
function find_first_best(para::Para, S::Integer, version::Integer)
    if version != 1 && version != 2
        throw(ArgumentError("version must be 1 or 2"))
    end
    beta, Theta, Uc, Un, G, Pi =
        para.beta, para.Theta, para.Uc, para.Un, para.G, para.Pi
    function res!(z, out)
        c = z[1:S]
        n = z[S+1:end]
        out[1:S] = Theta.*Uc(c, n)+Un(c, n)
        out[S+1:end] = Theta.*n - c - G
    end
    res = nlsolve(res!, 0.5*ones(2*S))

    if converged(res) == false
        error("Could not find first best")
    end

    if version == 1
        cFB = res.zero[1:S]
        nFB = res.zero[S+1:end]
        XiFB = Uc(cFB, nFB) #multiplier on the resource constraint.
        zFB = vcat(cFB, nFB, XiFB)
        return cFB, nFB, XiFB, zFB
    elseif version == 2
        cFB = res.zero[1:S]
        nFB = res.zero[S+1:end]
        IFB = Uc(cFB, nFB).*cFB + Un(cFB, nFB).*nFB
        xFB = \(eye(S) - beta*Pi, IFB)
        zFB = [vcat(cFB[s], xFB[s], xFB) for s in 1:S]
        return cFB, nFB, IFB, xFB, zFB
    end
end

"""
Computes optimal allocation for time ``t\geq 1`` for a given ``\mu``
"""
function time1_allocation(pas::Planners_Allocation_Sequential, mu::Real)
    para, S = pas.para, pas.S
    Theta, beta, Pi, G, Uc, Ucc, Un, Unn =
        para.Theta, para.beta, para.Pi, para.G,
        para.Uc, para.Ucc, para.Un, para.Unn
    function FOC!(z::Vector, out)
        c = z[1:S]
        n = z[S+1:2S]
        Xi = z[2S+1:end]
        out[1:S] = Uc(c,n) - mu*(Ucc(c,n).*c+Uc(c,n)) -Xi #foc c
        out[S+1:2S] = Un(c,n) - mu*(Unn(c,n).*n+Un(c,n)) + Theta.*Xi #foc n
        out[2S+1:end] = Theta.*n - c - G #resource constraint
        return out
    end
    #find the root of the FOC
    res = nlsolve(FOC!, pas.zFB)
    if res.f_converged == false
        error("Could not find LS allocation.")
    end
    z = res.zero
    c, n, Xi = z[1:S], z[S+1:2S], z[2S+1:end]
    #now compute x
    I  = Uc(c,n).*c +  Un(c,n).*n
    x = \(eye(S) - beta*para.Pi, I)
    return c, n, x, Xi
end

"""
Finds the optimal allocation given initial government debt `B_` and state `s_0`
"""
function time0_allocation(pas::Planners_Allocation_Sequential,
                          B_::AbstractFloat, s_0::Integer)
    para = pas.para
    Pi, Theta, G, beta =
        para.Pi, para.Theta, para.G, para.beta
    Uc, Ucc, Un, Unn =
        para.Uc, para.Ucc, para.Un, para.Unn
    #first order conditions of planner's problem
    function FOC!(z, out)
        mu, c, n, Xi = z[1], z[2], z[3], z[4]
        xprime = time1_allocation(pas, mu)[3]
        out .= vcat(
        Uc(c, n).*(c-B_) + Un(c, n).*n + beta*dot(Pi[s_0, :], xprime),
        Uc(c, n) - mu*(Ucc(c, n).*(c-B_) + Uc(c, n)) - Xi,
        Un(c, n) - mu*(Unn(c, n).*n+Un(c, n)) + Theta[s_0].*Xi,
        (Theta.*n - c - G)[s_0]
        )
    end
    #find root
    res = nlsolve(FOC!, [0.0, pas.cFB[s_0], pas.nFB[s_0], pas.XiFB[s_0]])
    if res.f_converged == false
        error("Could not find time 0 LS allocation.")
    end
    return (res.zero...)
end

"""
Find the value associated with multiplier `mu`
"""
function time1_value(pas::Planners_Allocation_Sequential, mu::Real)
    para = pas.para
    c, n, x, Xi = time1_allocation(pas, mu)
    U_val = para.U.(c, n)
    V = \(eye(pas.S) - para.beta*para.Pi, U_val)
    return c, n, x, V
end

"""
Computes Tau given `c`, `n`
"""
function Tau(para::Para, c::Union{Real,Vector}, n::Union{Real,Vector})
    Uc, Un = para.Uc.(c, n), para.Un.(c, n)
    return 1+Un./(para.Theta .* Uc)
end

"""
Simulates planners policies for `T` periods
"""
function simulate(pas::Planners_Allocation_Sequential,
                  B_::AbstractFloat, s_0::Integer, T::Integer,
                  sHist::Union{Vector, Void}=nothing)
    para = pas.para
    Pi, beta, Uc = para.Pi, para.beta, para.Uc
    if sHist == nothing
        sHist = QuantEcon.simulate(pas.mc, T, init=s_0)
    end
    cHist = zeros(T)
    nHist = zeros(T)
    Bhist = zeros(T)
    TauHist = zeros(T)
    muHist = zeros(T)
    RHist = zeros(T-1)
    #time0
    mu, cHist[1], nHist[1], _  = time0_allocation(pas, B_, s_0)
    TauHist[1] = Tau(pas.para, cHist[1], nHist[1])[s_0]
    Bhist[1] = B_
    muHist[1] = mu
    #time 1 onward
    for t in 2:T
        c, n, x, Xi = time1_allocation(pas,mu)
        u_c = Uc(c,n)
        s = sHist[t]
        TauHist[t] = Tau(pas.para, c, n)[s]
        Eu_c = dot(Pi[sHist[t-1],:], u_c)
        cHist[t], nHist[t], Bhist[t] = c[s], n[s], x[s]/u_c[s]
        RHist[t-1] = Uc(cHist[t-1], nHist[t-1])/(beta*Eu_c)
        muHist[t] = mu
    end
    return cHist, nHist, Bhist, TauHist, sHist, muHist, RHist
end

"""
Bellman equation for the continuation of the Lucas-Stokey Problem
"""
mutable struct BellmanEquation{TP <: Para,
                               TI <: Integer,
                               TV <: AbstractVector,
                               TM <: AbstractMatrix{TV},
                               TVV <: AbstractVector{TV}}
    para::TP
    S::TI
    xbar::TV
    time_0::Bool
    z0::TM
    cFB::TV
    nFB::TV
    xFB::TV
    zFB::TVV
end

"""
Initializes the class from the calibration `Para`
"""
function BellmanEquation(para::Para, xgrid::AbstractVector, policies0::Vector)
    S = size(para.Pi, 1) # number of states
    xbar = [minimum(xgrid), maximum(xgrid)]
    time_0 = false
    cf, nf, xprimef = policies0
    z0 = [vcat(cf[s](x), nf[s](x), [xprimef[s, sprime](x) for sprime in 1:S])
                        for x in xgrid, s in 1:S]
    cFB, nFB, IFB, xFB, zFB = find_first_best(para, S, 2)
    return BellmanEquation(para, S, xbar, time_0, z0, cFB, nFB, xFB, zFB)
end

"""
Finds the optimal policies
"""
function get_policies_time1(T::BellmanEquation,
                        i_x::Integer, x::AbstractFloat,
                        s::Integer, Vf::AbstractArray)
    para, S = T.para, T.S
    beta, Theta, G, Pi = para.beta, para.Theta, para.G, para.Pi
    U, Uc, Un = para.U, para.Uc, para.Un

    function objf(z::Vector, grad)
        c, xprime = z[1], z[2:end]
        n=c+G[s]
        Vprime = [Vf[sprime](xprime[sprime]) for sprime in 1:S]
        return -(U(c, n) + beta * dot(Pi[s, :], Vprime))
    end
    function cons(z::Vector, grad)
        c, xprime = z[1], z[2:end]
        n=c+G[s]
        return x - Uc(c, n)*c-Un(c, n)*n - beta*dot(Pi[s, :], xprime)
    end
    lb = vcat(0, T.xbar[1]*ones(S))
    ub = vcat(1-G[s], T.xbar[2]*ones(S))
    opt = Opt(:LN_COBYLA, length(T.z0[i_x, s])-1)
    min_objective!(opt, objf)
    equality_constraint!(opt, cons)
    lower_bounds!(opt, lb)
    upper_bounds!(opt, ub)
    maxeval!(opt, 300)
    maxtime!(opt, 10)
    init = vcat(T.z0[i_x, s][1], T.z0[i_x, s][3:end])
    for (i, val) in enumerate(init)
        if val > ub[i]
            init[i] = ub[i]
        elseif val < lb[i]
            init[i] = lb[i]
        end
    end
    (minf, minx, ret) = optimize(opt, init)
    T.z0[i_x, s] = vcat(minx[1], minx[1]+G[s], minx[2:end])
    return vcat(-minf, T.z0[i_x, s])
end
"""
Finds the optimal policies
"""
function get_policies_time0(T::BellmanEquation,
                        B_::AbstractFloat, s0::Integer, Vf::Array)
    para, S = T.para, T.S
    beta, Theta, G, Pi = para.beta, para.Theta, para.G, para.Pi
    U, Uc, Un = para.U, para.Uc, para.Un
    function objf(z, grad)
        c, xprime = z[1], z[2:end]
        n = c+G[s0]
        Vprime = [Vf[sprime](xprime[sprime]) for sprime in 1:S]
        return -(U(c, n) + beta*dot(Pi[s0, :], Vprime))
    end
    function cons(z::Vector, grad)
        c, xprime = z[1], z[2:end]
        n = c+G[s0]
        return -Uc(c, n)*(c-B_)-Un(c, n)*n - beta*dot(Pi[s0, :], xprime)
    end
    lb = vcat(0, T.xbar[1]*ones(S))
    ub = vcat(1-G[s0], T.xbar[2]*ones(S))
    opt = Opt(:LN_COBYLA, length(T.zFB[s0])-1)
    min_objective!(opt, objf)
    equality_constraint!(opt, cons)
    lower_bounds!(opt, lb)
    upper_bounds!(opt, ub)
    maxeval!(opt, 300)
    maxtime!(opt, 10)
    init = vcat(T.zFB[s0][1], T.zFB[s0][3:end])
    for (i, val) in enumerate(init)
        if val > ub[i]
            init[i] = ub[i]
        elseif val < lb[i]
            init[i] = lb[i]
        end
    end
    (minf, minx, ret) = optimize(opt, init)
    return vcat(-minf, vcat(minx[1], minx[1]+G[s0], minx[2:end]))
end


text:
function length(s::String)
    i = start(s)
    if done(s,i)
        return 0
    end
    n = 1
    while true
        c, j = next(s,i)
        if done(s,j)
            return n
        end
        n += 1
        i = j
    end
end

function isvalid(s::String, i::Integer)
    try
        next(s,i)
        true
    catch
        false
    end
end

max_stars_repo_name: dadarakt/jamon

max_issues_repo_name: dadarakt/jamon

max_forks_repo_name: dadarakt/jamon


text:
function prod_smap!(y::Vector{ComplexF64}, smaps::Matrix{ComplexF64}, x::Vector{ComplexF64}, numVox::Int64, numChan::Int64, numContr::Int=1)
    x = reshape(x,numVox,numContr)
    y = reshape(y,numVox,numContr,numChan)
    @inbounds for j=1:numChan
        sj = (@view smaps[:,j])
        @inbounds for i=1:numContr
            y[:,i,j] .= (@view x[:,i]).* sj
        end
    end
end

function ctprod_smap!(y::Vector{ComplexF64}, smapsC::Matrix{ComplexF64}, x::Vector{ComplexF64}, numVox::Int64, numChan::Int64, numContr::Int=1)
    x = reshape(x,numVox,numContr,numChan)
    y = reshape(y,numVox,numContr)
    y .= 0
    @inbounds for j=1:numChan
        sj = @view smapsC[:,j]
        @inbounds for i=1:numContr
            y[:,i] += (@view x[:,i,j]).* sj
        end
    end
end

function prod_smap(smaps::Matrix{ComplexF64}, x::Vector{ComplexF64}, numVox::Int64, numChan::Int64, numContr::Int=1)
    y = zeros(ComplexF64,numVox*numContr*numChan)
    prod_smap!(y,smaps,x,numVox,numChan,numContr)
    return y
end

function ctprod_smap(smapsC::Matrix{ComplexF64}, x::Vector{ComplexF64}, numVox::Int64, numChan::Int64, numContr::Int=1)
    y = zeros(ComplexF64,numVox*numContr)
    ctprod_smap!(y,smapsC,x,numVox,numChan,numContr)
    return y
end

"""
    SensitivityOp(sensMaps::Matrix{ComplexF64}, numEchoes::Int=1)

builds a `LinearOperator` which performs multiplication of a given image with
the coil sensitivities specified in `sensMaps`

# Arguments
* `sensMaps::Matrix{ComplexF64}`  - sensitivity maps ( 1. dim -> voxels, 2. dim-> coils)
* `numEchoes`                     - number of contrasts to which the opetaor will be applied
"""
function SensitivityOp(sensMaps::Matrix{ComplexF64}, numContr::Int=1)
    numVox, numChan = size(sensMaps)
    sensMapsC = conj.(sensMaps)
    return LinearOperator(numVox*numContr*numChan,numVox*numContr,false,false,
                         x->prod_smap(sensMaps,x,numVox,numChan,numContr),
                         nothing,
                         x->ctprod_smap(sensMapsC,x,numVox,numChan,numContr))
end

"""
    SensitivityOp(sensMaps::Array{T,4}, numContr::Int=1) where T

builds a `LinearOperator` which performs multiplication of a given image with
the coil sensitivities specified in `sensMaps`

# Arguments
* `sensMaps::Array{T,4}`  - sensitivity maps ( 1.-3. dim -> voxels, 4. dim-> coils)
* `numContr`             - number of contrasts to which the opetaor will be applied
"""
function SensitivityOp(sensMaps::Array{T,4}, numContr::Int=1) where T #{T,D}
  sensMaps_mat = reshape(sensMaps, div(length(sensMaps),size(sensMaps,4)),size(sensMaps,4))
  return SensitivityOp(sensMaps_mat, numContr)
end


text:
using Test
using LinearAlgebra
using Random
using StridedTensorXD
using StridedTensorXD: StridedView, UnsafeStridedView

Random.seed!(1234)

println("Base.Threads.nthreads() =  $(Base.Threads.nthreads())")

println("Running tests single-threaded:")
StridedTensorXD.disable_threads()
include("othertests.jl")
include("blasmultests.jl")

println("Running tests multi-threaded:")
StridedTensorXD.enable_threads()
include("othertests.jl")
include("blasmultests.jl")

StridedTensorXD.enable_threaded_mul()
include("blasmultests.jl")


text:
# Building Blocks
include("atomic_information.jl")

# %%codeblock
# default user input specifies 8 aluminum oxide building blocks, density 3.95 g cm^-3
D_INPUT = ["Al4:O6", "8", "3.95"]

# %%codeblock
function main()
    println("Building some blocks")
end

# %%codeblock
main()