text:
export ClassField

#TODO: create an interface using characters

mutable struct BadPrime <: Exception
  p
end

function Base.show(io::IO, E::BadPrime)
  if isdefined(E, :p)
    println("Bad prime $(E.p) encountered")
  else
    println("Unknown bad prime encountered")
  end
end

mutable struct ClassField_pp{S, T}
  rayclassgroupmap::S
  quotientmap::T
  a::FacElem{nf_elem, AnticNumberField}#Generator of the Kummer Extension

  sup::Array{NfOrdIdl, 1} # the support of a - if known
  sup_known::Bool

  K::NfRel{nf_elem} # the target with the roots of unity
  A::NfRel{nf_elem} # the target
  o::Int # the degree of K - note, in general this is a divisor of the degree of A
  defect::Int # div(degree(A), degree(K)) = div(degree(A), o)
  pe::NfRelElem{nf_elem} #The image of the generator of A in K
  AutG::Vector{NfRelToNfRelMor{nf_elem, nf_elem}}
  AutR::fmpz_mat
  bigK::KummerExt
  degree::Int # The degree of the relative extension we are searching for.
              # In other words, the order of the codomain of quotientmap

  function ClassField_pp{S, T}() where {S, T}
    z = new{S, T}()
    z.degree = -1
    return z
  end
end

function Base.show(io::IO, C::ClassField_pp)
  println(IOContext(io, :compact => true), "Cyclic class field of degree $(degree(C)) defined modulo $(defining_modulus(C))")
  if isdefined(C, :a)
    println(io, "Kummer generator ", C.a)
  end
  if isdefined(C, :K)
    println(io, "Kummer extension: ", C.K)
  end
  if isdefined(C, :A)
    println(io, "Defining  polynomial ", C.A.pol)
  end
end

mutable struct ClassField{S, T}
  rayclassgroupmap::S#Union{MapRayClassGrp{GrpAbFinGen}, MapClassGrp{GrpAbFinGen}}
  quotientmap::T#GrpAbFinGenMap

  conductor::Tuple{NfOrdIdl, Array{InfPlc, 1}}
  relative_discriminant::Dict{NfOrdIdl, Int}
  absolute_discriminant::Dict{fmpz,Int}
  cyc::Array{ClassField_pp{S, T}, 1}
  A::NfRel_ns{nf_elem}
  AbsAutGrpA::Vector{NfRel_nsToNfRel_nsMor{nf_elem}} #The generators for the absolute automorphism
                                                     #group of A
  degree::Int # The degree of the relative extension we are searching for.
              # In other words, the order of the codomain of quotientmap

  function ClassField{S, T}() where {S, T}
    z = new{S, T}()
    z.degree = -1
    return z
  end
end

function Base.show(io::IO, CF::ClassField)
  print(IOContext(io, :compact => true), "Class field defined mod ",
                   defining_modulus(CF), " of structure ",
                   codomain(CF.quotientmap))
end

###############################################################################
#
#  Modulus function
#
###############################################################################

@doc Markdown.doc"""
    defining_modulus(CF::ClassField)
The modulus, ie. an ideal the the set of real places, used to create the
class field.
"""
function defining_modulus(CF::ClassField)
  return _modulus(CF.rayclassgroupmap)
end

function defining_modulus(CF::ClassField_pp)
  return _modulus(CF.rayclassgroupmap)
end

function _modulus(mq::MapRayClassGrp)
  return mq.defining_modulus::Tuple{NfOrdIdl, Vector{InfPlc}}
end

function _modulus(mq::MapClassGrp)
  return (ideal(order(codomain(mq)), 1), InfPlc[])::Tuple{NfOrdIdl, Vector{InfPlc}}
end

###############################################################################
#
#  Base ring and Base field
#
###############################################################################

@doc Markdown.doc"""
    base_ring(A::ClassField)
The maximal order of the field that $A$ is defined over.
"""
function base_ring(A::ClassField)
  return order(codomain(A.rayclassgroupmap))
end

@doc Markdown.doc"""
    base_field(A::ClassField)
The number field that $A$ is defined over.
"""
function base_field(A::ClassField)
  return number_field(base_ring(A))
end

function base_ring(A::ClassField_pp)
  return order(codomain(A.rayclassgroupmap))
end

function base_field(A::ClassField_pp)
  return number_field(base_ring(A))
end

###############################################################################
#
#  Degree
#
###############################################################################

@doc Markdown.doc"""
    degree(A::ClassField)
The degree of $A$ over its base field, ie. the size of the defining ideal group.
"""
function degree(A::ClassField)
  if A.degree == -1
    A.degree = Int(order(codomain(A.quotientmap)))
  end
  return A.degree
end

function degree(A::ClassField_pp{S, T}) where {S, T}
  if A.degree == -1
    A.degree = Int(order(codomain(A.quotientmap)))
  end
  return A.degree
end

###############################################################################
#
#  Exponent
#
###############################################################################

@doc Markdown.doc"""
    exponent(A::ClassField)
The exponent of $A$ over its base field, i.e. the exponent of the Galois
group of the extension.
"""
function exponent(A::ClassField{S, T}) where {S, T}
  return Int(exponent(codomain(A.quotientmap)))
end
###############################################################################
#
#  Compositum
#
###############################################################################

@doc Markdown.doc"""
    compositum(a::ClassField, b::ClassField) -> ClassField
> The compositum of $a$ and $b$ as a (formal) class field.
"""
function compositum(a::ClassField, b::ClassField)
  @assert base_ring(a) == base_ring(b)
  c = lcm(defining_modulus(a)[1], defining_modulus(b)[1])
  d = lcm(degree(a), degree(b))
  c_inf = union(defining_modulus(a)[2], defining_modulus(b)[2])
  r, mr = ray_class_group(c, c_inf, n_quo = Int(d))
  C = ray_class_field(mr)
  @assert domain(C.rayclassgroupmap) == r
  h = norm_group_map(C, [a,b])
  U = intersect(kernel(h[1])[1], kernel(h[2])[1])
  q, mq = quo(codomain(C.quotientmap), U)
  return ray_class_field(mr, GrpAbFinGenMap(C.quotientmap * mq))
end

@doc Markdown.doc"""
  *(A::ClassField, B::ClassField) -> ClassField
The compositum of $a$ and $b$ as a (formal) class field.
"""
*(a::ClassField, b::ClassField) = compositum(a, b)

###############################################################################
#
#  Intersection
#
###############################################################################

@doc Markdown.doc"""
    intersect(a::ClassField, b::ClassField) -> ClassField
The intersect of $a$ and $b$ as a class field.
"""
function Base.intersect(a::ClassField, b::ClassField)
  @assert base_ring(a) == base_ring(b)
  c = lcm(defining_modulus(a)[1], defining_modulus(b)[1])
  c_inf = union(defining_modulus(a)[2], defining_modulus(b)[2])
  d = lcm(degree(a), degree(b))

  r, mr = ray_class_group(c, c_inf, n_quo = Int(d))
  C = ray_class_field(mr)
  h = norm_group_map(C, [a,b])
  U = kernel(h[1])[1] + kernel(h[2])[1]
  q, mq = quo(codomain(C.quotientmap), U)
  return ray_class_field(mr, GrpAbFinGenMap(C.quotientmap * mq))
end

###############################################################################
#
#  Subfield
#
###############################################################################

@doc Markdown.doc"""
    issubfield(a::ClassField, b::ClassField) -> Bool
Determines of $a$ is a subfield of $b$.
"""
function issubfield(a::ClassField, b::ClassField)
  @assert base_ring(a) == base_ring(b)
  c = lcm(defining_modulus(a)[1], defining_modulus(b)[1])
  c_inf = union(defining_modulus(a)[2], defining_modulus(b)[2])
  d = lcm(degree(a), degree(b))

  r, mr = ray_class_group(c, c_inf, n_quo = Int(d))
  C = ray_class_field(mr)
  h = norm_group_map(C, [a,b])
  return issubset(kernel(h[2])[1], kernel(h[1])[1])
end

###############################################################################
#
#  Equality
#
###############################################################################

@doc Markdown.doc"""
    ==(a::ClassField, b::ClassField)
Tests if $a$ and $b$ are equal.
"""
function ==(a::ClassField, b::ClassField)
  @assert base_ring(a) == base_ring(b)
  mq1 = a.quotientmap
  mq2 = b.quotientmap
  if !isisomorphic(codomain(mq1), codomain(mq2))
    return false
  end
  expo = Int(exponent(codomain(mq1)))
  c = lcm(defining_modulus(a)[1], defining_modulus(b)[1])
  c_inf = union(defining_modulus(a)[2], defining_modulus(b)[2])
  r, mr = ray_class_group(c, c_inf, n_quo = expo)
  C = ray_class_field(mr)
  @assert defining_modulus(C) == (c, c_inf)
  h = norm_group_map(C, [a,b])
  return iseq(kernel(h[2])[1], kernel(h[1])[1])
end

###############################################################################
#
#  IsCyclic
#
###############################################################################

@doc Markdown.doc"""
    iscyclic(C::ClassField)
Tests if the (relative) automorphism group of $C$ is cyclic (by checking
the defining ideal group).
"""
function iscyclic(C::ClassField)
  mp = C.quotientmap
  return iscyclic(codomain(mp))
end

###############################################################################
#
#  Maximal p-subfield
#
###############################################################################

@doc Markdown.doc"""
    maximal_p_subfield(C::ClassField, p::Int)
Returns the class field corresponding to the maximal subextension of
prime power degree
"""
function maximal_p_subfield(C::ClassField, p::Int)
  mg = C.quotientmap
  v = valuation(degree(C), p)
  q, mq = quo(codomain(mg), p^v, false)
  s, ms = snf(q)
  mg1 = GrpAbFinGenMap(mg*mq*ms)
  return ray_class_field(C.rayclassgroupmap, mg1)
end


text:
export isgrid, classify, quantize, quantize!, quantize_duration!
###############################################################################
# Grid
###############################################################################
function isgrid(grid)
    issorted(grid) || throw(ArgumentError("Grids must be sorted."))
    if grid[1] != 0 || grid[end] != 1
        throw(ArgumentError("Grids must start from 0 and end in 1."))
    end
    true
end

function closest_point(grid, x)
    best = 1
    dxbest = abs(x - grid[1])
    for i in 2:length(grid)
        dx = abs(x - grid[i])
        if dx < dxbest
            dxbest = dx
            best = i
        end
    end
    return best
end

function closest_realgrid(grid, x, tpq::Integer)
    best = 1
    dxbest = abs(x - grid[1]*tpq)
    for i in 2:length(grid)
        dx = abs(x - grid[i]*tpq)
        if dx < dxbest
            dxbest = dx
            best = i
        end
    end
    return best
end

###############################################################################
# Classifiers and quantizers
###############################################################################

"""
```julia
classify(notes::Notes, grid)
classify(note::AbstractNote, grid, tpq::Integer)
```
Classify given notes according to the given grid.

Returns an integer (or vector of integers) that corresponds to the index
of the closest grid point to the note position modulo the quarter note.
`1` means start of the grid and `length(grid)` means
end of the grid (i.e. *next* quarter note).
"""
function classify(note::AbstractNote, grid, tpq::Integer)
    posmod = mod(note.position, tpq)
    return closest_realgrid(grid, posmod, tpq)
end

function classify(notes::Notes, grid)
    isgrid(grid)
    r = zeros(Int, length(notes))
    for i in 1:length(notes)
        r[i] = classify(notes[i], grid, notes.tpq)
    end
    return r
end

function quantize!(note::AbstractNote, grid, tpq::Integer)

    number_of_quarters = div(note.position, tpq)
    b = classify(note, grid, tpq)
    note.position = round(Int, (number_of_quarters*tpq + grid[b]*tpq))
    return nothing
end

function quantize!(notes::Notes, grid)

    isgrid(grid)
    for note in notes
        quantize!(note, grid, notes.tpq)
    end
    return nothing
end

"""
```julia
quantize(notes::Notes, grid, duration = true)
quantize(note::AbstractNote, grid, tpq::Integer)
```
Return a quantized copy of the given notes on the given `grid`, which can be any
sorted iterable that starts on `0` and ends on `1`.

Each note is quantized (relocated) to its closest point of the `grid`, by first
identifying that point using [`classify`](@ref).
It is assumed that the grid is the same for all quarter notes of the track.

If `duration` is `true`, the function also quantizes the duration of the notes
on the same grid.

This function respects the notes' absolute position and quantizes in absolute position,
not relative.
"""
function quantize(notes::Notes, grid, duration = true)
    qnotes = deepcopy(notes)
    quantize!(qnotes, grid)
    duration && quantize_duration!(qnotes, grid)
    return qnotes
end

function quantize_duration!(note::AbstractNote, grid, tpq)
    durmod = mod(note.position + note.duration, tpq)
    best = closest_realgrid(grid, durmod, tpq)
    # Check if quantization removes duration
    final_dur = Int(grid[best]*tpq)
    dist = final_dur - durmod
    note.duration += final_dur - durmod
    if note.duration == 0
        if best != length(grid)
            note.duration += (grid[best+1] - grid[best])*tpq
        else
            note.duration += grid[2]*tpq # grid[1] is always 0 by definition
        end
    end
    return
end

"""
    quantize_duration!(notes::Notes, grid)
    quantize_duration!(note::Note, grid, tpq)
Quantize the duration of given notes on the `grid`.
"""
function quantize_duration!(notes::Notes, grid)
    for note in notes
        quantize_duration!(note, grid, notes.tpq)
    end
    return notes
end

ext: jl


text:
# BEGIN NamedTuples 5.0.0 deprecations

if VERSION < v"0.7.0-DEV.2738"
    # Create a nicer error message for the breaking change to the
    # `@NT( a::Int64, b::Float64 )( 1, 2.0 )` syntax. Note that we will avoid creating
    # this error on Julia 0.7 as NamedTuples.jl previously did not work and won't require
    # the deprecation.

    function (::Type{NT})(args...) where {NT <: NamedTuple}
        error(string(
            "`$NT(args...)` is no longer supported as it is ambiguous with the new " *
            "`$NT((args...,))` constructor."
        ))
    end
end

# END NamedTuples 5.0.0 deprecations


text:
export Avac, vacA, vacExpVal

