text:
export score, AverageEnergy, DifferentialEntropy, BetheFreeEnergy
export @average_energy

abstract type AbstractScoreObjective end

function score end

# Default version is differentiable
# Specialized versions like score(Float64, ...) are not differentiable, but could be faster
score(objective::AbstractScoreObjective, model)                                           = score(objective, model, AsapScheduler())
score(objective::AbstractScoreObjective, model, scheduler)                                = score(InfCountingReal, objective, model, scheduler)
score(::Type{T}, objective::AbstractScoreObjective, model) where { T <: Real }            = score(T, objective, model, AsapScheduler())
score(::Type{T}, objective::AbstractScoreObjective, model, scheduler) where { T <: Real } = score(InfCountingReal{T}, objective, model, scheduler)

# Bethe Free Energy objective

struct AverageEnergy       end
struct DifferentialEntropy end

struct BetheFreeEnergy{S} <: AbstractScoreObjective
    marginal_skip_strategy :: S
end

BetheFreeEnergy() = BetheFreeEnergy(SkipInitial())

marginal_skip_strategy(objective::BetheFreeEnergy) = objective.marginal_skip_strategy

function score(::Type{T}, objective::BetheFreeEnergy, model, scheduler) where { T <: InfCountingReal }

    node_bound_free_energies     = map((node) -> score(T, objective, FactorBoundFreeEnergy(), node, scheduler), getnodes(model))
    variable_bound_entropies     = map((v) -> score(T, objective, VariableBoundEntropy(), v, scheduler), getrandom(model))
    node_bound_free_energies_sum = collectLatest(T, T, node_bound_free_energies, reduce_with_sum)
    variable_bound_entropies_sum = collectLatest(T, T, variable_bound_entropies, reduce_with_sum)

    data_point_entropies_n     = mapreduce(degree, +, getdata(model), init = 0)
    constant_point_entropies_n = mapreduce(degree, +, getconstant(model), init = 0)
    form_point_entropies_n     = count(r -> is_point_mass_form_constraint(form_constraint(r)), getrandom(model))

    point_entropies = InfCountingReal(eltype(T), data_point_entropies_n + constant_point_entropies_n + form_point_entropies_n)

    return combineLatest((node_bound_free_energies_sum, variable_bound_entropies_sum), PushNew()) |> map(eltype(T), d -> float(d[1] + d[2] - point_entropies))
end

## Average energy function helpers

function score(::AverageEnergy, fform, ::Type{ <: Val }, marginals::Tuple{ <: Marginal{ <: NamedTuple{ N } } }, meta) where N
    joint = marginals[1]

    transform = let is_joint_clamped = is_clamped(joint), is_joint_initial = is_initial(joint)
        (data) -> Marginal(data, is_joint_clamped, is_joint_initial)
    end

    return score(AverageEnergy(), fform, Val{ N }, map(transform, values(getdata(joint))), meta)
end

## Differential entropy function helpers

score(::DifferentialEntropy, marginal::Marginal) = entropy(marginal)

function score(::DifferentialEntropy, marginal::Marginal{ <: NamedTuple })

    compute_score = let is_marginal_clamped = is_clamped(marginal), is_marginal_initial = is_initial(marginal)
        (data) -> score(DifferentialEntropy(), Marginal(data, is_marginal_clamped, is_marginal_initial))
    end

    return mapreduce(compute_score, +, values(getdata(marginal)))
end

## Average enery macro helper

import .MacroHelpers

macro average_energy(fformtype, lambda)
    @capture(lambda, (args_ where { whereargs__ } = body_) | (args_ = body_)) ||
        error("Error in macro. Lambda body specification is incorrect")

    @capture(args, (inputs__, meta::metatype_) | (inputs__, )) ||
        error("Error in macro. Lambda body arguments speicifcation is incorrect")

    fuppertype = MacroHelpers.upper_type(fformtype)
    whereargs  = whereargs === nothing ? [] : whereargs
    metatype   = metatype === nothing ? :Any : metatype

    inputs = map(inputs) do input
        @capture(input, iname_::itype_) || error("Error in macro. Input $(input) is incorrect")
        return (iname, itype)
    end

    q_names, q_types, q_init_block = rule_macro_parse_fn_args(inputs; specname = :marginals, prefix = :q_, proxy = :Marginal)

    result = quote
        function ReactiveMP.score(
            ::AverageEnergy,
            fform           :: $(fuppertype),
            marginals_names :: $(q_names),
            marginals       :: $(q_types),
            meta            :: $(metatype)
        ) where { $(whereargs...) }
            $(q_init_block...)
            $(body)
        end
    end

    return esc(result)
end


text:
module Day23
__revise_mode__ = :eval

my_input = 643719258
input_e1 = 389125467

# linkedlist as array: cups[i] is the label of the cup that comes after the cup labeled i

function make_cups(input, to = length(input))
    cups = Vector{Int}(undef, to)

    for i in 2:length(input)
        cups[input[i - 1]] = input[i]
    end
    for i in length(input) + 1:to - 1
        cups[i] = i + 1
    end

    if length(input) == to
        cups[input[end]] = input[1]
    else
        cups[input[end]] = length(input) + 1
        cups[end] = input[1]
    end

    cups
end

function simulate!(cups, curr)
    cup1 = cups[curr]
    cup2 = cups[cup1]
    cup3 = cups[cup2]
    cup4 = cups[cup3]

    dest = mod1(curr - 1, length(cups))
    while dest == cup1 || dest == cup2 || dest == cup3
        dest = mod1(dest - 1, length(cups))
    end
    cups[curr] = cup4
    cups[cup3] = cups[dest]
    cups[dest] = cup1

    cup4
end

function solve_a(input)
    inp = input |> digits |> reverse
    cups = make_cups(inp)
    curr = inp[1]
    for _ in 1:100
        curr = simulate!(cups, curr)
    end

    out = cups[1]
    prev = cups[1]
    while cups[prev] != 1
        out = muladd(out, 10, cups[prev])
        prev = cups[prev]
    end
    out
end

function solve_b(input)
    inp = input |> digits |> reverse
    cups = make_cups(inp, 1000000)
    curr = inp[1]
    for _ in 1:10000000
        curr = simulate!(cups, curr)
    end

    cup1 = cups[1]
    cup1 * cups[cup1]
end

solve(input = my_input) = solve_a(input), solve_b(input)

end

ext: jl

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
using Test
using CVChannel

"""
This script verifies that the entanglement-assisted communication value for
axisymmetric channels is bound as `eaCV(J_N) ≤ d * CV(J_N)`, where `d` is the
Hilbert space dimension of the channel's output. PPT axisymmetric states are also
separable, therefore, we can compute the communication value with `pptCVPrimal`.
"""

"""
    _verify_eaCV_bound_for_axisymmetric_channels(
        d::Int64,
        x_step::Float64,
        y_step::Float64
    )

Scans across the `d`-dimensional axisymmetric channel parameters evaluating both
the entanglement-assisted CV (`eaCVPrimal`) and the CV (`pptCVPrimal`).
The two axisymmetric parameters x and y are scanned across with step-sizes `x_step`
and `y_step` respectively.
For each channel, it is verified that the entanglement-assisted CV is bound by
`eaCV ≦ d * CV`.
"""
function _verify_eaCV_bound_for_axisymmetric_channels(d::Int64, x_step::Float64, y_step::Float64)
    y_bounds = CVChannel._axisymmetric_y_bounds(d)

    it = 0
    y_range = y_bounds[1]:y_step:y_bounds[2]
    for y in y_range
        it = it + 1
        println("(d = $d) verifying y-slice $it of $(length(y_range))")
        x_constraints = CVChannel._axisymmetric_x_constraints(d,y)
        for x in x_constraints[1]:x_step:x_constraints[2]
            @testset "(x,y) = ($x,$y)" begin
                ρ_axi = axisymmetricState(d,x,y)

                J_N = d*ρ_axi
                (ea_cv_N, ea_σAB_N) = eaCVPrimal(J_N, d, d)
                (ppt_cv_N, ppt_σAB_N) = pptCVPrimal(J_N, d, d)

                @test ea_cv_N ≤ d * ppt_cv_N || ea_cv_N ≈ d * ppt_cv_N
            end
        end
    end
end

println("Verifying bounds on entanglement-assisted CV")
@testset "Coarse-grained verification of entanglement-assisted CV bounds (d = $d)" for d in 2:8
    y_step = 0.1
    x_step = 0.2

    _verify_eaCV_bound_for_axisymmetric_channels(d,x_step,y_step)
end

@testset "Fine-grained verification of entanglement-assisted CV bounds (d = $d)" for d in 2:4
    y_step = 0.01
    x_step = 0.01

    _verify_eaCV_bound_for_axisymmetric_channels(d,x_step,y_step)
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
using CodecZlib

const VERSION_NUMBER = 0

"""
    decode_blueprint(io)

Decodes a blueprint string from Factorio.

Takes a base64 encoded string, removes the first byte (version header),
then decompresses it into JSON via zlib. Returns a JSON string.
"""
function decode_blueprint(in::IOBuffer)
    in_str = String(read(in))
    base64_str = base64decode(in_str[2:end])
    zlib_stream = ZlibDecompressorStream(IOBuffer(base64_str))
    return String(read(zlib_stream))
end

decode_blueprint(in::String) = decode_blueprint(IOBuffer(in))
decode_blueprint(in::Array{UInt8, 1}) = decode_blueprint(IOBuffer(in))

"""
    encode_blueprint(io)

Encodes a blueprint string from Factorio.

Takes a JSON string, compresses it via zlib, converts it to base64, then
adds a version header. Returns a base64 string
"""
function encode_blueprint(in::IOBuffer)
    zlib_stream = ZlibCompressorStream(in)
    base64_str = base64encode(String(read(zlib_stream)))
    return base64_str
end

encode_blueprint(in::String) = encode_blueprint(IOBuffer(in))
encode_blueprint(in::Array{UInt8, 1}) = encode_blueprint(IOBuffer(in))


text:
"""
    CSV.write(file::Union{String, IO}, file; kwargs...) => file
    table |> CSV.write(file::Union{String, IO}; kwargs...) => file

Write a [Tables.jl interface input](https://github.com/JuliaData/Tables.jl) to a csv file, given as an `IO` argument or String representing the file name to write to.

Keyword arguments include:
* `delim::Union{Char, String}=','`: a character or string to print out as the file's delimiter
* `quotechar::Char='"'`: character to use for quoting text fields that may contain delimiters or newlines
* `openquotechar::Char`: instead of `quotechar`, use `openquotechar` and `closequotechar` to support different starting and ending quote characters
* `escapechar::Char='\\'`: character used to escape quote characters in a text field
* `missingstring::String=""`: string to print
* `dateformat=Dates.default_format(T)`: the date format string to use for printing out Date & DateTime columns
* `append=false`: whether to append writing to an existing file/IO, if `true`, it will not write column names by default
* `writeheader=!append`: whether to write an initial row of delimited column names, not written by default if appending
* `header`: pass a list of column names (Symbols or Strings) to use instead of the column names of the input table
"""
function write end

_seekstart(p::Base.Process) = return
_seekstart(io::IO) = seekstart(io)

function with(f::Function, io::IO, append)
    !append && _seekstart(io)
    f(io)
end

function with(f::Function, file::String, append)
    open(file, append ? "a" : "w") do io
        f(io)
    end
end

const VALUE_BUFFERS = [IOBuffer()]

function _reset(io::IOBuffer)
    io.ptr = 1
    io.size = 0
end

function bufferedwrite(io, val::Number, df)
    print(io, val)
    return false
end

getvalue(x, df) = x
getvalue(x::T, df) where {T <: Dates.TimeType} = Dates.format(x, df === nothing ? Dates.default_format(T) : df)

function bufferedwrite(io, val, df)
    VALUE_BUFFER = VALUE_BUFFERS[Threads.threadid()]
    _reset(VALUE_BUFFER)
    print(VALUE_BUFFER, getvalue(val, df))
    return true
end

function bufferedescape(io, delim, oq, cq, e)
    VALUE_BUFFER = VALUE_BUFFERS[Threads.threadid()]
    n = position(VALUE_BUFFER)
    n == 0 && return
    needtoescape, needtoquote = check(n, delim, oq, cq)
    seekstart(VALUE_BUFFER)
    if needtoquote
        if needtoescape
            Base.write(io, oq)
            buf = VALUE_BUFFER.data
            for i = 1:n
                @inbounds b = buf[i]
                if b === cq || b === oq
                    Base.write(io, e, b)
                else
                    Base.write(io, b)
                end
            end
            Base.write(io, cq)
        else
            Base.write(io, oq, VALUE_BUFFER, cq)
        end
    else
        Base.write(io, VALUE_BUFFER)
    end
    return
end

function check(n, delim::Char, oq, cq)
    needtoescape = false
    buf = VALUE_BUFFERS[Threads.threadid()].data
    @inbounds needtoquote = buf[1] === oq
    d = delim % UInt8
    @simd for i = 1:n
        @inbounds b = buf[i]
        needtoquote |= (b === d) | (b === UInt8('\n')) | (b === UInt8('\r'))
        needtoescape |= b === cq
    end
    return needtoescape, needtoquote
end

function check(n, delim::String, oq, cq)
    needtoescape = false
    buf = VALUE_BUFFERS[Threads.threadid()].data
    @inbounds needtoquote = buf[1] === oq
    @simd for i = 1:n
        @inbounds b = buf[i]
        needtoquote |= (b === UInt8('\n')) | (b === UInt8('\r'))
        needtoescape |= b === cq
    end
    needtoquote |= occursin(delim, String(buf[1:n]))
    return needtoescape, needtoquote
end

