text:
using PlotlyJS
using JuMP

"""
    plot_inventory(supply_chain, storage, product)

Plots the amount of inventory of a product on-hand at a storage location at the beginning of each period.
"""
function plot_inventory(supply_chain, storage, product)
    return Plot(scatter(;x=1:supply_chain.horizon, y=[get_inventory_at_end(supply_chain, storage, product, t) for t in 1:supply_chain.horizon], mode="lines+markers"))
end

function plot_networks(supply_chain; geography="usa", showlegend=true)
    return hcat([plot_network(supply_chain, p) for p in 1:supply_chain.horizon])
end

"""
    plot_network(supply_chain, period=1; geography="usa", showlegend=true)

Plots the nodes of the supply chain on a map.

The geography must be one of: "world" | "usa" | "europe" | "asia" | "africa" | "north america" | "south america".
"""
function plot_network(supply_chain, period=1; geography="usa", showlegend=true)
    traces = AbstractTrace[]
    first = true
    for storage in supply_chain.storages
        if is_opened(supply_chain, storage, period)
            push!(traces,
                scattergeo(;lat=[storage.location.latitude],
                            lon=[storage.location.longitude],
                            legendgroup="storage",
                            showlegend=first,
                            hoverinfo="text",
                            text="$(storage.name)",
                            name="storages",
                            mode="makers",
                            marker_symbol="square",
                            marker_size=10,
                            marker_color="blue")
            )
            first = false
        end
    end

    first = true
    for plant in supply_chain.plants
        if is_opened(supply_chain, plant, period)
            push!(traces,
                scattergeo(;lat=[plant.location.latitude],
                            lon=[plant.location.longitude],
                            legendgroup="plant",
                            showlegend=first,
                            hoverinfo="text",
                            text="$(plant.name)",
                            name="plants",
                            mode="makers",
                            marker_symbol="triangle-up",
                            marker_size=10,
                            marker_color="red")
            )
            first = false
        end
    end

    first = true
    for supplier in supply_chain.suppliers
        if sum(get_shipments(supply_chain, supplier, p, period) for p in supply_chain.products) > 1e-10
            push!(traces,
                scattergeo(;lat=[supplier.location.latitude],
                            lon=[supplier.location.longitude],
                            legendgroup="supplier",
                            showlegend=first,
                            hoverinfo="text",
                            text="$(supplier.name)",
                            name="suppliers",
                            mode="makers",
                            marker_symbol="circle",
                            marker_size=10,
                            marker_color="yellow")
            )
            first = false
        end
    end

    for (i, customer) in enumerate(supply_chain.customers)
        push!(traces,
            scattergeo(;lat=[customer.location.latitude],
                        lon=[customer.location.longitude],
                        legendgroup="customer",
                        showlegend=(i==1),
                        hoverinfo="text",
                        text="$(customer.name)",
                        name="customers",
                        mode="makers",
                        marker_symbol="circle",
                        marker_color="green",
                        marker_opacity=0.35)
        )
    end

    geo = attr(scope=geography,
                showland=true,)

    layout = Layout(;title="Supply chain network", showlegend=showlegend, geo=geo)
    return Plot(traces, layout)
end

"""
    plot_costs(supply_chain)

Plots the costs of operating the supply chain.
"""
function plot_costs(supply_chain)
    trace1 = bar(;x=1:supply_chain.horizon,
                 y=[value(supply_chain.optimization_model[:total_fixed_costs_per_period][t]) for t in 1:supply_chain.horizon],
                 name="Fixed Costs")
    trace2 = bar(;x=1:supply_chain.horizon,
                 y=[value(supply_chain.optimization_model[:total_transportation_costs_per_period][t]) for t in 1:supply_chain.horizon],
                 name="Transportation Costs")
    trace3 = bar(;x=1:supply_chain.horizon,
                 y=[value(supply_chain.optimization_model[:total_costs_per_period][t])  - value(supply_chain.optimization_model[:total_fixed_costs_per_period][t]) - value(supply_chain.optimization_model[:total_transportation_costs_per_period][t]) for t in 1:supply_chain.horizon],
                 name="Other Costs")

    layout = Layout(;barmode="stack",
                    xaxis_title="Period", xaxis_tick0=1, xaxis_dtick=1, xaxis_rangemode="nonnegative",
                    yaxis_title="Cost (\$)")
    Plot([trace1, trace2, trace3], layout)
end

"""
    plot_flows(supply_chain, period=1; geography="usa", showlegend=true)

Plots the flows of products in the supply chain.
"""
function plot_flows(supply_chain, period=1; geography="usa", showlegend=true)

    colors = ["red", "blue", "green", "orange", "black", "purple"]

    origin_colors = Dict{Node, String}()
    color_index = 0

    traces = AbstractTrace[]
    for l in sort(collect(supply_chain.lanes), by=l->string(l.origin.location.name, l.destination.location.name))
        if sum(get_shipments(supply_chain, l, p, period) for p in supply_chain.products) > 1e-10
            if l.origin in keys(origin_colors)
                color = origin_colors[l.origin]
            else
                color_index = mod1(color_index + 1, length(colors))
                color = colors[color_index]
                push!(origin_colors, l.origin => color)
            end

            push!(traces,
                scattergeo(;lat=[l.origin.location.latitude, l.destination.location.latitude],
                            lon=[l.origin.location.longitude, l.destination.location.longitude],
                            hoverinfo="text",
                            text="$(l.origin.location.name) - $(l.destination.location.name)",
                            name="$(l.origin.location.name) - $(l.destination.location.name)",
                            mode="lines",
                            line_color=color)
            )
        end
    end

    geo = attr(scope=geography,
                showland=true,)

    layout = Layout(;title="Supply chain flows", showlegend=showlegend, geo=geo)
    Plot(traces, layout)
end


text:
"""
    DistributedODEIntegrator <: AbstractODEIntegrator

A simplified variant of [`ODEIntegrator`](https://github.com/SciML/OrdinaryDiffEq.jl/blob/6ec5a55bda26efae596bf99bea1a1d729636f412/src/integrators/type.jl#L77-L123).

"""
mutable struct DistributedODEIntegrator{algType,uType,tType,pType} <: DiffEqBase.AbstractODEIntegrator{algType,true,uType,tType}
    prob
    alg::algType
    u::uType
    p::pType
    dt::tType
    t::tType
    tstop::tType
    tdir::tType
    step::Int
    stepstop::Int # -1
    adjustfinal::Bool
    callback::DiffEqBase.CallbackSet
    u_modified::Bool
    cache
    sol
end

# called by DiffEqBase.init and solve (see below)
function DiffEqBase.__init(
    prob::DiffEqBase.AbstractODEProblem,
    alg::DistributedODEAlgorithm,
    args...;
    dt,  # required
    stepstop=-1,
    adjustfinal=false,
    callback=nothing,
    save_func=(u,t,integrator)->copy(u),
    saveat=nothing,
    save_everystep=false,
    kwargs...)

    u = prob.u0
    p = prob.p
    t = prob.tspan[1]
    tstop = prob.tspan[2]
    tdir = sign(tstop - t)

    sol = DiffEqBase.build_solution(prob, alg, typeof(t)[], typeof(u)[])
    if isnothing(saveat)
        saveat = [tstop]
    elseif saveat isa Number
        saveat = t:saveat:tstop
    end
    saving_callback = NonInterpolatingSavingCallback(save_func, DiffEqCallbacks.SavedValues(sol.t, sol.u); saveat, save_everystep)
    callbackset = DiffEqBase.CallbackSet(callback, saving_callback)
    isempty(callbackset.continuous_callbacks) || error("Continuous callbacks are not supported")
    integrator = DistributedODEIntegrator(prob, alg, u, p, dt, t, tstop, tdir, 0, stepstop, adjustfinal, callbackset, false, cache(prob, alg; dt=dt, kwargs...), sol)
    DiffEqBase.initialize!(callbackset,u,t,integrator)
    return integrator
end

# called by DiffEqBase.solve
function DiffEqBase.__solve(
    prob::DiffEqBase.AbstractODEProblem,
    alg::DistributedODEAlgorithm,
    args...;
    kwargs...)

    integrator = DiffEqBase.__init(prob, alg, args...; kwargs...)
    DiffEqBase.solve!(integrator)
end

# either called directly (after init), or by solve (via __solve)
function DiffEqBase.solve!(integrator::DistributedODEIntegrator)
    while integrator.t < integrator.tstop
        if integrator.adjustfinal && integrator.t + integrator.dt > integrator.tstop
            adjust_dt!(integrator, integrator.tstop - integrator.t)
        end
        DiffEqBase.step!(integrator)

        if integrator.step == integrator.stepstop
            break
        end
    end

    if isdefined(DiffEqBase, :finalize!)
        DiffEqBase.finalize!(integrator.callback, integrator.u, integrator.t, integrator)
    end
    if isempty(integrator.sol.t)
        push!(integrator.sol.t, integrator.t)
        push!(integrator.sol.u, integrator.u)
    end
    return integrator.sol
end

# either called directly, or via solve!
function DiffEqBase.step!(integrator::DistributedODEIntegrator)
    step_u!(integrator) # solvers need to define this interface
    integrator.t += integrator.dt
    integrator.step += 1

    # apply callbacks
    discrete_callbacks = integrator.callback.discrete_callbacks
    for callback in discrete_callbacks
        if callback.condition(integrator.u,integrator.t,integrator)
            callback.affect!(integrator)
        end
    end
end

# solvers need to define this interface
step_u!(integrator) = step_u!(integrator, integrator.cache)

function adjust_dt!(integrator::DistributedODEIntegrator, dt)
    # TODO: figure out interface for recomputing other objects (linear operators, etc)
    integrator.dt = dt
end

# not sure what this should do?
# defined as default initialize: https://github.com/SciML/DiffEqBase.jl/blob/master/src/callbacks.jl#L3
DiffEqBase.u_modified!(i::DistributedODEIntegrator,bool) = nothing

# this is roughly based on SavingCallback from DiffEqCallbacks, except that it
# doesn't interpolate; instead it will save the first step after
function NonInterpolatingSavingCallback(save_func, saved_values::DiffEqCallbacks.SavedValues;
    saveat=nothing,
    save_everystep=false,
)
    if isnothing(saveat) && save_everystep == false
        error("saveat or save_everystep must be defined")
    end
    if save_everystep
        condition = (u, t, integrator) -> true
    else
        saveat = collect(saveat)
        condition = (u, t, integrator) -> begin
            cond = false
            while !isempty(saveat) && first(saveat) <= t
                cond = true
                popfirst!(saveat)
            end
            return cond
        end
    end
    function affect!(integrator)
        push!(saved_values.t, integrator.t)
        push!(saved_values.saveval, save_func(integrator.u, integrator.t, integrator))
    end
    initialize(cb, u, t, integrator) = condition(u,t,integrator) && affect!(integrator)
    finalize(cb, u, t, integrator) = !save_everystep && !isempty(saveat) && affect!(integrator)
    DiffEqBase.DiscreteCallback(condition, affect!; initialize=initialize, finalize=finalize)
end

ext: jl