Avac(A::Union{BosonDestroy,FermionDestroy,σminus}) = scal(0)
Avac(A::Union{BosonCreate,FermionCreate,σplus}) = A
# vacuum is an eigenstate of σz
Avac(A::σ) = (A.a == z) ? scal(-1) : A
Avac(A::OpSum) = Avac(A.A) + Avac(A.B)
Avac(A::OpSumAnalytic) = OpSumAnalytic(A.ind,Avac(A.A))
Avac(A::OpProd) = begin
    Bv = Avac(A.B)
    Bv != A.B && return Avac(A.A*Bv)
    # if A.A and A.B commute, try what happens with Avac(A.B*A.A)
    if comm(A.A,A.B) == scal(0)
        Av = Avac(A.A)
        Av != A.A && return Avac(A.B*Av)
    end
    # did not find a way to simplify
    A
end
Avac(A::Scalar) = A

vacA(A::Union{BosonDestroy,FermionDestroy,σminus}) = A
vacA(A::Union{BosonCreate,FermionCreate,σplus}) = scal(0)

# vacuum is an eigenstate of σz
vacA(A::σ) = (A.a == z) ? scal(-1) : A
vacA(A::OpSum) = vacA(A.A) + vacA(A.B)
vacA(A::OpSumAnalytic) = OpSumAnalytic(A.ind,vacA(A.A))
vacA(A::OpProd) = begin
    if A.A isa Scalar
        A.A*vacA(A.B)
    else
        # for applying ⟨0|AB, see if ⟨0|A changes A, and if so,
        # keep going with the new operator
        vA = vacA(A.A)
        vA==A.A ?  A : vacA(vA * A.B)
    end
end
vacA(A::Scalar) = A

"""
    Avac(A::Operator), vacA(A::Operator)

Simplify operator by assuming it is applied to the vacuum from the left or
right, respectively. To be precise, `Avac(A)` returns ``A'`` such that ``A'|0⟩ =
A|0⟩``, while `vacA(A)` does the same for ``⟨0|A``."""
Avac, vacA