write(file::Union{String, IO}; kwargs...) = x->write(file, x; kwargs...)
function write(file::Union{String, IO}, itr; kwargs...)
    rows = Tables.rows(itr)
    sch = Tables.schema(rows)
    return write(sch, rows, file; kwargs...)
end

function printheader(io, header, delim, oq, cq, e, df)
    cols = length(header)
    for (col, nm) in enumerate(header)
        bufferedwrite(io, string(nm), nothing) && bufferedescape(io, delim, oq, cq, e)
        Base.write(io, ifelse(col == cols, UInt8('\n'), delim))
    end
    return
end

function write(sch::Tables.Schema{schema_names}, rows, file::Union{String, IO};
    delim::Union{Char, String}=',',
    quotechar::Char='"',
    openquotechar::Union{Char, Nothing}=nothing,
    closequotechar::Union{Char, Nothing}=nothing,
    escapechar::Char='\\',
    missingstring::AbstractString="",
    dateformat=nothing,
    append::Bool=false,
    writeheader::Bool=!append,
    header::Vector=String[],
    kwargs...) where {schema_names}
    oq, cq = openquotechar !== nothing ? (openquotechar % UInt8, closequotechar % UInt8) : (quotechar % UInt8, quotechar % UInt8)
    e = escapechar % UInt8
    names = isempty(header) ? schema_names : header
    cols = length(names)
    with(file, append) do io
        writeheader && printheader(io, names, delim, oq, cq, escapechar, dateformat)
        for row in rows
            Tables.eachcolumn(sch, row) do val, col, nm
                bufferedwrite(io, coalesce(val, missingstring), dateformat) && bufferedescape(io, delim, oq, cq, e)
                Base.write(io, ifelse(col == cols, UInt8('\n'), delim))
            end
        end
    end
    return file
end

# handle unknown schema case
function write(::Nothing, rows, file::Union{String, IO};
    delim::Union{Char, String}=',',
    quotechar::Char='"',
    openquotechar::Union{Char, Nothing}=nothing,
    closequotechar::Union{Char, Nothing}=nothing,
    escapechar::Char='\\',
    missingstring::AbstractString="",
    dateformat=nothing,
    append::Bool=false,
    writeheader::Bool=!append,
    header::Vector=String[],
    kwargs...)
    oq, cq = openquotechar !== nothing ? (openquotechar % UInt8, closequotechar % UInt8) : (quotechar % UInt8, quotechar % UInt8)
    e = escapechar % UInt8
    state = iterate(rows)
    if state === nothing
        if writeheader && !isempty(header)
            with(file, append) do io
                printheader(io, header, delim, oq, cq, escapechar, dateformat)
            end
        end
        return file
    end
    row, st = state
    names = isempty(header) ? propertynames(row) : header
    sch = Tables.Schema(names, nothing)
    cols = length(names)
    with(file, append) do io
        writeheader && printheader(io, names, delim, oq, cq, escapechar, dateformat)
        while true
            Tables.eachcolumn(sch, row) do val, col, nm
                bufferedwrite(io, coalesce(val, missingstring), dateformat) && bufferedescape(io, delim, oq, cq, e)
                Base.write(io, ifelse(col == cols, UInt8('\n'), delim))
            end
            state = iterate(rows, st)
            state === nothing && break
            row, st = state
        end
    end
    return file
end


text:
using stoosc
rc("font", family="arial")

function params_create_model!(params)

    parameter_names = params["parameter_names"]  # Ω must be last!
    parameter_values = params["parameter_values"]
    problem_type = params["problem_type"]
    model_name = params["model_name"]

    # Set up an SDE model
    model = load_model(model_name, problem_type)
    model = set_solver(model, saveat=0.01)
    model = set_timespan(model, (0.0, 200.0))
    model = set_parameter_value(model, parameter_names, parameter_values)
    if problem_type == "jump"
        outfun = sol -> hcat(sol[2, :], sol[3, :]) ./ parameter_values[end]
    else
        outfun = sol -> hcat(sol[2, :], sol[3, :])  # save only Y and Z
    end
    model = set_output(model, outfun)

    # Add events
    events = create_events([(:LD, 200 + 1, 0.5, 0.5)])
    model = set_input(model, events)

    params["model"] = model

end

function params_simulate_model!(params)

    model = params["model"]
    trajectories = params["trajectories"]

    solution = simulate_population(model, trajectories, show_progress=true)
    solution = select_time(solution, mintime=100.0)
    t = solution.time
    Y = reshape(solution.trajectories[:, 1, :],
        length(solution.trajectories[:, 1, 1])*trajectories, 1)[:, 1]
    Z = reshape(solution.trajectories[:, 2, :],
        length(solution.trajectories[:, 2, 1])*trajectories, 1)[:, 1]

    params["t"] = t
    params["Y"] = Y
    params["Z"] = Z

end

function params_create_figure(params)

    t = params["t"]
    Y = params["Y"]
    Z = params["Z"]
    plot_range = params["plot_range"]
    filename = params["filename"]
    histogram_bins = params["histogram_bins"]
    phase_plane_title = params["phase_plane_title"]

    fig, ax = subplots(figsize=(4.5, 2.5))
    ax.plot(t, Z[1:length(t)], color="black")
    ax.set_ylim(plot_range[2][1], plot_range[2][2])
    ax.set_xlim(0, 50)
    ax.set_title("Time series")
    ax.set_xlabel("Time")
    ax.set_ylabel("Z")
    fig.tight_layout(pad=0.1)

    if !isnothing(filename)
        time_filename = filename[1:end-4] * "_time" * filename[end-3:end]
        save_figure(fig, "sde_vs_jump", time_filename)
    end

    fig, ax = subplots(figsize=(3, 3))
    h = ax.hist2d(Y, Z, bins=histogram_bins, density=true, range=plot_range,
        rasterized=true)
    if isnothing(phase_plane_title)
        ax.set_title("Phase-plane histogram")
    else
        ax.set_title(phase_plane_title)
    end
    ax.set_xlabel("Y")
    ax.set_ylabel("Z")
    fig.tight_layout(pad=0.1)

    if !isnothing(filename)
        save_figure(fig, "sde_vs_jump", filename)
    end

end

params = Dict{String, Any}(
    "trajectories" => 100,
    "model_name" => "model01",
    "parameter_names" => ["I", "Ω"],
    "parameter_values" => [0, 1/0.01^2],
    "problem_type" => "jump",  # sde or jump
    "histogram_bins" => 25, # 25,
    "plot_range" => [[0.04, 0.16], [0.06, 0.14]],
    "histogram_bins" => 25,
    "filename" => "jump_σ=0.01_I=0.svg",
    "phase_plane_title" => "Gillespie, σ = 0.01, I = 0"
)

params_create_model!(params)
params_simulate_model!(params)
params_create_figure(params)

params["parameter_values"] = [0.3, 1/0.01^2]
params["filename"] = "jump_σ=0.01_I=0.3.svg"
params["phase_plane_title"] = "Gillespie, σ = 0.01, I = 0.3"

params_create_model!(params)
params_simulate_model!(params)
params_create_figure(params)

params["problem_type"] = "sde"
params["parameter_names"] = ["I", "σ"]
params["parameter_values"] = [0, 0.01]
params["filename"] = "sde_σ=0.01_I=0.svg"
params["phase_plane_title"] = "SDE model, σ = 0.01, I = 0"

params_create_model!(params)
params_simulate_model!(params)
params_create_figure(params)

params["parameter_values"] = [0.2, 0.01]
params["filename"] = "sde_σ=0.01_I=0.3.svg"
params["phase_plane_title"] = "SDE model, σ = 0.01, I = 0.3"

params_create_model!(params)
params_simulate_model!(params)
params_create_figure(params)

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
using Knet: Knet, AutoGrad, Data, adam, adam!, progress!, minibatch, save, relu,
    gpu, KnetArray, load, sigm, goldensection, hyperband
using Statistics: mean
using DataFrames
using Plots
using Dates

PARAMS = Dict("batchsize" => 32,
              "pdrop" => 0.1,
              "lr" => 0.01,
              "hidden" => 100,
              "optimizer" => adam!
              )

mutable struct Model
    config::Config;
    name::String
    model
    savepath::String
    datapath::String
    vocabulary
    params
    extractor
    fdict
end

Model(config::Config; name="model", savedir=SAVEDIR, datadir="",
      voc=nothing, params=PARAMS, extractor=Dict(), fdict=Dict()) = Model(
          config, name, nothing, savedir, datadir, voc, params, extractor,
          fdict)

Model(inputs::Array{Tuple{String, String},1},
      outputs::Array{Tuple{String, String},1};
      name="model", savedir=SAVEDIR, datadir="", voc=nothing, params=PARAMS,
      extractor=Dict(), fdict=Dict()) = Model(Config(inputs, outputs), name=name,
                                savedir=SAVEDIR, datadir=datadir,
                                voc=voc, params=params,
                                extractor=extractor, fdict=fdict)

getfeatures(m::Model; ftype="all") = getfeatures(m.config; ftype=ftype)
getfnames(m::Model; ftype="all") = getfnames(m.config; ftype=ftype)
getftypes(m::Model; ftype="all") = getftypes(m.config; ftype=ftype)

function iscategorical(m::Model)
    in(CATEGORY, getftypes(m; ftype="output")) ||
        in(BINARYCATEGORY, getftypes(m; ftype="output"))
end

isimagemodel(m::Model) = in(IMAGE, getftypes(m; ftype="input"))
istextmodel(m::Model) = in(TEXT, getftypes(m; ftype="input"))

function savemodel(m::Model, savepath=nothing)
    savepath = savepath == nothing ? joinpath(SAVEDIR, "$(m.name).jld2") :
        savepath
    save(savepath,
         "config", m.config,
         "name", m.name,
         "model", m.model,
         "savepath", m.savepath,
         "params", m.params,
         "extractor", m.extractor,
         "fdict", m.fdict)
end

function loadmodel(m::Model, loadpath::String)
    m.config = load(loadpath, "config")
    m.name = load(loadpath, "name")
    m.model = load(loadpath, "model")
    m.savepath = load(loadpath, "savepath")
    m.params = load(loadpath, "params")
    m.extractor = load(loadpath, "extractor")
    m.fdict = load(loadpath, "fdict")

    return m
end

function loadmodel(loadpath::String)
    m = Model(Config())
    loadmodel(m, loadpath)
    m
end

function predictdata(m::Model, example)
    data = Dict(fname => [value] for (fname, value) in example)
    data = preprocess(m, data)
    x = preparedata(m, data; output=false)

    iscategorical(m) ? predict(m.model, x) : m.model(x)
end

max_stars_repo_name: alperkesen/AutoML

max_issues_repo_name: alperkesen/AutoML


text:
using kRPC
using CSV
using LinearAlgebra
using DataFrames
using Plots
using Rotations
pyplot()

conn = kRPCConnect("Aerodata")
using kRPC.Remote.SpaceCenter
totable = 1
outp_a = 1
outp_b = 1
vessel = ActiveVessel()
flight = Flight(vessel)
orbit = Orbit(vessel)
body = Body(orbit)
body_frame = ReferenceFrame(body)

function compute_reference_frame(body, lat, lon, alt::Float64)
    landing_position = SurfacePosition(body, lat, lon, body_frame)
    q_long = (
        0.0,
        sin(-lon * 0.5 * pi / 180),
        0.0,
        cos(-lon * 0.5 * pi / 180))
    q_lat = (
        0.0,
        0.0,
        sin(lat * 0.5 * pi / 180),
        cos(lat * 0.5 * pi / 180))
    cr1 = CreateRelative(body_frame, position=landing_position, rotation=q_long)
    cr2 = CreateRelative(cr1,rotation=q_lat)
    return CreateRelative(cr2,position=(alt,0.0,0.0))
end

rf = compute_reference_frame(body, -0.09720792138114814, -74.55767985150328, 0.0)
ref = ReferenceFrame(vessel)
gfl = Flight(vessel, referenceFrame=rf)

trq_vis = kRPC.Remote.Drawing.AddLine((0.0,0.0,0.0),(0.0,0.0,0.0),rf)
strq_vis = kRPC.Remote.Drawing.AddLine((0.0,0.0,0.0),(0.0,0.0,0.0),rf)
sfrc_vis = kRPC.Remote.Drawing.AddLine((0.0,0.0,0.0),(0.0,0.0,0.0),rf)
kRPC.Remote.Drawing.Color!(trq_vis, (0.0,1.0,0.0))
kRPC.Remote.Drawing.Color!(strq_vis, (0.0,0.0,1.0))
global lut = -10000.0
global lav = Float64[]
global iter = 0
add_multiple_streams([
	kRPC.Remote.SpaceCenter.Delayed.AerodynamicTorque(gfl),
	kRPC.Remote.SpaceCenter.Delayed.Velocity(gfl),
	kRPC.Remote.SpaceCenter.Delayed.Position(vessel, rf)],
	function (at, vel, pos)
		if iter % 10 == 0
			sat,saf = kRPC.SendMessage(conn, [
				kRPC.Remote.SpaceCenter.Delayed.SimulateAerodynamicTorqueAt(gfl, body, pos, vel),
				kRPC.Remote.SpaceCenter.Delayed.SimulateAerodynamicForceAt(gfl, body, pos, vel)])
			kRPC.SendMessage(conn,
				[kRPC.Remote.Drawing.Delayed.End!(trq_vis, pos.+at),
			 	 kRPC.Remote.Drawing.Delayed.End!(strq_vis, pos.+sat),
			 	 kRPC.Remote.Drawing.Delayed.End!(sfrc_vis, pos.+saf),
			 	 kRPC.Remote.Drawing.Delayed.Start!(trq_vis, pos),
			 	 kRPC.Remote.Drawing.Delayed.Start!(strq_vis, pos),
			 	 kRPC.Remote.Drawing.Delayed.Start!(sfrc_vis, pos)])
		end
		global iter = iter + 1
	end)