text:
mutable struct QP{T <: AbstractFloat} <: AbstractLoss
    Q::Matrix{T}
    q::Vector{T}
    temp::Vector{T}

    function (::Type{QP})(d::Integer, t::Real, ::Type{T} = Float64) where T
        Q = randn(T,d,d)
        q = randn(T,d)
        F = qr(Q)
        temp = t .^ rand(d-2)
        D = Diagonal(vcat(1/t,temp,t))
        mul!(Q, D, Matrix(F.Q))
        lmul!(F.Q', Q)

        return new{T}(Q, q, copy(q))
    end
end

nfeatures(loss::QP) = size(loss.Q,2)

function loss!(loss::QP{T}, x::Vector{T}, g::Vector{T}) where T <: AbstractFloat
    mul!(loss.temp,loss.Q,x)
    result = convert(T, 0.5) * x⋅loss.temp + loss.q⋅x
    g .= loss.temp .+ loss.q
    return result
end


text:
export iauAb
"""
Apply aberration to transform natural direction into proper
direction.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   pnat    double[3]   natural direction to the source (unit vector)
   v       double[3]   observer barycentric velocity in units of c
   s       double      distance between the Sun and the observer (au)
   bm1     double      sqrt(1-|v|^2): reciprocal of Lorenz factor

Returned:
   ppr     double[3]   proper direction to source (unit vector)

Notes:

   1. The algorithm is based on Expr. (7.40) in the Explanatory
   Supplement (Urban & Seidelmann 2013), but with the following
   changes:
      - Rigorous rather than approximate normalization is applied.
      -The gravitational potential term from Expr. (7) in
      Klioner (2003) is added, taking into account only the Sun's
      contribution.  This has a maximum effect of about
      0.4 microarcsecond.

   2. In almost all cases, the maximum accuracy will be limited by the
   supplied velocity.  For example, if the SOFA iauEpv00 function is
   used, errors of up to 5 microarcseconds could occur.

  References:
   - Urban, S. & Seidelmann, P. K. (eds), Explanatory Supplement to
   the Astronomical Almanac, 3rd ed., University Science Books
   (2013).

   - Klioner, Sergei A., "A practical relativistic model for micro-
   arcsecond astrometry in space", Astr. J. 125, 1580-1597 (2003).

Called:
   iauPdp       scalar product of two p-vectors

This revision:   2013 October 9

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
function iauAb(pnat::AbstractVector{<:Real}, v::AbstractVector{<:Real}, s::Real, bm1::Float64)
   ppr = zeros(Float64, 3)

   ccall((:iauAb, libsofa_c), Cvoid,
         (Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, Ptr{Cdouble}),
         convert(Array{Float64, 1}, pnat),
         convert(Array{Float64, 1}, v),
         convert(Float64, s),
         convert(Float64, bm1), ppr)

   return SVector{3}(ppr)
end


text:
TOL = 0.5e-8

function newton(f::Function, ∂f::Function, x₀::Real;
                ε::Real = TOL, iₘₐₓ::Int = 100)
    xᵢ = x₀
    εᵢ = xᵢ
    i = 0
    while εᵢ >= ε && i < iₘₐₓ
        x₋, xᵢ = xᵢ, xᵢ - f(xᵢ) / ∂f(xᵢ)
        εᵢ = abs( xᵢ - x₋ )
        n += 1
    end
    xᵢ, εᵢ, i
end


text:
"""
    AbstractRiemannianDiffBackend

An abstract type for backends for differentiation.
"""
abstract type AbstractRiemannianDiffBackend end

@doc raw"""
    differential(M::AbstractManifold, f, t::Real, backend::AbstractDiffBackend)
    differential!(M::AbstractManifold, f, X, t::Real, backend::AbstractDiffBackend)

Compute the Riemannian differential of a curve $f: ℝ\to M$ on a manifold `M`
represented by function `f` at time `t` using the given backend.
It is calculated as the tangent vector equal to $\mathrm{d}f_t(t)[1]$.

The mutating variant computes the differential in place of `X`.
"""
differential(::AbstractManifold, ::Any, ::Real, ::AbstractRiemannianDiffBackend)

@doc raw"""
    gradient(M::AbstractManifold, f, p, backend::AbstractRiemannianDiffBackend)
    gradient!(M::AbstractManifold, f, X, p, backend::AbstractRiemannianDiffBackend)

Compute the Riemannian gradient ``∇f(p)`` of a real-valued function ``f:\mathcal M \to ℝ``
at point `p` on the manifold `M` using the specified [`AbstractRiemannianDiffBackend`](@ref).

The mutating variant computes the gradient in place of `X`.
"""
gradient(::AbstractManifold, ::Any, ::Any, ::AbstractRiemannianDiffBackend)

function differential!(
    M::AbstractManifold,
    f::Any,
    X,
    t,
    backend::AbstractRiemannianDiffBackend,
)
    return copyto!(X, differential(M, f, t, backend))
end

function gradient!(M::AbstractManifold, f, X, p, backend::AbstractRiemannianDiffBackend)
    return copyto!(X, gradient(M, f, p, backend))
end

@doc raw"""
    TangentDiffBackend <: AbstractRiemannianDiffBackend

A backend that uses a tangent space and a basis therein to derive an
intrinsic differentiation scheme.

Since it works in a tangent space, methods might require a retraction and an
inverse retraction as well as a basis.

In the tangent space itself, this backend then employs an (Euclidean)
[`AbstractDiffBackend`](@ref)

# Constructor

    TangentDiffBackend(diff_backend)

where `diff_backend` is an [`AbstractDiffBackend`](@ref) to be used on the tangent space.

With the keyword arguments

* `retraction` an [`AbstractRetractionMethod`](@ref) ([`ExponentialRetraction`](@ref) by default)
* `inverse_retraction` an [`AbstractInverseRetractionMethod`](@ref) ([`LogarithmicInverseRetraction`](@ref) by default)
* `basis` an [`AbstractBasis`](@ref) ([`DefaultOrthogonalBasis`](@ref) by default)
"""
struct TangentDiffBackend{
    TAD<:AbstractDiffBackend,
    TR<:AbstractRetractionMethod,
    TIR<:AbstractInverseRetractionMethod,
    TB<:AbstractBasis,
} <: AbstractRiemannianDiffBackend
    diff_backend::TAD
    retraction::TR
    inverse_retraction::TIR
    basis::TB
end
function TangentDiffBackend(
    diff_backend::TAD;
    retraction::TR=ExponentialRetraction(),
    inverse_retraction::TIR=LogarithmicInverseRetraction(),
    basis::TB=DefaultOrthonormalBasis(),
) where {
    TAD<:AbstractDiffBackend,
    TR<:AbstractRetractionMethod,
    TIR<:AbstractInverseRetractionMethod,
    TB<:AbstractBasis,
}
    return TangentDiffBackend{TAD,TR,TIR,TB}(
        diff_backend,
        retraction,
        inverse_retraction,
        basis,
    )
end

function differential(M::AbstractManifold, f, t::Real, backend::TangentDiffBackend)
    p = f(t)
    onb_coords = _derivative(zero(number_eltype(p)), backend.diff_backend) do h
        return get_coordinates(
            M,
            p,
            inverse_retract(M, p, f(t + h), backend.inverse_retraction),
            backend.basis,
        )
    end
    return get_vector(M, p, onb_coords, backend.basis)
end

function differential!(M::AbstractManifold, f, X, t::Real, backend::TangentDiffBackend)
    p = f(t)
    onb_coords = _derivative(zero(number_eltype(p)), backend.diff_backend) do h
        return get_coordinates(
            M,
            p,
            inverse_retract(M, p, f(t + h), backend.inverse_retraction),
            backend.basis,
        )
    end
    return get_vector!(M, X, p, onb_coords, backend.basis)
end

@doc raw"""
    gradient(M, f, p, backend::TangentDiffBackend)

This method uses the internal `backend.diff_backend` (Euclidean) on the function

```math
    f(\retr_p(\cdot))
```

which is given on the tangent space. In detail, the gradient can be written in
terms of the `backend.basis`. We illustrate it here for an [`AbstractOrthonormalBasis`](@ref),
since that simplifies notations:

```math
\operatorname{grad}f(p) = \operatorname{grad}f(p) = \sum_{i=1}^{d} g_p(\operatorname{grad}f(p),X_i)X_i
	= \sum_{i=1}^{d} Df(p)[X_i]X_i
```

where the last equality is due to the definition of the gradient as the Riesz representer of the differential.

If the backend is a forward (or backward) finite difference, these coefficients in the sum can be approximates as

```math
DF(p)[Y] ≈ \frac{1}{h}\bigl( f(\exp_p(hY)) - f(p) \bigr)
```
writing ``p=\exp_p(0)`` we see that this is a finite difference of ``f\circ\exp_p``, i.e. for
a function on the tangent space, so we can also use other (Euclidean) backends
"""
function gradient(M::AbstractManifold, f, p, backend::TangentDiffBackend)
    X = get_coordinates(M, p, zero_vector(M, p), backend.basis)
    onb_coords = _gradient(X, backend.diff_backend) do Y
        return f(retract(M, p, get_vector(M, p, Y, backend.basis), backend.retraction))
    end
    return get_vector(M, p, onb_coords, backend.basis)
end

function gradient!(M::AbstractManifold, f, X, p, backend::TangentDiffBackend)
    X2 = get_coordinates(M, p, zero_vector(M, p), backend.basis)
    onb_coords = _gradient(X2, backend.diff_backend) do Y
        return f(retract(M, p, get_vector(M, p, Y, backend.basis), backend.retraction))
    end
    return get_vector!(M, X, p, onb_coords, backend.basis)
end

@doc raw"""
    RiemannianProjectionBackend <: AbstractRiemannianDiffBackend

This backend computes the differentiation in the embedding, which is currently limited
to the gradient. Let ``mathcal M`` denote a manifold embedded in some ``R^m``, where ``m``
is usually (much) larger than the manifold dimension.
Then we require three tools

* A function ``f̃: ℝ^m → ℝ`` such that its restriction to the manifold yields the cost
  function ``f`` of interest.
* A [`project`](@ref) function to project tangent vectors from the embedding (at ``T_pℝ^m``)
  back onto the tangent space ``T_p\mathcal M``. This also includes possible changes
  of the representation of the tangent vector (e.g. in the Lie algebra or in a different data format).
* A [`change_representer`](@ref) for non-isometrically embedded manifolds,
  i.e. where the tangent space ``T_p\mathcal M`` of the manifold does not inherit
  the inner product from restriction of the inner product from the tangent space ``T_pℝ^m``
  of the embedding

For more details see [^AbsilMahonySepulchre2008], Section 3.6.1 for a derivation on submanifolds.

[^AbsilMahonySepulchre2008]:
    > Absil, P.-A., Mahony, R. and Sepulchre R.,
    > _Optimization Algorithms on Matrix Manifolds_
    > Princeton University Press, 2008,
    > doi: [10.1515/9781400830244](https://doi.org/10.1515/9781400830244)
    > [open access](http://press.princeton.edu/chapters/absil/)
"""
struct RiemannianProjectionBackend{TADBackend<:AbstractDiffBackend} <:
       AbstractRiemannianDiffBackend
    diff_backend::TADBackend
end

function gradient(M::AbstractManifold, f, p, backend::RiemannianProjectionBackend)
    amb_grad = _gradient(f, p, backend.diff_backend)
    return change_representer(M, EuclideanMetric(), p, project(M, p, amb_grad))
end

function gradient!(M::AbstractManifold, f, X, p, backend::RiemannianProjectionBackend)
    amb_grad = embed(M, p, X)
    _gradient!(f, amb_grad, p, backend.diff_backend)
    project!(M, X, p, amb_grad)
    return change_representer!(M, X, EuclideanMetric(), p, X)
end


text:
include("../base.jl")
import Random

f(x, y) = (y/5.5)+x^2

println("hola")
println("hola")
x = collect(-5:0.5:5)
y = collect(0:0.5:10)
nvalues = length(x)
objs = Array{Number}(undef, nvalues)
for i=1:nvalues
    objs[i] = f(x[i], y[i])
end

setIndividualType(GEPGenotype)
setRandomSeed(34)
setGPExperimentInfo()
setStopCondition(maxIterations=50)
setEvaluator([FitnessFunction(compareFunctions, objs, x, y, weight=-1)])
setGenerator(rampedHalfHalfGenerator, popSize = 50, generateOneByOne = false)
setSelector(tournamentSelector, 4)
setCrossoverOperator(kPointRecombinationCross, 2, probability=0.9)
setMutationOperator(geneMutation, probability=0.2)
setReplacementOperator(replaceAllPopulation, eliteSize=5)
setExperimentSummary(batchSize=5, printDuringExperiment=true)

@time runGenJ(verbose=false)

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
using FactorizationMachines.Methods: MethodParams, SGDMethod, sgd, sgd_train!
using FactorizationMachines.Evaluators: Evaluator
using FactorizationMachines.Tasks: TaskParams
using FactorizationMachines.Models: ModelParams
using FactorizationMachines.Predictors: FMPredictor

function train(X::FMMatrix, y::Vector{FMFloat};
        method::SGDMethod         = Methods.sgd(alpha = 0.01, num_epochs = 100, reg0 = .0, regv = .0, regw = .0),
        evaluator::Evaluator      = Evaluators.rmse(),
        task_params::TaskParams   = Tasks.regression(),
        model_params::ModelParams = Models.gauss(k0 = true, k1 = true, num_factors = 8, mean = .0, stddev = .01))

    model = @time Models.init(model_params, X, y)
    task = Tasks.init(task_params, X, y)
    predictor = FMPredictor(task, model)

    # Train the predictor using SGD
    sgd_train!(method, evaluator, predictor, X, y)

    predictor
end


text:
export CompactSet,
       NonCompactSet

"""
    CompactSet

An alias for compact set types.

### Notes

Most lazy operations are not captured by this alias because whether their result
is compact or not depends on the argument(s).
"""
const CompactSet = Union{
    AbstractCentrallySymmetric,
    AbstractPolytope,
    EmptySet
    }

"""
    NonCompactSet

An alias for non-compact set types.

### Notes

Most lazy operations are not captured by this alias because whether their result
is non-compact or not depends on the argument(s).
"""
const NonCompactSet = Union{
    HalfSpace,
    Hyperplane,
    Line2D,
    Universe
    }


text:
function TestHill()

close all
nx=20;
ny=1;
lx=2000;
ly=100;
x0=-1000;
y0=0;
Param.Grav=9.81;
Param.Coriolis=false;
Param.HyperDM=1.e15*(nx/30)^3.2;
Param.HyperD=1.e15*(nx/30)^3.2;
Param.Flat=false;

Boundary.WE='Period';
Boundary.BT='Period';
Param.Grid=CartGridCartGrid(nx,ny,lx,ly,x0,y0,OrientFaceCart,Boundary,Param);

OrdPoly=3;
CG.OrdPoly=OrdPoly;
[CG.w,CG.xw]=GaussLobattoQuad(CG.OrdPoly);
[CG.wX,CG.xwX]=GaussLobattoQuad(CG.OrdPoly);
[CG.wY,CG.xwY]=GaussLobattoQuad(CG.OrdPoly);
[CG.DW,CG.DS]=DerivativeMatrixSingle(CG);

Param.J=zeros(OrdPoly+1,OrdPoly+1,Param.Grid.NumFaces);
Param.X=zeros(OrdPoly+1,OrdPoly+1,3,Param.Grid.NumFaces);
Param.dXdx=zeros(OrdPoly+1,OrdPoly+1,2,2,Param.Grid.NumFaces);
Param.lat=zeros(OrdPoly+1,OrdPoly+1,Param.Grid.NumFaces);
% Param.J1=zeros(OrdPoly+1,OrdPoly+1,Param.Grid.NumFaces);
% Param.X1=zeros(OrdPoly+1,OrdPoly+1,3,Param.Grid.NumFaces);
% Param.dXdx1=zeros(OrdPoly+1,OrdPoly+1,2,2,Param.Grid.NumFaces);
% Param.lat1=zeros(OrdPoly+1,OrdPoly+1,Param.Grid.NumFaces);
for iF=1:Param.Grid.NumFaces
%   [Param.X1(:,:,:,iF),Param.J1(:,:,iF),Param.dXdx1(:,:,:,:,iF),Param.lat1(:,:,iF)]=...
%     JacobiSphereDG(CG,Param.Grid.Faces(iF),Param.Grid);
  for j=1:OrdPoly+1
    for i=1:OrdPoly+1
      [Param.X(i,j,:,iF),Param.J(i,j,iF),Param.dXdx(i,j,:,:,iF),Param.lat(i,j,iF)]...
        =JacobiSphere2(CG.xw(i),CG.xw(j),Param.Grid.Faces(iF),Param.Grid);
    end
  end
end

%Setting initial Values
Profh='Hill';
ProfVel='Hill';
Param.ProfVel=ProfVel;
Param.Profh=Profh;

[CG.Faces,CG.NumG,CG.NumI]=NumberingFemCG(Param.Grid,OrdPoly);
[CG.M]=MassCG(CG,Param);

U=zeros(CG.NumG,3);
U(:,1)=Project(@fh,CG,@JacobiSphere2,Param);
Param.hS=Project(@fhs,CG,@JacobiSphere2,Param);
[U(:,2),U(:,3)]=ProjectVec(@fVel,CG,@JacobiSphere2,Param);
fig=1;
fig=PlotCG(U(:,2),CG,@JacobiSphere2,Param,fig);

CFL=0.125;
dtau=100;
time=0;

IntMethod='RungeKutta';
nIter=5*3600*24/100;
v = VideoWriter ('Hill.avi');
open (v);
switch IntMethod
  case 'RungeKutta'
    for i=1:nIter
      F=FcnShallowCurl(U,CG,Param);
      UNeu=U+1/3*dtau*F;
      F=FcnShallowCurl(UNeu,CG,Param);
      UNeu=U+1/2*dtau*F;
      F=FcnShallowCurl(UNeu,CG,Param);
      U=U+dtau*F;
      time=time+dtau;
      if mod(i,50)==0
        PlotCG(U(:,2),CG,@JacobiSphere2,Param,fig);
        frame=getframe(gcf);
        writeVideo (v, frame);
      end
    end
end
%fig=PlotCG(U(:,1),CG,@JacobiSphere2,Param,fig);
close(v)
end


text:
module LogArrayFunctions

export LogArray
type LogArray <: AbstractArray
    content
    function LogArray(content)
        new(content)
    end
end

import Base.display
function display(L::LogArray)
    return display(L.content)
end

import Base.size
function size(A::LogArray,I::Int64)
    return size(A.content,I)
end

function size(A::LogArray)
    return size(A.content)
end

import Base.permutedims
function permutedims(A::LogArray,I::Array{Int64,1})
    return permutedims(A.content,I)
end

function *(A::LogArray,B::LogArray)
    logprefactor=maximum([maximum(A.content[:]) maximum(B.content[:])]);
    AA=exp.(A.content-logprefactor)
    BB=exp.(B.content-logprefactor)
    L=LogArray(1);
    L.content=2*logprefactor+log.(AA*BB)
    return L
end

function +(A::LogArray,B::LogArray)
    logprefactor=maximum([maximum(A.content[:]) maximum(B.content[:])]);
    AA=exp.(A.content-logprefactor)
    BB=exp.(B.content-logprefactor)
    L=LogArray(1);
    L.content=logprefactor+log.(AA+BB)
    return L
end

import Base.sum
function sum(A::LogArray,I...)
    logprefactor=maximum(A.content[:])
    AA=exp.(A.content-logprefactor)
    L=LogArray(1)
    L.content=logprefactor+log.(sum(AA,I...))
    return L
end

import Base.ndims
function ndims(A::LogArray)
    return ndims(A.content)
end

import Base.getindex
function getindex(A::LogArray, I...)
    tmp=getindex(A.content,I...)
    return LogArray(tmp)
end

import Base.transpose
function transpose(A::LogArray)
    return LogArray((A.content)')
end

function .*(A::LogArray,B::LogArray)
    logprefactor=maximum([maximum(A.content[:]) maximum(B.content[:])]);
    AA=exp.(A.content-logprefactor)
    BB=exp.(B.content-logprefactor)
    L=LogArray(1)
    L.content=2*logprefactor+log.(AA.*BB)
    return L
end

end #end module


text:
# using Base: AlwaysLockedST
function Benchmark.run(env::MuEnv, eval::AlphaZero.Benchmark.Evaluation)
  # net() = Network.copy(env.bestnns, on_gpu=eval.sim.use_gpu, test_mode=true)
  device = eval.sim.use_gpu ? Flux.gpu : Flux.cpu
  function net()
    nns = deepcopy(env.bestnns) |> device |> Flux.testmode!
    return (InitialOracle(nns), RecurrentOracle(nns))
  end
  if isa(eval, Benchmark.Single)
    simulator = Simulator(net, record_trace) do net
      Benchmark.instantiate(eval.player, env.gspec, net)
    end
  else
    @assert isa(eval, Benchmark.Duel)
    simulator = Simulator(net, record_trace) do net
      player = Benchmark.instantiate(eval.player, env.gspec, net)
      baseline = Benchmark.instantiate(eval.baseline, env.gspec, net)
      return TwoPlayers(player, baseline)
    end
  end
  samples, elapsed = @timed simulate(
    simulator, env.gspec, eval.sim)
  # gamma = env.params.self_play.mcts.gamma
  rewards, redundancy = rewards_and_redundancy(samples, gamma=1.0)
  opponent_won, draw, mu_won = count_wins(rewards) ./ length(rewards)
  # @info "Benchmark" white_lost draw white_won time=elapsed
  return (; opponent_won, draw, mu_won, time_benchmarked=elapsed)
#   return Report.Evaluation(
#     name(eval), mean(rewards), redundancy, rewards, nothing, elapsed)
end

struct Mu <: Benchmark.Player
  params :: MctsParams
  # make_oracles # Function make_oracles()
end

# name(::Mu) = "MuZero"

function Benchmark.instantiate(p::Mu, gspec::AbstractGameSpec, nns)
  return MuPlayer(nns, p.params)
end

function run_duel(env, benchmark)
  # report = []
  report = Dict()
  for (name,duel) in pairs(benchmark)
    outcome = Benchmark.run(env, duel)
    # push!(report, (;name=outcome))
    report[name] = outcome
  end
  return report
end
# benchmark
# env
# run_duel(env, benchmark)

# benchmark_sim = SimParams(
#     num_games=400,
#     num_workers=1,
#     batch_size=4,
#     use_gpu=false,
#     reset_every=1,
#     flip_probability=0.5,
#     alternate_colors=true)

# benchmark = [
#   Benchmark.Duel(
#     Mu(self_play.mcts),
#     Benchmark.MctsRollouts(self_play.mcts),
#     benchmark_sim)]

function count_wins(rewards)
  lost, draw, won = (count(==(i), rewards) for i in -1:1)
  @assert sum([lost, draw, won]) == length(rewards)
  # @info "Benchmark" lost draw won
  return lost, draw, won
end


text:
import GLFW
using Base.Test

window = GLFW.CreateWindow(800, 600, "InexactError")
@test !GLFW.WindowShouldClose(window)
# If you call the julia function GLFW.SetWindowShouldClose(window,
# truefalse), then only 0 or 1 can be passed. However, at least on
# some platforms (Linux), clicking on the window's close icon can send
# other values (like 189). Make sure such values don't cause trouble for
# GLFW.WindowShouldClose.
ccall( (:glfwSetWindowShouldClose, GLFW.lib), Void, (GLFW.WindowHandle, Cint), window, 189)
@test GLFW.WindowShouldClose(window)

# Test for a segfault triggered by a scheduled task
using ModernGL, GeometryTypes
window = GLFW.CreateWindow(1024, 768, "Segfault")
GLFW.MakeContextCurrent(window)

vao = Ref(GLuint(0))
glGenVertexArrays(1, vao)
glBindVertexArray(vao[])

vbo = Ref(GLuint(0))
glGenBuffers(1, vbo)
glBindBuffer(GL_ARRAY_BUFFER, vbo[])
vertices = Point3f0[(-1,-1,0), (1,-1,0), (0, 1, 1)]
glBufferData(GL_ARRAY_BUFFER, sizeof(vertices), vertices, GL_STATIC_DRAW)

glEnableVertexAttribArray(0)
glBindBuffer(GL_ARRAY_BUFFER, vbo[])
glVertexAttribPointer(0, length(vertices), GL_FLOAT, GL_FALSE, 0, C_NULL)

function render()
    glDrawArrays(GL_TRIANGLES, 0, length(vertices))
    GLFW.SwapBuffers(window)
    GLFW.PollEvents()
end

schedule(Task(render))

GLFW.DestroyWindow(window)

sleep(0.25)


text:
using Documenter
using Descartes

makedocs(
    sitename = "Descartes",
    format = Documenter.HTML(),
    modules = [Descartes],
    authors = "Steve Kelly",
    pages = ["Primitives" => "primitives.md",
             "Design" => "design.md"]
)

deploydocs(
    repo = "github.com/juliageometry/Descartes.jl.git"
)


text:
using ContinuousPOMDPTreeSearchExperiments
using ParticleFilters
using ARDESPOT
using BasicPOMCP
using POMCPOW
using POMDPs
using DiscreteValueIteration
using QMDP
using MCTS
using SubHunt
using POMDPToolbox
using DataFrames
using CSV

file_contents = readstring(@__FILE__())

pomdp = SubHuntPOMDP()

vs = ValueIterationSolver()
if !isdefined(:vp) || vp.mdp != pomdp
    vp = solve(vs, pomdp, verbose=true)
end
qs = QMDPSolver()
qp = QMDP.create_policy(qs, pomdp)
qp.alphas[:] = vp.qmat

@show max_depth = 20
@show N = 1000

alldata = DataFrame()
for max_time in logspace(-2, 1, 7)
    @show max_time
    solvers = Dict{String, Union{Solver,Policy}}(
        "qmdp" => qp,
        "ping_first" => PingFirst(qp),

        "pomcpow" => begin
            rng = MersenneTwister(13)
            solver = POMCPOWSolver(tree_queries=10_000_000_000,
                                   criterion=MaxUCB(17.0),
                                   final_criterion=MaxTries(),
                                   max_depth=max_depth,
                                   max_time=max_time,
                                   enable_action_pw=false,
                                   k_observation=6.0,
                                   alpha_observation=1/30,
                                   estimate_value=FOValue(vp),
                                   check_repeat_obs=false,
                                   default_action=ReportWhenUsed(qp),
                                   rng=rng
                                  )
        end,

        "pft_20" => begin
            rng = MersenneTwister(13)
            m = 20
            node_updater = ObsAdaptiveParticleFilter(deepcopy(pomdp),
                                               LowVarianceResampler(m),
                                               0.1, rng)
            solver = DPWSolver(n_iterations=typemax(Int),
                               exploration_constant=100.0,
                               depth=max_depth,
                               max_time=max_time,
                               k_state=2.0,
                               alpha_state=1/10,
                               check_repeat_state=false,
                               check_repeat_action=false,
                               estimate_value=RolloutEstimator(qp),
                               enable_action_pw=false,
                               # default_action=ReportWhenUsed(qp),
                               rng=rng
                              )
            belief_mdp = GenerativeBeliefMDP(deepcopy(pomdp), node_updater)
            solve(solver, belief_mdp)
        end,

        "pft_10" => begin
            rng = MersenneTwister(13)
            m = 10
            node_updater = ObsAdaptiveParticleFilter(deepcopy(pomdp),
                                               LowVarianceResampler(m),
                                               0.1, rng)
            solver = DPWSolver(n_iterations=typemax(Int),
                               exploration_constant=100.0,
                               depth=max_depth,
                               max_time=max_time,
                               k_state=2.0,
                               alpha_state=1/10,
                               check_repeat_state=false,
                               check_repeat_action=false,
                               estimate_value=RolloutEstimator(qp),
                               enable_action_pw=false,
                               # default_action=ReportWhenUsed(qp),
                               rng=rng
                              )
            belief_mdp = GenerativeBeliefMDP(deepcopy(pomdp), node_updater)
            solve(solver, belief_mdp)
        end,

    )

    for (k, solver) in solvers
    # test = ["qmdp", "pomcpow", "pft"]
    # for (k, solver) in [(s, solvers[s]) for s in test]
        @show k
        if isa(solver, Solver)
            planner = solve(solver, pomdp)
        else
            planner = solver
        end
        sims = []
        for i in 1:N
            srand(planner, i+50_000)
            filter = SIRParticleFilter(deepcopy(pomdp), 100_000, rng=MersenneTwister(i+90_000))

            md = Dict(:solver=>k, :i=>i, :max_time=>max_time)
            sim = Sim(deepcopy(pomdp),
                      planner,
                      filter,
                      rng=MersenneTwister(i+70_000),
                      max_steps=100,
                      metadata=md
                     )

            push!(sims, sim)
        end

        data = run_parallel(sims)
        # data = run(sims)

        rs = data[:reward]
        println(@sprintf("reward: %6.3f ± %6.3f", mean(rs), std(rs)/sqrt(length(rs))))
        if isempty(alldata)
            alldata = data
        else
            alldata = vcat(alldata, data)
        end
    end
end

datestring = Dates.format(now(), "E_d_u_HH_MM")
copyname = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "icaps_2018", "data", "subhunt_time_$(datestring).jl")
write(copyname, file_contents)
filename = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "icaps_2018", "data", "subhunt_time_$(datestring).csv")
println("saving to $filename...")
CSV.write(filename, alldata)
println("done.")


text:
using TSne #tSNE
using MultivariateStats #PCA

data_path = ARGS[1] # ARGS[1] must be a directory

usernames = [];
validFiles = [];

reduce_dim = 40

function getUsername(filename)
    m = match(r"^(.*)\.txt$", filename)
    (s,) = size(m.captures)
    if s < 1
        return ""
    else
        return String(m.captures[1])
    end
end

# n^2
function closestPoints(X::Array{Float64, 2}, i, numPoints)

    n = size(X, 1)

    distances=[]
    for j in 1:n
        push!(distances, norm(dataPoints[j, :] - dataPoints[i, :]))
    end

    return X[sortperm(distances)[1:numPoints], :], distances, sortperm(distances)
end

function densityCluster(dataPoints)
    n = size(dataPoints, 1)

    MINPTS = 5
    RADIUS = sqrt(n * 3 / 10)

    y = zeros(n)

    currentClusterNum = 1

    for i in 1:n
        if y[i] != 0
            continue
        end

        distances=[]
        for j in 1:n
            if i !== j
                push!(distances, norm(dataPoints[j, :] - dataPoints[i, :]))
            end
        end

        if sort(distances)[MINPTS] <= RADIUS
            println("calling expandCluster on point ", i)
            currentClusterNum = expandCluster!(y, dataPoints, i, currentClusterNum, RADIUS, MINPTS)
        end
    end

    return y
end

function expandCluster!(y, dataPoints, idx, currentClusterNum, radius, minPts)

    queue = [idx]

    visited = Set()

    n = size(dataPoints, 1)

    count = 0

    for o in queue
        if y[o] == 0
            y[o] = currentClusterNum
            count += 1
        else
            continue
        end

        (sortedPoints, distances, perm) = closestPoints(dataPoints, idx, n-1)

        for j in 1:n
            if distances[j] <= radius && y[perm[j]] == 0
                push!(queue, perm[j])
                push!(visited, perm[j])
            end
        end
    end

    if count < minPts
        for pt in visited
            y[pt] = 0
        end
        return currentClusterNum
    end

    return currentClusterNum + 1
end

for (_, _, fs) in walkdir(data_path)
    usernames = filter(u -> u != "", map(getUsername, fs))
    validFiles = map(u -> joinpath(data_path, "$(u).txt"), usernames);
end

followed = Dict()

mutable struct FollowedAccount
    numFollowing::Int
    following::Array{String}
end

for (n, file) in enumerate(validFiles)
    strm = open(file, "r")
    followedAccs = readlines(strm)
    for acc in followedAccs
        if haskey(followed, acc)
            followed[acc].numFollowing += 1
            push!(followed[acc].following, usernames[n])
        else
            followed[acc] = FollowedAccount(1, [usernames[n]])
        end
    end
    close(strm)
end

THRESHOLD = 1 # with THRESHOLD=1, there is no threshold. Experiments have been done setting this to as high as 15.

followingToCheck = []

for (accname, accdata) in followed
    if accdata.numFollowing >= THRESHOLD
        push!(followingToCheck, accname)
    end
end

t_noneinthr = []
for username in usernames
    inAcc = map(acc -> in(username, followed[acc].following), followingToCheck)
    if !any(inAcc)
        push!(t_noneinthr, username)
    end
end

usernamesFiltered = filter(name -> !in(name, t_noneinthr), usernames)

(dim,) = size(followingToCheck)

(n,) = size(usernamesFiltered)
dataPoints = zeros(n,dim)

for (usernameIdx, username) in enumerate(usernamesFiltered)

    for (followedIdx, followedAcc) in enumerate(followingToCheck)
        if in(username, followed[followedAcc].following)
            dataPoints[usernameIdx, followedIdx] = 1
        end
    end
end

function rescale(A, dim::Integer=1)
    res = A .- mean(A, dim)
    res ./= map!(x -> x > 0.0 ? x : 1.0, std(A, dim))
    return res
end

reshapedData = convert(Matrix{Float64}, dataPoints)
X = rescale(reshapedData, 1)
M = fit(PCA, X'; maxoutdim = reduce_dim)

y = densityCluster(transform(M, X')')
@show(y)

for value in unique(y)
    @printf("Cluster %d: %s\n", value, usernamesFiltered[find(y.==value)])
end


text:
module VortexFilaments

using Reexport
using DocStringExtensions

include("vortexfilament.jl")

end # module


text:
#
# Data for the simulation
#

using Parameters

@with_kw struct Data

    ε :: Float64 = 5.0
    σ :: Float64 = 0.5
    T :: Float64 = 298.15

    cutoff :: Float64 = 2.0
    bdim :: Float64 = 100.0

    N :: Int64 = 10_000

    # Data for computation
    nc :: Int64 = (bdim % cutoff == 0) ? trunc(Int64, bdim/cutoff) : trunc(Int64, bdim/cutoff + 1)
    ε4 :: Float64 = 4*ε
    ε24 :: Float64 = 24*ε
    σ6 :: Float64 = σ^6
    σ12 :: Float64 = σ^12
    cutoff2 :: Float64 = cutoff^2

end


text:
# module Hangul

# 첫가끝
const Initial = :Initial
const Medial  = :Medial
const Final   = :Final

# Hangul Jamo
primitive type YetJamo{T} <: AbstractChar 32 end

struct YetJamoError <: Exception
    message::String
end

struct YetSyllable
    initial::Union{Nothing,YetJamo{Initial}}
    medial::Union{Nothing,YetJamo{Medial}}
    final::Union{Nothing,YetJamo{Final}}
end

YetJamo{T}(c::UInt32) where T = reinterpret(YetJamo{T}, c)
Base.codepoint(c::YetJamo{T}) where T = reinterpret(UInt32, c)

YetSyllable(i::YetJamo{Initial}) = YetSyllable(i, nothing, nothing)
YetSyllable(m::YetJamo{Medial}) = YetSyllable(nothing, m, nothing)
YetSyllable(f::YetJamo{Final}) = YetSyllable(nothing, nothing, f)

YetSyllable(i::YetJamo{Initial}, m::YetJamo{Medial}) = YetSyllable(i, m, nothing)
YetSyllable(i::YetJamo{Initial}, f::YetJamo{Final}) = YetSyllable(i, nothing, f)
YetSyllable(m::YetJamo{Medial}, f::YetJamo{Final}) = YetSyllable(nothing, m, f)

const 초성채움 = Char(0x115F)
const 중성채움 = Char(0x1160)

function to_string(ysyl::YetSyllable)::String
    초성, 중성, 종성 = ysyl.initial, ysyl.medial, ysyl.final
    if 초성 isa Nothing
        if 중성 isa Nothing
            if 종성 isa Nothing
                string()
            else
                string(초성채움, 중성채움, 종성)
            end
        else # 초성x, 중성o
            if 종성 isa Nothing
                string(초성채움, 중성)
            else
                string(초성채움, 중성, 종성)
            end
        end
    else # 초성o
        if 중성 isa Nothing
            if 종성 isa Nothing
                string(초성, 중성채움)
            else
                string(초성, 중성채움, 종성)
            end
        else # 초성o, 중성o
            if 종성 isa Nothing
                string(초성, 중성)
            else
                string(초성, 중성, 종성)
            end
        end
    end
end

const 초성ㄱ = 0x1100
const 중성ㅏ = 0x1161
const 종성ㄱ = 0x11a8

function is_yetjamo(c::Union{UInt16,UInt32})::Bool
    초성ㄱ <= c <= 0x115f || # ㄱ   ... HCF
    0xa960 <= c <= 0xa97c || # ㄷㅁ ... ㆆㆆ # Hangul Jamo Extended-A
    0x1160 <= c <= 0x11a7 || # HJF ... ㅗㅒ
    0xd7b0 <= c <= 0xd7c6 || # ㅗㅕ... .ㅔ  # Hangul Jamo Extended-B
    종성ㄱ <= c <= 0x11ff || # ㄱ   ... ㄴㄴ
    0xd7cb <= c <= 0xd7fb    # ㄴㄹ ... ㅍㅌ  # Hangul Jamo Extended-B
end

function YetJamo(c::UInt32)::YetJamo # YetJamoError
    if 초성ㄱ <= c <= 0x115f || # ㄱ   ... HCF
       0xa960 <= c <= 0xa97c    # ㄷㅁ ... ㆆㆆ # Hangul Jamo Extended-A
        YetJamo{Initial}(c)
    elseif 0x1160 <= c <= 0x11a7 || # HJF ... ㅗㅒ
           0xd7b0 <= c <= 0xd7c6    # ㅗㅕ... .ㅔ  # Hangul Jamo Extended-B
        YetJamo{Medial}(c)
    elseif 종성ㄱ <= c <= 0x11ff || # ㄱ   ... ㄴㄴ
           0xd7cb <= c <= 0xd7fb    # ㄴㄹ ... ㅍㅌ  # Hangul Jamo Extended-B
        YetJamo{Final}(c)
    else
        throw(YetJamoError(string(c)))
    end
end

function YetSyllable(syl::Syllable)::YetSyllable
    if syl.initial isa Nothing
        초성 = nothing
    else
        초성 = YetJamo{Initial}(초성ㄱ + indexof(초성표, syl.initial) - 1)
    end
    if syl.medial isa Nothing
        중성 = nothing
    else
        중성 = YetJamo{Medial}(중성ㅏ + indexof(중성표, syl.medial) - 1)
    end
    if syl.final isa Nothing
        종성 = nothing
    else
        종성 = YetJamo{Final}(종성ㄱ + indexof(종성표, syl.final) - 1 - 1)
    end
    YetSyllable(초성, 중성, 종성)
end

# 매핑 좀 해줘요
#
#   Hangul Jamo - Range: 1100–11FF
#   http://www.unicode.org/charts/PDF/U1100.pdf
#
#   Hangul Compatibility Jamo - Range: 3130–318F
#   http://www.unicode.org/charts/PDF/U3130.pdf

function throw_compose_error(ysyl::YetSyllable) # ComposeError
    throw(ComposeError(string("매핑 좀 해줘요 ", ysyl)))
end

function Syllable(ysyl::YetSyllable)::Syllable # ComposeError
    if ysyl.initial isa Nothing
        초성 = nothing
    else
        인덱스 = Int(ysyl.initial - 초성ㄱ) + 1
        if 1 <= 인덱스 <= 19 # length(초성표)
            초성 = Letter{Consonant}(초성표[인덱스])
        else
            throw_compose_error(ysyl) # 매핑 좀 해줘요
        end
    end
    if ysyl.medial isa Nothing
        중성 = nothing
    else
        인덱스 = Int(ysyl.medial - 중성ㅏ) + 1
        if 1 <= 인덱스 <= 21 # length(중성표)
            중성 = Letter{Vowel}(중성표[인덱스])
        else
            throw_compose_error(ysyl) # 매핑 좀 해줘요
        end
    end
    if ysyl.final isa Nothing
        종성 = nothing
    else
        인덱스 = Int(ysyl.final - 종성ㄱ + 1) + 1
        if 1 <= 인덱스 <= 28 # length(종성표)
            종성 = Letter{Consonant}(종성표[인덱스])
        else
            throw_compose_error(ysyl) # 매핑 좀 해줘요
        end
    end
    Syllable(초성, 중성, 종성)
end

# module Hangul


text:
using ProgressMeter

include("./interface.jl")

## fiber cut scenarios and their affected IP links
function GetAllScenarios(IPTopo, OpticalTopo, cutoff, is_ffc_scenario, ffc_fiber_cut_num)
    if is_ffc_scenario
        # simutaneous fiber cut
        scenarios = []
        probs = []
        # healthy state
        n_unifibers = length(OpticalTopo["bidirect_links"])
        s0 = ones(n_unifibers)
        push!(scenarios, s0)
        push!(probs, 1)

        # failure states
        for bits in collect(combinations(1:n_unifibers,ffc_fiber_cut_num))
            s = ones(n_unifibers)
            current_prob = 1
            for bit in bits
                s[bit] = 0
                current_prob = current_prob * OpticalTopo["fiber_probs"][bit]
            end
            push!(scenarios, s)
            push!(probs, current_prob)  # because of uni-directional fibers
        end
        probs[1] = 1 - sum(probs[2:end])
    else
        # probabilistic fiber cut with cutoff
        scenarios, probs = subScenarios(OpticalTopo, IPTopo, cutoff, first=true, last=false)
    end
    # Note that here the scenario is on directional graph
    # println(scenarios, probs)

    oscenarios = []
    oprobs = probs
    # translate optical scenario to bidirectional because optical fibers are bidirectional
    for s in 1:size(scenarios, 1)
        current_scenario = ones(length(OpticalTopo["links"]))
        for i in 1:length(scenarios[s])
            if scenarios[s][i] == 0
                xx = findfirst(x -> x == (OpticalTopo["bidirect_links"][i][1], OpticalTopo["bidirect_links"][i][2]), OpticalTopo["links"])
                yy = findfirst(x -> x == (OpticalTopo["bidirect_links"][i][2], OpticalTopo["bidirect_links"][i][1]), OpticalTopo["links"])
                current_scenario[xx] = 0
                current_scenario[yy] = 0
            end
        end
        push!(oscenarios, current_scenario)
    end
    # println(oscenarios, oprobs)

    OpticalScenarios = Dict()
    OpticalScenarios["code"] = oscenarios
    OpticalScenarios["prob"] = oprobs

    ## map optical failures to IP failures (IP links are bidirectional)
    iscenarios_all = []
    iprobs_all = []
    for s in 1:length(oscenarios)
        iscenarios = ones(length(IPTopo["links"]))
        iprobs = oprobs[s]
        for i in 1:length(oscenarios[s])
            if oscenarios[s][i] == 0
                # this fiber fails
                for k in 1:length(IPTopo["links"])
                    # if i in IPTopo["link_fiberroute"][k]
                    if in(i, IPTopo["link_fiberroute"][k])
                        iscenarios[k] = 0
                    end
                end
            end
        end
        push!(iscenarios_all, iscenarios)
        push!(iprobs_all, Float64(iprobs))
    end
    IPScenarios = Dict()
    IPScenarios["code"] = iscenarios_all
    IPScenarios["prob"] = iprobs_all
    return IPScenarios, OpticalScenarios
end

## generate subscenarios
function subScenarios(optical_topo, IPTopo, cutoff; first=true, last=true, progress=true)
    original = optical_topo["bidirect_fiber_probs"]
    p = ProgressMeter.ProgressUnknown("Computing scenarios cutoff=$(cutoff)...")
    if progress
        scenarios, probabilities = subScenariosRecursion(optical_topo, IPTopo, original, cutoff, progress=p)
    else
        scenarios, probabilities = subScenariosRecursion(optical_topo, IPTopo, original, cutoff)
    end
    if first == false
        scenarios = scenarios[2:end]
        probabilities = probabilities[2:end]
    end
    if last
        push!(scenarios, zeros(length(scenarios[1])))
        push!(probabilities, 1 - sum(probabilities))
    end
    if sum(probabilities) < 1
        probabilities = probabilities ./ sum(probabilities)
    end
    ProgressMeter.finish!(p)
    return scenarios, probabilities
end

## recursively generate failure scenarios with probability
function subScenariosRecursion(optical_topo, IPTopo, original, cutoff, remaining=[], offset=0, partial=[], scenarios=[], probabilities=[]; progress=nothing)
    if (size(partial,1) == 0)   #first
        push!(scenarios, ones(size(original,1)))
        push!(probabilities, prod(1 .- original))
        remaining = original
    else (size(partial,1) > 0)
        probs = 1 .- original
        bitmap = ones(size(original, 1))   #create bitmap
        for i in 1:length(partial)
            probs[partial[i]] = original[partial[i]]
            bitmap[partial[i]] = 0  # this is bidirectioanl fiber failure bitmap
        end
        ## translate bidirectional fiber failure to directional fiber failure for connectivity check
        directional_bitmap = ones(length(optical_topo["links"]))
        for i in 1:length(bitmap)
            if bitmap[i] == 0
                xx = findfirst(x -> x == (optical_topo["bidirect_links"][i][1], optical_topo["bidirect_links"][i][2]), optical_topo["links"])
                yy = findfirst(x -> x == (optical_topo["bidirect_links"][i][2], optical_topo["bidirect_links"][i][1]), optical_topo["links"])
                directional_bitmap[xx] = 0
                directional_bitmap[yy] = 0
            end
        end
        product = prod(probs)
        if progress !== nothing
            ProgressMeter.next!(progress, showvalues = [(:cutoff,cutoff), (:scenarios_added,length(probabilities)), (:last, product)])
        end
        if product >= cutoff
            if OpticalGraphConnectivity(optical_topo, directional_bitmap)  # remove fiber cut events that disconnect fiber topo
                ## map optical failures to IP failures
                iscenarios = ones(length(IPTopo["links"]))
                for i in 1:length(directional_bitmap)
                    if directional_bitmap[i] == 0
                        # this fiber fails
                        for k in 1:length(IPTopo["links"])
                            if in(i, IPTopo["link_fiberroute"][k])
                                iscenarios[k] = 0
                            end
                        end
                    end
                end

                ## check if this IP failure will disconnecct IP topology
                if IPGraphConnectivity(IPTopo,iscenarios)
                    push!(scenarios, bitmap)  # return bidirectional fiber failure
                    push!(probabilities, product)
                end
            end
        else
            return
        end
    end

    for i in 1:size(remaining,1)
        offset = size(original,1) - size(remaining,1)
        n = offset + i
        subScenariosRecursion(optical_topo, IPTopo, original, cutoff, remaining[i+1:end], offset, vcat(partial, [n]), scenarios, probabilities, progress=progress)
    end
    return scenarios, probabilities
end

## make sure fiber cut does not make the fiber topology disconnected (in real life no fiber cut should not separate a graph)
function OpticalGraphConnectivity(optical_topo, optical_failure_bitmap)
    num_nodes = length(optical_topo["nodes"])
    graph = LightGraphs.SimpleDiGraph(num_nodes)
    for f in 1:length(optical_topo["bidirect_links"])
        if optical_failure_bitmap[f] > 0
            LightGraphs.add_edge!(graph, optical_topo["bidirect_links"][f][1], optical_topo["bidirect_links"][f][2])
        end
    end
    if is_connected(graph)
        return true
    else
        return false
    end
end

## make sure fiber cut does not make the IP topology disconnected (in real life no fiber cut should not separate a graph)
function IPGraphConnectivity(IP_topo, IP_failure_bitmap)
    num_nodes = length(IP_topo["nodes"])
    graph = LightGraphs.SimpleDiGraph(num_nodes)
    for f in 1:length(IP_topo["links"])
        if IP_failure_bitmap[f] > 0
            LightGraphs.add_edge!(graph, IP_topo["links"][f][1], IP_topo["links"][f][2])
        end
    end
    if is_connected(graph)
        return true
    else
        return false
    end
end

max_stars_repo_name: hipersys-team/arrow

max_issues_repo_name: arrow-bow/arrow

max_forks_repo_name: arrow-bow/arrow


text:
using Test, FinancialMonteCarlo,Statistics;
@show "Black Scholes Model"
S0=100.0;
K=100.0;
r=[0.00,0.02];
T=1.0;
d=0.01;
D=90.0;

Nsim=100000;
Nstep=30;
sigma=0.2;
mc=MonteCarloConfiguration(Nsim,Nstep);
toll=0.8

rfCurve=ZeroRate(r,T);

FwdData=Forward(T)
EUData=EuropeanOption(T,K)
EUBin=BinaryEuropeanOption(T,K)
AMData=AmericanOption(T,K)
AmBin=BinaryEuropeanOption(T,K)
BarrierData=BarrierOptionDownOut(T,K,D)
AsianFloatingStrikeData=AsianFloatingStrikeOption(T)
AsianFixedStrikeData=AsianFixedStrikeOption(T,K)
Model=BlackScholesProcess(sigma,Underlying(S0,d));

display(Model)

@show FwdPrice=pricer(Model,rfCurve,mc,FwdData);
@show EuPrice=pricer(Model,rfCurve,mc,EUData);
@show EuBinPrice=pricer(Model,rfCurve,mc,EUBin);
@show AmPrice=pricer(Model,rfCurve,mc,AMData);
@show AmBinPrice=pricer(Model,rfCurve,mc,AmBin);
@show BarrierPrice=pricer(Model,rfCurve,mc,BarrierData);
@show AsianPrice1=pricer(Model,rfCurve,mc,AsianFloatingStrikeData);
@show AsianPrice2=pricer(Model,rfCurve,mc,AsianFixedStrikeData);


text:
using PDBTools
using ComplexMixtures

# Load a PDB file of the system
system = readPDB("./simulation/equilibrated.pdb")

# Trajectory file
# The trajectory file is available at:
# https://drive.google.com/file/d/1ug43ncCLsBATaJrT9zlbaqK6AORVvhhx/view?usp=sharing
trajectory_file = "../trajectories/traj_Polyacry.dcd"

# Select the atoms corresponding DMF molecules
dmf = select(system,"resname DMF")

# Select the atoms corresponding to the Poly-acrylamide
acr = select(system,"resname FACR or resname ACR or resname LACR")

# Set the solute and the solvent selections for ComplexMixtures
solute = Selection(acr,nmols=1)
solvent = Selection(dmf,natomspermol=12)

# Set the trajectory structure
traj = Trajectory(trajectory_file,solute,solvent)

# Use a large dbulk distance for better KB convergence
opt = ComplexMixtures.Options(dbulk=20.)

# Compute the mddf and associated properties
mddf_dmf_acr = mddf(traj,opt)

# Save results to file for later use
save(mddf_dmf_acr,"./results/mddf_dmf_acr.json")


text:
export Layout, LForwardReference, LFixedVector, LVarVector
export is_template_fixed, is_template_variable, simple_size, total_size, create_bytes

"""
Layout

All structs used to describe the memory layout (of a C-data structure) need to be
subtypes of this.
Some controlling objects used in such templates to describe vectors and pointers
have also this type.
A `Layout` structure and a memory pointer are needed to construct an `CAccessor` object.
"""
abstract type Layout end

abstract type LVector{T} <: Layout end

Base.eltype(::Type{<:LVector{T}}) where T = T

# Layout Elements
"""
    LFixedVector{T,N}

Denote a fixed size vector with element type `T` and size `N`.
"""
struct LFixedVector{T,N} <: LVector{T}
    p::NTuple{N,T}
end
Base.length(::Type{LFixedVector{T,N}}, ::Any) where {T,N} = N
Base.eltype(::Type{LFixedVector{T,N}}) where {T,N} = T

"""
    LVarVector{T,F}

Denote a variable length vector with element type `T` in a template.
`F` is a function, which calculates the length of the vector, given the
accessor object containing the vector.

Example:
    struct A <: Layout
        len::Int
        vec::NVarVector{Float64, (x) -> x.len}
    end

"""
struct LVarVector{T,F}  <: LVector{T}
    p::NTuple{0,T}
end
Base.length(::Type{LVarVector{T,F}}, x) where {T,F} = F(x)
Base.eltype(::Type{LVarVector{T,F}}) where {T,F} = T

struct LForwardReference{M,L} <: Layout
    p::Ptr{Nothing}
end
Base.eltype(::Type{LForwardReference{M,L}}) where {M,L} = M.name.module.eval(L)

const TEMPLATE_FIXED = true
const TEMPLATE_VAR = false
"""
    is_template_variable(type)

Has the layout described by `type` a variable size
(for example variable sized vector in last field of a struct)?
"""
is_template_variable(T::Type, deep::Bool=false) = !is_template_fixed(T, deep)

"""
    is_template_fixed(type)

Has the layout described by `type` a fixed size.
"""
is_template_fixed(T::Type, deep::Bool=false) = is_template_fixed(T, deep, Dict())
function is_template_fixed(::Type{T}, deep::Bool, dup) where T
    isprimitivetype(T) || throw(ArgumentError("$T is not a supported layout type"))
    TEMPLATE_FIXED
end
function is_template_fixed(::Type{S}, deep::Bool, dup) where {T,S<:Ptr{T}}
    T <: Ptr && throw(ArgumentError("$S is not a supported layout type"))
    get!(dup, S) do
        dup[S] = TEMPLATE_FIXED
        d = is_template_fixed(T, deep, dup)
        deep ? d : TEMPLATE_FIXED
    end
end
function is_template_fixed(::Type{S}, deep::Bool, dup) where {S<:LForwardReference}
    is_template_fixed(Ptr{eltype(S)}, deep, dup)
end
function is_template_fixed(::Type{S}, deep::Bool, dup) where {T,N,S<:LFixedVector{T,N}}
    get!(dup, S) do
        dup[S] = TEMPLATE_FIXED
        k = is_template_fixed(T, deep, dup)
        if N > 1 && k == TEMPLATE_VAR
            throw(ArgumentError("$S with variable length elements"))
        end
        N == 0 ? TEMPLATE_FIXED : k
    end
end
function is_template_fixed(::Type{S}, deep::Bool, dup) where {T,S<:LVarVector{T}}
    get!(dup, S) do
        dup[S] = TEMPLATE_VAR
        is_template_fixed(T, deep, dup)
        TEMPLATE_VAR
    end
end
function is_template_fixed(::Type{T}, deep::Bool, dup) where {T<:Layout}
    get!(dup, T) do
        k = dup[T] = TEMPLATE_FIXED
        if !isbitstype(T)
            text = isconcretetype(T) ? "bits" : "concrete"
            throw(ArgumentError("$T is not a $text type struct"))
        end
        fields = fieldnames(T)
        n = length(fields)
        for i = 1:n
            f = fields[i]
            F = fieldtype(T, f)
            k = is_template_fixed(F, deep, dup)
            if i < n && k == TEMPLATE_VAR
                throw(ArgumentError("$F has variable length in '$T.$f' - not last field"))
            end
        end
        k
    end
end

"""
    create_bytes(T)

Return a `Vector{UInt8}` capable of keeping all data of an object of type `T`
including structure subfields. Ptr-fields are populated with correct pointers into
the same vector recursively.
"""
function create_bytes(::Type{T}, veclens=()) where T
    n = total_size(T, veclens) * 2
    buf = fill(UInt8(0x0), n)
    m = create_bytes!(buf, T, 0, veclens)
    resize!(buf, m)
    buf
end

function create_bytes!(bytes::Vector{UInt8}, ::Type{T}, offset, veclens) where T <: Layout
    len = simple_size(T, veclens)
    off = offset
    ptr = pointer_from_vector(bytes) + off
    noff = len
    j = 0
    for i = 1:fieldcount(T)
        f = fieldoffset(T, i)
        p = ptr + f
        off = offset + f
        F = fieldtype(T, i)
        vl = is_template_variable(F, true) ? veclens[j += 1] : ()

        if F <: Union{Ptr,LForwardReference}
            noff = align(noff)
            unsafe_store!(Ptr{Ptr{Nothing}}(p), ptr + noff)
            len = create_bytes!(bytes, eltype(F), noff + offset, vl)
            noff += len
        else
            create_bytes!(bytes, F, off, vl)
        end
    end
    noff
end

function create_bytes!(bytes::Vector{UInt8}, ::Type{T}, offset, veclens) where {N,F,T<:LFixedVector{F,N}}
    len = simple_size(T, veclens)
    off = offset
    ptr = pointer_from_vector(bytes) + off
    noff = len
    j = 0
    for i = 1:N
        vl = is_template_variable(F, true) ? veclens[j += 1] : ()
        f = simple_size(F, vl) * (i - 1)

        if F <: Union{Ptr,LForwardReference}
            noff = align(noff)
            unsafe_store!(ptr + f, ptr + noff)
            len = create_bytes!(bytes, eltype(F), noff + off, vl)
            noff += len
        else
            create_bytes!(bytes, F, offset + f, vl)
        end
    end
    noff
end
function create_bytes!(bytes::Vector{UInt8}, ::Type{T}, offset, veclens) where {F,T<:LVarVector{F}}
    len = simple_size(T, veclens)
    off = offset
    ptr = pointer_from_vector(bytes) + off
    noff = len
    vlen(i) = i <= length(veclens) ? veclens[i] : ()
    j = 0
    N = vlen(j += 1)
    for i = 1:N
        vl = is_template_variable(F, true) ? vlen(j += 1) : ()
        f = simple_size(F, vl) * (i - 1)

        if F <: Union{Ptr,LForwardReference}
            noff = align(noff)
            unsafe_store!(ptr + f, ptr + noff)
            len = create_bytes!(bytes, eltype(F), noff + off, vl)
            noff += len
        else
            create_bytes!(bytes, F, off + f, vl)
        end
    end
    noff
end

function create_bytes!(::Vector{UInt8}, ::Type{T}, offset, veclens) where T
    simple_size(T, ())
end

function align(p::Integer, s::Integer=sizeof(Ptr)) # s must be 2^n
    t = s - 1
    (p + t )  & ~t
end

simple_size(T::Type, veclens) = blength(T, veclens, Val(false))
total_size(T::Type, veclens) = blength(T, veclens, Val(true))

function blength(::Type{T}, veclens, v::Val{P}) where {P,F,N,T<:LFixedVector{F,N}}
    s = sizeof(T)
    j = 0
    for _ = 1:N
        j, s = blength_helper(F, veclens, j, s, v, T)
    end
    if j < length(veclens)
        throw(ArgumentError("too many variable length specifiers for $T only $j are needed"))
    end
    s
end

function blength(::Type{T}, veclens, v::Val{P}) where {P,S,T<:LVarVector{S}}
    isempty(veclens) && return 0
    n = first(veclens)
    n == 0 && return 0
    vl(i) = i < length(veclens) ? veclens[i+1] : ()
    sum(blength(S, vl(i), v) for i = 1:n)
end

blength(::Type{Ptr{T}}, veclens, v::Val{P}) where {P,T} = P ? blength(T, veclens, v) : 0

function blength(::Type{T}, veclens, v::Val{P}) where {P,T<:Layout}
    s = sizeof(T)
    j = 0
    for i = 1:fieldcount(T)
        F = fieldtype(T, i)
        j, s = blength_helper(F, veclens, j, s, v, T)
    end
    if j < length(veclens)
        throw(ArgumentError("too many variable length specifiers for $T- only $j are needed"))
    end
    s
end

blength(::Type{T}, veclens, ::Val) where T = sizeof(T)

function blength_helper(::Type{F}, veclens, j, s, v::Val{P}, T) where {P,F}
    if is_template_variable(F, true)
        j += 1
        if j > length(veclens)
            throw(ArgumentError("not enough variable length specifiers for $T"))
        end
        vl = veclens[j]
    else
        vl = ()
    end
    al = Base.datatype_alignment(F)
    s = align(s, al)
    s += F <: LVarVector ? blength(F, vl, v) :
        P && F <: Union{Ptr,LForwardReference} ? blength(eltype(F), vl, v) : 0

    j, s
end


text:
using Test

import TreeViews: hastreeview, numberofnodes, treelabel, nodelabel, treenode

# test for sane defaults
struct TVT_default
    a
    b
    c
    d
end

teststruct = TVT_default(1, "asd", UInt8(2), 22.22)

hastreeview(x::TVT_default) = true

@test hastreeview(teststruct) == true
@test numberofnodes(teststruct) == 4
@test sprint(io -> treelabel(io, teststruct)) == "TVT_default"
@test sprint(io -> nodelabel(io, teststruct, 3)) == "c"
@test treenode(teststruct, 3) == teststruct.c

# customization
struct TVT_customized
    a
    b
    c
    d
end

teststruct = TVT_customized(1, "asd", UInt8(2), 22.22)

hastreeview(::TVT_customized) = true

numberofnodes(x::TVT_customized) = 2
treelabel(io::IO, x::TVT_customized, ::MIME"text/plain") = print(io, "customized")
treelabel(io::IO, x::TVT_customized, ::MIME"text/html") = print(io, "HTML")

function nodelabel(io::IO, x::TVT_customized, i::Integer, ::MIME"text/plain")
    i <= 2 || throw(BoundsError(x, i))
    print(io, "customized$i")
end
function nodelabel(io::IO, x::TVT_customized, i::Integer, ::MIME"text/html")
    i <= 2 || throw(BoundsError(x, i))
    print(io, "HTML$i")
end

function treenode(x::TVT_customized, i::Integer)
    i == 1 && return x.a
    i == 2 && return TVT_default(x.b, x.c, x.c, x.d)
end

@test numberofnodes(teststruct) == 2
@test sprint(io -> treelabel(io, teststruct)) == "customized"
@test sprint(io -> treelabel(io, teststruct, "text/html")) == "HTML"
@test sprint(io -> treelabel(io, teststruct, MIME"text/html"())) == "HTML"
@test_throws BoundsError sprint(io -> nodelabel(io, teststruct, 3))
@test sprint(io -> nodelabel(io, teststruct, 2)) == "customized2"
@test treenode(teststruct, 1) == teststruct.a
@test treenode(teststruct, 2) == TVT_default(teststruct.b, teststruct.c, teststruct.c, teststruct.d)

@test_deprecated(@test sprint(io -> treelabel(io, teststruct, 1)) == "customized1")
@test_deprecated(@test sprint(io -> treelabel(io, teststruct, 1, "text/html")) == "HTML1")


text:
using SpikingNetworks
@unit
parseParameters"""
Cm = 1*uF
Iapp = 2*uA
gL = 0.1*msiemens
EL = -65*mV
ENa = 55*mV
EK = -90*mV
gNa = 35*msiemens
gK = 9*msiemens"""

parseEquations"""
dv/dt=(-gNa*m(v)^3*h*(v-ENa)-gK*n^4*(v-EK)-gL*(v-EL)+Iapp)/Cm
dh/dt=5*(α_h(v)*(1-h)-β_h(v)*h)
dn/dt=5*(α_n(v)*(1-n)-β_n(v)*n)
m(v)=α_m(v)/(α_m(v)+β_m(v))
α_m(v)=-0.1*(v+35)/(exp(-0.1*(v+35))-1)
β_m(v)=4*exp(-(v+60)/18)
α_h(v)=0.07*exp(-(v+58)/20)
β_h(v)=1./(exp(-0.1*(v+28))+1)
α_n(v)=-0.01*(v+34)/(exp(-0.1*(v+34))-1)
β_n(v)=0.125*exp(-(v+44)/80)"""

net=init"""
v=-70*ones(N)
h=1*ones(N)
n=0*ones(N)"""

stateMonitor"[1]"

net.dt=0.1
run!(net,100/net.dt)
using MatlabPlot
mplot(net.vRecord)

# workspace()


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

mutable struct _CacheModel
    cache::Vector{UInt8}
    variable_type::Vector{_VariableType}
    variable_primal::Vector{Float64}
    variable_lower::Vector{Float64}
    variable_upper::Vector{Float64}
    constraints::Vector{Expr}
    constraint_lower::Vector{Float64}
    constraint_upper::Vector{Float64}
    objective::Expr
    sense::MOI.OptimizationSense
    function _CacheModel()
        return new(
            zeros(UInt8, 64),
            _VariableType[],
            Float64[],
            Float64[],
            Float64[],
            Expr[],
            Float64[],
            Float64[],
            :(),
            MOI.FEASIBILITY_SENSE,
        )
    end
end

function Base.read!(io::IO, model::Model)
    cache = _CacheModel()
    _parse_header(io, cache)
    while !eof(io)
        _parse_section(io, cache)
    end
    model.model = _to_model(cache)
    return
end

function _resize_variables(model::_CacheModel, n::Int)
    resize!(model.variable_type, n)
    fill!(model.variable_type, _CONTINUOUS)
    resize!(model.variable_lower, n)
    fill!(model.variable_lower, -Inf)
    resize!(model.variable_upper, n)
    fill!(model.variable_upper, Inf)
    resize!(model.variable_primal, n)
    fill!(model.variable_primal, 0.0)
    return
end

function _resize_constraints(model::_CacheModel, n::Int)
    resize!(model.constraint_lower, n)
    fill!(model.constraint_lower, -Inf)
    resize!(model.constraint_upper, n)
    fill!(model.constraint_upper, Inf)
    resize!(model.constraints, n)
    for i in 1:n
        model.constraints[i] = :()
    end
    return
end

_is_valid_number(::Type{Int}, x::UInt8) = UInt8('0') <= x <= UInt8('9')

function _is_valid_number(::Type{Float64}, x::UInt8)
    if _is_valid_number(Int, x)
        return true
    elseif x == UInt8('+') || x == UInt8('-')
        return true
    elseif x == UInt8('.') || x == UInt8('e') || x == UInt8('E')
        return true
    end
    return false
end

function _next_token(::Type{T}, io::IO, cache::Vector{UInt8}) where {T}
    # Skip all spaces
    byte = UInt8(' ')
    while byte == UInt8(' ')
        byte = read(io, UInt8)
    end
    @assert _is_valid_number(T, byte)
    cache[1] = byte
    i = 1
    while _is_valid_number(T, peek(io, UInt8))
        i += 1
        cache[i] = read(io, UInt8)
    end
    return i
end

function _next(::Type{Float64}, io::IO, model::_CacheModel)
    nnz = _next_token(Float64, io, model.cache)
    @assert nnz > 0
    return parse(Float64, String(model.cache[1:nnz]))
end

function _next(::Type{Int}, io::IO, model::_CacheModel)
    nnz = _next_token(Int, io, model.cache)
    @assert nnz > 0
    y = 0
    mult = 1
    for i in nnz:-1:1
        y += mult * (model.cache[i] - UInt8('0'))
        mult *= 10
    end
    return y
end

"""
    _read_til_newline(io::IO)

This function reads until it finds a new line character. This is useful for
skipping comments.
"""
function _read_til_newline(io::IO)
    while read(io, UInt8) != UInt8('\n')
    end
    return
end

_force_expr(expr::Expr) = expr
_force_expr(expr) = Expr(:call, :+, expr)

function _parse_expr(io::IO, model::_CacheModel)
    char = Char(read(io, UInt8))
    if char == 'o'
        opcode = _next(Int, io, model)
        _read_til_newline(io)
        arity, op_func = _AMPL_TO_JULIA[opcode]
        op_sym = Symbol(op_func)
        if arity == -1
            arity = _next(Int, io, model)
            _read_til_newline(io)
            if op_sym == :sum
                op_sym = :+
            elseif op_sym == :minimum
                op_sym = :min
            elseif op_sym == :maximum
                op_sym = :max
            end
        end
        parent = Expr(:call, op_sym)
        for _ in 1:arity
            child = _parse_expr(io, model)
            push!(parent.args, child)
        end
        return parent
    elseif char == 'v'
        index = _next(Int, io, model)
        _read_til_newline(io)
        return MOI.VariableIndex(index + 1)
    else
        @assert char == 'n'
        ret = _next(Float64, io, model)
        _read_til_newline(io)
        return ret
    end
end

function _to_model(data::_CacheModel)
    model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
    x = MOI.add_variables(model, length(data.variable_primal))
    for (xi, lb, ub) in zip(x, data.variable_lower, data.variable_upper)
        if lb > -Inf
            MOI.add_constraint(model, xi, MOI.GreaterThan(lb))
        end
        if ub < Inf
            MOI.add_constraint(model, xi, MOI.LessThan(ub))
        end
    end
    for (xi, type) in zip(x, data.variable_type)
        if type == _INTEGER
            MOI.add_constraint(model, xi, MOI.Integer())
        elseif type == _BINARY
            MOI.add_constraint(model, xi, MOI.ZeroOne())
        end
    end
    MOI.set.(model, MOI.VariablePrimalStart(), x, data.variable_primal)
    MOI.set(model, MOI.ObjectiveSense(), data.sense)
    nlp = MOI.Nonlinear.Model()
    MOI.Nonlinear.set_objective(nlp, data.objective)
    for (i, expr) in enumerate(data.constraints)
        lb, ub = data.constraint_lower[i], data.constraint_upper[i]
        if lb == ub
            MOI.Nonlinear.add_constraint(nlp, expr, MOI.EqualTo(lb))
        elseif -Inf < lb < ub < Inf
            MOI.Nonlinear.add_constraint(nlp, expr, MOI.Interval(lb, ub))
        elseif -Inf == lb && ub < Inf
            MOI.Nonlinear.add_constraint(nlp, expr, MOI.LessThan(ub))
        else
            @assert -Inf < lb && ub == Inf
            MOI.Nonlinear.add_constraint(nlp, expr, MOI.GreaterThan(lb))
        end
    end
    evaluator =
        MOI.Nonlinear.Evaluator(nlp, MOI.Nonlinear.SparseReverseMode(), x)
    block = MOI.NLPBlockData(evaluator)
    MOI.set(model, MOI.NLPBlock(), block)
    return model
end

function _parse_header(io::IO, model::_CacheModel)
    # Line 1
    # We don't support the binary format.
    byte = read(io, UInt8)
    if byte != UInt8('g')
        error("Unable to parse NL file : unsupported mode $(Char(byte))")
    end
    # L1 has some magic bytes for AMPL internals (to quote David, "The numbers
    # on the first line matter to AMPL; for other uses, it is best simply to
    # supply the ones shown above.")
    _read_til_newline(io)
    # Line 2
    # The number of variables
    n_var = _next(Int, io, model)
    _resize_variables(model, n_var)
    # The number of constraints
    _resize_constraints(model, _next(Int, io, model))
    # The number of linear objectives
    @assert 0 <= _next(Int, io, model) <= 1
    # The number of range constraints
    @assert _next(Int, io, model) >= 0
    # The number of equations
    @assert _next(Int, io, model) >= 0
    # The number of logical constraints. This one is optional, so just read til
    # the end of the line.
    # @assert _next(Int, io, model) == 0
    _read_til_newline(io)
    # Line 3
    # The number of nonlinear constraints
    @assert _next(Int, io, model) >= 0
    # The number of nonlinear objectives
    @assert 0 <= _next(Int, io, model) <= 1
    _read_til_newline(io)
    # Line 4
    # The number of nonlinear network constraints
    @assert _next(Int, io, model) == 0
    # The number of linear network constraints
    @assert _next(Int, io, model) == 0
    _read_til_newline(io)
    # Line 5
    # The number of nonlienar variables in constraints
    nlvc = _next(Int, io, model)
    # The number of nonlienar variables in objectives
    nlvo = _next(Int, io, model)
    # The number of nonlienar variables in constraints and objectives (both)
    nlvb = _next(Int, io, model)
    _read_til_newline(io)
    # Line 6
    # The number of linear network variables
    @assert _next(Int, io, model) == 0
    # The number of user-defined functions
    @assert _next(Int, io, model) == 0
    # The number of "arith"
    # TODO(odow): I don't know what this is.
    @assert _next(Int, io, model) == 0
    # The "flags" entry. This is mainly used for specifying that we want duals.
    # Ignore when reading.
    _next(Int, io, model)
    _read_til_newline(io)
    # Line 7
    # Number of binary variables
    nbv = _next(Int, io, model)
    # Number of integer variables
    niv = _next(Int, io, model)
    # Number of integer variables in nonlinear constraints AND objectives
    nlvbi = _next(Int, io, model)
    # Number of integer variables in nonlinear constraints
    nlvci = _next(Int, io, model)
    # Number of integer variables in nonlinear objectives
    nlvoi = _next(Int, io, model)
    _read_til_newline(io)
    # Line 8
    # Read the number of nonzeros in Jacobian and gradient, but don't do
    # anything with that information.
    @assert _next(Int, io, model) >= 0
    @assert _next(Int, io, model) >= 0
    _read_til_newline(io)
    # Line 9
    # We don't support reading variable and constraint names
    @assert _next(Int, io, model) == 0
    @assert _next(Int, io, model) == 0
    _read_til_newline(io)
    # Line 10
    # We don't support reading common subexpressions
    for _ in 1:5
        if _next(Int, io, model) > 0
            error("Unable to parse NL file : we don't support common exprs")
        end
    end
    _read_til_newline(io)
    # ==========================================================================
    # Deal with the integrality of variables. This is quite complicated, so go
    # read the README in this folder.
    offsets = [
        nlvb - nlvbi,
        nlvbi,
        nlvc - (nlvb + nlvci),
        nlvci,
        max(0, nlvo - nlvc - nlvoi),
        nlvoi,
        n_var - (max(nlvc, nlvo) + nbv + niv),
        nbv,
        niv,
    ]
    types = [
        _CONTINUOUS,
        _INTEGER,
        _CONTINUOUS,
        _INTEGER,
        _CONTINUOUS,
        _INTEGER,
        _CONTINUOUS,
        _BINARY,
        _INTEGER,
    ]
    offset = 0
    for i in 1:9
        for _ in 1:offsets[i]
            offset += 1
            model.variable_type[offset] = types[i]
        end
    end
    return
end

function _parse_section(io::IO, model::_CacheModel)
    char = Char(read(io, UInt8))
    _parse_section(io, Val(char), model)
    return
end

function _parse_section(::IO, ::Val{T}, ::_CacheModel) where {T}
    return error("Unable to parse NL file: unhandled header $T")
end

function _parse_section(io::IO, ::Val{'C'}, model::_CacheModel)
    index = _next(Int, io, model) + 1
    _read_til_newline(io)
    expr = _force_expr(_parse_expr(io, model))
    current = model.constraints[index]
    if current == :()
        model.constraints[index] = expr
    else
        model.constraints[index] = Expr(:call, :+, current, expr)
    end
    return
end

function _parse_section(io::IO, ::Val{'O'}, model::_CacheModel)
    @assert _next(Int, io, model) == 0
    sense = _next(Int, io, model)
    if sense == 1
        model.sense = MOI.MAX_SENSE
    else
        @assert sense == 0
        model.sense = MOI.MIN_SENSE
    end
    _read_til_newline(io)
    expr = _force_expr(_parse_expr(io, model))
    if model.objective == :()
        model.objective = expr
    else
        model.objective = Expr(:call, :+, model.objective, expr)
    end
    return
end

function _parse_section(io::IO, ::Val{'x'}, model::_CacheModel)
    index = _next(Int, io, model)
    _read_til_newline(io)
    for _ in 1:index
        xi = _next(Int, io, model) + 1
        v = _next(Float64, io, model)
        model.variable_primal[xi] = v
        _read_til_newline(io)
    end
    return
end

# TODO(odow): we don't read in dual starts.
function _parse_section(io::IO, ::Val{'d'}, model::_CacheModel)
    index = _next(Int, io, model)
    _read_til_newline(io)
    for _ in 1:index
        _read_til_newline(io)
    end
    return
end

function _parse_section(io::IO, ::Val{'r'}, model::_CacheModel)
    _read_til_newline(io)
    for i in 1:length(model.constraint_lower)
        type = _next(Int, io, model)
        if type == 0
            model.constraint_lower[i] = _next(Float64, io, model)
            model.constraint_upper[i] = _next(Float64, io, model)
        elseif type == 1
            model.constraint_upper[i] = _next(Float64, io, model)
        elseif type == 2
            model.constraint_lower[i] = _next(Float64, io, model)
        elseif type == 3
            # Free constraint
        else
            @assert type == 4
            value = _next(Float64, io, model)
            model.constraint_lower[i] = value
            model.constraint_upper[i] = value
        end
        _read_til_newline(io)
    end
    return
end

function _parse_section(io::IO, ::Val{'b'}, model::_CacheModel)
    _read_til_newline(io)
    for i in 1:length(model.variable_lower)
        type = _next(Int, io, model)
        if type == 0
            model.variable_lower[i] = _next(Float64, io, model)
            model.variable_upper[i] = _next(Float64, io, model)
        elseif type == 1
            model.variable_upper[i] = _next(Float64, io, model)
        elseif type == 2
            model.variable_lower[i] = _next(Float64, io, model)
        elseif type == 3
            # Free variable
        else
            @assert type == 4
            value = _next(Float64, io, model)
            model.variable_lower[i] = value
            model.variable_upper[i] = value
        end
        _read_til_newline(io)
    end
    return
end

# We ignore jacobian counts for now
function _parse_section(io::IO, ::Val{'k'}, model::_CacheModel)
    _read_til_newline(io)
    for _ in 2:length(model.variable_lower)
        _read_til_newline(io)
    end
    return
end

function _parse_section(io::IO, ::Val{'J'}, model::_CacheModel)
    i = _next(Int, io, model) + 1
    nnz = _next(Int, io, model)
    _read_til_newline(io)
    expr = Expr(:call, :+)
    for _ in 1:nnz
        x = _next(Int, io, model)
        c = _next(Float64, io, model)
        if !iszero(c)
            push!(expr.args, Expr(:call, :*, c, MOI.VariableIndex(x + 1)))
        end
        _read_til_newline(io)
    end
    if length(expr.args) == 1
        # Linear part is just zeros
    elseif model.constraints[i] == :()
        model.constraints[i] = expr
    else
        model.constraints[i] = Expr(:call, :+, expr, model.constraints[i])
    end
    return
end

function _parse_section(io::IO, ::Val{'G'}, model::_CacheModel)
    i = _next(Int, io, model) + 1
    nnz = _next(Int, io, model)
    _read_til_newline(io)
    expr = Expr(:call, :+)
    for _ in 1:nnz
        x = _next(Int, io, model)
        c = _next(Float64, io, model)
        if !iszero(c)
            push!(expr.args, Expr(:call, :*, c, MOI.VariableIndex(x + 1)))
        end
        _read_til_newline(io)
    end
    if length(expr.args) == 1
        # Linear part is just zeros
    elseif model.objective == :()
        model.objective = expr
    else
        model.objective = Expr(:call, :+, expr, model.objective)
    end
    return
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
using PyTensorFlow
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

if Sys.islinux()
py"""
import tensorflow as tf
libSatOp = tf.load_op_library('../Ops/Saturation/build/libSatOp.so')
@tf.custom_gradient
def sat_op(s0,pt,permi,poro,qw,qo,muw,muo,sref,dt,h):
    sat = libSatOp.sat_op(s0,pt,permi,poro,qw,qo,muw,muo,sref,dt,h)
    def grad(dy):
        return libSatOp.sat_op_grad(dy, sat, s0,pt,permi,poro,qw,qo,muw,muo,sref,dt,h)
    return sat, grad
"""
elseif Sys.isapple()
py"""
import tensorflow as tf
libSatOp = tf.load_op_library('../Ops/Saturation/build/libSatOp.dylib')
@tf.custom_gradient
def sat_op(s0,pt,permi,poro,qw,qo,muw,muo,sref,dt,h):
    sat = libSatOp.sat_op(s0,pt,permi,poro,qw,qo,muw,muo,sref,dt,h)
    def grad(dy):
        return libSatOp.sat_op_grad(dy, sat, s0,pt,permi,poro,qw,qo,muw,muo,sref,dt,h)
    return sat, grad
"""
elseif Sys.iswindows()
py"""
import tensorflow as tf
libSatOp = tf.load_op_library('../Ops/Saturation/build/libSatOp.dll')
@tf.custom_gradient
def sat_op(s0,pt,permi,poro,qw,qo,muw,muo,sref,dt,h):
    sat = libSatOp.sat_op(s0,pt,permi,poro,qw,qo,muw,muo,sref,dt,h)
    def grad(dy):
        return libSatOp.sat_op_grad(dy, sat, s0,pt,permi,poro,qw,qo,muw,muo,sref,dt,h)
    return sat, grad
"""
end

sat_op = py"sat_op"


text:
module FLPprocess

using DataFrames
using Dates
using OrderedCollections
using ShiftedArrays
using CSV

#=
Remember to:
Pokes doesn't process anymore the file name to extract name etc, move this to process session
=#

include("utilities.jl")
include("create_DataIndex.jl")
include("process_pokes.jl")
include("process_bouts.jl")
include("process_streaks.jl")
include("process_session.jl")
include("process_dataset.jl")

export get_DataIndex, process_pokes, process_bouts, process_streaks, process_session, process_dataset

end # module


text:
# An immutable implementation of AbstractMOP for internal use.

# TODO make num_vars a type parameter?
# TODO make vectors (`full_lb` etc.) statically sized?
@with_kw struct StaticMOP{
	F <: AbstractFloat,
	OT <: Tuple{Vararg{<:AbstractObjective}},
	} <: AbstractMOP{false}

	num_vars :: Int
	num_objectives :: Int # no of scalar valued outputs

	tuple_of_objectives :: OT
	objf_output_mapping :: Base.ImmutableDict{UInt64,Vector{Int}}

	full_lb :: Vector{F}
	full_ub :: Vector{F}

	full_width :: Vector{F} = full_ub - full_lb

	full_lb_int :: Vector{F} = [isinf(l) ? l : 0.0 for l ∈ full_lb ]
	full_ub_int :: Vector{F} = [ isinf(u) ? u : 1.0 for u ∈ full_ub ]

	internal_indices :: Vector{Int} = vcat( (collect(objf_output_mapping[hash(objf)]) for objf in tuple_of_objectives )... )
	reverse_indices :: Vector{Int} = sortperm(internal_indices)
end

function StaticMOP( mop :: AbstractMOP, out_el_type :: Type = Nothing )
	num_out = num_objectives(mop)

	@assert num_out > 0 "There must at least be 1 objective."

	mop_objectives = list_of_objectives(mop)
	if out_el_type <: AbstractFloat
		tuple_of_objectives = Tuple( [ OutTypeWrapper(objf, Vector{out_el_type}) for objf in mop_objectives ] )
	else
		tuple_of_objectives = Tuple( mop_objectives )
	end

	objf_output_mapping = Base.ImmutableDict(
		( hash( tuple_of_objectives[i]) => output_indices( mop_objectives[i], mop) for i = eachindex(mop_objectives))...
	)
	return StaticMOP(;
		num_vars = num_vars(mop),
		num_objectives = num_out,
		tuple_of_objectives,
		objf_output_mapping,
		full_lb = full_lower_bounds(mop),
		full_ub = full_upper_bounds(mop)
	)
end

full_lower_bounds( mop :: StaticMOP ) = mop.full_lb
full_upper_bounds( mop :: StaticMOP ) = mop.full_ub

list_of_objectives( mop :: StaticMOP ) = mop.tuple_of_objectives

num_vars( mop :: StaticMOP ) = mop.num_vars

num_objectives( mop :: StaticMOP ) = mop.num_objectives

output_indices( mop :: StaticMOP ) = mop.internal_indices
output_indices( objf :: AbstractObjective, mop :: StaticMOP ) = mop.objf_output_mapping[ hash(objf) ]

full_lower_bounds_internal( mop :: StaticMOP ) = mop.full_lb_int
full_upper_bounds_internal( mop :: StaticMOP ) = mop.full_ub_int


text:
### w-theta form of the non-convex AC equations

export
    ACTPowerModel, StandardACTForm

""
abstract type AbstractACTForm <: AbstractPowerFormulation end

""
abstract type StandardACTForm <: AbstractACTForm end

""
const ACTPowerModel = GenericPowerModel{StandardACTForm}

"default AC constructor"
ACTPowerModel(data::Dict{String,Any}; kwargs...) =
    GenericPowerModel(data, StandardACTForm; kwargs...)

"`t[ref_bus] == 0`"
function constraint_theta_ref(pm::GenericPowerModel{T}, n::Int, i::Int) where T <: AbstractACTForm
    @constraint(pm.model, pm.var[:nw][n][:va][i] == 0)
end

""
function variable_voltage(pm::GenericPowerModel{T}, n::Int=pm.cnw; kwargs...) where T <: AbstractACTForm
    variable_voltage_angle(pm, n; kwargs...)
    variable_voltage_magnitude_sqr(pm, n; kwargs...)
    variable_voltage_product(pm, n; kwargs...)
end

function constraint_voltage(pm::GenericPowerModel{T}, n::Int) where T <: StandardACTForm
    t = pm.var[:nw][n][:va]
    w = pm.var[:nw][n][:w]
    wr = pm.var[:nw][n][:wr]
    wi = pm.var[:nw][n][:wi]

    for (i,j) in keys(pm.ref[:nw][n][:buspairs])
        @constraint(pm.model, wr[(i,j)]^2 + wi[(i,j)]^2 == w[i]*w[j])
        @NLconstraint(pm.model, wi[(i,j)] == tan(t[i] - t[j])*wr[(i,j)])
    end
end

"""
```
t[f_bus] - t[t_bus] <= angmax
t[f_bus] - t[t_bus] >= angmin
```
"""
function constraint_voltage_angle_difference(pm::GenericPowerModel{T}, n::Int, f_bus, t_bus, angmin, angmax) where T <: StandardACTForm
    va_fr = pm.var[:nw][n][:va][f_bus]
    va_to = pm.var[:nw][n][:va][t_bus]

    @constraint(pm.model, va_fr - va_to <= angmax)
    @constraint(pm.model, va_fr - va_to >= angmin)
end

""
function add_bus_voltage_setpoint(sol, pm::GenericPowerModel{T}) where T <: AbstractACTForm
    add_setpoint(sol, pm, "bus", "vm", :w; scale = (x,item) -> sqrt(x))
    add_setpoint(sol, pm, "bus", "va", :va)
end


text:
@testset "Hilbert symbols" begin
  v = [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, 1, 1,
       -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, 1, 1, 1, -1, 1, -1, 1, -1, 1, 1,
       -1, -1, -1, -1, 1, 1, -1, 1, 1, -1, -1, 1, 1, 1, -1, 1, -1, -1, 1, 1 ]

  for i in 1:8
    for j in 1:8
      @test hilbert_symbol(i, j, 2) == v[(i - 1) * 8 + j]
      @test hilbert_symbol(fmpz(i), fmpz(j), 2) == v[(i - 1) * 8 + j]
      @test hilbert_symbol(fmpz(i), fmpz(j), fmpz(2)) == v[(i - 1) * 8 + j]
      @test hilbert_symbol(fmpq(i), fmpq(j), 2) == v[(i - 1) * 8 + j]
      @test hilbert_symbol(fmpq(i), fmpq(j), fmpz(2)) == v[(i - 1) * 8 + j]
    end
  end

  for p in PrimesSet(3, 100)
    for a in 1:100
      for b in 1:00
        h = hilbert_symbol(a, b, p)
        a = fmpz(a)
        b = fmpz(b)
        r = (-1)^(valuation(a, p) * valuation(b, p)) * a^(valuation(b, p)) * b^(valuation(a, p))
        @test h == jacobi(r, fmpz(p))
      end
    end
  end

  Qx, x = PolynomialRing(FlintQQ, "x")
  K, b = NumberField(x^3-3*x-1, "a")
  OK = maximal_order(K)
  for P in prime_ideals_up_to(OK, 200)
    @test hilbert_symbol(b, -3, P) == 1
  end

  # Test where Magma div(x, y) differs from julia div(x, y) (internally)
  K, a = MaximalRealSubfield(8, "a") # x^2 - 2
  z = 9278908160780559301//4*a+6561375391013480455//2
  w = K(-2)
  p = prime_decomposition(maximal_order(K), 2)[1][1]
  @test hilbert_symbol(z, w, p) == 1
end


text:
function setuptaxonomy()
    #set up taxonomy example
    #observations
    w_strings = ["Laptop","Monitor","Gamepad", #COMPUTERS
                 "Coffee machine","Vacuum cleaner","Electric toothbrush", #small APPLIANCES
                 "Grapes","Strawberries","Limes", #FRUIT
                 "Pancake mix","Baking soda","Baker's yeast","Muffin cups", #BAKING
                ]

    num_obs = size(w_strings,1)
    w_vec = collect(1:num_obs)

    #actions
    a_strings = ["Laptop sleeve","Monitor cable","Video game",
                 "Coffee","Vacuum cleaner bags","Brush heads",
                 "Cheese","Cream","Cane sugar",
                 "Maple syrup","Vinegar","Flour","Chocolate chips",
                 "COMPUTERS","APPLIANCES","FRUIT","BAKING","Electronics","Food"]
    num_acts = size(a_strings,1)
    a_vec = collect(1:num_acts)

    #set up uniform p(w)
    p_w = ones(num_obs)/num_obs

    #define utility function
    #(everything is hardcoded here, which is a bit hacky but will do the job)
    #the function expects integer indices
    function U(a::Integer,w::Integer)
        u_correct_w = 3;
        u_correct_category = 2.2;
        u_correct_supercategory = 1.6;

        #correct item
        if a<14 && a==w
            return u_correct_w
        end

        #flour is also fine for o=muffin cups
        if w==13 && a==12
            return u_correct_w
        end

        #For pancake mix both FRUIT and BAKING is fine
        if w==10 && a==16
            return u_correct_category
        end

        #extra if-clause is required for muffin cups
        if w==13 && a==17
            return u_correct_category
        end

        #correct category
        if a<18
            cat = ceil(w/3)
            if (a-13) == cat
                return u_correct_category
            end
        end

        #correct supercategory
        supcat = ceil(w/6)
        if (a-17) == supcat
            return u_correct_supercategory
        end

        #separate case for a==19
        if a==19 && w==13
           return u_correct_supercategory
        end

        #incorrect action
        return 0

    end

    return w_vec, w_strings, a_vec, a_strings, p_w, U
end

function setuptaxonomy_animals_plants()
    #set up taxonomy example
    #observations
    w_strings = ["Persian","Siamese","British Shorthair", #CATS
                 "German Shepherd","Rottweiler","Dachshund", #DOGS
                 "Oak","Birch","Pine", #TREES
                 "Rose","Dandelion","Sunflower", #FLOWERS
                 "Hibiscus" #Tree and flower (to break symmetry)
                ]

    num_obs = size(w_strings,1)
    w_vec = collect(1:num_obs)

    #actions are observations + categories + supercategories
    a_strings = [w_strings; "CAT";"DOG";"TREE";"FLOWER";"Animal";"Plant"]
    num_acts = size(a_strings,1)
    a_vec = collect(1:num_acts)

    #set up uniform p(w)
    p_w = ones(num_obs)/num_obs

    #define utility function
    #(everything is hardcoded here, which is a bit hacky but will do the job)
    function U(a,w)
        u_correct_w = 3;
        u_correct_category = 2.2;
        u_correct_supercategory = 1.6;

        #correct observation
        if w==a
            return u_correct_w
        end

        #correct category
        if w < 13
            cat = ceil(w/3)
            if (a-13) == cat
                return u_correct_category
            end
        end

        if w == 13  #hibiscus case
            if (a==16 || a==17)
                return u_correct_category
            end
        end

        #correct supercategory
        supcat = ceil(w/6)
        if (a-17) == supcat
            return u_correct_supercategory
        end

        if w == 13 #hibiscus case
            if a==19
                return u_correct_supercategory
            end
        end

        #incorrect action
        return 0

    end

    return w_vec, w_strings, a_vec, a_strings, p_w, U
end


text:
import RowEchelon

#=
Homework 3
Section 1.7
Problem 3
=#

# Create the matrix
A = [2 -4 ; -3 6]

# Row Reduce the matrix
A_rref = RowEchelon.rref(A)

# Determine whether there is a row of all zeroes
hasAllZeroRow = false
index = 1
currentZeros = 0
while index <= size(A_rref)[2]
	global index
	global hasAllZeroRow
	global currentZeros

	# Parse through each entry of the rows, searching for an all 0 row.
	for entry in A_rref[index, :]
		global hasAllZeroRow
		global currentZeros

		if entry == 0
			# If here then:
			# Current entry is a zero, increment zero counter
			currentZeros += 1
		end
	end
	if currentZeros == size(A_rref)[2]
		# If here then:
		# A row with all zeroes was found
		hasAllZeroRow = true
	end
	currentZeros = 0
	index += 1
end

if hasAllZeroRow == true
	# If here then:
	# The vectors are not linearly independent
	println("These two vectors are linearly dependent")
else
	# If here then:
	# The vectors are not linearly independent
	println("These two vectors are linearly independent")
end


text:
using Winston
using Compat
using Compat.Test
using Compat.Printf
using Colors
import Cairo
using Random
using Printf

include("canvas.jl")
include("crashes.jl")
include("compare.jl")


text:
### A Pluto.jl notebook ###
# v0.17.5

using Markdown
using InteractiveUtils

# ╔═╡ 587f852c-f357-4260-ae87-a60eaf101d36
# hideall
let
    docs_dir = dirname(dirname(@__DIR__))
    pkg_dir = dirname(docs_dir)

    using Pkg: Pkg
    Pkg.activate(docs_dir)
    Pkg.develop(; path=pkg_dir)
    Pkg.instantiate()

    # Putting the include here to avoid Pluto getting confused about cell order.
    include(joinpath(docs_dir, "src", "tutorials_utils.jl"))
end;

# ╔═╡ 206214e5-ab7a-4cae-bc85-7a77ed909c39
using CSV

# ╔═╡ 7f497bba-b0cc-40ef-a32f-4f69fa1cc355
using DataFrames

# ╔═╡ 178595bd-d701-47c6-a33c-4a33bad90a4f
using TuringGLM

# ╔═╡ 27c478b2-730f-11ec-1fbd-0b63384eaaf4
md"""
Let's cover the **Linear Regression** example with the `kidiq` dataset (Gelman & Hill, 2007), which is data from a survey of adult American women and their respective children.
Dated from 2007, it has 434 observations and 4 variables:

* `kid_score`: child's IQ
* `mom_hs`: binary/dummy (0 or 1) if the child's mother has a high school diploma
* `mom_iq`: mother's IQ
* `mom_age`: mother's age
"""

# ╔═╡ 041e7dca-b554-450c-9f43-79494e54e48e
url = "https://github.com/TuringLang/TuringGLM.jl/raw/main/data/kidiq.csv"

# ╔═╡ c0c65854-57af-488d-b80b-d37601e8024b
kidiq = CSV.read(download(url), DataFrame)

# ╔═╡ 31287120-a04c-46b2-ba06-6168c3435a78
md"""
Using `kid_score` as dependent variable and `mom_hs` along with `mom_iq` as independent variables with a moderation (interaction) effect:
"""

# ╔═╡ 1e4c7e9e-b32a-4eca-8621-d3b748d45d9e
fm = @formula(kid_score ~ mom_hs * mom_iq)

# ╔═╡ 0cd3dfce-3001-4c6e-b550-ca6964f19e35
md"""
Let's create our CustomPrior object.
No need for the third (auxiliary) prior for this model so we leave it as `nothing`:
"""

# ╔═╡ 9766b37c-55a0-4f74-924e-66b92eab7429
priors = CustomPrior(Normal(0, 2.5), Normal(10, 20), nothing);

# ╔═╡ 56498ac7-3476-42eb-9c12-078562fff51d
md"""
We instantiate our model with `turing_model` without specifying any model, thus the default model will be used: `Gaussian()`.
Notice that we are specifying the `priors` keyword argument:
"""

# ╔═╡ 99cbb309-393f-4a13-9454-1dee747c88a6
model = turing_model(fm, kidiq; priors);

# ╔═╡ 3e227e48-259c-40db-90e4-cd62afc307b2
chn = sample(model, NUTS(), 2_000);

# ╔═╡ 954f928a-5fcb-4c14-a74c-09bfe92f50bd
# hide
plot_chains(chn)

# ╔═╡ 9bb7a28f-a980-4cbc-9931-b38df79ae94b
md"""
## References

Gelman, A., & Hill, J. (2007). Data analysis using regression and multilevel/hierarchical models. Cambridge university press.
"""

# ╔═╡ Cell order:
# ╠═587f852c-f357-4260-ae87-a60eaf101d36
# ╠═27c478b2-730f-11ec-1fbd-0b63384eaaf4
# ╠═206214e5-ab7a-4cae-bc85-7a77ed909c39
# ╠═7f497bba-b0cc-40ef-a32f-4f69fa1cc355
# ╠═041e7dca-b554-450c-9f43-79494e54e48e
# ╠═c0c65854-57af-488d-b80b-d37601e8024b
# ╠═178595bd-d701-47c6-a33c-4a33bad90a4f
# ╠═31287120-a04c-46b2-ba06-6168c3435a78
# ╠═1e4c7e9e-b32a-4eca-8621-d3b748d45d9e
# ╠═0cd3dfce-3001-4c6e-b550-ca6964f19e35
# ╠═9766b37c-55a0-4f74-924e-66b92eab7429
# ╠═56498ac7-3476-42eb-9c12-078562fff51d
# ╠═99cbb309-393f-4a13-9454-1dee747c88a6
# ╠═3e227e48-259c-40db-90e4-cd62afc307b2
# ╠═954f928a-5fcb-4c14-a74c-09bfe92f50bd
# ╠═9bb7a28f-a980-4cbc-9931-b38df79ae94b