"""
    vacExpVal(A::Operator,S::Operator=scal(1))

Calculate the vacuum expectation value ``⟨0|S^\\dagger A S|0⟩``, i.e., the
expectation value ``⟨ψ|A|ψ⟩`` for the state defined by ``|ψ⟩= S|0⟩```.
"""
function vacExpVal(A::Operator,stateop::Operator=scal(1))
    # simplify down as much as possible by applying vacuum from left and right
    vsAsv = vacA(Avac(stateop' * A * stateop))
    # only operators that should survive here as operators are σs
    _vacExpVal(vsAsv)
end
_vacExpVal(A::Scalar) = A
# the terms that survive until here have at most a single σ of any particle
# so it does not matter if we are inside a product
_vacExpVal(A::σ) = A.a == z ? scal(-1) : scal(0)
_vacExpVal(A::OpSum) = _vacExpVal(A.A) + _vacExpVal(A.B)
_vacExpVal(A::OpSumAnalytic) = OpSumAnalytic(A.ind,_vacExpVal(A.A))
# we know that the operators here commute (all a and a† have disappeared, and at most a single σ remaining for each particle)
_vacExpVal(A::OpProd) = _vacExpVal(A.A) * _vacExpVal(A.B)


text:
# Add Max Trust Angle Constraints

using SparseArrays, LinearAlgebra

include("../constraints/constraints.jl")

function getAlpha(theta, deg = true)
    if deg
        return -tand(theta)
    end
    return -tan(theta)
end

"""
    makeMaxThrustConstraint(NSteps::Int64, nDim::Int64, thetaMax::Float64)

Creates an [`AL_Multiple_simpleCone`](@ref) constraint such that the trust
(or equivalently control) vector angle is capped by `thetaMax`.

Usage: `NSteps` states the number of timesteps in the trajectory. `nDim` states
the dimensionality of the problem (e.g. 1D, 2D, or 3D). `thetaMax` is the max
thrust angle.
"""
function makeMaxAngleConstraint(NSteps::Int64, nDim::Int64, thetaMax::Float64,
                                    deg = true)
    # sizeTraj = 3 * NSteps * nDim + 2 * nDim

    # This helper function gets the index of ux for any k
    startUInd(k) = (2 * nDim + 1) + (3 * nDim * (k - 1))

    indicatorList = Int64[startUInd(k) for k in 1:NSteps]

    sInd = Diagonal([ones(Int64, nDim - 1); 0]) # e.g. Diagonal([1; 1; 0])
    tInd = [zeros(Int64, nDim - 1); 1]          # e.g. [0; 0; 1]

    alpha = getAlpha(thetaMax, deg)

    return makeAL_Multiple_AngleCone(alpha, sInd, tInd,
                                        nDim, indicatorList)
end

function makeMaxAngleConstraint(NSteps::Int64, nDim::Int64, thetaMax::Float64,
                                    sIndicate::Diagonal{Int64,Array{Int64,1}},
                                    tIndicate::Array{Int64,1}, deg = true)
    # sizeTraj = 3 * NSteps * nDim + 2 * nDim

    # This helper function gets the index of ux for any k
    startUInd(k) = (2 * nDim + 1) + (3 * nDim * (k - 1))

    indicatorList = Int64[startUInd(k) for k in 1:NSteps]

    alpha = getAlpha(thetaMax, deg)

    return makeAL_Multiple_simpleCone(alpha, sIndicate, tIndicate,
                                        nDim, indicatorList)
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
const DEFAULT_SHRINKAGE_SHAPE = 100.0
const DEFAULT_SHRINKAGE_SCALE = 1.0
const DEFAULT_STARTING_SHAPE = 2.0

mutable struct MatrixShrinkageLikelihoods <: MyXMLElement
    mult_els::Vector{XMLOrNothing}
    gp_els::Vector{XMLOrNothing}
    local_scale_els::Vector{XMLOrNothing}
    global_scale_els::Vector{XMLOrNothing}
    bb_els::Vector{XMLOrNothing}
    ms_el::XMLOrNothing
    global_prior_els::Vector{XMLOrNothing}
    loadings::MatrixParameter
    shapes::Vector{Float64}
    scales::Vector{Float64}
    local_scales::Matrix{Float64}

    function MatrixShrinkageLikelihoods(loadings::MatrixParameter)

        k, p = size(loadings.mat)
        shapes = fill(DEFAULT_SHRINKAGE_SHAPE, k)
        shapes[1] = DEFAULT_STARTING_SHAPE
        scales = fill(DEFAULT_SHRINKAGE_SCALE, k)

        return new(xml_vec(k - 1),
                    xml_vec(k),
                    xml_vec(k),
                    xml_vec(k),
                    xml_vec(k),
                    nothing,
                    xml_vec(k),
                    loadings,
                    shapes,
                    scales,
                    ones(k, p)
                    )
    end
end

function get_fac_dim(msl::MatrixShrinkageLikelihoods)
    return length(msl.gp_els)
end

function make_xml(msl::MatrixShrinkageLikelihoods)

    loadings = make_xml(msl.loadings)
    load_rows = loadings[bn.PARAMETER]
    k = length(load_rows)
    p = length(split(attribute(load_rows[1], bn.VALUE)))

    l_ids = Vector{String}(undef, k)
    for i = 1:k
        l_ids[i] = attribute(load_rows[i], bn.ID)
    end

    for i = 2:k
        msl.mult_els[i - 1] =
                make_parameter(id="rowMult$i",
                                value=[msl.shapes[i] * msl.scales[i]],
                                lower=0.0)
    end

    for i = 1:k
        msl.local_scale_els[i] = make_parameter(id="localScale$i",
                                    value=vec(msl.local_scales[i, :]),
                                    lower=0.0)
    end

    msl.gp_els[1] = make_parameter(id="globalPrecision1",
                                    value=[msl.shapes[1] * msl.scales[1]],
                                    lower=0.0)

    for i = 2:k
        prod_param = new_element(bn.PRODUCT_PARAMETER)
        set_id!(prod_param, "globalPrecision$i")
        add_ref_el(prod_param, msl.gp_els[i - 1])
        add_ref_el(prod_param, msl.mult_els[i - 1])
        msl.gp_els[i] = prod_param
    end

    for i = 1:k
        trans_param = new_element(bn.TRANSFORMED_PARAMETER)
        set_id!(trans_param, "globalScale$i")
        set_attribute(trans_param, bn.TYPE, bn.POWER)
        add_ref_el(trans_param, msl.gp_els[i])
        pt_el = new_child(trans_param, bn.POWER_TRANSFORM)
        set_attribute(pt_el, bn.POWER, "-0.5")
        t_el = new_child(pt_el, bn.TRANSFORM)
        set_attribute(t_el, bn.TYPE, bn.POWER)
        msl.global_scale_els[i] = trans_param
    end

    for i = 1:k
        bb_el = new_element(bn.BAYESIAN_BRIDGE)
        set_id!(bb_el, "bb$i")
        add_ref_el(bb_el, load_rows[i])
        gs_el = new_child(bb_el, bn.GLOBAL_SCALE)
        add_ref_el(gs_el, msl.global_scale_els[i])
        ls_el = new_child(bb_el, bn.LOCAL_SCALE)
        add_ref_el(ls_el, msl.local_scale_els[i])
        exp_el = new_child(bb_el, bn.EXPONENT)
        add_parameter(exp_el, value=[0.5])
        msl.bb_els[i] = bb_el
    end

    ms_el = new_element(bn.MATRIX_SHRINKAGE_LIKELIHOOD)
    set_id!(ms_el, "$(get_id(loadings)).prior")
    add_ref_el(ms_el, loadings)
    rows_el = new_child(ms_el, bn.ROW_PRIORS)
    for i = 1:k
        add_ref_el(rows_el, msl.bb_els[i])
    end
    msl.ms_el = ms_el

    for i = 1:k

        gp_prior = new_element(bn.GAMMA_PRIOR)
        if i == 1
            set_id!(gp_prior, "globalPrecisionPrior1")
        else
            set_id!(gp_prior, "rowMultPrior$i")
        end
        set_attribute(gp_prior, bn.SHAPE, msl.shapes[i])
        set_attribute(gp_prior, bn.SCALE, msl.scales[i])
        if i == 1
            add_ref_el(gp_prior, msl.gp_els[1])
        else
            add_ref_el(gp_prior, msl.mult_els[i - 1])
        end

        msl.global_prior_els[i] = gp_prior
    end

    all_els = [msl.mult_els; msl.gp_els; msl.local_scale_els;
                msl.global_scale_els; msl.bb_els; msl.ms_el; msl.global_prior_els]

    return all_els

end

function get_priors(xml::MatrixShrinkageLikelihoods)
    return XMLElement[xml.ms_el; xml.global_prior_els]
end

function get_loggables(xml::MatrixShrinkageLikelihoods)
    return [xml.global_scale_els; xml.mult_els; xml.local_scale_els]
end

function get_precision_prior(xml::MatrixShrinkageLikelihoods, ind::Int; make_xml::Bool = true)
    if make_xml
        make_xml(xml)
    end
    return xml.global_prior_els[ind]
end

function make_loadings_gradient(::MyXMLElement, msl::MatrixShrinkageLikelihoods)
    make_xml(msl)
    return msl.ms_el
end

function get_normal_prior(msl::MatrixShrinkageLikelihoods)
    make_xml(msl)
    return msl.ms_el
end


text:
#!/usr/bin/env julia

using Luxor, Colors

function test_circular_arrows_1(pos)
    gsave()
    froma = rescale(rand(1:100), 1, 100, 0, 2pi)
    toa =   rescale(rand(1:100), 1, 100, 0, 2pi)
    sethue("black")
    arrow(pos, 100, froma, toa, linewidth=rand(1:6), arrowheadlength=rand(10:30))
    text(string("from: ", round(rad2deg(froma), 1)), pos)
    text(string("to: ", round(rad2deg(toa), 1)), pos.x, pos.y+10)
    sethue("magenta")
    arrow(pos, 100, toa, froma; linewidth=rand(1:6), arrowheadlength=rand(10:30))
    text(string("from: ", round(rad2deg(toa), 1)), pos.x, pos.y+20)
    text(string("to: ", round(rad2deg(froma), 1)), pos.x, pos.y+30)
    grestore()
end

function test_circular_arrows_2(pos, w)
  sethue("black")
  for a in 50:10:w
    randomhue()
    starta = rand(0:pi/12:2pi)
    finisha = rand(0:pi/12:2pi)
    arrow(pos, a, starta, finisha, linewidth=rand(1:6), arrowheadlength=rand(1:20))
  end
end

function arrow_test(fname)
  pagewidth, pageheight = 2000, 2000
  Drawing(pagewidth, pageheight, fname)
  origin() # move 0/0 to center
  background("ivory")
  setopacity(0.5)
  setline(2)
  pagetiles = Tiler(pagewidth, pageheight, 4, 4, margin=50)
  for (pos, n) in pagetiles
    if isodd(n)
      test_circular_arrows_1(pos)
    else
      test_circular_arrows_2(pos, pagetiles.tilewidth/2)
    end
  end
  finish()
  println("arrow-test saved in $(fname)")
end

srand(42)
fname = "/tmp/arrow-arctest.pdf"
arrow_test(fname)


text:
export merge_graphs;

# Change a symbol name using prefix & postfix
function translate_keyname(key, prefix, postfix, skip_basic, input)
    newkey = Symbol(prefix * string(key) * postfix)
    if (skip_basic & (key == :I || key == input))
        newkey = key
    end
    return newkey
end

function change_keywords!(
    graph;
    prefix = "G",
    postfix = "",
    skip_basic = true,
    input = :A,
)
    translate =
        key -> translate_keyname(key, prefix, postfix, skip_basic, input)

    # Go through all dicts and update symbols according to translate
    for dict_org in (graph.coeffs, graph.operations, graph.parents)
        dict = deepcopy(dict_org)
        for (key, value) in dict
            delete!(dict_org, key)
            new_key = translate(key)
            dict_org[new_key] = value
            if (dict_org == graph.parents) # Also translate parents
                val = [value...] # Turn into vector
                for s = 1:2
                    val[s] = translate(val[s])
                end
                dict_org[new_key] = Tuple(val)
            end
        end
    end
    # Update outputs
    return graph.outputs[:] = translate.(graph.outputs)
end

"""
    graph=merge_graphs(
    graph1,
    graph2;
    prefix1 = "",
    prefix2 = "G2",
    skip_basic1 = true,
    skip_basic2 = true,
    cref1 = Vector(),
    cref2 = Vector(),
    input1 = :A,
    input2 = :A,

)

Takes all the nodes and edges in `graph1` and `graph2` and generates a new
graph. The node names are in `graph1` are changed by adding a prefix `prefix1`
and `graph2` correspondingly. The nodes `:I` and `:A` are unchanged if
`skip_basicX=true`. All coefficient references `crefX` are modified accordingly.
"""
function merge_graphs(
    graph1,
    graph2;
    prefix1 = "",
    prefix2 = "G2",
    skip_basic1 = true,
    skip_basic2 = true,
    cref1 = Vector(),
    cref2 = Vector(),
    input1 = :A,
    input2 = :A,
)
    T1 = eltype(graph1)
    T2 = eltype(graph2)

    T = promote_type(T1, T2)

    g1 = deepcopy(graph1)
    change_keywords!(
        g1,
        prefix = prefix1,
        skip_basic = skip_basic1,
        input = input1,
    )
    g2 = deepcopy(graph2)
    change_keywords!(
        g2,
        prefix = prefix2,
        skip_basic = skip_basic2,
        input = input2,
    )

    operations = merge(g1.operations, g2.operations)
    parents = merge(g1.parents, g2.parents)

    # Type logic for the coefficients
    coeffs = Dict{Symbol,Tuple{T,T}}()
    for dict in (g1.coeffs, g2.coeffs)
        for (key, value) in dict
            val = (convert(T, value[1]), convert(T, value[2]))
            coeffs[key] = (convert(T, value[1]), convert(T, value[2]))
        end
    end
    outputs = vcat(g1.outputs, g2.outputs)

    cref1_copy = copy(cref1)
    cref2_copy = copy(cref2)

    empty!(cref1)
    empty!(cref2)

    for c in cref1_copy
        push!(cref1, (Symbol(prefix1 * String(c[1])), c[2]))
    end
    for c in cref2_copy
        push!(cref2, (Symbol(prefix2 * String(c[1])), c[2]))
    end

    return Compgraph(operations, parents, coeffs, outputs)
end


text:
using DashBootstrapComponents

progress = dbc_progress("25%", value=25);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
const IMGUI_DEMO_COLOR_EDIT_STATES = Dict(
    "basic_color_edit_no_alpha" =>
        ColorEdits.State("color 1", Cfloat[1.0, 0.0, 0.2, 0.0]),
    "basic_color_edit_with_alpha" =>
        ColorEdits.State("color 2", Cfloat[0.4, 0.7, 0.0, 0.5]),
)


text:
function fakedata(name::AbstractString, args...)

    # gaussians with separation
    name == "gaussian" && return fakedata_gaussian(args...)

    # no dataset found
    error("dataset name not recognized.")
end

function fakedata_gaussian(
        pts_per_clust::Int=10,
        nclust::Int=2,
        xmin=0.,
        xmax=nclust*7.,
        nx=round(Int,(xmax-xmin)*10),
        σ=1.0,
        amin=1.,
        amax=2.
    )
    x = linspace(xmin,xmax,nx)
    μ = linspace(xmin,xmax,nclust+2)[2:end-1]
    a = linspace(xmin,xmax,pts_per_clust+2)[2:end-1]

    npts = pts_per_clust*nclust
    data = zeros(nx,npts)
    labels = zeros(Int,npts)

    i = 1
    for c = 1:nclust
        for n = 1:pts_per_clust
            data[:,i] = gauss_func(a[n],μ[c],σ,x)
            labels[i] = c
            i += 1
        end
    end

    return data, labels
end

gauss_func(a,μ,σ,x) = a*exp( (-(x-μ).^2) / (2σ^2) )


text:
mutable struct Counter
    value
    count::Int

    Counter(value) = new(value, 1)
end

function str(counter::Counter)
    (counter.count > 1 ? string(counter.count) : "") * counter.value
end

function encode(s)
    if s == ""
        return ""
    end

    encoded = [Counter(s[1])]
    for ch in s[2:end]
        if ch == encoded[end].value
            encoded[end].count += 1
        else
            push!(encoded, Counter(ch))
        end
    end
    return join(map(str, encoded))
end

function decode(s)
    out = ""
    n = 0
    for ch in s
        maybe_number = Int(ch) - 48
        if 0 <= maybe_number <= 9
            n *= 10
            n += maybe_number
        else
            out *= repeat(ch, max(n, 1))
            n = 0
        end
    end
    return out
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
function read_data(data_folder::String = "./data_input/")
	"""
	Read data from https://open.toronto.ca/dataset/ttc-routes-and-schedules/.
	Based on these files TTC routes and graph will be created.
	"""
    trips =  DataFrame(CSV.File(data_folder * "trips.txt"));
    routes = DataFrame(CSV.File(data_folder * "routes.txt"));
    stop_times = DataFrame(CSV.File(data_folder * "stop_times.txt"));
    shapes = DataFrame(CSV.File(data_folder * "shapes.txt"));
    stops = DataFrame(CSV.File(data_folder * "stops.txt"));
    return trips, routes, stop_times, shapes, stops
end

function get_longest_trip_data(stop_times::DataFrame, trips_filtered::DataFrame, stops::DataFrame)
    """
    Description:
        It takes the data and then find the longest trip per each route and leave only this trip in the dataset.
    Returns:
        Dataframe, where each route is associated with the longest trip only.
    """
    # Find stops for each trip
    routes_stops = innerjoin(stop_times, trips_filtered, on=:trip_id)[:,[:route_id,:trip_id,:stop_id,:stop_sequence,:arrival_time,:departure_time]]
    # Find the longest trip for each route
    routes_max_trip = combine(find_max_index -> find_max_index[argmax(find_max_index[!,"stop_sequence"]), [:trip_id]], groupby(routes_stops, :route_id))
    # Leave only those trips, which are the longest one.
    routes_stops_longest = innerjoin(routes_stops, routes_max_trip, on=:trip_id, makeunique=true)
    # Join stop LLA info
    route_stops_geo = innerjoin(routes_stops_longest, stops, on = :stop_id)[:,[:route_id,:stop_id,:stop_sequence,:stop_lat,:stop_lon,:arrival_time,:departure_time]]

    # Direction A->B
    route_stops_geo_dir_A = deepcopy(route_stops_geo)
    route_stops_geo_dir_A[!,"route_id"] = Symbol.(string.(route_stops_geo_dir_A[!, "route_id"])  .* "_a")

    # Direction B->A
    route_stops_geo_dir_B = deepcopy(route_stops_geo)
    new_stop_sequence = Int64[]
    for r in unique(route_stops_geo_dir_B[!,"route_id"])
        stops_reversed = sort(route_stops_geo_dir_B[route_stops_geo_dir_B.route_id .== r, :stop_sequence],rev=true)
        for i in stops_reversed
           push!(new_stop_sequence,i)
        end
    end
    route_stops_geo_dir_B[!,"stop_sequence"] = new_stop_sequence
    route_stops_geo_dir_B[!,"route_id"] = Symbol.(string.(route_stops_geo_dir_B[!, "route_id"]) .* "_b")

    # UNION
    both_direction_routes = vcat(route_stops_geo_dir_A,route_stops_geo_dir_B)

    # Global id for each stop, to ensure that even two identical points will have different ids.
    both_direction_routes[!,"my_stop_id"] = [i for i in 1:nrow(both_direction_routes)]

    return both_direction_routes
end

function get_routes_dict(data::DataFrame)
    """
    Description:
        Creates structure to represent route as a dictionary of stops, which are located on the map.
    Returns:
        1) Dict(route_id => Dict(stop_id => (lat,lon), stop_id => (lat,lon), ...))
        2) Number_of_stops::Int64
    """
    # Initialize final dictionary
    routes_d = OrderedDict{Symbol, OrderedDict{Int,Tuple{Float64,Float64}}}()
    # For each shape we find stops in sequence and give each stop unique id.
    # Once we have all stops gathered we can assign them to a particular route.
    for r in unique(data[!,"route_id"])
        stops_selected = sort(data[data.route_id .== r, :], :stop_sequence)
        route_points = OrderedDict{Int,Tuple{Float64,Float64}}()
        for srow in eachrow(stops_selected)
            geo_stop_point = (srow.stop_lat, srow.stop_lon)
            route_points[srow.my_stop_id] = geo_stop_point
        end
        routes_d[r] = route_points
    end

    routes_d, maximum(data[!,"my_stop_id"])
end

function convert_my_date(str_time::String)
    """
    Description:
        Takes string time like 22:55:33 and calculate number of seconds from a day start.
        This function is used instead of date - date, because there are times like 26:58:03 in the dataset.
    Returns:
        Number of seconds from a day start.
    """
    hours = parse(Int64,split(str_time,":")[1])
    minutes = parse(Int64,split(str_time,":")[2])
    seconds = parse(Int64,split(str_time,":")[3])
    sum_sec = hours*60*60 + minutes*60 + seconds
    return sum_sec
end

function trips_time(data::DataFrame)
    """
    Description:
        Takes arrival time and departure time for stops in the data and calculate number of seconds between two stops on the route.
    Returns:
        Dictionary, where keys are (stop_origin, stop_destination) and values are seconds between stops.
    """
    mydata = deepcopy(data)
    mydata[!, "arrival_time"] = map(x -> convert_my_date(x),mydata[!,"arrival_time"])
    mydata[!, "departure_time"] = map(x -> convert_my_date(x),mydata[!,"departure_time"])
    stops_time = Dict{Tuple{Int,Int},Int64}()
    unq_routes = unique(mydata[!,"route_id"])
    for route in unq_routes
       per_route = mydata[mydata[!,"route_id"].==route,:]
       for i in 2:nrow(per_route)
            stop_a, dep_a = per_route[i-1,:my_stop_id], per_route[i-1,:departure_time]
            stop_b, arr_b = per_route[i,:my_stop_id], per_route[i,:arrival_time]
            stops_time[(stop_a,stop_b)] = (arr_b - dep_a)
            stops_time[(stop_b,stop_a)] = (arr_b - dep_a)
       end
    end
    return stops_time
end

function assign_route_type_to_route(routes_dict::OrderedDict{Symbol, OrderedDict{Int,Tuple{Float64,Float64}}}, routes_filtered::DataFrame)
	"""
	Description:
		Assigns route type for each node on the route.
	Returns:
		Route type of each node.
	"""
    route_types = Dict{Int,Int}()
    for route in keys(routes_dict)

        if length(routes_filtered[routes_filtered[!,"route_id"] .== parse(Int64,split(string(route),"_")[1]),:][!,"route_type"]) > 1
            println("More than 1 route type")
            break
        end

        temp = Dict(zip(keys(routes_dict[route]),
                repeat(routes_filtered[routes_filtered[!,"route_id"] .== parse(Int64,split(string(route),"_")[1]),:][!,"route_type"],
                    length(keys(routes_dict[route])))))
        merge!(route_types, temp)
    end # each node has a type of route attached to it - may be useful in simulations
    return route_types
end

function routes_types_both_dir_func(routes_filtered::DataFrame)
	route_both_dir_type = Dict{Symbol,Symbol}()
	for i in eachrow(routes_filtered)
		for dire in ["a","b"]
			dir_key = string(i.route_id,"_",dire)
			route_both_dir_type[Symbol(dir_key)] = i.route_type == 1 ? :subway : :streetcar
		end
	end
	return route_both_dir_type
end

function read_and_prep_TTC_data(;data_folder::String = "./data_input/")
	"""
	Description:
		Gather all functions which should be used to prepare TTC data into a single function.
	Returns:
		TTC graph, number of seconds between stops, Dictionary with each route types and associated stops, and the DataFrame with filtered routes to be used in TTC routes visualization.
	"""
    trips, routes, stop_times, shapes, stops = read_data(data_folder);
    route_ids_without_buses = routes[coalesce.(routes.route_type .!= 3, false), :][!,"route_id"];
    routes_filtered = routes[in.(routes[!,"route_id"], [route_ids_without_buses]), :];
    trips_filtered = trips[in.(trips[!,"route_id"], [route_ids_without_buses]), :];
    trip_longest_df = get_longest_trip_data(stop_times, trips_filtered, stops);
    routes_dict, n_stops = get_routes_dict(trip_longest_df);
    secs_between_stops = trips_time(trip_longest_df);
    route_types = assign_route_type_to_route(routes_dict, routes_filtered);
	routes_types_both_dir = routes_types_both_dir_func(routes_filtered);
    gg = create_graph_from_TTC(routes_dict,n_stops,secs_between_stops);
    println("TTC data was read")
    return gg,secs_between_stops,routes_dict,route_types,routes_filtered,routes_types_both_dir
end
println("Read and prep TTC functions were read")

max_stars_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_issues_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_forks_repo_name: NykPol/EpidemicInUrbanNetworkToronto


text:
######################
# Logit and Logistic #
######################
using StatsFuns: logit, logistic

struct Logit{T<:Real} <: Bijector{0}
    a::T
    b::T
end

(b::Logit)(x::Real) = logit((x - b.a) / (b.b - b.a))
(b::Logit)(x) = mapvcat(b, x)

(ib::Inverse{<:Logit})(y::Real) = (ib.orig.b - ib.orig.a) * logistic(y) + ib.orig.a
(ib::Inverse{<:Logit})(y) = mapvcat(ib, y)

logabsdetjac(b::Logit, x::Real) = -log((x - b.a) * (b.b - x) / (b.b - b.a))
logabsdetjac(b::Logit, x) = mapvcat(x) do x
    logabsdetjac(b, x)
end


text:
// review  https://replit.com/languages/julia  there are issues in Test2.jl and Test3.jl

//Test 1 in Julia
function a()
	x = 1
		print(x)

end

max_stars_repo_name: tylerholmes/Interpreter

max_issues_repo_name: tylerholmes/Interpreter

max_forks_repo_name: tylerholmes/Interpreter


text:
module SimpleNeRF

using Random

using ArgCheck
using Distributions
using Flux
using Functors
using Optimisers
using Configurations
using ChainRulesCore
using ChainRulesCore: NoTangent
using MLUtils: stack, unsqueeze, unbatch
# using PyCall
# using PyCallChainRules.Jax: JaxFunctionWrapper

# const jax = PyNULL()

# function __init__()
#   copy!(jax, pyimport("jax"))
# end

# # From https://github.com/FluxML/Functors.jl/issues/35
# # and https://github.com/FluxML/Optimisers.jl/pull/57
# const INIT = Base._InitialValue();
# function total(f, x; init = INIT)
#     fmap(x; exclude = Optimisers.isnumeric, walk = (f, z) -> foreach(f, Optimisers._trainable(z))) do y
#       val = f(y)
#       init = init===INIT ? val : (init+val)
#     end
#     init
# end

# function fmapreduce(f, op, x; init = INIT, walk = (f, x) -> foreach(f, Functors.children(x)), kw...)
#     fmap(x; walk, kw...) do y
#       init = init===INIT ? f(y) : op(init, f(y))
#     end
#     init===INIT ? Base.mapreduce_empty(f, op) : init
# end

include("repeat.jl")
include("interpolations.jl")

include("data.jl")
include("model.jl")
include("render.jl")

include("train.jl")

end


text:
function mk_score(x)
    s = 0
    n = length(x)

    for i = 1:n-1
        for j = i : n
            s = s + sign(x[j] - x[i])
        end
    end

    return s
end

function variance_s(x, n)
    # calculate the unique data
    unique_x = unique(x)
    g = length(unique_x)

    # calculate the var(s)
    if n == g            # there is no tie
        var_s = (n*(n-1)*(2*n+5))/18

    else                 # there are some ties in data
        tp = zeros(length(unique_x))
        demo = ones(n)

        for i in 1:g
            tp[i] = sum(demo[x == unique_x[i]])
        end
        var_s = (n*(n-1)*(2*n+5) - sum(tp*(tp-1)*(2*tp+5)))/18
    end

    return var_s
end

function z_score(s, var_s)
    if s > 0
        z = (s - 1)/sqrt(var_s)
    elseif s == 0
        z = 0
    elseif s < 0
        z = (s + 1)/sqrt(var_s)
    end

    return z
end

function p_value(z, α)
    # two tail test
    p = 2*(1-normcdf(abs(z)))
    h = abs(z) > norminvcdf(1-α/2)

    if (z < 0) && h
        trend = "decreasing"
    elseif (z > 0) && h
        trend = "increasing"
    else
        trend = "no trend"
    end

    return p, h, trend
end


text:
using ModelingToolkit, Test
@parameters t a
@variables x y
@derivatives D'~t
eqs = [D(x) ~ a*x - x*y,
       D(y) ~ -3y + x*y]
f = build_function(eqs,[x,y],[a],t,expression=Val{false},target=ModelingToolkit.CTarget())
f2 = eval(build_function([x.rhs for x in eqs],[x,y],[a],t)[2])
du = rand(2); du2 = rand(2)
u = rand(2)
p = rand(1)
_t = rand()
f(du,u,p,_t)
f2(du2,u,p,_t)
@test du == du2


text:
using Test
using VsRecorderBase

@testset "VsRecorderBase.jl" begin
    @test true
end


text:
include("learning_rate/alrsvi.jl")
include("learning_rate/inversedecay.jl")
include("vi_optimizers.jl")
include("analytic.jl")
include("analyticVI.jl")
include("numericalVI.jl")
include("sampling.jl")

function post_process!(model::AbstractGP{T,<:Likelihood,<:Inference}) where {T}
    nothing
end

Base.length(::Inference) = 1

Base.iterate(l::Inference) = (l,nothing)
Base.iterate(l::Inference, ::Any) = nothing

isstochastic(l::Inference) = l.Stochastic

const GibbsorVI = Union{<:GibbsSampling,<:AnalyticVI}

#Conversion from natural to standard distribution parameters
function global_update!(gp::Abstract_GP) where {T,L}
    gp.Σ .= -0.5*inv(gp.η₂)
    gp.μ .= gp.Σ*gp.η₁
end

## Default function for getting a view on y
@inline view_y(l::Likelihood,y::AbstractVector,i::AbstractVector) = view(y,i)

## Default function for getting gradient ##
function grad_logpdf(l::Likelihood,y::Real,f::Real)
    ForwardDiff.gradient(x->logpdf(l,y,x),[f])[1]
end

function grad_logpdf(l::Likelihood,y::Real,f::AbstractVector)
    ForwardDiff.gradient(x->logpdf(l,y,x),f)
end


text:
using SearchLight

SearchLight.Configuration.load() |> SearchLight.connect


text:
module Integrators

export rk4

using LinearAlgebra

function rk4(f::Function, y0::Array{Float64, 1}, t0::Float64,
             t1::Float64, h::Float64; inplace::Bool=true)
    y = y0
    n = round(Int, (t1 - t0)/h)
    t = t0
    if ~inplace
        hist = zeros(n, length(y0))
    end
    for i in 1:n
        k1 = h * f(t, y)
        k2 = h * f(t + 0.5*h, y + 0.5*k1)
        k3 = h * f(t + 0.5*h, y + 0.5*k2)
        k4 = h * f(t + h, y + k3)
        y = y + (k1 + 2*k2 + 2*k3 + k4)/6
        if ~inplace
            hist[i, :] = y
        end
        t = t0 + i*h
    end
    if ~inplace
        return hist
    else
        return y
    end
end

function rk4_prop(jac::Function, y0::Array{Float64, 1}, t0::Float64,
                  t1::Float64, h::Float64)
    y = y0
    n = round(Int, (t1 - t0)/h)
    t = t0
    A_acc = I
    for i in 1:n
        J = jac(t, y)
        A = I + h*J + h^2/factorial(2)*J^2 + h^3/factorial(3)*J^3 + h^4/factorial(4)*J^4
        A_acc = A*A_acc
        y = A*y
        t = t0 + i*h
    end
    return A_acc
end

end

max_stars_repo_name: eviatarbach/mmda

max_issues_repo_name: eviatarbach/mmda

max_forks_repo_name: eviatarbach/mmda


text:
InitialValues.@def push!! [x]

@static if isdefined(InitialValues, Symbol("@def_monoid"))
    InitialValues.@def_monoid append!!
else
    InitialValues.@def append!!
end


text:
#!/usr/bin/julia

# Copyright (C) 2022 Matheus Fernandes Bigolin <mfrdrbigolin@disroot.org>
# SPDX-License-Identifier: MIT

"Day Twenty-Five, Sea Cucumber."

include("Utils.jl")

using .Utils

"Find the first step in which no sea `cucumbers` move."
function solve(cucumbers :: Matrix{Complex{Int}})
    old = copy(cucumbers)
    new = copy(cucumbers)

    sz = size(cucumbers)

    step = 0
    changed = true
    while changed
        step += 1
        changed = false

        for dir in [(1, 0), (0, 1)]
            for ind in CartesianIndices(new)
                if old[ind] != complex(dir...)
                    continue
                end

                moveto = CartesianIndex(@. mod1((ind.I + dir), sz))

                if old[moveto] == 0
                    new[ind] = 0
                    new[moveto] = complex(dir...)
                    changed = true
                end
            end

            old = copy(new)
        end
    end

    return step
end

function main()
    usage_and_exit(length(ARGS) != 1)

    data = readlines(first(ARGS))
    numdata = [[[0, 1, im][findfirst(s, ".>v")] for s in r] for r in data]
    matdata = hcat(numdata...)

    println(solve(matdata))
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
function expandgrid(levels...)
    lengths = length.(levels)
    inner = 1
    outer = prod(lengths)
    grid = []
    for i in 1:length(levels)
        outer = div(outer, lengths[i])
        push!(grid, repeat(levels[i], inner=inner, outer=outer))
        inner *= lengths[i]
    end
    Tuple(grid)
end


text:
using DataFrames
using Gadfly

function savetableimg(A::BenchmarkTable, title::String, filename::String,
                        relproc::Int, indp = 1:size(A, 2))
    #require("DataFrames")
    #require("Gadfly")
    x = A.cfgs
    nn = size(A, 1)
    np = size(A, 2)
    dfa = Array(DataFrames.DataFrame, length(indp))

    for i = 1:length(indp)
        proc = indp[i]
        y = Array(Float64, nn)
        for j = 1:nn
            rel = get(A[j, relproc], BenchmarkLite.Gps())
            y[j] = rel/get(A[j, proc], BenchmarkLite.Gps())
        end
        dfa[i] = DataFrames.DataFrame(x=x, y=y, method=string(rtable.procs[proc]))
    end
    df = vcat(dfa...)

    p = Gadfly.plot(df, x="x", y="y", color="method",
            Scale.x_log10, Scale.y_log10, Theme(line_width=0.7mm, default_point_size=0.5mm), Geom.line, Geom.point,
            Guide.XLabel("Array length [n]"), Guide.YLabel("Rel. [s / G ops.] to FAO"),
            Guide.title(title))
    Gadfly.draw(SVG(string(filename,".svg"), 16cm, 10cm), p)
    Gadfly.draw(PNG(string(filename,".png"), 16cm, 10cm), p)
    return nothing
end


text:
using BCTRNN
using DiffEqSensitivity
using OrdinaryDiffEq
import DiffEqFlux: FastChain, FastDense
import Flux: ClipValue, ADAM

# Not in Project.toml
using Plots
gr()

include("sine_wave_dataloader.jl")

function train_sine_ss_fc(epochs, solver=Tsit5();
  sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)),
  T=Float32, model_size=5,
  mtkize=false, gen_jac=false, lr=0.02, kwargs...)

  train_dl = generate_2d_data(T)

  f_in = 2
  f_out = 1

  im = BCTRNN.InputAllToAll()
  sm = BCTRNN.SynsFullyConnected()
  #om = BCTRNN.OutputAll()
  om = BCTRNN.OutputIdxs(collect(1:model_size))
  wiring = BCTRNN.WiringConfig(f_in, model_size, im,sm,om)

  model = FastChain(BCTRNN.Mapper(f_in),
                    BCTRNN.LTCSynState(wiring, solver, sensealg; T, mtkize, gen_jac, kwargs...),
                    FastDense(wiring.n_out, f_out))

  hs = []
  for (k,v) in wiring.matrices
    push!(hs, heatmap(v, title=k))
  end
  display(plot(hs..., layout=length(hs)))

  cb = BCTRNN.MyCallback(T; ecb=mycb, nepochs=epochs, nsamples=length(train_dl))
  #opt = GalacticOptim.Flux.Optimiser(ClipValue(0.5), ADAM(0.02))
  opt = BCTRNN.ClampBoundOptim(BCTRNN.get_bounds(model,T)..., ClipValue(T(1.0)), ADAM(T(lr)))
  BCTRNN.optimize(model, BCTRNN.loss_seq, cb, opt, train_dl, epochs, T), model