max_stars_repo_name: BenChung/SuccessiveConvexification

max_issues_repo_name: BenChung/SuccessiveConvexification

max_forks_repo_name: BenChung/SuccessiveConvexification


text:
using Test

include("exercism-matrix.jl")

@testset "construct E" begin
    @test E == [
        0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
        0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0;
        0 1 0 0 1 0 1 0 0 0 0 1 0 1 0 0 1 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
        1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
        0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0;
        0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0;
        0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
        0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0;
    ]
    @test isconst(@__MODULE__, :E)
end

@testset "Sadxercism" begin
    @test frown(E) == [
        0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
        0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0;
        0 1 0 0 1 0 1 0 0 0 0 1 0 1 0 0 1 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
        1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
        0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0;
        0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0;
        0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
        0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0;
    ]

    Ẽ = copy(E)
    frown!(Ẽ)
    @test Ẽ == [
        0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
        0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0;
        0 1 0 0 1 0 1 0 0 0 0 1 0 1 0 0 1 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
        1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
        0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0;
        0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0;
        0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
        0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0;
    ]

    @test Ẽ != E
end

@testset "rot -90° / +270°" begin
    @test rot270(E) == [
        0 0 0 0 0 1 0 0 0 0 0;
        0 1 1 1 1 0 1 1 1 1 0;
        1 0 0 0 0 0 0 0 0 0 1;
        1 0 0 0 0 0 0 0 0 0 1;
        0 0 0 1 0 0 0 0 0 0 0;
        0 0 1 0 0 0 0 0 0 0 0;
        0 0 0 1 0 0 1 0 0 0 0;
        0 0 0 0 0 0 0 1 0 0 0;
        0 0 0 0 0 0 0 0 1 0 0;
        0 0 0 0 0 0 0 0 1 0 0;
        0 0 0 0 0 0 0 1 0 0 0;
        0 0 0 1 0 0 1 0 0 0 0;
        0 0 1 0 0 0 0 0 0 0 0;
        0 0 0 1 0 0 0 0 0 0 0;
        1 0 0 0 0 0 0 0 0 0 1;
        1 0 0 0 0 0 0 0 0 0 1;
        0 1 1 1 1 0 1 1 1 1 0;
        0 0 0 0 0 1 0 0 0 0 0;
    ]

    @test rot270(frown(E)) == [
        0 0 0 0 0 1 0 0 0 0 0;
        0 1 1 1 1 0 1 1 1 1 0;
        1 0 0 0 0 0 0 0 0 0 1;
        1 0 0 0 0 0 0 0 0 0 1;
        0 0 0 1 0 0 0 0 0 0 0;
        0 0 1 0 0 0 0 0 0 0 0;
        0 0 0 1 0 0 0 0 1 0 0;
        0 0 0 0 0 0 0 1 0 0 0;
        0 0 0 0 0 0 1 0 0 0 0;
        0 0 0 0 0 0 1 0 0 0 0;
        0 0 0 0 0 0 0 1 0 0 0;
        0 0 0 1 0 0 0 0 1 0 0;
        0 0 1 0 0 0 0 0 0 0 0;
        0 0 0 1 0 0 0 0 0 0 0;
        1 0 0 0 0 0 0 0 0 0 1;
        1 0 0 0 0 0 0 0 0 0 1;
        0 1 1 1 1 0 1 1 1 1 0;
        0 0 0 0 0 1 0 0 0 0 0;
    ]
end

@testset "rot +90°" begin
    @test rot90(E) == [
        0 0 0 0 0 1 0 0 0 0 0;
        0 1 1 1 1 0 1 1 1 1 0;
        1 0 0 0 0 0 0 0 0 0 1;
        1 0 0 0 0 0 0 0 0 0 1;
        0 0 0 0 0 0 0 1 0 0 0;
        0 0 0 0 0 0 0 0 1 0 0;
        0 0 0 0 1 0 0 1 0 0 0;
        0 0 0 1 0 0 0 0 0 0 0;
        0 0 1 0 0 0 0 0 0 0 0;
        0 0 1 0 0 0 0 0 0 0 0;
        0 0 0 1 0 0 0 0 0 0 0;
        0 0 0 0 1 0 0 1 0 0 0;
        0 0 0 0 0 0 0 0 1 0 0;
        0 0 0 0 0 0 0 1 0 0 0;
        1 0 0 0 0 0 0 0 0 0 1;
        1 0 0 0 0 0 0 0 0 0 1;
        0 1 1 1 1 0 1 1 1 1 0;
        0 0 0 0 0 1 0 0 0 0 0;
    ]

    @test rot90(frown(E)) == [
        0 0 0 0 0 1 0 0 0 0 0;
        0 1 1 1 1 0 1 1 1 1 0;
        1 0 0 0 0 0 0 0 0 0 1;
        1 0 0 0 0 0 0 0 0 0 1;
        0 0 0 0 0 0 0 1 0 0 0;
        0 0 0 0 0 0 0 0 1 0 0;
        0 0 1 0 0 0 0 1 0 0 0;
        0 0 0 1 0 0 0 0 0 0 0;
        0 0 0 0 1 0 0 0 0 0 0;
        0 0 0 0 1 0 0 0 0 0 0;
        0 0 0 1 0 0 0 0 0 0 0;
        0 0 1 0 0 0 0 1 0 0 0;
        0 0 0 0 0 0 0 0 1 0 0;
        0 0 0 0 0 0 0 1 0 0 0;
        1 0 0 0 0 0 0 0 0 0 1;
        1 0 0 0 0 0 0 0 0 0 1;
        0 1 1 1 1 0 1 1 1 1 0;
        0 0 0 0 0 1 0 0 0 0 0;
    ]
end

@testset "stickerwall" begin
    @test stickerwall(E) == [
        0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
        0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0;
        0 1 0 0 1 0 1 0 0 0 0 1 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 1 0 0 0 0 1 0 1 0 0 1 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
        1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
        0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0;
        0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0;
        0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
        0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
        1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
        0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0;
        0 1 0 0 1 0 1 0 0 0 0 1 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 1 0 0 0 0 1 0 1 0 0 1 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
        1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
        0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0;
        0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0;
        0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
        0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0;
    ]
end


text:
function post(𝒜::GLGM06,
              𝑃::InitialValueProblem{<:AbstractContinuousSystem},
              𝑂::Options)

    # ==================================
    # Initialization and discretization
    # ==================================

    𝑂 = merge(𝒜.options.defaults, 𝑂, 𝒜.options.specified)
    max_order = 𝑂[:max_order]
    δ, T = 𝑂[:δ], 𝑂[:T]
    N = round(Int, T / δ)

    # compute and unrwap discretized system
    𝑃_discrete = discretize(𝑃, δ, algorithm=𝑂[:discretization],
                                  sih_method=𝑂[:sih_method],
                                  exp_method=𝑂[:exp_method],
                                  set_operations="zonotope")
    Ω0, Φ = 𝑃_discrete.x0, 𝑃_discrete.s.A
    Ω0red = reduce_order(Ω0, max_order)

    # =====================
    # Flowpipe computation
    # =====================

    # preallocate output
    Rsets = Vector{ReachSet{<:Zonotope{Float64}}}(undef, N)

    info("Reachable States Computation...")
    @timing begin
    if inputdim(𝑃_discrete) == 0
        reach_homog!(Rsets, Ω0red, Φ, N, δ)

    else
        U = inputset(𝑃_discrete)
        reach_inhomog!(Rsets, Ω0red, U, Φ, N, δ, max_order)
    end
    end # timing

    Rsol = ReachSolution(Flowpipe(Rsets), 𝑂)

    # ===========
    # Projection
    # ===========

    if 𝑂[:project_reachset]
        info("Projection...")
        Rsol = @timing project(Rsol)
    end

    return Rsol
end


text:
## ================================================================================================
## FlatVariables - used for packing variables for optimization
## ================================================================================================

struct FlatVariables{T<:Real}
  X::Vector{T}
  idx::Dict{Symbol, UnitRange{Int}}
end

function FlatVariables(fg::AbstractDFG, varIds::Vector{Symbol})

  index = 1
  idx = Dict{Symbol, UnitRange{Int}}()
  for vid = varIds
    v = getVariable(fg, vid)
    dims = getDimension(v)
    idx[vid] = index:(index+dims-1)
    index += dims
  end
  return FlatVariables(Vector{Float64}(undef, index-1), idx)
end

function Base.setindex!(flatVar::FlatVariables{T}, val::Vector{T}, vId::Symbol) where T<:Real
  if length(val) == length(flatVar.idx[vId])
    flatVar.X[flatVar.idx[vId]] .= val
  else
    error("array could not be broadcast to match destination")
  end
end

function Base.getindex(flatVar::FlatVariables{T}, vId::Symbol) where T<:Real
  return flatVar.X[flatVar.idx[vId]]
end

## ================================================================================================
## Parametric Factors
## ================================================================================================

"""
    $SIGNATURES

Returns the parametric measurement for a factor as a tuple (measurement, inverse covariance) for parametric inference (assumign Gaussian).
Defaults to find the parametric measurement at field `Z` followed by `z`.

Notes
- Users should overload this method should their factor not default to `.Z<:ParametricType` or `.z<:ParametricType`
"""
function getParametricMeasurement end

function getParametricMeasurement(Z)
  error("$Z not supported, please use non-parametric or open an issue if it should be")
end

function getParametricMeasurement(Z::Normal)
  meas = mean(Z)
  iσ = 1/std(Z)^2
  return [meas], reshape([iσ],1,1)
end

function getParametricMeasurement(Z::MvNormal)
  meas = mean(Z)
  iΣ = invcov(Z)
  return meas, iΣ
end

function getParametricMeasurement(s::FunctorInferenceType)
  if hasfield(typeof(s), :Z)
    Z = s.Z
    @info "getParametricMeasurement falls back to using field `.Z` by default. Extend it for more complex factors." maxlog=1
  elseif hasfield(typeof(s), :z)
    Z = s.z
    @info "getParametricMeasurement falls back to using field `.z` by default. Extend it for more complex factors." maxlog=1
  else
    error("$s not supported, please use non-parametric or open an issue if it should be")
  end
  return getParametricMeasurement(Z)
end

## ================================================================================================
## Parametric solve with Mahalanobis distance - CalcFactor
## ================================================================================================

"""
    $TYPEDEF

Internal parametric extension to [`CalcFactor`](@ref) used for buffering measurement and calculating Mahalanobis distance
"""
struct CalcFactorMahalanobis{CF}
  calcfactor!::CF
  varOrder::Vector{Symbol}
  meas::Vector{Float64}
  iΣ::Matrix{Float64}
end

function CalcFactorMahalanobis(fct::DFGFactor)
  cf = getFactorType(fct)
  varOrder = getVariableOrder(fct)
  meas, iΣ = getParametricMeasurement(cf)
  calcf = CalcFactor(cf, _getFMdThread(fct), 0, 0, (), [])
  return CalcFactorMahalanobis(calcf, varOrder, meas, iΣ)
end

# This is where the actual parametric calculation happens, CalcFactor equivalent for parametric
function (cfp::CalcFactorMahalanobis)(variables...)
  # call the user function (be careful to call the new CalcFactor version only!!!)
  res = cfp.calcfactor!(cfp.meas, variables...)
  # 1/2*log(1/(  sqrt(det(Σ)*(2pi)^k) ))  ## k = dim(μ)
  return res' * cfp.iΣ * res
end

function calcFactorMahalanobisDict(fg)
  calcFactors = Dict{Symbol, CalcFactorMahalanobis}()
  for fct in getFactors(fg)
    calcFactors[fct.label] = CalcFactorMahalanobis(fct)
  end
  return calcFactors
end

function _totalCost(cfdict::Dict{Symbol, CalcFactorMahalanobis},
                    flatvar,
                    X )
  #
  obj = 0
  for (fid, cfp) in cfdict

    varOrder = cfp.varOrder

    Xparams = [view(X, flatvar.idx[varId]) for varId in varOrder]

    # call the user function
    retval = cfp(Xparams...)

    # 1/2*log(1/(  sqrt(det(Σ)*(2pi)^k) ))  ## k = dim(μ)
    obj += 1/2*retval
  end

  return obj
end

# build the cost function
function _totalCost(fg::AbstractDFG,
                    flatvar,
                    X )
  #
  obj = 0
  for fct in getFactors(fg)

    cf = getFactorType(fct)
    varOrder = getVariableOrder(fct)

    Xparams = [view(X, flatvar.idx[varId]) for varId in varOrder]

    retval = cf(Xparams...)

    # 1/2*log(1/(  sqrt(det(Σ)*(2pi)^k) ))  ## k = dim(μ)
    obj += 1/2*retval
  end

  return obj
end