end

@time train_sine_ss_fc(1000, AutoTsit5(Rosenbrock23(autodiff=false)); model_size=5, mtkize=true, lr=0.004)


text:
export project_l2!

function project_l2!(
                     x     ::Union{Vector{TF},Vector{Complex{TF}}},
                     sigma ::TF          #maximum 2-norm ||x||_2
                     ) where {TF<:Real}

nl2 = norm(x,2)
if nl2 <= sigma
  return x
else
  rmul!(x,sigma/nl2)
end

return x
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
#!/usr/bin/env julia

using StructuralDDP
using Test
mytol = 1e-4

# load models
include("models/CapitalAdjustModel.jl")
include("models/CapitalAdjustModel2.jl")
createmodel(model::Symbol; kwargs...) = eval(model)(; kwargs...)

include("test_transitionmatrix.jl")
include("test_rewardmatrix_options.jl")
include("test_intdim.jl")
include("test_monotonicity_concavity.jl")
include("test_simulate.jl")
# include("test_speed.jl")


text:
include("test_utils.jl")
include("test_softmax.jl")
include("test_embedding.jl")
#include("test_parallel.jl")


text:
struct EliminationTree
    par::Vector{Int}
    child::Vector{Vector{Int}}
    peo::Vector{Int}
    function EliminationTree(L::SparseMatrixCSC; peo::Vector{Int})
        n = size(L, 1)
        if !is_chordal(L; peo=1:n)
            error(ArgumentError("L must be chordal"))
        end
        par = get_etree(L)
        child = get_children_from_par(par)
        new(par, child, peo)
    end
end

struct CliqueTree
    etree::EliminationTree
    vreps::Vector{Int}
    snd_par::Vector{Int}
    snd_child::Vector{Vector{Int}}
    snds::Vector{Vector{Int}}
    seps::Vector{Vector{Int}}
    postordering::Vector{Int}
    perm::Vector{Int}
    function CliqueTree(L::SparseMatrixCSC; peo=nothing)
        n = size(L, 1)
        if isnothing(peo)
            peo = collect(1:n)
        end
        etree = EliminationTree(L; peo=peo)
        vreps, snd_par, snd_membership = max_supernode_etree(L, etree.par)
        snd_children = get_children_from_par(snd_par)
        postordering = get_postordering(snd_par, snd_children)

        n_snds = length(vreps)
        snds = [Vector{Int}(undef, 0) for _ in 1:n_snds]
        for v in 1:n
            push!(snds[snd_membership[v]], v)
        end

        seps = [Vector{Int}(undef, 0) for _ in 1:n_snds]
        for (ind, vrep) in enumerate(vreps)
            seps[ind] = filter(x->!(x in snds[ind]), rowvals(L)[nzrange(L, vrep)])
        end

        new(etree, vreps, snd_par, snd_children, snds, seps, postordering, collect(1:n))
    end
end

# construct a permutation as in VA section 4.6
# st the supernodes have consecutive vertices
"""
    order_snds!(ct::CliqueTree)

Orders supernodes such that elements of each supernode are numbered consecutively
and the order is a topological ordering of the reprsentative vertices in the
supernodal elimination tree. (See [VA15, 4.6])

# Reference
[VA15] [Chordal Graphs and Semidefinite Optimization](https://www.seas.ucla.edu/~vandenbe/publications/chordalsdp.pdf)
by Lieven Vandenberghe and Martin Andersen
"""
function order_snds!(ct::CliqueTree)
    n = length(ct.perm)
    perm = zeros(Int, n)
    ind = 1
    for j in 1:length(ct.snds)
        snd_ind = ct.postordering[j]
        snd = ct.snds[snd_ind]
        len_snd = length(snd)
        perm[ind:ind+len_snd-1] .= snd
        snd .= ind:ind+len_snd-1
        ind += len_snd
    end
    ct.perm .= perm

    iperm = invperm(perm)
    for sep in ct.seps, i in 1:length(sep)
        sep[i] = iperm[sep[i]]
    end
    ct.vreps .= first.(ct.snds)
end

"""
    generate_clique_tree(A::SparseMatrixCSC{Tv, Ti}) where {Tv <: AbstractFloat, Ti <: Integer}

Generates a CliqueTree from chordally sparse matrix `A`.

# Reference
[VA15] [Chordal Graphs and Semidefinite Optimization](https://www.seas.ucla.edu/~vandenbe/publications/chordalsdp.pdf) by Lieven Vandenberghe and Martin Andersen
"""
function generate_clique_tree(A::SparseMatrixCSC{Tv, Ti}) where {Tv <: AbstractFloat, Ti <: Integer}
    size(A, 1) != size(A, 2) && error(ArgumentError("A must be square"))
    if !is_chordal(A)
        error(ArgumentError("Matrix must be chordal"))
    end
    peo = !is_chordal(A; peo=1:n) ? collect(1:n) : maximum_cardinality_search(A)
    # NOTE: this sets the diagonal equal to 0
    L = tril(A[peo, peo], -1)

    return CliqueTree(L, peo)
end

"""
    get_postordering(par, child)

Gets a postordering of the tree represented by vector of parents `par` and vector
of children `child`.
"""
function get_postordering(par, child)
    n = length(par)
    root = findfirst(x->par[x] == 0, 1:n)
    isnothing(root) && error(ArgumentError("Can't find root; set par[root] := 0"))
    visit_order = zeros(Int, n)

    # post order via DFS (https://en.wikipedia.org/wiki/Depth-first_search)
    curr = n
    S = [root]
    while !isempty(S)
        v = pop!(S)
        visit_order[v] = curr
        curr -= 1
        push!(S, child[v]...)
    end
    postordering = sort(1:n, by=x->visit_order[x])

    return postordering
end