export solveFactorGraphParametric
"""
    $SIGNATURES

Batch solve a Gaussian factor graph using Optim.jl. Parameters can be passed directly to optim.
Notes:
  - Only :Euclid and :Circular manifolds are currently supported, own manifold are supported with `algorithmkwargs` (code may need updating though)
"""
function solveFactorGraphParametric(fg::AbstractDFG;
                                    useCalcFactor::Bool=true, #TODO dev param will be removed
                                    solvekey::Symbol=:parametric,
                                    autodiff = :forward,
                                    algorithm=Optim.Optim.BFGS,
                                    algorithmkwargs=(), # add manifold to overwrite computed one
                                    options = Optim.Options(allow_f_increases=true,
                                                            time_limit = 100,
                                                            # show_trace = true,
                                                            # show_every = 1,
                                                            ))

  #Other options
  # options = Optim.Options(time_limit = 100,
  #                     iterations = 1000,
  #                     show_trace = true,
  #                     show_every = 1,
  #                     allow_f_increases=true,
  #                     g_tol = 1e-6,
  #                     )
  varIds = listVariables(fg)

  #TODO mabye remove sorting, just for convenience
  sort!(varIds, lt=natural_lt)

  flatvar = FlatVariables(fg, varIds)

  for vId in varIds
    flatvar[vId] = getVariableSolverData(fg, vId, solvekey).val[:,1]
  end

  initValues = flatvar.X

  mc_mani = IIF.MixedCircular(fg, varIds)
  alg = algorithm(;manifold=mc_mani, algorithmkwargs...)
  # alg = algorithm(; algorithmkwargs...)

  if useCalcFactor
    cfd = calcFactorMahalanobisDict(fg)
    tdtotalCost = Optim.TwiceDifferentiable((x)->_totalCost(cfd, flatvar, x), initValues, autodiff = autodiff)
  else
    tdtotalCost = Optim.TwiceDifferentiable((x)->_totalCost(fg, flatvar, x), initValues, autodiff = autodiff)
  end

  result = Optim.optimize(tdtotalCost, initValues, alg, options)
  rv = Optim.minimizer(result)

  H = Optim.hessian!(tdtotalCost, rv)

  Σ = pinv(H)

  d = Dict{Symbol,NamedTuple{(:val, :cov),Tuple{Vector{Float64},Matrix{Float64}}}}()

  for key in varIds
    r = flatvar.idx[key]
    push!(d,key=>(val=rv[r],cov=Σ[r,r]))
  end

  return d, result, flatvar.idx, Σ
end

#TODO maybe consolidate with solveFactorGraphParametric
#TODO WIP
```
    $SIGNATURES
Solve for frontal values only with values in seprarators fixed
```
function solveConditionalsParametric(fg::AbstractDFG,
                                    frontals::Vector{Symbol};
                                    solvekey::Symbol=:parametric,
                                    autodiff = :forward,
                                    algorithm=Optim.BFGS,
                                    algorithmkwargs=(), # add manifold to overwrite computed one
                                    options = Optim.Options(allow_f_increases=true,
                                                            time_limit = 100,
                                                            # show_trace = true,
                                                            # show_every = 1,
                                                            ))

  varIds = listVariables(fg)

  #TODO mabye remove sorting, just for convenience
  sort!(varIds, lt=natural_lt)
  separators = setdiff(varIds, frontals)

  varIds = [frontals; separators]

  flatvar = FlatVariables(fg, varIds)

  for vId in varIds
    flatvar[vId] = getVariableSolverData(fg, vId, solvekey).val[:,1]
  end
  initValues = flatvar.X

  frontalsLength = sum(map(v->getDimension(getVariable(fg, v)), frontals))

  # build variables for frontals and seperators
  # fX = view(initValues, 1:frontalsLength)
  fX = initValues[1:frontalsLength]
  # sX = view(initValues, (frontalsLength+1):length(initValues))
  sX = initValues[frontalsLength+1:end]

  mc_mani = MixedCircular(fg, varIds)
  alg = algorithm(;manifold=mc_mani, algorithmkwargs...)

  tdtotalCost = Optim.TwiceDifferentiable((x)->_totalCost(fg, flatvar,x), initValues, autodiff = autodiff)

  result = Optim.optimize((x)->_totalCost(fg, flatvar, [x;sX]), fX, alg, options)
  # result = optimize(x->totalCost([x;sX]), fX, alg, options)

  rv = Optim.minimizer(result)

  H = Optim.hessian!(tdtotalCost, [rv; sX])

  Σ = pinv(H)

  d = Dict{Symbol,NamedTuple{(:val, :cov),Tuple{Vector{Float64},Matrix{Float64}}}}()

  for key in frontals
    r = flatvar.idx[key]
    push!(d,key=>(val=rv[r],cov=Σ[r,r]))
  end

  return d, result, flatvar.idx, Σ
end

## ================================================================================================
## MixedCircular Manifold for Optim.jl
## ================================================================================================

"""
    MixedCircular
Mixed Circular Manifold. Simple manifold for circular and cartesian mixed for use in optim
"""
struct MixedCircular <: Optim.Manifold
  isCircular::BitArray
end

function MixedCircular(fg::AbstractDFG, varIds::Vector{Symbol})
  circMask = Bool[]
  for k = varIds
    append!(circMask, getVariableType(fg, k) |> getManifolds .== :Circular)
  end
  MixedCircular(circMask)
end

function Optim.retract!(c::MixedCircular, x)
  for (i,v) = enumerate(x)
    c.isCircular[i] && (x[i] = rem2pi(v, RoundNearest))
  end
  return x
end
Optim.project_tangent!(S::MixedCircular,g,x) = g

## ================================================================================================
## UNDER DEVELOPMENT Parametric solveTree utils
## ================================================================================================

"""
    $SIGNATURES
Get the indexes for labels in FlatVariables
"""
function collectIdx(varinds, labels)
  idx = Int[]
  for lbl in labels
    append!(idx, varinds[lbl])
  end
  return idx
end

"""
    $SIGNATURES
Calculate the marginal distribution for a clique over subsetVarIds.
"""
function calculateMarginalCliqueLikelihood(vardict, Σ, varindxs, subsetVarIds)

  μₘ = Float64[]
  for lbl in subsetVarIds
    append!(μₘ, vardict[lbl].val)
  end

  Aidx = collectIdx(varindxs, subsetVarIds)
  Σₘ = Σ[Aidx, Aidx]

  return createMvNormal(μₘ, Σₘ)
end

"""
    $SIGNATURES

"""
function calculateCoBeliefMessage(soldict, Σ, flatvars, separators, frontals)
  Aidx = IIF.collectIdx(flatvars,separators)
  Cidx = IIF.collectIdx(flatvars,frontals)

  #marginalize separators
  A = Σ[Aidx, Aidx]
  #marginalize frontals
  C = Σ[Cidx, Cidx]
  # cross
  B = Σ[Aidx, Cidx]

  Σₘ = deepcopy(A)
  if length(separators) == 0

    return (varlbl=Symbol[], μ=Float64[], Σ=Matrix{Float64}(undef,0,0))

  elseif length(separators) == 1

    # create messages
    return (varlbl = deepcopy(separators), μ = soldict[separators[1]].val, Σ = A)

  elseif length(separators) == 2
    A = Σₘ[1, 1]
    C = Σₘ[2, 2]
    B = Σₘ[1, 2]

    #calculate covariance between separators
    ΣA_B = A - B*inv(C)*B'
    # create messages
    m2lbl = deepcopy(separators)
    m2cov = isa(ΣA_B, Matrix) ? ΣA_B : fill(ΣA_B,1,1)
    m2val = soldict[m2lbl[2]].val - soldict[m2lbl[1]].val
    return (varlbl = m2lbl, μ = m2val, Σ = m2cov)

  else
    error("Messages with more than 2 seperators are not supported yet")
  end
end

## ================================================================================================
## Parametric utils
## ================================================================================================

## SANDBOX of usefull development functions to be cleaned up

"""
    $SIGNATURES
Add parametric solver to fg, batch solve using [`solveFactorGraphParametric`](@ref) and update fg.
"""
function solveFactorGraphParametric!(fg::AbstractDFG; init::Bool=true, kwargs...)
  if !(:parametric in fg.solverParams.algorithms)
    addParametricSolver!(fg; init=init)
  elseif init
    initParametricFrom!(fg)
  end

  vardict, result, varIds, Σ = solveFactorGraphParametric(fg; kwargs...)

  updateParametricSolution!(fg, vardict)

  return vardict, result, varIds, Σ
end

"""
    $SIGNATURES
Initialize the parametric solver data from a different solution in `fromkey`.
"""
function initParametricFrom!(fg::AbstractDFG, fromkey::Symbol = :default; parkey::Symbol = :parametric)
  for var in getVariables(fg)
      fromvnd = getSolverData(var, fromkey)
      if fromvnd.dims == 1
        nf = fit(Normal, fromvnd.val)
        getSolverData(var, parkey).val[1,1] = nf.μ
        getSolverData(var, parkey).bw[1,1] = nf.σ
        # m = var.estimateDict[:default].mean
      else
        #FIXME circular will not work correctly with MvNormal
        nf = fit(MvNormal, fromvnd.val)
        getSolverData(var, parkey).val[1:fromvnd.dims] .= mean(nf)[1:fromvnd.dims]
        getSolverData(var, parkey).bw = cov(nf)
      end
  end
end

"""
    $SIGNATURES
Add the parametric solveKey to all the variables in fg if it doesn't exists.
"""
function addParametricSolver!(fg; init=true)
  if !(:parametric in fg.solverParams.algorithms)
      push!(fg.solverParams.algorithms, :parametric)
      foreach(v->IIF.setDefaultNodeDataParametric!(v, getVariableType(v), initialized=false), getVariables(fg))
      if init
        initParametricFrom!(fg)
      end
  else
      error("parametric solvekey already exists")
  end
  nothing
end

function updateVariablesFromParametricSolution!(fg::AbstractDFG, vardict)
  for (v,val) in vardict
    vnd = getVariableSolverData(fg, v, :parametric)
    vnd.val .= val.val
    if size(vnd.bw) != size(val.cov)
      vnd.bw = val.cov
    else
      vnd.bw .= val.cov
    end
  end
end

"""
    $SIGNATURES
Update the fg from solution in vardict and add MeanMaxPPE (all just mean). Usefull for plotting
"""
function updateParametricSolution!(sfg, vardict)

  for (v,val) in vardict
      vnd = getSolverData(getVariable(sfg, v), :parametric)
      # fill in the variable node data value
      vnd.val .= val.val
      #calculate and fill in covariance
      vnd.bw = val.cov
      #fill in ppe as mean
      ppe = MeanMaxPPE(:parametric, val.val, val.val, val.val)
      getPPEDict(getVariable(sfg, v))[:parametric] = ppe
  end
end

function createMvNormal(val,cov)
    #TODO do something better for properly formed covariance, but for now just a hack...FIXME
    if all(diag(cov) .> 0.001) && isapprox(cov, transpose(cov), rtol=1e-4)
        return MvNormal(val,Symmetric(cov))
    else
        @error("Covariance matrix error", cov)
        # return nothing # FIXME, blanking nothing during #459 consolidation
        return MvNormal(val, ones(length(val)))
    end
end

function createMvNormal(v::DFGVariable, key=:parametric)
    if key == :parametric
        vnd = getSolverData(v, :parametric)
        dims = vnd.dims
        return createMvNormal(vnd.val[1:dims,1], vnd.bw[1:dims, 1:dims])
    else
        @warn "Trying MvNormal Fit, replace with PPE fits in future"
        return fit(MvNormal,getSolverData(v, key).val)
    end
end


text:
using LinearAlgebra

function solve()
    n = parse(Int, readline())
    a = [parse(Int, x) for x in split(readline())]

    sort!(a, rev=true)

    c = zeros(Int, n)
    for i=1:n
        c[i] = n + 1 - 2i
    end

    a⋅c
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
# RELOG: Reverse Logistics Optimization
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using DataStructures
using JSON
using JSONSchema
using Printf
using Statistics

function validate(json, schema)
    result = JSONSchema.validate(json, schema)
    if result !== nothing
        if result isa JSONSchema.SingleIssue
            msg = "$(result.reason) in $(result.path)"
        else
            msg = convert(String, result)
        end
        throw("Error parsing input file: $(msg)")
    end
end

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


text:
using HDF5

Base.@kwdef mutable struct CalInfo <: Info
    refant::Union{AbstractString,Nothing} = nothing
    cal_K::AbstractArray{Float32,2} = Array{Float32,2}(undef, 0, 0)
    cal_B::AbstractArray{ComplexF32,3} = Array{ComplexF32,3}(undef, 0, 0, 0)
    cal_G::AbstractArray{ComplexF32,2} = Array{ComplexF32,2}(undef, 0, 0)
    cal_all::AbstractArray{ComplexF32,3} = Array{ComplexF32,3}(undef, 0, 0, 0)
end # mutable struct CalInfo

function to_hdf5(group::HDF5.Group, info::CalInfo; chunks=size(info.cal_B))
    invoke(to_hdf5, Tuple{HDF5.Group, Info}, group, info;
           cal_K=(chunks=chunks[1:2],),
           cal_B=(chunks=chunks,),
           cal_G=(chunks=chunks[1:2],),
           cal_all=(chunks=chunks,)
          )
end


text:
using Test
using Vcov: completecases, materialize, simple, robust, cluster, names, nclusters
using GroupedArrays
const tests = [
    "estimators"
]

printstyled("Running tests:\n", color=:blue, bold=true)

for test in tests
    include("$test.jl")
    println("\033[1m\033[32mPASSED\033[0m: $(test)")
end


text:
f(z::Int) = begin z = 0 end
f(z::Array{Int}) = begin z[1] = 0 end

x = 1
@show typeof(x)
@show isimmutable(x)
println("Before call by value: ", x)
f(x)
println("After call by value: ", x,"\n")

x = [1]
@show typeof(x)
@show isimmutable(x)
println("Before call by reference: ", x)
f(x)
println("After call by reference: ", x)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
@doc raw"""

"""
module MathewsLakshmananOscillator

    export hamiltonian

end


text:
function RoadTrip(N, c, sol)
cities_visited = 1
city = 1
println(c[city])
while cities_visited ≠ N
    city = findfirst(sol[city,:] .== 1.0)
    cities_visited += 1
    println(c[city])
end
end


text:
# set up