"""
    get_etree(A)

Compute the parent function of the elimination tree of a symmetric sparse matrix `A`
"""
function get_etree(A::SparseMatrixCSC{Tv, Ti}) where {Tv <: AbstractFloat, Ti <: Integer}
    (!istril(A) && !issymmetric(A)) && error(ArgumentError("A must be symmetric or lower triangular"))
    n = size(A, 1)
    par = zeros(Ti, n)
    ancestor = zeros(Ti, n)

    U = istril(A) ? sparse(A') : A
    _etree(U, par, ancestor)
    count(x->x==0, par) > 1 && error(ArgumentError("A must be connected"))
    return par
end

# (almost) nonallocating verion of above
# Computes using the upper triangular part
# TODO: change to lower triangular??
# Davis. Direct Methods for Sparse Linear Systems, pg 42
function _etree(U::SparseMatrixCSC{Tv, Ti}, par::Vector{Ti}, ancestor::Vector{Ti}) where {Tv <: AbstractFloat, Ti <: Integer}
    n = size(U, 1)
    for col::Ti in 1:n, p in nzrange(U, col) #A.colptr[col]:(A.colptr[col+1]-1)
        i = rowvals(U)[p]
        while !iszero(i) && i < col
            i_next = ancestor[i]
            ancestor[i] = col
            if iszero(i_next)
                par[i] = col
            end
            i = i_next
        end
    end
    return nothing
end

# Algorithm from Pothen and Sun
#   restated in VA, Algorithm 4.1
"""
    max_supernode_etree(L::SparseMatrixCSC, etree_par::Vector{Int})

Constructs the supernodal elimination tree of the chordal graph represented by
lower triangular matrix `L` and with elimination tree `etree_par`. Returns the
reprsentation vertices, the supernodal elimination tree, and a supernode memership
indicator.

Implemented using [VA15, Algorithm 4.1], which was originally formulated by [PS89].

# References
* [PS89] [Compact clique tree data structures in sparse matrix factorizations](https://minds.wisconsin.edu/bitstream/handle/1793/59224/TR897.pdf?sequence=1&isAllowed=y) by Alex Pothen and Chunguang Sun

* [VA15] [Chordal Graphs and Semidefinite Optimization](https://www.seas.ucla.edu/~vandenbe/publications/chordalsdp.pdf) by Lieven Vandenberghe and Martin Andersen
"""
function max_supernode_etree(L::SparseMatrixCSC, etree_par::Vector{Int})
    n = size(L, 1)
    etree_children = get_children_from_par(etree_par)
    post_ord = get_postordering(etree_par, etree_children)
    deg⁺ = get_higher_deg(L)

    vreps = Vector{Int}(undef, 0)
    sizehint!(vreps, n)
    snd_membership = zeros(Int, n)
    snd_par = zeros(Int, n)

    for v = post_ord
        ŵ = -1
        for w in etree_children[v]
            deg⁺[v] > deg⁺[w] - 1 && continue
            ŵ = w
            break
        end

        # v is a rep vertex
        if ŵ == -1
            push!(vreps, v)
            snd_membership[v] = v
            u = v
        # otherwise, add v to snd(ŵ)
        else
            u = snd_membership[ŵ]
            snd_membership[v] = u
        end

        for w in etree_children[v]
            w == ŵ && continue
            # q(z) = u, a(z) = v; z is the vertex in vreps that satisfies w ∈ snd(z)
            z = snd_membership[w]
            snd_par[z] = u
            # a[z] = v
        end
    end

    # numbers the supernodes
    par = zeros(Int, length(vreps))
    @views for (ind, v) in enumerate(snd_par[vreps])
        vrep_ind = findfirst(i->vreps[i] == v, 1:length(vreps))
        if (!isnothing(vrep_ind)) par[ind] = vrep_ind end
    end

    # updates membership similarly
    # TODO: make this efficient
    for (ind, v) in enumerate(snd_membership)
        snd_membership[ind] = findfirst(i->vreps[i] == v, 1:length(vreps))
    end

    return vreps, par, snd_membership

end

# input: parent vector for tree
# returns vector v st v[i] = {k | k is a child of i}
"""
    get_children_from_par(par)

Compute the children of each vertex `v` in a tree specified by `par`
"""
function get_children_from_par(par::Vector{Int})
    n = length(par)
    children = [Vector{Int}(undef, 0) for _ in 1:n]

    for i in 1:n
        par_i = par[i]
        if par_i != 0
            push!(children[par_i], i)
        end
    end

    return children
end

# NOTE: Assumes that L has order σ = 1:n & is lower triangular w zeros on diagonal
"""
    get_higher_deg(L::SparseMatrixCSC)

Returns the higher degrees of each node in the graph represented by lower triangular
matrix `L`.

NOTE: the algorithm assums that `L` has order `σ = 1:n` and zeros on the diagonal.
"""
function get_higher_deg(L::SparseMatrixCSC)
    n = size(L, 1)
    deg⁺ = zeros(Int, n)

    for i in 1:n
        deg⁺[i] = L.colptr[i+1] - L.colptr[i]
    end

    return deg⁺
end

# TODO: Iterator for cliques
"""
    get_cliques(ct::CliqueTree)

Returns the cliques in the graph represented by CliqueTree `ct`.
"""
function get_cliques(ct::CliqueTree)
    return [vcat(ct.snds[i], ct.seps[i]) for i in 1:length(ct.snds)]
end


text:
module Noise

using LinearAlgebra

export Worley, Perlin

struct Worley
    fs
    Worley(d, n) = new(rand(d, n))
end

Worley(n) = Worley(2, n)

function (G::Worley)(n, p)
    Δ = map(f -> norm(p - f), eachcol(G.fs)) |> sort
    return Δ[n]
end

struct Perlin
    field
end

Perlin(xs, ys) = [rand(2) for _ in 1:xs, _ in 1:ys] |> Perlin

function (G::Perlin)(p)
    M = size(G.field) |> collect |> Diagonal
    p = M * p
    A = floor.(p)
    B = A .+ 1
end

# 1,1 => [(1,1), (1, 2), (2, 1), (2, 2)]
# 1,1,1 => [(1,1,1), (1,1,2), (1,2,1), (1,2,2), (2,1,1), (2, 1, 2), (2, 2, 2)]

end # module


text:
module MiniPL

import Base: *,+,-,÷,<,!,&,==,showerror

include("Interpreter.jl")

DEBUG = false

end # module

max_stars_repo_name: juhamyllari/mini-pl-julia

max_issues_repo_name: juhamyllari/mini-pl-julia

max_forks_repo_name: juhamyllari/mini-pl-julia


text:
# Legendre-Gauss quadrature weights and nodes
# This code was originally written for Matlab
#  by Greg von Winckel - 02/25/2004
# posted on Mathworks website and licensed under BSD License
# http://fr.mathworks.com/matlabcentral/fileexchange/4540-legendre-gauss-quadrature-weights-and-nodes
#
#

function lgwt(N,a,b)

	N=N-1;
	N1=N+1; N2=N+2;
	xu=linspace(-1,1,N1)';

	y=(cos((2*(0:N)'+1)*pi/(2*N+2))+(0.27/N1)*sin(pi*xu*N/N2))';

	L=zeros(N1,N2);

	Lp=zeros(N1,N2);
	y0=2;

	while maximum(abs(y-y0))>eps()
		L[:,1]=1;
		Lp[:,1]=0;
		L[:,2]=y;
		for k=2:N1
			L[:,k+1]=( (2*k-1)*y.*L[:,k]-(k-1)*L[:,k-1] )/k;
		end
		Lp=(N2)*( L[:,N1]-y.*L[:,N2] )./(1-y.^2);
		y0=y;
		y=y0-L[:,N2]./Lp;
	end

	x=(a*(1-y)+b*(1+y))/2;

	w=(b-a)./((1-y.^2).*Lp.^2)*(N2/N1)^2;

	return x,w
end

max_stars_repo_name: flavioluiz/PortHamiltonian

max_issues_repo_name: flavioluiz/PortHamiltonian

max_forks_repo_name: flavioluiz/PortHamiltonian


text:
module TestCaretWrapper

include(joinpath("..", "fixture_learners.jl"))
using .FixtureLearners
nfcp = NumericFeatureClassification()

using FactCheck

using MLBase
importall Orchestra.Types
importall Orchestra.Transformers.CaretWrapper
CW = CaretWrapper
using PyCall
@pyimport rpy2.robjects as RO
@pyimport rpy2.robjects.packages as RP
@pyimport rpy2.robjects.numpy2ri as N2R
N2R.activate()
RP.importr("caret")

function behavior_check(caret_learner::String, impl_options=Dict())
  # Predict with Orchestra learner
  srand(1)
  pycall(RO.r["set.seed"], PyObject, 1)
  learner = CRTLearner({
    :learner => caret_learner,
    :impl_options => impl_options
  })
  orchestra_predictions = fit_and_transform!(learner, nfcp)

  # Predict with backend learner
  srand(1)
  pycall(RO.r["set.seed"], PyObject, 1)
  (r_dataset_df, label_factors) = CW.dataset_to_r_dataframe(
    nfcp.train_instances, nfcp.train_labels
  )
  label_factors = collect(label_factors)
  caret_formula = RO.Formula("Y ~ .")
  r_fit_control = pycall(RO.r[:trainControl], PyObject,
    method = "none"
  )
  if isempty(impl_options)
    r_model = pycall(RO.r[:train], PyObject,
      caret_formula,
      method = caret_learner,
      data = r_dataset_df,
      trControl = r_fit_control,
      tuneLength = 1
    )
  else
    r_model = pycall(RO.r[:train], PyObject,
      caret_formula,
      method = caret_learner,
      data = r_dataset_df,
      trControl = r_fit_control,
      tuneGrid = RO.DataFrame(impl_options)
    )
  end
  (r_instance_df, _) = CW.dataset_to_r_dataframe(nfcp.test_instances)
  original_predictions = collect(pycall(RO.r[:predict], PyObject,
    r_model,
    newdata = r_instance_df
  ))
  original_predictions = map(x -> label_factors[x], original_predictions)

  # Verify same predictions
  @fact orchestra_predictions => original_predictions
end

facts("CARET learners") do
  context("CRTLearner gives same results as its backend") do
    caret_learners = ["svmLinear", "nnet", "earth"]
    for caret_learner in caret_learners
      behavior_check(caret_learner)
    end
  end
  context("CRTLearner with options gives same results as its backend") do
    behavior_check("svmLinear", {:C => 5.0})
  end

  context("CRTLearner throws on incompatible feature") do
    instances = {
      1 "a";
      2 3;
    }
    labels = [
      "a";
      "b";
    ]

    learner = CRTLearner()
    @fact_throws fit!(learner, instances, labels)
  end
end

end # module


text:
module OmegaViz
using Omega
import Plots
import Makie
import AbstractPlotting
import Omega.Space
using Omega.Space: getdim

include("contour.jl")

export ωcontour

end # module


text:
using Test
using NeXLCore

@testset "Materials" begin
    @testset "K411" begin
        k411 = material(
            "K411",
            n"O" => 0.423667,
            n"Mg" => 0.0884654,
            n"Si" => 0.253817,
            n"Ca" => 0.110563,
            n"Fe" => 0.112087,
            properties = Dict{Symbol,Any}(:Density => 3.5),
        )

        @test k411[n"Ca"] == 0.110563
        @test k411[n"Ce"] == 0.0
        @test massfraction(k411)[n"Ca"] == k411[n"Ca"]
        @test k411[n"Ca"] == 0.110563

        nk411 = normalizedmassfraction(k411)

        @test isapprox(nk411[n"O"], 0.428553, atol = 1.0e-6)
        @test isapprox(nk411[n"Mg"], 0.0894855, atol = 1.0e-6)
        @test isapprox(nk411[n"Si"], 0.256744, atol = 1.0e-6)
        @test isapprox(nk411[n"Ca"], 0.111838, atol = 1.0e-6)
        @test isapprox(nk411[n"Fe"], 0.11338, atol = 1.0e-6)

        @test analyticaltotal(k411) == 0.9885994

        ak411 = atomicfraction(k411)
        @test isapprox(ak411[n"O"], 0.602876, atol = 1.0e-5)
        @test isapprox(ak411[n"Mg"], 0.0828675, atol = 1.0e-5)
        @test isapprox(ak411[n"Si"], 0.205753, atol = 1.0e-5)
        @test isapprox(ak411[n"Ca"], 0.0628072, atol = 1.0e-5)
        @test isapprox(ak411[n"Fe"], 0.0456961, atol = 1.0e-5)

        k411_2 = material(
            "K411",
            n"O" => 0.423667,
            n"Mg" => 0.0884654,
            n"Si" => 0.253817,
            n"Ca" => 0.110563,
            n"Fe" => 0.112087,
            atomicweights = Dict(n"Ca" => 41.0, n"O" => 16.1),
            properties = Dict{Symbol,Any}(:Density => 3.5),
        )

        @test a(n"Ca", k411) == a(n"Ca")
        @test a(n"Ca", k411_2) == 41.0
        @test a(n"O", k411_2) == 16.1
        @test a(n"Si", k411) == a(n"Si")
        @test a(n"Fe", k411) == a(n"Fe")
        @test a(n"Si", k411_2) == a(n"Si")
        @test a(n"Fe", k411_2) == a(n"Fe")

        d411 = delete(k411, n"O")
        @test d411[n"O"]==0.0
        @test all(el->k411[el]==d411[el], keys(d411))
        d411 = delete(k411, [ n"O", n"Fe" ])
        @test d411[n"O"]==0.0
        @test d411[n"Fe"]==0.0
        @test all(el->k411[el]==d411[el], keys(d411))
    end
    @testset "Parse" begin
        albite = parse(Material, "AlNaSi3O8")
        @test isapprox(analyticaltotal(albite), 1.0, atol = 1.0e-8)
        @test isapprox(albite[n"O"], 0.488116, atol = 0.00001)
        @test isapprox(albite[n"Al"], 0.102895, atol = 0.00001)
        @test isapprox(albite[n"Si"], 0.321316, atol = 0.00001)
        @test isapprox(albite[n"Na"], 0.0876726, atol = 0.00001)

        apatite = parse(Material, "Ca5(PO4)3F")
        @test isapprox(apatite[n"O"], 0.38071, atol = 0.00001)
        @test isapprox(apatite[n"F"], 0.0376719, atol = 0.00001)
        @test isapprox(apatite[n"P"], 0.184257, atol = 0.00001)
        @test isapprox(apatite[n"Ca"], 0.397361, atol = 0.00001)

        apatite = parse(Material, "Ca5(PO4)3(OH)")
        @test isapprox(apatite[n"O"], 0.41407, atol = 0.00001)
        @test isapprox(apatite[n"H"], 0.0020066, atol = 0.00001)
        @test isapprox(apatite[n"P"], 0.184987, atol = 0.00001)
        @test isapprox(apatite[n"Ca"], 0.398936, atol = 0.00001)

        apatite = parse(Material, "Ca5(PO4)3OH")
        @test isapprox(apatite[n"O"], 0.41407, atol = 0.00001)
        @test isapprox(apatite[n"H"], 0.0020066, atol = 0.00001)
        @test isapprox(apatite[n"P"], 0.184987, atol = 0.00001)
        @test isapprox(apatite[n"Ca"], 0.398936, atol = 0.00001)

        apatite = parse(Material, "Ca5(PO4)3OH")
        @test isapprox(apatite[n"O"], 0.41407, atol = 0.00001)
        @test isapprox(apatite[n"H"], 0.0020066, atol = 0.00001)
        @test isapprox(apatite[n"P"], 0.184987, atol = 0.00001)
        @test isapprox(apatite[n"Ca"], 0.398936, atol = 0.00001)

        quartz1 = parse(Material, "Si23O46")
        quartz2 = parse(Material, "SiO2")
        @test all(e -> isapprox(quartz1[e], quartz2[e], atol = 0.00001), keys(quartz1))
        @test isapprox(quartz1, quartz2, atol = 1.0e-6)

        ss = parse(Material, "0.8*Fe+0.15*Ni+0.04*Cr")
        @test ss[n"Fe"] == 0.80
        @test ss[n"Ni"] == 0.15
        @test ss[n"Cr"] == 0.04

        m = parse(Material, "0.6*SiO2+0.4*Al2O3")
        @test isapprox(m[n"Si"], 0.6 * mat"SiO2"[n"Si"], atol = 1.0e-6)
        @test isapprox(m[n"Al"], 0.4 * mat"Al2O3"[n"Al"], atol = 1.0e-6)
        @test isapprox(
            m[n"O"],
            0.6 * mat"SiO2"[n"O"] + 0.4 * mat"Al2O3"[n"O"],
            atol = 1.0e-6,
        )

        @test isapprox(
            parse(Material, "Ba(Al2Si2O8)", name = "Paracelsian"),
            mat"Ba(Al2Si2O8)",
            atol = 1.0e-8,
        )
        @test isapprox(mat"Ba(Al2Si2O8)"[n"Ba"], 0.36576, atol = 0.00001)

        k458=parse(Material, "(0.4179±0.0020)*Ba+(0.0301±0.0006)*Zn+(0.2305±0.0034)*Si+(0.3186±0.0040)*O",
                name="SRM-1873 K-458", description="https://www-s.nist.gov/srmors/certificates/1873.pdf")
        @test isapprox(k458,
                material("Test", Dict(n"Ba"=>uv(0.4179,0.0020), n"Zn"=>uv(0.0301,0.0006), n"Si"=>uv(0.2305,0.0034), n"O"=>uv(0.3186,0.0040))))
        @test value(k458[n"Zn"])==0.0301
        @test σ(k458[n"Zn"])==0.0006

        para = material(
            "Paracelsian",
            n"O" => 0.340906,
            n"Ba" => 0.36576,
            n"Si" => 0.149607,
            n"Al" => 0.143727,
            properties = Dict{Symbol,Any}(:Density => 3.5),
        )
        @test isapprox(para, mat"Ba(Al2Si2O8)", atol = 1.0e-5)
        para = material(
            "Paracelsian",
            n"O" => 0.340906,
            n"Ba" => 0.36576,
            n"Si" => 0.149607,
            n"Al" => 0.143727,
            properties = Dict{Symbol,Any}(:Density => 3.5),
        )
        @test isapprox(para, mat"Ba(Al2Si2O8)", atol = 1.0e-5)

        para = material(
            "Paracelsian",
            n"O" => 0.340906,
            n"Ba" => 0.36576,
            n"Si" => 0.149607,
            n"Al" => 0.143727,
            properties = Dict{Symbol,Any}(
                :Density => 3.5,
                :Description => "Barium-rich feldspar",
            ),
        )
        @test isapprox(para, mat"Ba(Al2Si2O8)", atol = 1.0e-5)
        # Define a couple of additional properties...
        para[:Color] = "Colorless to white, light yellow tint"
        para[:System] = :Monoclinic
        @test para[:Description] == "Barium-rich feldspar"
        @test para[:Density] == 3.5
        @test density(para) == 3.5
        @test para[:Color] == "Colorless to white, light yellow tint"
        @test para[:System] == :Monoclinic
    end
    @testset "K412 w. uncertainty" begin
        k412u = material(
            "K412",
            n"Fe" => uv(0.0774, 0.0014),
            n"Al" => uv(0.0491, 0.0010),
            n"O" => uv(0.4276, 0.0006),
            n"Ca" => parse(UncertainValue, "0.1090-+0.0012"),
            n"Si" => parse(UncertainValue, "0.2120+-0.0005"),
            n"Mg" => parse(UncertainValue, "0.1166±0.0010"),
            density = 3.45,
            pedigree = "NIST SRM-470",
            conductivity = :Insulator,
            description = "NIST SRM glass - olivine simulant",
        )
        @test isapprox(value(k412u[n"O"]), 0.4276, atol = 0.00001)
        @test isapprox(σ(k412u[n"O"]), 0.0006, atol = 0.00001)
        @test isapprox(k412u[n"O"], uv(0.4276, 0.0006), atol = 0.00001)
        @test isapprox(value(k412u[n"Si"]), 0.2120, atol = 0.00001)
        @test isapprox(σ(k412u[n"Si"]), 0.0005, atol = 0.00001)
        @test isapprox(analyticaltotal(k412u), 0.9917, atol = 0.00005)
        @test k412u[:Density] == 3.45
        @test startswith(k412u[:Pedigree], "NIST SRM")
        @test k412u[:Conductivity] == :Insulator
        @test startswith(k412u[:Description], "NIST SRM glass")

        k412x = mat"0.4535⋅SiO2+0.0996⋅FeO+0.1933⋅MgO+0.1525⋅CaO+0.0927⋅Al2O3"
        k412y = parse(
            Material,
            "0.4535*SiO2+0.0996*FeO+0.1933*MgO+0.1525*CaO+0.0927*Al2O3",
            name = "K412",
        )
        @test all(isapprox(value(k412x[el]), value(k412u[el]), atol = 0.0001) for el in keys(k412x))
        @test isapprox(k412x, k412y, atol = 0.0001)
        @test name(k412y) == "K412"
        q = parse(
            Material,
            "SiO2",
            name = "Quartz",
            density = 2.32,
            pedigree = "Stoichiometric",
        )
        @test q[:Density] == 2.32
        @test q[:Pedigree] == "Stoichiometric"
        @test name(q) == "Quartz"
    end
end


text:
"""
Actor-Critic with a surprise-modulated learning rate for the actor (core)

See also actorcriticlearner.jl

Critic implemented by td-lambda
This implementation follows Sutton & Barto,
Chapter 7.7 "Eligibility Traces for Actor-Critic Methods"
The policy update follows the scheme given in chapter 6.6 where
the policy changes inversely proportional to the action selection probability
"""
mutable struct ActorCriticModulated
    nr_states::Int
    nr_actions::Int
    α_critic::Float64
    α_actor::Float64
    α_actor_eff::Float64
    γ::Float64
    λ_critic::Float64
    λ_actor::Float64
    initvalue::Float64
    V::Array{Float64, 1}
    policyparameters::Array{Float64, 2}
    eligTrace_critic::Array{Float64, 1}
    eligTrace_actor::Array{Float64, 2}
    td_error::Array{Float64, 1}
    use_replacing_trace::Bool
    do_debug_print::Bool
    actioncost::Float64
end
function ActorCriticModulated(; nr_states = 10, nr_actions = 4,
    α_critic = 0.1, α_actor = 0.1, γ = .9,
    λ_critic = 0.9, λ_actor = 0.9, initvalue = 0.,
    V = zeros(Float64, nr_states) .+ initvalue,
    policyparameters = zeros(Float64, nr_states, nr_actions),
    eligTrace_critic = zeros(Float64, nr_states),
    eligTrace_actor = zeros(Float64, nr_states, nr_actions),
    use_replacing_trace = true,
    do_debug_print = false,
    actioncost = 0.)

    ActorCriticModulated(nr_states, nr_actions, α_critic, α_actor, α_actor, γ,
        λ_critic, λ_actor, initvalue, V,
        policyparameters, eligTrace_critic, eligTrace_actor, [0.],
        use_replacing_trace, do_debug_print, actioncost)
end
export ActorCriticModulated
function getlearningrate!(learner::ActorCriticModulated,
                        modulatingfactor)
    learner.α_actor_eff = learner.α_actor * modulatingfactor
end


text:
mutable struct Simulation
	w0        :: Matrix{Float64}
	w1        :: Matrix{Float64}
	cons      :: Matrix{Float64}
	util      :: Matrix{Float64}
	ystate    :: Matrix{Int}
	inc       :: Matrix{Float64}
	prob_work :: Matrix{Float64}
	ret_age   :: Vector{Int}
	worker    :: Matrix{Int}
	p :: Param
	function Simulation(p::Param)
		this = new()
		this.w0        = fill(NaN,p.nsims,p.nT)  # beginning of period wealth
		this.w1        = fill(NaN,p.nsims,p.nT)  # end of period wealth
		this.cons      = fill(NaN,p.nsims,p.nT)
		this.util      = fill(NaN,p.nsims,p.nT)
		this.ystate    = zeros(Int,p.nsims,p.nT)
		this.inc       = fill(NaN,p.nsims,p.nT)
		this.worker    = zeros(Int,p.nsims,p.nT)
		this.prob_work = fill(NaN,p.nsims,p.nT)
		this.ret_age   = zeros(Int,p.nsims)
		this.p = deepcopy(p)
		this
	end
end

mutable struct BSimulation
	ainit        :: Matrix{Float64}
	w0        :: Matrix{Float64}
	w1        :: Matrix{Float64}
	cons      :: Matrix{Float64}
	ystate    :: Matrix{Int}
	inc       :: Matrix{Float64}
	prob_nofile :: Matrix{Float64}
	file_age   :: Vector{Int}
	bkflag    :: Matrix{Bool}  # on/off
	vfile :: Matrix{Float64}
	vnofile :: Matrix{Float64}
	p :: Param
	function BSimulation(p::Param)
		this = new()
		this.ainit        = fill(NaN,p.nsims,p.nT)  # beginning of period wealth
		this.w0        = fill(NaN,p.nsims,p.nT)  # beginning of period wealth
		this.w1        = fill(NaN,p.nsims,p.nT)  # end of period wealth
		this.cons      = fill(NaN,p.nsims,p.nT)
		this.vfile      = fill(NaN,p.nsims,p.nT)
		this.vnofile      = fill(NaN,p.nsims,p.nT)
		this.ystate    = zeros(Int,p.nsims,p.nT)
		this.inc       = fill(NaN,p.nsims,p.nT)
		this.bkflag    = falses(p.nsims,p.nT)
		this.prob_nofile = fill(NaN,p.nsims,p.nT)
		this.file_age   = zeros(Int,p.nsims)
		this.p = deepcopy(p)
		this
	end
end

function sim(m::FModel,p::Param)

	s = Simulation(p)

	working   = trues(p.nsims)  # indicator of currentlcy working or not
	vmat      = zeros(2,p.nsims)

	N = Normal(0,1)  # normal dist N(0,σ)

	# random setup
	Random.seed!(p.rseed)
	wshocks = quantile.(N, rand(p.nsims,p.nT)) .* p.sigma
	dshocks = rand(p.nsims,p.nT) # discrete choice shocks
	rshocks = rand(p.nsims,p.nT) # returning to work shocs
	w0shocks = rand(p.nsims)

	for it in 1:p.nT
		if it == 1
			# draw from initial distributions
			s.w0[:, it] = p.initw0 .+ w0shocks *(p.initw1 - p.initw0)
			s.worker[:, it] .= 1   # all work in frist period
		else

			# working state
			s.worker[ working, it] .= 1 # chose working, so stay worker
			s.worker[.!(working), it] .= 2 # chose retire
			s.ret_age[(s.worker[:,it-1] .== 1) .& (s.worker[:,it] .== 2)] .= it

			# income
			s.inc[.!(working), it] .= p.pension
			s.inc[ working, it] .= income(it,p,wshocks[working,it])

			# end of period wealth in period it
			s.w0[ :     , it] .= s.w1[:  , it-1]*p.R .+ s.inc[ : , it]
			# s.w0[working, it] .= s.w1[working, it-1]*p.R .+ s.inc[ working, it]
		end
		# consumption
		s.cons[ working,it]    = gety(interp(m.c[1,it].env,s.w0[ working, it]))
		s.util[ working,it]    = u(s.cons[ working,it], true, s.p)

		s.cons[.!(working),it] = gety(interp(m.c[2,it].env,s.w0[.!(working), it]))
		s.util[.!(working),it] = u(s.cons[.!(working),it], false, s.p)
		# end of period wealth
		s.w1[:,it] =  s.w0[:,it] - s.cons[:,it]
		clamp!(s.w1,p.a_low, p.a_high)

		if any(s.cons[ :,it] .< 0)
			println("neg cons")
			println(s.cons[s.cons[ :,it] .< 0, it])
			println(s.w0[.!(working), it])
			error()
		end

		# CCP to remain worker
		vmat[1,:] = vfun(1,it, s.cons[:,it] ,s.w0[:,it], m.v[1,it],p )
		vmat[2,:] = vfun(2,it, s.cons[:,it] ,s.w0[:,it], m.v[2,it],p )
		s.prob_work[:,it] = (s.worker[:,it] .== 1) .* ccp(vmat,p)  # prob to remain worker. 0 for current non-workers

		# discrete choice of workers
		working[:] = s.prob_work[:,it] .> dshocks[:,it]

		# random transition back into work state for retirees
		working[(s.worker[:,it] .== 2)] = p.delta .> rshocks[(s.worker[:,it] .== 2),it]
	end
	s
end

function sim(m::GModel,p::Param)

	s = Simulation(p)

	# random
	Random.seed!(p.rseed)
	wshocks = rand(p.nsims,p.nT)
	dshocks = rand(p.nsims,p.nT)
	w0shocks = rand(p.nsims)
	s.ystate[:,1] .= rand(1:p.ny,p.nsims)  # randomly allocate to an income state

	iy = 0 # y-state index
	working   = trues(p.nsims)  # indicator of currentlcy working or not
	vmat      = zeros(2,p.nsims)
	# s.shocks[:] = rand(p.nT * p.nsims)
	GG = cumsum(m.ywgt,dims = 2)
	for it in 2:p.nT
		for i in 1:p.nsims
			s.ystate[i,it] = searchsortedfirst(GG[s.ystate[i,it-1],:], wshocks[i,it])
		end
	end

	for it in 1:p.nT
		if it == 1
			# draw from initial distributions
			s.w0[:, it] = p.initw0 .+ w0shocks * (p.initw1 - p.initw0)
			s.worker[:, it] .= 1   # all work in frist period
			s.inc[ working, it] .= income(it,p,m.yvec[s.ystate[working,it]])
		else

			# working state
			s.worker[ working, it] .= 1 # chose working, so stay worker
			s.worker[.!(working), it] .= 2 # chose retire
			s.ret_age[(s.worker[:,it-1] .== 1) .& (s.worker[:,it] .== 2)] .= it

			# income
			s.inc[.!(working), it] .= p.pension #* income(it,p,m.yvec[s.ystate[.!(working),it]])
			s.inc[ working, it] .= income(it,p,m.yvec[s.ystate[working,it]])

			# start of period wealth in period it
			s.w0[ :     , it] .= s.w1[:      , it-1]*p.R.+ s.inc[ : ,it]
			# s.w0[working, it] .= s.w1[working, it-1]*p.R .+ s.inc[ working, it]
		end
		# consumption and value
		# can be made much faster by grouping individuals by `iy` state
		for i in 1:p.nsims
			iy = s.ystate[i,it]
			if working[i]
				s.cons[i,it]    = gety(interp(m.c[1,iy,it].env,[s.w0[ i, it]]))[1]
			else
				s.cons[i,it]    = gety(interp(m.c[2,iy,it].env,[s.w0[ i, it]]))[1]
			end
			vmat[1,i] = vfun(1,it, [s.cons[i,it] ],[s.w0[i,it]], m.v[1,iy,it],p )[1]
			vmat[2,i] = vfun(2,it, [s.cons[i,it] ],[s.w0[i,it]], m.v[2,iy,it],p )[1]
		end
		# end of period wealth
		s.util[ working,it]    = u(s.cons[ working,it], true, s.p)
		s.util[ .!(working),it]    = u(s.cons[  .!(working),it], false, s.p)

		s.w1[:,it] = s.w0[:,it] - s.cons[:,it]

		# CCP to remain worker
		s.prob_work[:,it] = (s.worker[:,it] .== 1) .* ccp(vmat,p)  # prob to remain worker

		# discrete choice
		working[:] = s.prob_work[:,it] .> dshocks[:,it]
	end
	s
end

function sim(m::BModel,p::Param)

	s = BSimulation(p)

	# s.bkflag: 1 = flag off, 2 = flag on

	# random
	Random.seed!(p.rseed)
	wshocks = rand(p.nsims,p.nT)
	dshocks = rand(p.nsims,p.nT)
	w0shocks = rand(p.nsims)
	s.ystate[:,1] .= rand(1:p.ny,p.nsims)  # randomly allocate to an income state

	iy = 0 # y-state index

	# everybody starts as not
	filing    = falses(p.nsims)  # indicator of whether is not filing in current period
	notfiling = trues(p.nsims)
	vmat      = zeros(2,p.nsims)

	# s.shocks[:] = rand(p.nT * p.nsims)
	GG = cumsum(m.ywgt,dims = 2)
	for it in 2:p.nT
		for i in 1:p.nsims
			s.ystate[i,it] = searchsortedfirst(GG[s.ystate[i,it-1],:], wshocks[i,it])
		end
	end

	for it in 1:p.nT
		# println(it)
		notfiling[:] .= .!(filing)
		if it == 1
			# draw from initial distributions
			s.w0[:, it] = p.initw0 .+ w0shocks * (p.initw1 - p.initw0)
			s.ainit[:,it] = s.w0[:, it]

			# everybody starts with bkflag off
			s.inc[: , it] .= income(it,p,m.yvec[s.ystate[:,it]])  # all get their income
		else

			# set bk flag at start of period depending on last periods choices:
			# if not filing, flag does not change
			s.bkflag[notfiling , it] = s.bkflag[notfiling , it-1]
			# else, flag goes from false to true
			s.bkflag[filing , it] .= true  # have chosen file in previous period, so flag is on
			s.file_age[filing] .= it-1  # age at which filed

			# income
			s.inc[: , it] .= income(it,p,m.yvec[s.ystate[:,it]])  # all get their income

			# s.inc[.!(notfiling), it] .= 0.0
			# s.inc[ notfiling, it] .= income(it,p,m.yvec[s.ystate[notfiling,it]])

			# start of period wealth in period it
			# s.w0[ :     , it] .= s.w1[:      , it-1]*p.R
			s.ainit[:,it] = s.w1[:, it-1]*p.R
			s.w0[notfiling, it] .= s.ainit[notfiling,it] .+ s.inc[notfiling, it]
			s.w0[filing    , it] .= 0.0 .+ s.inc[filing, it]
		end

		# end of recording current state
		fill!(filing, false)

		# figure out discrete choice in non bk state first
		inoflag = findall(.!(s.bkflag[:,it]))
		  iflag = findall(   s.bkflag[:,it])
		for i in inoflag
			iy = s.ystate[i,it]
			vmat[1,i] = gety(interp(m.v[1,iy,it].env,[s.w0[ i, it]]))[1]
			vmat[2,i] = s.ainit[i,it] < 0 ? m.v[2,iy,it].env[1].y : -Inf   # file only with neg assets
			s.vnofile[i,it] = vmat[1,i]
			s.vfile[i,it] = vmat[2,i]
		end

		s.prob_nofile[inoflag,it] = ccp(vmat[:,inoflag],p)  # prob to not file
		s.prob_nofile[iflag  ,it] .= NaN   # no discrete choice with a flag on
		# discrete choice
		filing[inoflag] = s.prob_nofile[inoflag,it] .< dshocks[inoflag,it]
		filing[iflag]  .= false

		# consumption depends on whehter flag and if not, which dchoice
		inoflag_nofile = findall( .!(s.bkflag[:,it]) .& .!(filing) )
		inoflag_file   = findall( .!(s.bkflag[:,it]) .& filing )
		for i in inoflag_nofile
			s.cons[i,it] = gety(interp(m.c[1,s.ystate[i,it],it].env,[s.w0[ i, it]]))[1]
			# end of period wealth
			s.w1[i,it] = s.w0[i,it] - s.cons[i,it]
		end
		for i in inoflag_file
			s.cons[i,it]    = m.c[2,s.ystate[i,it],it].env[1].y
			# end of period wealth
			s.w1[i,it] = 0.0
		end
		for i in iflag
			s.cons[i,it] = gety(interp(m.cbk[s.ystate[i,it],it].env,[s.w0[ i, it]]))[1]
			s.w1[i,it] = s.w0[i,it] - s.cons[i,it]

		end

	end
	s
end

function rbs()
	m,p = DCEGM.runbk()
	sim(m,p)
end

function rs(;par = Dict())
	m,p = DCEGM.runf(par = par)
	sim(m,p)
end
function rsp(;par = Dict())
	plot_s(rs(par=par))
end

function rsg(;par = Dict())
	m,p = DCEGM.rung(par = par)
	sim(m,p)
end
function rsgp(;par = Dict())
	plot_s(rsg(par=par))
end


text:
"""
    result = pslq(x::Vector; criteria)

Perform PSLQ iterations for the vector x until termination criteria are fulfilled.
`result` always contains the fields: `code`, `iter`, `A`, `B`, `H`.
The returned result is a named tuple and may contain different result variants:

`result.code == 1`:

 A solution `m` was found `|m' * x| <= cr.atol + cr.rtol*|x|`
 with `m = result.B[:,result.col]`

`result.code == 2`:
    No solution `m` possible with `|m| <= cr.solutionmax`

`result.code == 3`:
    No decision after number of iterations exceeds `cr.itermax`

`result.code == 4`:
    Inaccuracy in calculation - retry with higher precision input data

The input criteria are transformed to object `cr` by method `make_criteria(criteria, x)`.
"""
function pslq(x::Vector{T}; criteria=nothing) where T<:Real
    n = length(x)
    H = make_H(x)
    # criteria = make_criteria(criteria, x)
    A = Matrix(I*integertype(T)(1), n, n)
    B = copy(A)
    D = UnitLowerTriangular(similar(B)) # working area
    iter = 0
    result = (code=0, result=similar(B, n, 0))
    while iter < criteria.itermax && result.code == 0
        iter += 1
        pslq_step!(H, A, B, D)
        result = convergence(x, H, A, B, iter, criteria)
    end
    result
end

"""
    result = convergence(x, H, A, B, iter, criteria)

Termination check and provide result for `pslq`.
The arguments `H`, `A`, `B` are as described in the algoritm paper.
`result.code == 0` means that iterations should continue.
"""
function convergence(x::Vector{T}, H::AbstractMatrix, A::Matrix, B::Matrix,
                     iter::Integer, criteria::NamedTuple) where T
    n = size(x, 1)
    function xtb(j)
        s1 = s2 = T(0)
        for k = 1:n
            bkjx = B[k,j] * x[k]
            s1 += bkjx
            s2 += abs(bkjx)
        end
        s1, s2
    end
    resvec = []
    for j = 1:n
        s1, s2 = xtb(j)
        if abs(s1) <= s2 * criteria.rtol + criteria.atol
            push!(resvec, (code=1, col=j, xb=s1, xc=s2))
        end
    end
    if length(resvec) >= 1
        r = resvec[1]
        return (code=r.code, iter=iter, col=r.col, xb=r.xb, xc=r.xc, A=A, B=B, H=H)
    end
    (code=0, iter=iter, A=A, B=B, H=H)
end

function make_criteria(criteria, x)
    return (rtol=10*eps(eltype(x)), atol=10*eps(norm(x)), itermax=10000)
end


text:
using CrossMappings
using Test
using StatsBase

include("crossmapping.jl")
include("convergentcrossmapping.jl")


text:
## Huber, M., 2016. Nearly optimal Bernoulli factories for linear functions.
## Combinatorics, Probability and Computing, 25(4), pp.577-591.
## input: f() returns a Bernoulli(p) variate, ϵ < 1-C*p
## output: a Bernoulli(C*p) variate, and the number of calls to f
function _huber2016(f::F, C::Float64, ϵ::Float64, rng::RNG) where
  {F<:Function, RNG<:AbstractRNG}
  γ::Float64 = 0.5
  k::Float64 = 2.3/(γ*ϵ)
  i::Int64 = 1
  ϵ = min(ϵ, 0.644)
  R::Bool = true
  flips::Int64 = 0
  while i != 0 && R
    α::Float64 = (C-1.0)/C
    while i != 0 && i < k
      B::Bool = f()
      flips += 1
      G::Int64 = _sampleGeometric(α, rng)
      i += (1-B)*G - 1
    end
    if i >= k
      R = rand(rng) < (1.0+γ*ϵ)^(-i)
      C *= 1.0+γ*ϵ
      ϵ *= (1.0-γ)
      k /= (1.0-γ)
    end
  end
  coin::Bool = i == 0
  return coin, flips
end


text:
function plot_marginal(r::SamplingResults{I}, a) where I<:InferenceType
    l = @layout([a{0.01h}; [b c]])
    title = plot(
        title = "Marginal $a posterior, computed using $(string(typeof(r.interpretation)))",
        grid = false,
        showaxis = false,
        bottom_margin = -5Plots.px
    )
    p1 = histogram(
        r[a], normed=true, color=:black,
        xlabel=a, ylabel="p($a | data)", label=""
    )
    p2 = plot(
        r[a], color=:black,
        xlabel="Trace index", ylabel=a, label=""
    )
    plot(title, p1, p2, layout=l)
end

function plot_marginal(r::SamplingResults{I}, a, savepath::S, savename::S) where{I<:InferenceType, S<:AbstractString}
    mkpath(savepath)
    fname = joinpath(savepath, savename)
    p = plot_marginal(r, a)
    plot!(p, size=(600,400))
    savefig(p, fname)
end

export plot_marginal


text:
# """
# idris *.idr -o out.qb --codegen qb --cg-opt "--javaName" --cg-opt "--symemu"
#
# Requirements:
# - (v1.3) pkg> add ArgParse
# - (v1.3) pkg> add MLStyle
using MLStyle
using ArgParse

literal_map(kind, x) =
    @match String(kind) begin
        "float" => parse(Float64, x)
        "int" => parse(Int64, x)
        "bigInt" => parse(BigInt, x)
        "char" => x[1]
        "string" => x
        "bool" => parse(Bool, x)
        "unit" => nothing
        "symbol" => QuoteNode(Symbol(x))
        _ => error(kind)
    end

function ExternalCall(name, args)
    name = Symbol(name)
    :(__RTS.$name($(args...)))
end

function ExternalVar(name)
    name = Symbol(name)
    return :(__RTS.$name)
end

function Var(name)
    name = Symbol(name)
    return name
end

function Call(name, args)
    name = Symbol(name)
    return :($name($(args...)))
end

function  Defun(name, args, body)
    name = Symbol(name)
    args = Symbol.(args)
    :(function $name($(args...)); $(body...) end)
end

function Introduction(name)
    name = Symbol(name)
    :($name = nothing)
end

function Update(name, exp)
    name = Symbol(name)
    :($name = $exp)
end

function Constant(n)
    n
end

function Switch(var, xs, body)
    isempty(xs) && return Expr(:block, body...)
    body = Expr(:block, body...)
    for (cc, stmts) in reverse(xs)
        body = Expr(:elseif, Expr(:block, :($var == $cc)), Expr(:block, stmts...), body)
    end

    body.head = :if
    return body
end

function If(cond, t, e)
    return Expr(:if, cond, Expr(:block, t...), Expr(:block, e...))
end

function EffectExpr(exp)
   exp
end

function Return(exp)
    :(return $exp)
end

function read_and_gen(io)
    ctor_stack = []
    obj_stack = []
    left_stack = []
    left = 1
    while true
        while (left !== 0)
            s = readuntil(io, '\n'; keep=false)
            if isempty(s)
                continue
            end
            pats = split(s)
            dispatch = pats[1]
            left -= 1
            if dispatch == "constructor"
                cons = pats[2]
                n = parse(Int, pats[3])
                push!(left_stack, left)
                left = n
                ctor = getproperty(@__MODULE__, Symbol(cons))
                push!(ctor_stack, (ctor, n))
            elseif dispatch == "literal"
                kind = pats[2]
                length = parse(Int, pats[3])
                buf = String(read(io, length))
                readuntil(io, '\n'; keep=false)
                buf = literal_map(kind, buf)
                push!(obj_stack, buf)

            elseif dispatch == "list"
                n = parse(Int, pats[2])
                push!(left_stack, left)
                left = n
                push!(ctor_stack, (nothing, n))
            else
                error("malformed qb format")
            end
        end

        isempty(ctor_stack) && begin
            @assert length(obj_stack) === 1
            return obj_stack[1]
        end
        (ctor, n) = pop!(ctor_stack)
        args = []
        for _ in 1:n
            push!(args, pop!(obj_stack))
        end
        reverse!(args)

        if ctor === nothing
            v = args
        else
            v = ctor(args...)
        end
        push!(obj_stack, v)

        left = pop!(left_stack)
    end
end

s = ArgParseSettings()
@add_arg_table! s begin
    "filename"
        help = "input QB file"
        required = true
    "out"
        help = "output Julia file"
        required = true
end

parsed_args = parse_args(ARGS, s)

f = open(parsed_args["filename"])
big = read_and_gen(f)
close(f)

out = parsed_args["out"]
if out == "std"
    println(Expr(:block, :(include("rts.jl")), big...))
else
    open(out, "w") do f
        println(f, Expr(:block, :(include("rts.jl")), big...))
    end
end


text:
using BenchmarkTools
using StaticArrays
# this is a slower but pretty solution. took inspiration from a reddit comment

const start = SVector{300}(map(str -> [parse(Int,str[l]) for l in findall(r"[[:digit:]]+", str)], readlines("inputs/6-10/06.txt"))[1])

function day6(start, days)
  c = zeros(Int, 9)
  for x in start
    c[x] += 1
  end
  for _ in 1:days-1
    c = circshift(c, -1)
    c[7] += c[end]
  end
  sum(c)
end

println("Day 6 Part 1: ",day6(start, 80))
println("Day 6 Part 2: ",day6(start, 256))

const A = SMatrix{9,9}([
 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  1  0  0  0  0
 0  0  0  0  0  1  0  0  0
 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
 1  0  0  0  0  0  0  0  0
])

function day6transition( days)
  c = zeros(Int, 9)
  for x in start
    c[x] += 1
  end
  StC = SVector{9}(c)
  (A^days) * StC |> sum
end

@btime day6transition( 256)


text:
using ImageFiltering: imfilter, reflect, FIR
using OffsetArrays: OffsetArray

function synthesize(syn::TransformSystem{NS}, y::AbstractArray) where {TF,D,NS<:AbstractNsolt{TF,D}}
    y = reshape_coefs(syn.shape, syn.operator, y)
    pvx = synthesize(syn.operator, y; syn.options...)
    polyphase2mdarray(pvx, decimations(syn.operator))
end

synthesize(fb::PolyphaseFB{TF,D}, pvy::PolyphaseVector{TY,D}; kwargs...) where {TF,TY,D} = PolyphaseVector(synthesize(fb, pvy.data, pvy.nBlocks; kwargs...), pvy.nBlocks)

function synthesize(nsolt::AbstractNsolt, py::AbstractMatrix, nBlocks::NTuple; kwargs...) where {TF,TX,D}
    sy = ishiftFilterSymmetry(nsolt, py)

    nShifts = fld.(size(sy, 2), nBlocks)
    irotatedimsfcns = ([ t->irotatedimspv(t, blk) for blk in nBlocks ]...,)
    ty = concatenateAtoms(nsolt, sy, nShifts, irotatedimsfcns; kwargs...)

    finalStep(nsolt, ty; kwargs...)
end

function finalStep(nsolt::RnsoltTypeI, py::AbstractMatrix; kwargs...)
    M = prod(nsolt.decimationFactor)
    fM, cM = fld(M, 2), cld(M, 2)

    tyup = nsolt.W0' * py[1:nsolt.nChannels[1], :]
    tylw = nsolt.U0' * py[(1:nsolt.nChannels[2]) .+ nsolt.nChannels[1], :]

    return nsolt.CJ' * [ tyup[1:cM, :]; tylw[1:fM, :] ]
end

function finalStep(nsolt::RnsoltTypeII, py::AbstractMatrix; kwargs...)
    M = prod(nsolt.decimationFactor)
    fM, cM = fld(M, 2), cld(M, 2)

    tyup = nsolt.W0' * py[1:nsolt.nChannels[1], :]
    tylw = nsolt.U0' * py[(1:nsolt.nChannels[2]) .+ nsolt.nChannels[1], :]

    return nsolt.CJ' * [ tyup[1:cM, :]; tylw[1:fM, :] ]
end

function concatenateAtoms(nsolt::RnsoltTypeI, px::AbstractMatrix, nShifts::NTuple, irotatedimsfcns::NTuple; border=:circular)
    px = UniformScaling{eltype(nsolt)}(1) * px
    for d in reverse(nsolt.perm)
        xu = @view px[1:nsolt.nChannels[1], :]
        xl = @view px[(1:nsolt.nChannels[2]) .+ nsolt.nChannels[1], :]

        params_d = (1:nsolt.nStages[d], nsolt.Udks[d])
        foreach(reverse.(params_d)...) do k, U
            xl .= U' * xl
            unnormalized_butterfly!(xu, xl)
            adjshiftcoefs!(Val(border), k, xu, xl, nShifts[d])
            half_butterfly!(xu, xl)
        end
        px = irotatedimsfcns[d](px)
    end
    return px
end

function concatenateAtoms(nsolt::RnsoltTypeII, px::AbstractMatrix, nShifts::NTuple, irotatedimsfcns; border=:circular)
    px = UniformScaling{eltype(nsolt)}(1) * px
    mnP, mxP = minmax(nsolt.nChannels...)
    for d in reverse(nsolt.perm)
        xu = @view px[1:mnP, :]
        xl = @view px[end-mnP+1:end, :]
        xum = @view px[1:mxP, :]
        xml = @view px[end-mxP+1:end, :]

        params_d = (1:nsolt.nStages[d], nsolt.Wdks[d], nsolt.Udks[d])
        foreach(reverse.(params_d)...) do k, W, U
            if nsolt.nChannels[1] < nsolt.nChannels[2]
                xml .= U' * xml
            else
                xum .= W' * xum
            end
            unnormalized_butterfly!(xu, xl)
            adjshiftcoefs!(Val(border), 2k, xum, xl, nShifts[d])
            half_butterfly!(xu, xl)
            if nsolt.nChannels[1] < nsolt.nChannels[2]
                xu .= W' * xu
            else
                xl .= U' * xl
            end
            unnormalized_butterfly!(xu, xl)
            adjshiftcoefs!(Val(border), 2k-1, xu, xml, nShifts[d])
            half_butterfly!(xu, xl)
        end
        px = irotatedimsfcns[d](px)
    end
    return px
end

ishiftFilterSymmetry(::Rnsolt, x::AbstractMatrix) = x

function finalStep(nsolt::CnsoltTypeI, px::AbstractMatrix; kwargs...)
    nsolt.FJ' * (nsolt.V0' * px)[1:prod(nsolt.decimationFactor), :]
end

function finalStep(nsolt::CnsoltTypeII, px::AbstractMatrix; kwargs...)
    nsolt.FJ' * (nsolt.V0' * px)[1:prod(nsolt.decimationFactor), :]
end

function concatenateAtoms(nsolt::CnsoltTypeI, px::AbstractMatrix, nShifts::NTuple, irotatedimsfcns::NTuple; border=:circular)
    px = UniformScaling{eltype(nsolt)}(1) * px
    hP = fld(nsolt.nChannels, 2)
    for d in reverse(nsolt.perm)
        xu = @view px[1:hP, :]
        xl = @view px[(1:hP) .+ hP, :]

        params_d = (1:nsolt.nStages[d], nsolt.Wdks[d], nsolt.Udks[d], nsolt.θdks[d])
        foreach(reverse.(params_d)...) do k, W, U, θ
            xu .= W' * xu
            xl .= U' * xl

            B = getMatrixB(nsolt.nChannels, θ)
            px .= B' * px
            adjshiftcoefs!(Val(border), k, xu, xl, nShifts[d])
            px .= B * px
        end
        px = irotatedimsfcns[d](px)
    end
    return px
end

function concatenateAtoms(nsolt::CnsoltTypeII, px::AbstractMatrix, nShifts::NTuple, irotatedimsfcns; border=:circular)
    px = UniformScaling{eltype(nsolt)}(1) * px
    fP, cP = fld(nsolt.nChannels, 2), cld(nsolt.nChannels, 2)
    for d in reverse(nsolt.perm)
        xW = @view px[1:fP, :]
        xU = @view px[(1:fP) .+ fP, :]
        xb = @view px[(fP+1):end, :]
        xŴ = @view px[1:cP, :]
        xÛ = @view px[(1:cP) .+ fP, :]
        xe  = @view px[1:end-1, :]

        params_d = (1:nsolt.nStages[d], nsolt.Wdks[d], nsolt.Udks[d], nsolt.θ1dks[d], nsolt.Ŵdks[d], nsolt.Ûdks[d], nsolt.θ2dks[d])
        foreach(reverse.(params_d)...) do k, W, U, θ1, Ŵ, Û, θ2
            xŴ .= Ŵ' * xŴ
            xÛ .= Û' * xÛ
            B2 = getMatrixB(nsolt.nChannels, θ2)
            xe .= B2' * xe
            adjshiftcoefs!(Val(border), 2k, xW, xb, nShifts[d])
            xe .= B2 * xe

            xW .= W' * xW
            xU .= U' * xU
            B1 = getMatrixB(nsolt.nChannels, θ1)
            xe .= B1' * xe
            adjshiftcoefs!(Val(border), 2k-1, xW, xU, nShifts[d])
            xe .= B1 * xe
        end
        px = irotatedimsfcns[d](px)
    end
    return px
end

ishiftFilterSymmetry(nsolt::Cnsolt, x::AbstractMatrix) = nsolt.Φ' * x

function synthesize(jts::JoinedTransformSystems{MS}, y::AbstractArray) where {MS<:Multiscale}
    subsynthesize(jts.shape, y, jts.transforms...)
end

# shape == Shapes.Separated
function subsynthesize(shape::Shapes.Separated, sy::AbstractArray, abop::AbstractOperator, args...)
    ya = [ subsynthesize(shape, sy[2:end], args...), sy[1]... ]
    synthesize(abop, ya)
end

subsynthesize(::Shapes.Separated, sy::AbstractArray, abop::AbstractOperator) = synthesize(abop, sy[1])

# shape == Shapes.Combined
function subsynthesize(shape::Shapes.Combined, sy::AbstractArray, abop::AbstractOperator, args...)
    rx = subsynthesize(shape, sy[2:end], args...)
    synthesize(abop, cat(rx, sy[1]; dims=ndims(sy[1]) ))
end

subsynthesize(::Shapes.Combined, sy::AbstractArray, abop::AbstractOperator) = synthesize(abop, sy[1])

# shape == Shapes.Vec
function subsynthesize(shape::Shapes.Vec, sy::AbstractArray, abop::AbstractOperator, args...)
    lny = nchannels(abop) * prod(fld.(abop.shape.insize, decimations(abop))) - prod(args[1].shape.insize)

    ya = [ vec(subsynthesize(shape, sy[lny+1:end], args...)); sy[1:lny] ]
    synthesize(abop, ya)
end

subsynthesize(::Shapes.Vec, sy::AbstractArray, abop::AbstractOperator) = synthesize(abop, sy)

function synthesize(pfs::ParallelFilters{TF,D}, y::AbstractArray; resource=CPU1(FIR())) where {TF,D}
    df = pfs.decimationFactor
    ord = pfs.polyphaseOrder

    nShift = df .* cld.(ord, 2) .+ 1
    region = UnitRange.(1 .- nShift, df .* (ord .+ 1) .- nShift)

    sxs = map(y, pfs.kernelspair[2]) do yp, sfp
        upimg = upsample(yp, df)
        ker = reflect(OffsetArray(sfp, region...))
        imfilter(resource, upimg, ker, "circular")
    end
    sum(sxs)
end

function synthesize(cs::TransformSystem{PF}, y::AbstractVector) where {TF,D,PF<:ParallelFilters{TF,D}}
    ty = reshape_coefs(cs.shape, cs.operator, y)
    synthesize(cs.operator, ty; cs.options...)
end

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
x, y = generators(MPoly{Float64}, :x, :y)

@test string((x+y)^3) == "x^3 + 3.0x^2*y + 3.0x*y^2 + y^3"
@test string((x*y)^(-2)) == "x^(-2)*y^(-2)"