m = matread(quad_data_file_name)
qnwfuncs = [qnwlege, qnwcheb, qnwsimp, qnwtrap, qnwbeta, qnwgamma,
            qnwequi, qnwnorm, qnwunif, qnwlogn]

a = -2.0
b = 3.0
n = 11

# 3-d parameters -- just some random numbers
a_3 = [-1.0, -2.0, 1.0]
b_3 = [1.0, 12.0, 1.5]
n_3 = [7, 5, 9]

mu_3d = [1.0, 2.0, 2.5]
sigma2_3d = [1.0 0.1 0.0; 0.1 1.0 0.0; 0.0 0.0 1.2]

# 1-d nodes and weights
x_cheb_1,  w_cheb_1 = qnwcheb(n, a, b)
x_equiN_1, w_equiN_1 = qnwequi(n, a, b, "N")
x_equiW_1, w_equiW_1 = qnwequi(n, a, b, "W")
x_equiH_1, w_equiH_1 = qnwequi(n, a, b, "H")
x_lege_1, w_lege_1 = qnwlege(n, a, b)
x_norm_1, w_norm_1 = qnwnorm(n, a, b)
x_logn_1, w_logn_1 = qnwlogn(n, a, b)
x_simp_1, w_simp_1 = qnwsimp(n, a, b)
x_trap_1, w_trap_1 = qnwtrap(n, a, b)
x_unif_1, w_unif_1 = qnwunif(n, a, b)
x_beta_1, w_beta_1 = qnwbeta(n, b, b+1)
x_gamm_1, w_gamm_1 = qnwgamma(n, b)

@testset "Testing quad.jl" begin

    @testset "Testing method resolution" begin

        for f in qnwfuncs
            m1 = f(11, 1, 3)
            m2 = f([11], 1, 3)
            m3 = f(11, [1], 3)
            m4 = f(11, 1, [3])
            m5 = f([11], [1], 3)
            m6 = f([11], 1, [3])
            m7 = f([11], [1], [3])

            # Stack nodes/weights in columns
            @test isapprox([m1[1] m1[2]], [m2[1] m2[2]])
            @test isapprox([m1[1] m1[2]], [m3[1] m3[2]])
            @test isapprox([m1[1] m1[2]], [m4[1] m4[2]])
            @test isapprox([m1[1] m1[2]], [m5[1] m5[2]])
            @test isapprox([m1[1] m1[2]], [m6[1] m6[2]])
            @test isapprox([m1[1] m1[2]], [m7[1] m7[2]])
        end
    end

    @testset "testing nodes/weights against Matlab" begin
        # generate names
        for name in ["cheb", "equiN", "equiW", "equiH", "lege", "norm",
                     "logn", "simp", "trap", "unif", "beta", "gamm"]
            x_str_name = "x_$(name)_1"
            w_str_name = "w_$(name)_1"
            jl_x, jl_w = eval(symbol(x_str_name)), eval(symbol(w_str_name))
            ml_x, ml_w = m[x_str_name],  m[w_str_name]
            ml_x = squeeze(ml_x, 2)
            @test isapprox(jl_x, ml_x; atol=1e-5)
            @test isapprox(jl_w, squeeze(ml_w, 2); atol=1e-5)
        end
    end

    @testset "testing quadrect 1d against Matlab" begin
        f1(x) = exp(-x)
        f2(x) = 1.0 ./ (1.0 + 25.0 .* x.^2.0)
        f3(x) = abs(x).^0.5

        # dim 1: num nodes, dim2: method, dim3:func
        data1d = Array(Float64, 6, 6, 3)
        kinds = ["trap", "simp", "lege", "N", "W", "H"]
        n_nodes = [5, 11, 21, 51, 101, 401]  # number of nodes
        a, b = -1, 1

        for (k_i, k) in enumerate(kinds)
            for (n_i, n) in enumerate(n_nodes)
                num_in = length(k) == 1 ? n^2 : n
                for (f_i, f) in enumerate([f1, f2, f3])
                    data1d[n_i, k_i, f_i] = quadrect(f, num_in, a, b, k)
                end
            end
        end

        # NOTE: drop last column -- corresponds to "R" and we have different
        #       random numbers than Matlab.
        ml_data_1d = m["int_1d"][:, 1:6, :]

        @test isapprox(data1d[:, 1, :], ml_data_1d[:, 1, :])  # trap
        @test isapprox(data1d[:, 2, :], ml_data_1d[:, 2, :])  # simp
        @test isapprox(data1d[:, 3, :], ml_data_1d[:, 3, :])  # lege
        @test isapprox(data1d[:, 4, :], ml_data_1d[:, 4, :])  # N
        @test isapprox(data1d[:, 5, :], ml_data_1d[:, 5, :])  # W
        @test isapprox(data1d[:, 6, :], ml_data_1d[:, 6, :])  # H
    end

    @testset "testing quadrect 2d against Matlab" begin
        f1(x) = exp(x[:, 1] + x[:, 2])
        f2(x) = exp(- x[:, 1] .* cos(x[:, 2].^2))

        a = ([0.0, 0.0], [-1.0, -1.0])
        b = ([1.0, 2.0], [1.0, 1.0])

        # dim 1: num nodes, dim2: method
        data2d1 = Array(Float64, 6, 6)
        kinds = ["lege", "trap", "simp", "N", "W", "H"]
        n_nodes = [5, 11, 21, 51, 101, 401]  # number of nodes

        for (k_i, k) in enumerate(kinds)
            for (n_i, n) in enumerate(n_nodes)
                num_in = length(k) == 1 ? n^2 : n
                data2d1[n_i, k_i] = quadrect(f1, num_in, a[1], b[1], k)
            end
        end

        # NOTE: drop last column -- corresponds to "R" and we have different
        #       random numbers than Matlab.
        ml_data_2d1 = m["int_2d1"][:, 1:6]

        @test isapprox(data2d1[:, 1],ml_data_2d1[:, 1])  # trap
        @test isapprox(data2d1[:, 2], ml_data_2d1[:, 2])  # simp
        @test isapprox(data2d1[:, 3], ml_data_2d1[:, 3])  # lege
        @test isapprox(data2d1[:, 4], ml_data_2d1[:, 4])  # N
        @test isapprox(data2d1[:, 5], ml_data_2d1[:, 5])  # W
        @test isapprox(data2d1[:, 6], ml_data_2d1[:, 6])  # H
    end

end  # @testset


text:
# sort by column
function sortbyc(x,col)
    x = sortrows(x,by=x->x[col])
end


text:
function update_Xβ_dist!(rr, ϕβ)
    ϕ = ϕβ[1:rr.M]
    β = ϕβ[rr.M+1:end]
    A_mul_B!(rr.Xβ, rr.X, β)
    rr.Θ .= exp.(-rr.Xβ)
    rr.dist.params[:] = ϕ
    return
end

function compute_loglik!(rr)
    for i in rr.I1
        rr.τs[i,1] = rr.W[i,1]*rr.Θ[i]
        rr.loglik[i] = logpdf(rr.dist, rr.τs[i,1])-rr.Xβ[i]
    end
    for i in rr.I2
        for k in 1:(rr.W[i,2] < Inf ? 2 : 1)
            rr.τs[i, k]    = rr.W[i,k]*rr.Θ[i]
            rr.cdfs[i, k]  = cdf(rr.dist, rr.τs[i,k])
        end
        rr.Δcdf[i] = rr.cdfs[i,2] - rr.cdfs[i,1]
        rr.loglik[i] = log(rr.Δcdf[i])
    end
    return sum(rr.loglik)
end

function diff_ders!(gl, hl, gli, hli, fin, Δcdf)
    @inbounds for j1 in eachindex(gl)
        gl[j1] = (fin*gli[2][j1]  - gli[1][j1])/Δcdf
    end
    @inbounds for j2 in eachindex(gl), j1 in eachindex(gl)
        hl[j1,j2] = (fin*hli[2][j1,j2]- hli[1][j1,j2])/Δcdf - gl[j1]*gl[j2]
    end
end

function subtract_ders!(grad, hes, rr, int_coefs)
    grad[:] = 0.
    hes[:,:] = 0.
    gl, hl = zeros(rr.M+1), zeros(rr.M+1,rr.M+1)
    gli, hli = [zeros(rr.M+1) for i in 1:2], [zeros(rr.M+1,rr.M+1) for i in 1:2]
    for i in rr.I1
        d²l!(gl, hl, rr.dist, rr.τs[i,1], int_coefs)
        subtract_ders!(grad, hes, gl, hl, @view(rr.X[i,:]), rr.M, rr.N)
    end
    for i in rr.I2
        fin = rr.W[i,2] < Inf
        @inbounds for k in 1:(fin ? 2 : 1)
            d²l_int!(gli[k], hli[k], rr.dist, rr.τs[i,k], rr.cdfs[i,k], int_coefs)
        end
        diff_ders!(gl, hl, gli, hli, fin, rr.Δcdf[i])
        subtract_ders!(grad, hes, gl, hl, @view(rr.X[i,:]), rr.M, rr.N)
    end
end

function subtract_ders!(grad, hes, gl, hl, x, M , N)
    @inbounds for j1 = 1:(M+N)
        grad[j1] -= gl[min(j1,M+1)]*x[max(j1-M,1)]
    end
    @inbounds for j2 = 1:(M+N), j1 = 1:(M+N)
        hes[j1,j2] -= hl[min(j1,M+1),min(j2,M+1)]*x[max(j1-M,1)]*x[max(j2-M,1)]
    end
end

function d²l_int!(gl, hl, pdist::Distribution, t, x, int_coefs)
    @inbounds for j1 = 1:length(pdist.params)
        gl[j1] = clenshaw_asin(x,int_coefs.g[j1])
    end
    gl[end] = -t*pdf(pdist, t)
    last_line = @view hl[end, :]
    dl!(last_line, pdist, t, int_coefs)
    scale!(last_line, gl[end])
    @inbounds for j2 in 1:length(pdist.params), j1 in j2:length(pdist.params)
        hl[j1,j2] = clenshaw_asin(x,int_coefs.hggt[j1,j2])
    end
    @inbounds for j2 in 1:length(pdist.params)+1, j1 in 1:j2-1
        hl[j1,j2] = hl[j2,j1]
    end
end

auxvec(pdist::Distribution) = Array{Float64,1}(0)

initialize_aft(S::AbstractVector,X::AbstractArray, pdist::Distribution, N) =
vcat(pdist.params,zeros(N))


text:
## QUESTION 2

include("md.jl");
qstart = loadmdtrial("../../swp/1069_met_start.mdv")

δpe = 0.1                     # interaction distance
δke = 0.1                     # interaction distance
N = 20000                   # simulation steps
M = 10                      # steps between saved measurements
P = int(N/M)

println("Running ",N," steps for ",P," measurements")

q = addsteps(qstart, N)     # N new steps
pe = Float64[]
ke = Float64[]
v  = Float64[]

# Run for 10 steps at a time, saving values along the way
while !isfinished(q)
    metropolis!(q, M, δpe, δke)    # M steps between saves
    push!(pe, q.pet[q.currentStep])
    push!(ke, q.ket[q.currentStep])
    push!(v,  q.vir[q.currentStep])
end

# save our results
writedlm("../../swp/set04/q1/pe.dat", pe)
writedlm("../../swp/set04/q1/ke.dat", ke)
writedlm("../../swp/set04/q1/v.dat",  v)


text:
module TestUtils

using Distributions
using LinearAlgebra
using Random
using Test

__rand(::Nothing, args...) = rand(args...)
__rand(rng::AbstractRNG, args...) = rand(rng, args...)

__rand!(::Nothing, args...) = rand!(args...)
__rand!(rng::AbstractRNG, args...) = rand!(rng, args...)

"""
    test_mvnormal(
        g::AbstractMvNormal, n_tsamples::Int=10^6, rng::AbstractRNG=Random.GLOBAL_RNG
    )

Test that `AbstractMvNormal` implements the expected API.
"""
function test_mvnormal(
    g::AbstractMvNormal, n_tsamples::Int=10^6, rng::Union{AbstractRNG, Nothing}=nothing
)
    d = length(g)
    μ = mean(g)
    Σ = cov(g)
    @test length(μ) == d
    @test size(Σ) == (d, d)
    @test var(g) ≈ diag(Σ)
    @test entropy(g) ≈ 0.5 * logdet(2π * ℯ * Σ)
    ldcov = logdetcov(g)
    @test ldcov ≈ logdet(Σ)
    vs = diag(Σ)
    @test g == typeof(g)(params(g)...)
    @test g == deepcopy(g)
    @test minimum(g) == fill(-Inf, d)
    @test maximum(g) == fill(Inf, d)
    @test extrema(g) == (minimum(g), maximum(g))
    @test isless(extrema(g)...)

    # test sampling for AbstractMatrix (here, a SubArray):
    subX = view(__rand(rng, d, 2d), :, 1:d)
    @test isa(__rand!(rng, g, subX), SubArray)

    # sampling
    @test isa(__rand(rng, g), Vector{Float64})
    X = __rand(rng, g, n_tsamples)
    emp_mu = vec(mean(X, dims=2))
    Z = X .- emp_mu
    emp_cov = (Z * Z') * inv(n_tsamples)

    mean_atols = 8 .* sqrt.(vs ./ n_tsamples)
    cov_atols = 10 .* sqrt.(vs .* vs') ./ sqrt.(n_tsamples)
    for i = 1:d
        @test isapprox(emp_mu[i], μ[i], atol=mean_atols[i])
    end
    for i = 1:d, j = 1:d
        @test isapprox(emp_cov[i,j], Σ[i,j], atol=cov_atols[i,j])
    end

    X = rand(MersenneTwister(14), g, n_tsamples)
    Y = rand(MersenneTwister(14), g, n_tsamples)
    @test X == Y
    emp_mu = vec(mean(X, dims=2))
    Z = X .- emp_mu
    emp_cov = (Z * Z') * inv(n_tsamples)
    for i = 1:d
        @test isapprox(emp_mu[i]   , μ[i]  , atol=mean_atols[i])
    end
    for i = 1:d, j = 1:d
        @test isapprox(emp_cov[i,j], Σ[i,j], atol=cov_atols[i,j])
    end

    # evaluation of sqmahal & logpdf
    U = X .- μ
    sqm = vec(sum(U .* (Σ \ U), dims=1))
    for i = 1:min(100, n_tsamples)
        @test sqmahal(g, X[:,i]) ≈ sqm[i]
    end
    @test sqmahal(g, X) ≈ sqm

    lp = -0.5 .* sqm .- 0.5 * (d * log(2.0 * pi) + ldcov)
    for i = 1:min(100, n_tsamples)
        @test logpdf(g, X[:,i]) ≈ lp[i]
    end
    @test logpdf(g, X) ≈ lp

    # log likelihood
    @test loglikelihood(g, X) ≈ sum(i -> Distributions._logpdf(g, X[:,i]), 1:n_tsamples)
    @test loglikelihood(g, X[:, 1]) ≈ logpdf(g, X[:, 1])
    @test loglikelihood(g, [X[:, i] for i in axes(X, 2)]) ≈ loglikelihood(g, X)
end

end


text:
export mc_integral

function mc_integral(f::Function, range::StepRangeLen)
    start = Float64(range.ref)
    δ = Float64(range.step)
    len = Int(range.len)
    stop = start + δ * (len-1)

    xs = rand(len) .* (stop-start) .+ start
    ys = f.(xs)

    return sum(ys) * δ
end


text:
using Distributions, StatsBase

function cross_entropy(f, P; k=1, m=50, m_extra=20, m_elite=m, base_dist=Uniform(-1,1))
    for _=1:k
        # Sample from the current distribution
        samples = clamp.(rand(P, m), -1f0, 1f0)
        dim = size(samples, 1)

        # Add some noise to avoid the distribution collapsing entirely
        # samples .+= rand(Uniform(-1f-3, 1f-3), size(samples)...)

        # Add some extra samples from an exploration policy to account for dramatic shifts from previous iterations
        samples = hcat(samples, rand(base_dist, dim, m_extra))

        vals = [f(samples[:,i]) for i=1:size(samples,2)]
        weights = softmax(-vals)

        order = sortperm(vals)
        elite_samples = Float64.(samples[:, order[1:m_elite]])
        try
            P = fit(MvNormal, elite_samples, weights[order[1:m_elite]])
        catch
            println("weights: ", weights)
            @error "Error in fitting MVNormal distribution"
        end
    end
    P
end

function mcmc(f, P; Q=(sz...)->0.01*randn(sz...), k=1, m=size(P,2), m_extra=20, base_dist = Uniform(-1,1))
    dim, Npart=size(P)
    @assert Npart == m
    for _=1:k
        # Perturb particles and add some additional ones
        # P = Float32.(clamp.(P .+ Q(dim, Npart), -1, 1))
        P = hcat(P, rand(base_dist, dim, m_extra))

        # compute the particle values
        vals = [f(P[:,i]) for i=1:size(P,2)]
        weights = exp.(-vals)

        # resample
        samples = sample(1:Npart+m_extra, Weights(weights), Npart)
        P = P[:, samples]
    end
    P
end

best_estimate(d::MvNormal) = d.μ
uncertainty(d::MvNormal) = det(d.Σ)
best_estimate(P::Array) = mean(P, dims=2)
uncertainty(P::Array) = mean(std(P, dims=2))

# # Demo the inference
# target_dist = MvNormal([-0.5,-0.5], [0.1, 0.1])
# foo(z) = -logpdf(target_dist, z)
# x = -2:0.01:2
# y = x
#
# ## Particle filter
#
# contour(x, y, (x,y)->exp(logpdf(target_dist, [x,y])))
# P = randn(2,100) # Initial particle set
# scatter!(P[1,:], P[2,:], label="iteration 0")
# P = mcmc(foo, P)
# scatter!(P[1,:], P[2,:], label="iteration 1")
# P = mcmc(foo, P)
# scatter!(P[1,:], P[2,:], label="iteration 2")
#
# println("best estimate: ", best_estimate(P))
#
# ## cross entropy
# p1 = contour(x, y, (x,y)->exp(logpdf(target_dist, [x,y])), cbar=false)
# P = MvNormal([0,0], Diagonal([1.,1.]))
# p2 = contour(x, y, (x,y)->exp(logpdf(P, [x,y])), seriescolor=:blues)
# P = cross_entropy(foo, P)
# p3 = contour(x, y, (x,y)->exp(logpdf(P, [x,y])), seriescolor=:blues)
# P = cross_entropy(foo, P)
# p4 = contour(x, y, (x,y)->exp(logpdf(P, [x,y])), seriescolor=:blues)
#
# plot(p1, p2, p3, p4)
#
# println("best estimate: ", best_estimate(P))


text:
using Test, SpecialFunctions
import TruncatedNormal as TruncNorm

@test TruncNorm.tnmom2(-Inf, Inf) == 1
@test TruncNorm.tnmom2(+Inf, +Inf) == Inf
@test TruncNorm.tnmom2(-Inf, -Inf) == Inf
@test TruncNorm.tnmom2(-Inf, 0) == 1
@test TruncNorm.tnmom2(0, +Inf) == 1

@test isnan(TruncNorm.tnmom2(1, 0))
@test isnan(TruncNorm.tnmom2(+Inf, 0))
@test isnan(TruncNorm.tnmom2(0, -Inf))
@test isnan(TruncNorm.tnmom2(+Inf, -Inf))
@test isnan(TruncNorm.tnmom2(NaN, NaN))
@test isnan(TruncNorm.tnmom2(0, NaN))
@test isnan(TruncNorm.tnmom2(NaN, 0))

for x = -10:10
    @test TruncNorm.tnmom2(x, +Inf) ≈ 1 + x * √(2/π) / erfcx(+x/√2)
    @test TruncNorm.tnmom2(-Inf, x) ≈ 1 - x * √(2/π) / erfcx(-x/√2)
end

@test TruncNorm.tnmom2(-Inf, Inf, 0, 1) == 1

for x = -10:10
    @test TruncNorm.tnmom2(x, +Inf, 0, 1) ≈ TruncNorm.tnmom2(x, +Inf)
    @test TruncNorm.tnmom2(-Inf, x, 0, 1) ≈ TruncNorm.tnmom2(-Inf, x)
end

@test TruncNorm.tnmom2(-1e200, 1e200) ≈ 1
@test TruncNorm.tnmom2(-1e6, -999000) ≈ 9.9800100000199999999999799599399200002e11
@test TruncNorm.tnmom2(+1e6, +Inf) ≈ 1.00000000000199999999999800000e12
@test TruncNorm.tnmom2(-Inf, -1e6) ≈ 1.00000000000199999999999800000e12

@test TruncNorm.tnmom2(-Inf, Inf, 2, 0) ≈ 4
@test TruncNorm.tnmom2(0, 20, 3, 3//2) ≈ 11.4986153820554548159603543747
@test TruncNorm.tnmom2(10, 20, 3, 3//2) ≈ 106.111726253622178129044932133
@test TruncNorm.tnmom2(-10, 20, 3, 3//2) ≈ 11.2499999999999997949173964131


text:
using Documenter, SpinGlassDynamics
# makedocs(
#     modules=[SpinGlassDynamics],
#     sitename="SpinGlassTensors.jl",
#     format=Documenter.LaTeX()
# )
makedocs(
    modules=[SpinGlassDynamics],
    sitename="SpinGlassDynamics.jl",
    format=Documenter.HTML(
        prettyurls=get(ENV, "CI", nothing) == "true"
    )
)
deploydocs(
    repo="github.com/euro-hpc-pl/SpinGlassDynamics.jl.git",
)


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# deployment (allocation) of ambulances to stations

# make a single random deployment
# assumes ambulances are homogeneous
function makeRandDeployment(numAmbs::Int, numStations::Int;
	stationCapacities::Union{Vector{Int},Nothing} = nothing,
	rng::AbstractRNG = GLOBAL_RNG)::Deployment

	@assert(numStations > 0)
	if stationCapacities == nothing
		return sort(rand(rng, 1:numStations, numAmbs))
	else
		@assert(numStations == length(stationCapacities))
		@assert(numAmbs <= sum(stationCapacities))
		remainingCapacity = copy(stationCapacities)
		unfilledStations = Set(findall(remainingCapacity))
		deployment = zeros(Int, numAmbs)
		for i = 1:numAmbs
			j = rand(rng, unfilledStations) # station index
			deployment[i] = j
			remainingCapacity[j] -= 1
			if remainingCapacity[j] == 0
				delete!(unfilledStations, j)
			end
		end
		return sort(deployment)
	end
end

# generate a number of unique random deployments
# ignores station capacities, assumes ambulances are homogeneous
function makeRandDeployments(numAmbs::Int, numStations::Int, numDeployments::Int;
	rng::AbstractRNG = GLOBAL_RNG)
	@assert(numAmbs >= 1)
	@assert(numStations >= 1)
	@assert(numDeployments >= 1)
	@assert(numDeployments <= binomial(numAmbs + numStations - 1, numStations - 1))

	deployments = Set{Deployment}()
	while length(deployments) < numDeployments
		push!(deployments, makeRandDeployment(numAmbs, numStations; rng = rng))
	end
	deployments = collect(deployments)

	return deployments
end
function makeRandDeployments(sim::Simulation, numDeployments::Int;
	rng::AbstractRNG = GLOBAL_RNG)
	return makeRandDeployments(sim.numAmbs, sim.numStations, numDeployments; rng = rng)
end

function deploymentToStationsNumAmbs(deployment::Deployment, numStations::Int)
	# deployment[i] gives the station index for ambulance i
	stationsNumAmbs = zeros(Int, numStations)
	for stationIndex in deployment
		stationsNumAmbs[stationIndex] += 1
	end
	return stationsNumAmbs
end

function stationsNumAmbsToDeployment(stationsNumAmbs::Vector{Int})
	# stationsNumAmbs[i] gives the number of ambulances for station i
	# will assume that all ambulances are equivalent
	deployment = Deployment() # deployment[i] gives the station index for ambulance i
	for (stationIndex, numAmbs) in enumerate(stationsNumAmbs), i = 1:numAmbs
		push!(deployment, stationIndex)
	end
	return deployment
end

function getDeployment(sim::Simulation)::Deployment
	return map(a -> a.stationIndex, sim.ambulances)
end

function getStationsNumAmbs(sim::Simulation)
	stationsNumAmbs = zeros(Int, sim.numStations)
	for ambulance in sim.ambulances
		stationsNumAmbs[ambulance.stationIndex] += 1
	end
	return stationsNumAmbs
	# return deploymentToStationsNumAmbs(getDeployment(sim), sim.numStations)
end

# set ambulance to be located at given station
# this is hacky, expects ambulance will only need changing, this may not always be true
# mutates: ambulance
function setAmbStation!(sim::Simulation, ambulance::Ambulance, station::Station)
	ambulance.stationIndex = station.index
	ambulance.route = Route()
	initRoute!(sim, ambulance.route; startLoc = station.location, startFNode = station.nearestNodeIndex, startFNodeDist = station.nearestNodeDist)
end

# apply deployment 'deployment' for sim.ambulances and sim.stations
# mutates: sim.ambulances
function applyDeployment!(sim::Simulation, deployment::Deployment)
	n = length(deployment)
	@assert(n == sim.numAmbs)
	@assert(all(d -> in(d, 1:sim.numStations), deployment))
	# @assert(all(1 .<= deployment .<= sim.numStations))
	for i = 1:n
		setAmbStation!(sim, sim.ambulances[i], sim.stations[deployment[i]])
	end
end

function applyStationsNumAmbs!(sim::Simulation, stationsNumAmbs::Vector{Int})
	applyDeployment!(sim, stationsNumAmbsToDeployment(stationsNumAmbs))
end

# runs simulation for the deployment
# mutates: sim
function simulateDeployment!(sim::Simulation, deployment::Deployment)
	reset!(sim)
	applyDeployment!(sim, deployment)
	simulateToEnd!(sim)
end

# runs simulation for each deployment
# returns vector of function applied to end of each simulation run
# mutates: sim
function simulateDeployments!(sim::Simulation, deployments::Vector{Deployment}, f::Function;
	showEta::Bool = false)
	numDeployments = length(deployments)
	results = Vector{Any}(undef, numDeployments)
	t = time() # for estimating expected time remaining
	for i = 1:numDeployments
		simulateDeployment!(sim, deployments[i])
		results[i] = f(sim)
		if showEta
			remTime = (time() - t) / i * (numDeployments - i) / 60 # minutes
			print(rpad(string(" remaining run time (minutes): ", ceil(remTime, digits = 1)), 100, " "), "\r")
		end
	end
	reset!(sim)
	return results
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

function exec_map!(f::Base.Callable, executor::Transducers.Executor, Y::AbstractVector, X::AbstractVector)
    Folds.map!(f, Y, X, executor)
    return Y
end

function exec_map!(f::Base.Callable, executor::Transducers.DistributedEx, Y::AbstractVector, X::AbstractVector)
    Y .= Folds.map(f, X, executor)
    return Y
end


text:
using CSDP
using Base.Test
using Base.LinAlg.BlasInt

@testset "Interact with BLAS" begin
    vec = Cdouble[1.0, 2.0, 0.0, -1.0]
    l = length(vec)
    inc = 1
    n1 = ccall((BLAS.@blasfunc(dasum_), LinAlg.BLAS.libblas),
               Cdouble,
               (Ptr{BlasInt}, Ptr{Cdouble}, Ptr{BlasInt}),
               &l, vec, &inc)
    @test abs(n1 - 4) < 1e-15
end

@testset "Call libcsdp.norm1" begin
    vec = Cdouble[1.0, 2.0, 0.0, -1.0]
    try
        n1 = ccall((:norm1, CSDP.csdp),
                   Cdouble,
                   (Cint, Ptr{Cdouble}),
                   length(vec), vec)
        @test abs(n1 - 4) < 1e-15
    catch
        println("n1 = $n1, vec=$vec, length(vec)=$(length(vec))")
        rethrow()
    end
end

@testset "Example" begin
    cd("../examples/") do
        include(joinpath(pwd(), "example.jl"))
        @test size(X) == (7, 7)
        @test length(y) == 2
        @test size(Z) == (7, 7)
        X✓ = [3 3  0 0 0 0 0;
              3 3  0 0 0 0 0;
              0 0 16 0 0 0 0;
              0 0  0 0 0 0 0;
              0 0  0 0 0 0 0;
              0 0  0 0 0 0 0;
              0 0  0 0 0 0 0] / 24
        @test norm(AbstractArray(X) - X✓) < 1e-6
        y✓ = [3, 4] / 4
        @test norm(AbstractArray(y) - y✓) < 1e-6
        Z✓ = [1 -1 0 0 0 0 0;
             -1  1 0 0 0 0 0;
              0  0 0 0 0 0 0;
              0  0 0 8 0 0 0;
              0  0 0 0 8 0 0;
              0  0 0 0 0 3 0;
              0  0 0 0 0 0 4] / 4
        @test norm(AbstractArray(Z) - Z✓) < 1e-6
    end
end

@testset "Linear tests" begin
    include(joinpath(Pkg.dir("MathProgBase"),"test","conicinterface.jl"))
end

# FIXME fails on Windows 32 bits... Maybe I should put linear vars/cons
# in a diagonal matrix in SemidefiniteModels.jl instead of many 1x1 blocks
@static if !is_windows() || Sys.WORD_SIZE != 32
    @testset "Conic linear tests" begin
        coniclineartest(CSDP.CSDPSolver(), duals=true, tol=1e-6)
    end

    @testset "Conic SOC tests" begin
        conicSOCtest(CSDP.CSDPSolver(write_prob="soc.prob"), duals=true, tol=1e-6)
    end

    @testset "Conic SOC rotated tests" begin
        conicSOCRotatedtest(CSDP.CSDPSolver(), duals=true, tol=1e-6)
    end
end

@testset "Conic SDP tests" begin
    conicSDPtest(CSDP.CSDPSolver(), duals=false, tol=1e-6)
end


text:
convert(::Type{Mag}, x::Int64) = Mag(x)

convert(::Type{Mag}, x::UInt64) = Mag(x)

convert(::Type{Mag}, x::Float64) = Mag(x)

convert(::Type{ArbFloat{P}}, x::Int64) where {P} = ArbFloat{P}(x)

convert(::Type{ArbFloat{P}}, x::UInt64) where {P} = ArbFloat{P}(x)

convert(::Type{ArbFloat{P}}, x::Float64) where {P} = ArbFloat{P}(x)

convert(::Type{ArbReal{P}}, x::Int64) where {P} = ArbReal{P}(x)

convert(::Type{ArbReal{P}}, x::UInt64) where {P} = ArbReal{P}(x)

convert(::Type{ArbReal{P}}, x::Float64) where {P} = ArbReal{P}(x)

convert(::Type{ArbComplex{P}}, x::Int64) where {P} = ArbComplex{P}(x)

convert(::Type{ArbComplex{P}}, x::UInt64) where {P} = ArbComplex{P}(x)

convert(::Type{ArbComplex{P}}, x::Float64) where {P} = ArbComplex{P}(x)


text:
"""
    overwrite_netcdf2d(fname, var, varname)

Overwrite variable `varname` with field `var` in netcdf file `fname`
"""
function overwrite_netcdf2d(fname, var, varname)
    ncin = NetCDF.open(fname)
    ncwrite(var, fname, varname)
end


text:
using Flux
using Images
using CUDA
using Zygote
using Metalhead
using Statistics
using LinearAlgebra
using Plots
using IJulia
using Tar
using BSON:@save

if has_cuda()
    CUDA.allowscalar(false)
    device = gpu
    @info "CUDA is on"
else
    device = cpu
    @info "CUDA is off"
end

mutable struct vgg19
    style1
    style2
    style3
    style4
    style5
    content1
end

Flux.@functor vgg19

function vgg19()
    vgg = VGG19().layers
    style1 = Chain(vgg[1:1]...)
    style2 = Chain(vgg[2:4]...)
    style3 = Chain(vgg[5:7]...)
    style4 = Chain(vgg[8:12]...)
    style5 = Chain(vgg[13:17]...)
    content1 = Chain(vgg[18:18]...)
    vgg19(style1, style2, style3, style4, style5, content1)
end

function content(m::vgg19, x)
    x |> m.style1 |> m.style2 |> m.style3 |> m.style4 |> m.style5 |> m.content1
end

function style(m::vgg19, x)
    s1 = m.style1(x)
    s2 = m.style2(s1)
    s3 = m.style3(s2)
    s4 = m.style4(s3)
    s5 = m.style5(s4)
    (s1, s2, s3, s4, s5)
end

function gram_matrix(mat)
    h, w, c, n = size(mat)
    nm = reshape(mat, h * w, c)
    gram = transpose(nm) * nm
    gram
end

function content_loss(current, target)
    1 / 2 * sum((current - target).^2)
end

function style_loss(current_features, target_features)
    num = length(current_features)
    ls = 0f0
    weights = [i for i in 1:num]
    for i in 1:num
        h, w, c, n = size(current_features[i])
        lloss = 1 / (4 * c^2 * (h * w)^2) * sum((gram_matrix(current_features[i]) .- gram_matrix(target_features[i]))^2)
        ls += weights[i] * lloss
    end
    return ls
end

function tv_loss(target_image)
    img = target_image
    h, w, c, q = size(img)
    ver_comp = sum((img[2:end,:,:,:] - img[1:end - 1,:,:,:]).^2)
    hor_comp = sum((img[:,2:end,:,:] - img[:,1:end - 1,:,:]).^2)
    1 / (4 * h * w * c) * (ver_comp + hor_comp).^2
end

function total_loss(model::vgg19, target_image, contentf, stylef; α=1e-3, β=1, γ=1, splat::Bool=false)
    target_style = style(model, target_image)
    target_content = content(model, target_image)
    c_loss = content_loss(target_content, contentf)
    s_loss = style_loss(target_style, stylef) / 10
    t_loss = tv_loss(target_image)
    total = α * c_loss + β * s_loss + γ * t_loss

    if splat
        (total, α * c_loss, β * s_loss, γ * t_loss)
    else
        total
    end
end

function img_to_array(img)
    r, g, b = red.(img), green.(img), blue.(img)
    cat(r, g, b, dims=3)
end

function norm(x)
    means = [123.68, 116.779, 103.939 ]
    stds = [58.393, 57.12, 57.375]
    x .*= 255.0f0
    for i in 1:3
        x[:,:,i] .-= means[i]
        x[:,:,i] ./= stds[i]
    end
    x
end

function preprocess_array(img_arr, dsize::Tuple{Integer,Integer}=(226, 226))
    Flux.unsqueeze(imresize((2 * map(clamp01nan, Float32.(img_arr)) .- 1), dsize), 4)
end

function deprocess_array(arr)
    # Last step towards image transformation
    nim = clamp.(arr, -1f0, 1f0) |> cpu
    nim .-= minimum(nim)
    nim ./= maximum(nim)
    if length(size(nim)) == 4
        return nim[:,:,:,1]
    end
    map(clamp01nan, nim)
end

function dep_array(arr)
    nim = arr |> cpu

    means = [123.68, 116.779, 103.939 ]
    stds = [58.393, 57.12, 57.375]
    for i in 1:3
        nim[:,:,i] .*= stds[i]
        nim[:,:,i] .+= means[i]
    end
    nim ./= 255.0f0
end

function arr_to_img(iar)
    arr = Array(iar)
    return RGB.([arr[:,:,i] for i in 1:3]...)
end

styles = Dict(
    :starry => "https://upload.wikimedia.org/wikipedia/commons/thumb/e/ea/Van_Gogh_-_Starry_Night_-_Google_Art_Project.jpg/300px-Van_Gogh_-_Starry_Night_-_Google_Art_Project.jpg",
    :scream => "https://upload.wikimedia.org/wikipedia/commons/thumb/c/c5/Edvard_Munch%2C_1893%2C_The_Scream%2C_oil%2C_tempera_and_pastel_on_cardboard%2C_91_x_73_cm%2C_National_Gallery_of_Norway.jpg/1200px-Edvard_Munch%2C_1893%2C_The_Scream%2C_oil%2C_tempera_and_pastel_on_cardboard%2C_91_x_73_cm%2C_National_Gallery_of_Norway.jpg",
    :wave => "https://upload.wikimedia.org/wikipedia/commons/thumb/a/a5/Tsunami_by_hokusai_19th_century.jpg/800px-Tsunami_by_hokusai_19th_century.jpg",
    :giorgio => "https://upload.wikimedia.org/wikipedia/commons/thumb/d/da/Claude_Monet%2C_Saint-Georges_majeur_au_cr%C3%A9puscule.jpg/1200px-Claude_Monet%2C_Saint-Georges_majeur_au_cr%C3%A9puscule.jpg",
    :naif => "https://instagram.fbcn11-1.fna.fbcdn.net/v/t51.2885-15/e35/s1080x1080/116798277_581494632757362_4406532804858518791_n.jpg?_nc_ht=instagram.fbcn11-1.fna.fbcdn.net&_nc_cat=106&_nc_ohc=7zx1nibyyZ8AX8Os3uD&tp=1&oh=b9301bf53db67ea59518f2d1ddec5024&oe=604ED2A8",
    :wheat => "https://d16kd6gzalkogb.cloudfront.net/magazine_images/Vincent-van-Gogh-Whaet-Field-with-Cypresses.-Image-via-wikimedia.org_.jpg",
    :sunrise => "https://www.artblr.com/upload/userfiles/images/Sunrise.jpg",
    :nature => "https://sothebys-com.brightspotcdn.com/dims4/default/cd1ca73/2147483647/strip/true/crop/859x859+71+0/resize/1200x1200!/quality/90/?url=http%3A%2F%2Fsothebys-brightspot.s3.amazonaws.com%2Fdotcom%2Fa4%2Fb2%2F1f464df248a8bbfe8466856b32bc%2Fcover-pics.jpg",
    :curly => "https://thumbs.dreamstime.com/b/colorful-linear-wavy-texture-endless-abstract-pattern-template-design-decoration-88890825.jpg",
    :brick => "https://external-content.duckduckgo.com/iu/?u=https%3A%2F%2Fimages.freecreatives.com%2Fwp-content%2Fuploads%2F2016%2F12%2FFree-Brick-Wall-Texture.jpg&f=1&nofb=1",
    :oil => "https://encrypted-tbn0.gstatic.com/images?q=tbn:ANd9GcSXq4Jf5qCXKxJTOLqbULdFo8nsmN6qW_F0ZQ&usqp=CAU",
    :cubism => "https://www.artyfactory.com/art_appreciation/art_movements/art-movements/cubism/picasso_cubism.jpg",
    :composition => "https://upload.wikimedia.org/wikipedia/commons/b/b4/Vassily_Kandinsky%2C_1913_-_Composition_7.jpg",
)

function transfer_color(source::Array{RGB{Normed{UInt8,8}},2}, target::Array{RGB{Normed{UInt8,8}},2})::Array{RGB{Normed{UInt8,8}},2}
    px = x -> [x.r,x.g,x.b]
    s̄ = mean(px, source)
    t̄ = mean(px, target)

    covₛ = mean(x -> (px(x) - s̄) * transpose(px(x) - s̄), source)
    covₜ = mean(x -> (px(x) - t̄) * transpose(px(x) - t̄), target)
    vaₛ, veₛ = eigen(covₛ)
    vaₜ, veₜ = eigen(covₜ)

    covsqₛ = veₛ * diagm(vaₛ).^(1 / 2) * veₛ'
    covsqₜ = veₜ * inv(diagm(vaₜ).^(1 / 2)) * veₜ'
    A = covsqₛ * covsqₜ
    b = s̄ - A * t̄
    map(target) do x
        r = A * px(x) + b
        RGB(map(clamp01nan, r)...)
    end
end

function get_images(style::Symbol=:starry, transfer::Bool=true; content_url=nothing, dsize=(226, 226))
    if content_url === nothing
        cim = Images.load("img.jpg")
    else
        cim = mktemp() do fn, f
            download(content_url, fn)
            load(fn)
        end
    end
    if dsize === nothing
        dsize = size(cim)
    end
    sim = mktemp() do fn, f
        download(styles[style], fn)
        load(fn)
    end
    if transfer
        colim = transfer_color(sim, cim)
    else
        colim = copy(cim)
    end
    return preprocess_array.(map(img_to_array, (cim, colim, sim)), (dsize,))
end

function save_image(fn::String, image)
    im = arr_to_img(deprocess_array(image))
    Images.save(fn, im)
    return nothing
end

function make_composition(o, s, t)::Array{Float32,3}
    i1, i2, i3 = map(deprocess_array, (o, s, t))
    hcat(i1, i2, i3)
end

function save_composition(fn::String, original, style, final)
    narr = make_composition(copy.((original, style, final))...)
    if !isdir("outputs")
        mkdir("outputs")
    end

    save_image("outputs/" * fn, narr)
    return nothing
end

function save_run_and_clean()
    if !isdir("runs")
        mkdir("runs")
    end
    files = map(x -> parse(Int, split(split(x, "_")[end], ".")[begin]), readdir("runs"))
    last = length(files) > 0 ? maximum(files) + 1 : 1

    Tar.create("outputs", "runs/run_$last.tar")
    rm("outputs", recursive=true)
return nothing
end

function train(sty::Symbol=:starry, n::Int=50;plt::Bool=false,α=1e-2, β=10, γ=1e-1, init::Symbol=:content,lr::Float32=0.02f0,transfer::Bool=true,content_url=nothing,dsize=(226, 226),save_every=nothing)
    rm("outputs", force=true, recursive=true)
    GC.gc()
    CUDA.reclaim()
    @info "Loading model"
    vgg = @time vgg19() |> device
    @info "Loading images"
    ocontent_image, content_image, style_image = @time get_images(sty, transfer, content_url=content_url, dsize=dsize)
    content_image = content_image |> device
    style_image = style_image |> device
    if init == :rand
        target_image = (2 .* rand(Float32, size(content_image)) .- 1) |> device
    elseif init == :content
        target_image = copy(content_image)
    elseif init == :style
        target_image = copy(style_image)
    elseif init == :diff
        target_image = content_image - style_image
    elseif init == :sum
        target_image = (content_image + style_image) ./ 2
    end
    opt = ADAM(lr, (0.9, 0.999))

    stylef = style(vgg, style_image)
    contentf = content(vgg, content_image)
    ipa = Flux.params(target_image)

    losses = Array{Float32,1}()
    c_loss = Array{Float32,1}()
    s_loss = Array{Float32,1}()
    t_loss = Array{Float32,1}()

    loss(x1; kwargs...) = total_loss(vgg, x1, contentf, stylef; α=α, β=β, γ=γ, kwargs...)
    last_iteration = 0
    save_composition("output_0.jpg", ocontent_image, style_image, target_image)
    for i = 1:n
        @info "Iteration $i"
        gs = gradient(ipa) do
            loss(target_image)
        end
        Flux.update!(opt, ipa, gs)
        ls, cl, sl, tl = loss(target_image, splat=true)
        push!(losses, ls)
        push!(c_loss, cl)
        push!(s_loss, sl)
        push!(t_loss, tl)
        if plt
            try
                IJulia.clear_output(true)
            catch e
            end
            Plots.display(arr_to_img(make_composition(ocontent_image, style_image, target_image)))
        end
        if save_every !== nothing
            if i % save_every == 0
                save_composition("output_$i.jpg", ocontent_image, style_image, target_image)
            end
        end
        last_iteration = i
    end
    save_composition("output_$last_iteration.jpg", ocontent_image, style_image, target_image)
    @info "Finished"
        if plt
        plot(losses, title="Losses", label="Total loss", )
        plot!(c_loss, label="Content loss")
        plot!(s_loss, label="Style loss")
        f1 = plot!(t_loss, label="Total variation loss", yaxis=(:log, (1, Inf)))
        png(f1, "outputs/losses_plot.png")
        Plots.display(f1)

        h, w, c, n = size(content_image)
        histogram(reshape(deprocess_array(ocontent_image), h * w, c), α=0.3, layout=3, label="content")
        histogram!(reshape(clamp.(deprocess_array(target_image), 0f0, 1f0), h * w, c), α=0.3, layout=3, label="target")
        h1 = histogram!(reshape(deprocess_array(style_image), h * w, c), α=0.3, layout=3, label="style")
        png(h1, "outputs/histograms.png")
        Plots.display(h1)
    end
    @save "outputs/losses.bson" losses c_loss s_loss t_loss α β γ
    save_run_and_clean()
    return target_image
end

max_stars_repo_name: ramonus/deep-style-transfer

max_issues_repo_name: ramonus/deep-style-transfer

max_forks_repo_name: ramonus/deep-style-transfer


text:
using ApproxFun, MultivariateOrthogonalPolynomials, BlockArrays, FillArrays, Plots
import PyPlot
x,y = Fun(Triangle())
Δ = Laplacian() : TriangleWeight(1.0,1.0,1.0,JacobiTriangle(1.0,1.0,1.0))
V = x*y^2
L = Δ + 200^2*V
u = L \ ones(Triangle())

M = L[Block.(1:50),Block.(1:50)];

PyPlot.spy(M)

@time cfs = M \ [1; Zeros(size(M,1)-1)];
u = Fun(domainspace(Δ), cfs)
p = MultivariateTriangle.contourf(u)

u(0.1,0.2)

cfs1000 = cfs
cfs1001 = cfs


text:
using DataFrames
using Schrute
using Test

@testset "Schrute.jl" begin

    @test DataFrames._ncol(theOffice()) == 12
    @test DataFrames._nrow(theOffice()) == 55130

end


text:
struct VarianceComponents{Q <: AbstractMatrix{Float64},
                          R <: AbstractMatrix{Float64},
                          S <: AbstractMatrix{Float64},
                          T <: AbstractMatrix{Float64}}
    Ψ::Q
    Σ::R
    L::S
    Λ::T
end

function treeDimension(tsm::TraitSimulationModel)
    return get_dimension(tsm.treeModel)
end

function dataDimension(tsm::TraitSimulationModel)
    return isnothing(tsm.extensionModel) ?
            treeDimension(tsm) : dataDimension(tsm.extensionModel)
end

function get_loadings(tsm::TraitSimulationModel)
    if isnothing(tsm.extensionModel)
        return I
    end
    return get_loadings(tsm.extensionModel)
end

function get_residual_variance(tsm::TraitSimulationModel)
    if isnothing(tsm.extensionModel)
        return Diagonal(zeros(dataDimension(tsm)))
    end
    return get_residual_variance(tsm.extensionModel)
end

function tree_variance(tsm::TraitSimulationModel; pss::Float64 = Inf)
    Ψ = vcv(tsm.treeModel.tree, tsm.taxa)
    n = size(Ψ, 1)
    Ψ .+= 1.0 / pss * ones(n, n)
    return Hermitian(Ψ)
end

function data_distribution(tsm::TraitSimulationModel; pss::Float64 = Inf)
    comps = varianceComponents(tsm)
    Σ = data_variance(comps)
    n = length(tsm.taxa)
    μ = vec(ones(n) * (tsm.treeModel.μ' * comps.L))

    return MvNormal(μ, Symmetric(Σ))
end

function factor_marginal_distribution(tsm::TraitSimulationModel; pss::Float64 = Inf)
    Ψ = tree_variance(tsm, pss = pss)
    Σ = tsm.treeModel.Σ
    return Hermitian(kron(Σ, Ψ))
end

function varianceComponents(tsm::TraitSimulationModel; pss::Float64 = Inf)
    Ψ = tree_variance(tsm, pss = pss)
    Σ = tsm.treeModel.Σ
    L = get_loadings(tsm)
    Λ = get_residual_variance(tsm)
    return VarianceComponents(Ψ, Σ, L, Λ)
end

function factor_variance(vc::VarianceComponents)
    return Hermitian(kron(vc.Σ, vc.Ψ))
end

function factor_data_covariance(vc::VarianceComponents)
    return kron(vc.Σ * vc.L, vc.Ψ)
end

function data_variance(vc::VarianceComponents)
    return Hermitian(kron(vc.L' * vc.Σ * vc.L, vc.Ψ) .+
        kron(vc.Λ, Diagonal(ones(size(vc.Ψ, 1)))))
end

function dimensions(tsm::TraitSimulationModel)
    n = length(tsm.taxa)
    k, p = size(get_loadings(tsm))
    dim = n * (k + p)
    f_inds = 1:(n * k)
    y_inds = 1:(n * p)
    y_inds = y_inds .+ (n * k)
    return (factor_indices = f_inds, data_indices = y_inds, dimension = dim)
end

function factor_and_data_joint_distribution(tsm::TraitSimulationModel; pss::Float64 = Inf)
    comps = varianceComponents(tsm)
    n = length(tsm.taxa)

    f_inds, y_inds, dim = dimensions(tsm)

    V = zeros(dim, dim)
    V[f_inds, f_inds] .= factor_variance(comps)
    V[f_inds, y_inds] .= factor_data_covariance(comps)
    V[y_inds, f_inds] .= (@view V[f_inds, y_inds])'
    V[y_inds, y_inds] .= data_variance(comps)

    μ_fac = ones(n) * tsm.treeModel.μ'
    μ_data = μ_fac * comps.L

    μ = vcat(vec(μ_fac), vec(μ_data))
    return MvNormal(μ, Symmetric(V))
end

function factor_conditional_distribution(tsm::TraitSimulationModel,
        data::AbstractArray{Float64}; pss::Float64 = Inf)
    @unpack μ, Σ = factor_and_data_joint_distribution(tsm, pss = pss)
    f_inds, y_inds, dim = dimensions(tsm)
    Vff = @view Σ[f_inds, f_inds]
    Vyf = @view Σ[y_inds, f_inds]
    Vyy = @view Σ[y_inds, y_inds]
    μf = @view μ[f_inds]
    μy = @view μ[y_inds]
    y = vec(data)

    Pyy = inv(Hermitian(Vyy))

    cμ = μf + Vyf' * (Pyy * (y - μy))
    cΣ = Symmetric(Vff) - Symmetric(Vyf' * Pyy * Vyf)

    return MvNormal(cμ, cΣ)
end


text:
using Documenter, StateSpaceModels

makedocs(
    modules = [StateSpaceModels],
    doctest  = false,
    clean    = true,
    format   = Documenter.HTML(),
    sitename = "StateSpaceModels.jl",
    authors = "Raphael Saavedra, Guilherme Bodin, and Mario Souto",
    pages = [
        "Home" => "index.md",
        "manual.md",
        "examples.md",
        "reference.md"
    ]
)

deploydocs(
    repo = "github.com/LAMPSPUC/StateSpaceModels.jl.git",
)


text:
using Revise
using SnpArrays
using DelimitedFiles
using Random
using Knockoffs
using LinearAlgebra
using Distributions
using DataFrames
using CSV
using ProgressMeter
using UnicodePlots

function recombination_segments(breakpoints::Vector{Int}, snps::Int)
    start = 1
    result = UnitRange{Int}[]
    for bkpt in breakpoints
        push!(result, start:bkpt)
        start = bkpt + 1
    end
    push!(result, breakpoints[end]+1:snps)
    return result
end

"""
    simulate_IBD(parent_plinkfile, offsprings)

Simulate genotypes where half the samples are siblings with the other half.
This is done by first randomly sampling a male and female from real genotype data
to represent parent. Assume they have 2 children. Then generate offspring
individuals by copying segments of the parents genotype directly to the
offspring to represent IBD segments. The number of segments is between 2 and 5
per chromosome and is chosen uniformly across the chromosome.

# References
https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1003520
"""
function simulate_IBD(parent_plinkfile::AbstractString, offsprings::Int)
    destin = "/scratch/users/bbchu/ukb/siblings/"
    iseven(offsprings) || error("number of offsprings should be even")
    # read original data
    xdata = SnpData(parent_plinkfile)
    x = xdata.snparray
    p = size(x, 2)
    male_idx = findall(x -> x == "1", xdata.person_info[!, :sex])
    female_idx = findall(x -> x == "2", xdata.person_info[!, :sex])
    # simulate new samples
    xsib = SnpArray(destin * "ukb.10k.sibs.chr10.bed", offsprings, p)
    fathers = Int[]
    mothers = Int[]
    pmeter = Progress(offsprings >> 1)
    for i in 1:(offsprings >> 1)
        # assign parents
        dad = rand(male_idx)
        mom = rand(female_idx)
        push!(fathers, dad)
        push!(mothers, mom)
        # child 1
        recombinations = rand(1:4)
        breakpoints = sort!(rand(1:p, recombinations))
        parent1, parent2 = rand() < 0.5 ? (dad, mom) : (mom, dad)
        segments = recombination_segments(breakpoints, p)
        for j in 1:length(segments)
            parent = isodd(j) ? parent1 : parent2
            segment = segments[j]
            copyto!(@view(xsib[2i - 1, segment]), @view(x[parent, segment]))
        end
        # child 2
        recombinations = rand(1:4)
        breakpoints = sort!(rand(1:p, recombinations))
        parent1, parent2 = rand() < 0.5 ? (dad, mom) : (mom, dad)
        segments = recombination_segments(breakpoints, p)
        for j in 1:length(segments)
            parent = isodd(j) ? parent1 : parent2
            segment = segments[j]
            copyto!(@view(xsib[2i, segment]), @view(x[parent, segment]))
        end
        # update progress
        next!(pmeter)
    end
    # also generate bim and fam file
    cp(parent_plinkfile * ".bim", destin * "ukb.10k.sibs.chr10.bim")
    open(destin * "ukb.10k.sibs.chr10.fam", "w") do io
        for i in 1:(offsprings >> 1)
            println(io, 2i - 1, ' ', 2i - 1, ' ', fathers[i], ' ', mothers[i], " -9")
            println(io, 2i,     ' ', 2i,     ' ', fathers[i], ' ', mothers[i], " -9")
        end
    end
    return xsib
end

Random.seed!(2021)
parent_plinkfile = "/scratch/users/bbchu/ukb/subset/ukb.10k.chr10"
offsprings = 10000
xsib = simulate_IBD(parent_plinkfile, offsprings)

# check GRM
Φ = grm(xsib, method=:Robust)
histogram(vec(Φ))