text:
function OrdinaryDiffEq.compute_step!(fpsolver::FPSolver{<:NLFunctional},
                                      integrator::DDEIntegrator)
  # update ODE integrator to next time interval together with correct interpolation
  if fpsolver.iter == 1
    advance_ode_integrator!(integrator)
  else
    update_ode_integrator!(integrator)
  end

  compute_step_fixedpoint!(fpsolver, integrator)
end

function OrdinaryDiffEq.compute_step!(fpsolver::FPSolver{<:NLAnderson,false},
                                      integrator::DDEIntegrator)
  @unpack cache, iter = fpsolver
  @unpack aa_start = cache

  # perform Anderson acceleration
  previter = iter - 1
  if previter == aa_start
    # update cached values for next step of Anderson acceleration
    cache.dzold = cache.dz
    cache.z₊old = integrator.u
  elseif previter > aa_start
    # actually perform Anderson acceleration
    integrator.u = OrdinaryDiffEq.anderson(integrator.u, cache)
    integrator.destats.nsolve += 1
  end

  # update ODE integrator to next time interval together with correct interpolation
  if iter == 1
    advance_ode_integrator!(integrator)
  else
    # force recomputation of all interpolation data if Anderson acceleration was performed
    update_ode_integrator!(integrator, previter > aa_start)
  end

  # compute next step
  compute_step_fixedpoint!(fpsolver, integrator)
end

function OrdinaryDiffEq.compute_step!(fpsolver::FPSolver{<:NLAnderson,true},
                                      integrator::DDEIntegrator)
  @unpack cache, iter = fpsolver
  @unpack aa_start = cache

  # perform Anderson acceleration
  previter = iter - 1
  if previter == aa_start
    # update cached values for next step of Anderson acceleration
    @.. cache.dzold = cache.dz
    @.. cache.z₊old = integrator.u
  elseif previter > aa_start
    # actually perform Anderson acceleration
    OrdinaryDiffEq.anderson!(integrator.u, cache)
    integrator.destats.nsolve += 1
  end

  # update ODE integrator to next time interval together with correct interpolation
  if iter == 1
    advance_ode_integrator!(integrator)
  else
    # force recomputation of all interpolation data if Anderson acceleration was performed
    update_ode_integrator!(integrator, previter > aa_start)
  end

  # compute next step
  compute_step_fixedpoint!(fpsolver, integrator)
end

function compute_step_fixedpoint!(fpsolver::FPSolver{<:Union{NLFunctional,NLAnderson},
                                                     false},
                                  integrator::DDEIntegrator)
  @unpack t,opts = integrator
  @unpack cache = fpsolver
  ode_integrator = integrator.integrator

  # recompute next integration step
  OrdinaryDiffEq.perform_step!(integrator, integrator.cache, true)

  # compute residuals
  dz = integrator.u .- ode_integrator.u
  atmp = OrdinaryDiffEq.calculate_residuals(dz, ode_integrator.u, integrator.u,
                                            opts.abstol, opts.reltol, opts.internalnorm, t)

  # cache results
  if isdefined(cache, :dz)
    cache.dz = dz
  end

  opts.internalnorm(atmp, t)
end

function compute_step_fixedpoint!(fpsolver::FPSolver{<:Union{NLFunctional,NLAnderson},
                                                     true},
                                  integrator::DDEIntegrator)
  @unpack t,opts = integrator
  @unpack cache = fpsolver
  @unpack dz,atmp = cache
  ode_integrator = integrator.integrator

  # recompute next integration step
  OrdinaryDiffEq.perform_step!(integrator, integrator.cache, true)

  # compute residuals
  @.. dz = integrator.u - ode_integrator.u
  OrdinaryDiffEq.calculate_residuals!(atmp, dz, ode_integrator.u, integrator.u,
                                      opts.abstol, opts.reltol, opts.internalnorm, t)

  opts.internalnorm(atmp, t)
end

## resize!

function Base.resize!(fpcache::FPFunctionalCache, i::Int)
  resize!(fpcache.atmp, i)
  resize!(fpcache.dz, i)
  nothing
end

function Base.resize!(fpcache::FPAndersonCache, fpsolver::FPSolver{<:NLAnderson},
                      integrator::DDEIntegrator, i::Int)
  resize!(fpcache, fpsolver.alg, i)
end

function Base.resize!(fpcache::FPAndersonCache, fpalg::NLAnderson, i::Int)
  @unpack z₊old,Δz₊s = fpcache

  resize!(fpcache.atmp, i)
  resize!(fpcache.dz, i)
  resize!(fpcache.dzold, i)
  resize!(z₊old, i)

  # update history of Anderson cache
  max_history_old = length(Δz₊s)
  max_history = min(fpalg.max_history, fpalg.max_iter, i)

  resize!(fpcache.γs, max_history)
  resize!(fpcache.Δz₊s, max_history)

  if max_history != max_history_old
    fpcache.Q = typeof(fpcache.Q)(undef, i, max_history)
    fpcache.R = typeof(fpcache.R)(undef, max_history, max_history)
  end

  max_history = length(Δz₊s)
  if max_history > max_history_old
    for i in (max_history_old + 1):max_history
      Δz₊s[i] = zero(z₊old)
    end
  end

  nothing
end

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
using Documenter
using StatsBase
using MixedModels

makedocs(
    root = joinpath(dirname(pathof(MixedModels)), "..", "docs"),
    sitename = "MixedModels",
    pages = [
        "index.md",
        "constructors.md",
        "optimization.md",
        "GaussHermite.md",
        "bootstrap.md",
        # "SimpleLMM.md",
        # "MultipleTerms.md",
        # "SingularCovariance.md",
        # "SubjectItem.md",
        # "benchmarks.md"
    ],
)

deploydocs(repo = "github.com/JuliaStats/MixedModels.jl.git", push_preview = true)

ext: jl


text:
using CUDA
export randomMPS
function randomMPS(L::Int,S::Int,D::Int;FloatType=Float64,ArrayType=Array)
    mps = MPS(L,S,D;FloatType=Float64,ArrayType=Array)
    for i = 1:L
        mps[i] .= rand(size(mps[i])...)
    end
    if ArrayType==CuArray
        mps = cu(mps)
    end
    return mps
end

export randomMPO
function randomMPO(L::Int,S::Int,D::Int;FloatType=Float64,ArrayType=Array)
    mpo = MPO(L,S,D;FloatType=Float64,ArrayType=Array)
    for i = 1:L
        mpo[i] .= rand(size(mpo[i])...)
    end
    if ArrayType==CuArray
        mpo = cu(mpo)
    end
    return mpo
end

max_stars_repo_name: qiyang-ustc/JMPS

max_issues_repo_name: qiyang-ustc/JMPS

max_forks_repo_name: qiyang-ustc/JMPS


text:
module PolyUtils

using ForwardDiff

# taken from Jacobi.lagrange
function lagrange(i, x, z)
    nz = length(z)

    l = one(z[1])

    for k = 1:(i-1)
        l = l * (x-z[k]) / (z[i]-z[k])
    end

    for k = (i+1):nz
        l = l * (x-z[k]) / (z[i]-z[k])
    end

    return l
end

dlagrange(i, x, z) = ForwardDiff.derivative(x -> lagrange(i, x, z), x)

end  # module


text:
struct PMmap
    bus::Dict{Int, PSY.Bus}
    arcs::Dict{
        NamedTuple{(:from_to, :to_from), Tuple{Tuple{Int, Int, Int}, Tuple{Int, Int, Int}}},
        t where t <: PSY.ACBranch,
    }
    arcs_dc::Dict{
        NamedTuple{(:from_to, :to_from), Tuple{Tuple{Int, Int, Int}, Tuple{Int, Int, Int}}},
        t where t <: PSY.DCBranch,
    }
end

function get_branch_to_pm(ix::Int, branch::PSY.PhaseShiftingTransformer)
    PM_branch = Dict{String, Any}(
        "br_r" => PSY.get_r(branch),
        "rate_a" => PSY.get_rate(branch),
        "shift" => PSY.get_α(branch),
        "rate_b" => PSY.get_rate(branch),
        "br_x" => PSY.get_x(branch),
        "rate_c" => PSY.get_rate(branch),
        "g_to" => 0.0,
        "g_fr" => 0.0,
        "b_fr" => PSY.get_primary_shunt(branch) / 2,
        "f_bus" => PSY.get_number(PSY.get_arc(branch).from),
        "br_status" => Float64(PSY.get_available(branch)),
        "t_bus" => PSY.get_number(PSY.get_arc(branch).to),
        "b_to" => PSY.get_primary_shunt(branch) / 2,
        "index" => ix,
        "angmin" => -π / 2,
        "angmax" => π / 2,
        "transformer" => true,
        "tap" => PSY.get_tap(branch),
    )
    return PM_branch
end

function get_branch_to_pm(ix::Int, branch::PSY.Transformer2W)
    PM_branch = Dict{String, Any}(
        "br_r" => PSY.get_r(branch),
        "rate_a" => PSY.get_rate(branch),
        "shift" => 0.0,
        "rate_b" => PSY.get_rate(branch),
        "br_x" => PSY.get_x(branch),
        "rate_c" => PSY.get_rate(branch),
        "g_to" => 0.0,
        "g_fr" => 0.0,
        "b_fr" => PSY.get_primary_shunt(branch) / 2,
        "f_bus" => PSY.get_number(PSY.get_arc(branch).from),
        "br_status" => Float64(PSY.get_available(branch)),
        "t_bus" => PSY.get_number(PSY.get_arc(branch).to),
        "b_to" => PSY.get_primary_shunt(branch) / 2,
        "index" => ix,
        "angmin" => -π / 2,
        "angmax" => π / 2,
        "transformer" => true,
        "tap" => 1.0,
    )
    return PM_branch
end

function get_branch_to_pm(ix::Int, branch::PSY.TapTransformer)
    PM_branch = Dict{String, Any}(
        "br_r" => PSY.get_r(branch),
        "rate_a" => PSY.get_rate(branch),
        "shift" => 0.0,
        "rate_b" => PSY.get_rate(branch),
        "br_x" => PSY.get_x(branch),
        "rate_c" => PSY.get_rate(branch),
        "g_to" => 0.0,
        "g_fr" => 0.0,
        "b_fr" => PSY.get_primary_shunt(branch) / 2,
        "f_bus" => PSY.get_number(PSY.get_arc(branch).from),
        "br_status" => Float64(PSY.get_available(branch)),
        "t_bus" => PSY.get_number(PSY.get_arc(branch).to),
        "b_to" => PSY.get_primary_shunt(branch) / 2,
        "index" => ix,
        "angmin" => -π / 2,
        "angmax" => π / 2,
        "transformer" => true,
        "tap" => PSY.get_tap(branch),
    )
    return PM_branch
end

function get_branch_to_pm(ix::Int, branch::PSY.ACBranch)
    PM_branch = Dict{String, Any}(
        "br_r" => PSY.get_r(branch),
        "rate_a" => PSY.get_rate(branch),
        "shift" => 0.0,
        "rate_b" => PSY.get_rate(branch),
        "br_x" => PSY.get_x(branch),
        "rate_c" => PSY.get_rate(branch),
        "g_to" => 0.0,
        "g_fr" => 0.0,
        "b_fr" => PSY.get_b(branch).from,
        "f_bus" => PSY.get_number(PSY.get_arc(branch).from),
        "br_status" => Float64(PSY.get_available(branch)),
        "t_bus" => PSY.get_number(PSY.get_arc(branch).to),
        "b_to" => PSY.get_b(branch).to,
        "index" => ix,
        "angmin" => PSY.get_angle_limits(branch).min,
        "angmax" => PSY.get_angle_limits(branch).max,
        "transformer" => false,
        "tap" => 1.0,
    )
    return PM_branch
end

function get_branch_to_pm(ix::Int, branch::PSY.HVDCLine)
    PM_branch = Dict{String, Any}(
        "loss1" => PSY.get_loss(branch).l1,
        "mp_pmax" => PSY.get_reactive_power_limits_from(branch).max,
        "model" => 2,
        "shutdown" => 0.0,
        "pmaxt" => PSY.get_active_power_limits_to(branch).max,
        "pmaxf" => PSY.get_active_power_limits_from(branch).max,
        "startup" => 0.0,
        "loss0" => PSY.get_loss(branch).l1,
        "pt" => 0.0,
        "vt" => PSY.get_magnitude(PSY.get_arc(branch).to),
        "qmaxf" => PSY.get_reactive_power_limits_from(branch).max,
        "pmint" => PSY.get_active_power_limits_to(branch).min,
        "f_bus" => PSY.get_number(PSY.get_arc(branch).from),
        "mp_pmin" => PSY.get_reactive_power_limits_from(branch).min,
        "br_status" => Float64(PSY.get_available(branch)),
        "t_bus" => PSY.get_number(PSY.get_arc(branch).to),
        "index" => ix,
        "qmint" => PSY.get_reactive_power_limits_to(branch).min,
        "qf" => 0.0,
        "cost" => 0.0,
        "pminf" => PSY.get_active_power_limits_from(branch).min,
        "qt" => 0.0,
        "qminf" => PSY.get_reactive_power_limits_from(branch).min,
        "vf" => PSY.get_magnitude(PSY.get_arc(branch).from),
        "qmaxt" => PSY.get_reactive_power_limits_to(branch).max,
        "ncost" => 0,
        "pf" => 0.0,
    )
    return PM_branch
end

function get_branches_to_pm(sys::PSY.System)
    PM_ac_branches = Dict{String, Any}()
    PM_dc_branches = Dict{String, Any}()
    PMmap_ac = Dict{
        NamedTuple{(:from_to, :to_from), Tuple{Tuple{Int, Int, Int}, Tuple{Int, Int, Int}}},
        t where t <: PSY.ACBranch,
    }()
    PMmap_dc = Dict{
        NamedTuple{(:from_to, :to_from), Tuple{Tuple{Int, Int, Int}, Tuple{Int, Int, Int}}},
        t where t <: PSY.DCBranch,
    }()

    for (ix, branch) in enumerate(get_available_components(PSY.Branch, sys))
        if isa(branch, PSY.DCBranch)
            PM_dc_branches["$(ix)"] = get_branch_to_pm(ix, branch)
            if PM_dc_branches["$(ix)"]["br_status"] == true
                f = PM_dc_branches["$(ix)"]["f_bus"]
                t = PM_dc_branches["$(ix)"]["t_bus"]
                PMmap_dc[(from_to = (ix, f, t), to_from = (ix, t, f))] = branch
            end
        else
            PM_ac_branches["$(ix)"] = get_branch_to_pm(ix, branch)
            if PM_ac_branches["$(ix)"]["br_status"] == true
                f = PM_ac_branches["$(ix)"]["f_bus"]
                t = PM_ac_branches["$(ix)"]["t_bus"]
                PMmap_ac[(from_to = (ix, f, t), to_from = (ix, t, f))] = branch
            end
        end
    end

    return PM_ac_branches, PM_dc_branches, PMmap_ac, PMmap_dc
end

function get_buses_to_pm(buses::IS.FlattenIteratorWrapper{PSY.Bus})
    PM_buses = Dict{String, Any}()
    PMmap_buses = Dict{Int, PSY.Bus}()

    pm_bustypes = Dict{PSY.BusTypes.BusType, Int}(
        PSY.BusTypes.ISOLATED => 4,
        PSY.BusTypes.PQ => 1,
        PSY.BusTypes.PV => 2,
        PSY.BusTypes.REF => 3,
        PSY.BusTypes.SLACK => 3,
    )

    for bus in buses
        number = PSY.get_number(bus)
        PM_bus = Dict{String, Any}(
            "zone" => 1,
            "bus_i" => number,
            "bus_type" => pm_bustypes[PSY.get_bustype(bus)],
            "vmax" => PSY.get_voltage_limits(bus).max,
            "area" => 1,
            "vmin" => PSY.get_voltage_limits(bus).min,
            "index" => PSY.get_number(bus),
            "va" => PSY.get_angle(bus),
            "vm" => PSY.get_magnitude(bus),
            "base_kv" => PSY.get_base_voltage(bus),
            "pni" => 0.0,
            "qni" => 0.0,
            "name" => PSY.get_name(bus),
        )
        PM_buses["$(number)"] = PM_bus
        if PSY.get_bustype(bus) != PSY.BusTypes.ISOLATED::PSY.BusTypes.BusType
            PMmap_buses[number] = bus
        end
    end
    return PM_buses, PMmap_buses
end

function pass_to_pm(sys::PSY.System, time_periods::Int)

    ac_lines, dc_lines, PMmap_ac, PMmap_dc = get_branches_to_pm(sys)
    buses = PSY.get_components(PSY.Bus, sys)
    pm_buses, PMmap_buses = get_buses_to_pm(buses)
    PM_translation = Dict{String, Any}(
        "bus" => pm_buses,
        "branch" => ac_lines,
        "baseMVA" => PSY.get_base_power(sys),
        "per_unit" => true,
        "storage" => Dict{String, Any}(),
        "dcline" => dc_lines,
        "gen" => Dict{String, Any}(),
        "switch" => Dict{String, Any}(),
        "shunt" => Dict{String, Any}(),
        "load" => Dict{String, Any}(),
    )

    # TODO: this function adds overhead in large number of time_steps
    # We can do better later.

    PM_translation = PM.replicate(PM_translation, time_periods)

    PM_map = PMmap(PMmap_buses, PMmap_ac, PMmap_dc)

    return PM_translation, PM_map
end


text:
#!/usr/bin/env julia

using Base.Test
using MathProgBase
using Logging
using Drs

A = Float64[3 2 1; 2 5 3]
b = Float64[10, 15]
c = Float64[-2, -3, -4]

s = linprog(c, A, '<', b, -Inf, Inf, DrsMathProgSolver())

@test s.status == :Optimal
@test s.objval == -20
@test s.sol == [0, 0, 5]

max_stars_repo_name: abogushov/lplib

max_issues_repo_name: abogushov/lplib

max_forks_repo_name: abogushov/lplib


text:
"""
`SolutionPDAE`: Solution of a partitioned differential algebraic equation

Contains all fields necessary to store the solution of an PDAE.

### Fields

* `nd`: dimension of the dynamical variable ``q``
* `nm`: dimension of the constraint submanifold
* `nt`: number of time steps to store
* `ni`: number of initial conditions
* `t`:  time steps
* `q`:  solution `q[nd, nt+1, ni]` with `q[:,0,:]` the initial conditions
* `p`:  solution `p[nd, nt+1, ni]` with `p[:,0,:]` the initial conditions
* `λ`:  Lagrange multiplier `λ[nd, nt+1, ni]`
* `ntime`: number of time steps to compute
* `nsave`: store every nsave'th time step (default: 1)
* `nwrite`: save data to disk after every nwrite'th time step (default: ntime)
* `counter`: counter for copied solution entries
* `woffset`: counter for file offset
* `h5`: HDF5 file for storage
"""
abstract type SolutionPDAE{dType, tType, N} <: DeterministicSolution{dType, tType, N} end

# Create SolutionPDAEs with serial and parallel data structures.
for (TSolution, TDataSeries, Tdocstring) in
    ((:SSolutionPDAE, :SDataSeries, "Serial Solution of a partitioned differential algebraic equation."),
     (:PSolutionPDAE, :PDataSeries, "Parallel Solution of a partitioned differential algebraic equation."))
    @eval begin
        $Tdocstring
        mutable struct $TSolution{dType, tType, N} <: SolutionPDAE{dType, tType, N}
            nd::Int
            nm::Int
            nt::Int
            ni::Int
            t::TimeSeries{tType}
            q::$TDataSeries{dType,N}
            p::$TDataSeries{dType,N}
            λ::$TDataSeries{dType,N}
            ntime::Int
            nsave::Int
            nwrite::Int
            counter::Vector{Int}
            woffset::Int
            h5::HDF5File

            function $TSolution{dType, tType, N}(nd, nm, nt, ni, t, q, p, λ, ntime, nsave, nwrite) where {dType <: Number, tType <: Real, N}
                new(nd, nm, nt, ni, t, q, p, λ, ntime, nsave, nwrite, zeros(Int, ni), 0)
            end
        end

        function $TSolution(equation::Union{IODE{DT,TT},VODE{DT,TT},PDAE{DT,TT},IDAE{DT,TT},VDAE{DT,TT}}, Δt::TT, ntime::Int; nsave::Int=DEFAULT_NSAVE, nwrite::Int=DEFAULT_NWRITE, filename=nothing) where {DT,TT}
            @assert nsave > 0
            @assert ntime == 0 || ntime ≥ nsave
            @assert nwrite == 0 || nwrite ≥ nsave
            @assert mod(ntime, nsave) == 0

            if nwrite > 0
                @assert mod(nwrite, nsave) == 0
                @assert mod(ntime, nwrite) == 0
            end

            N  = equation.n > 1 ? 3 : 2
            nd = equation.d
            nm = equation.m
            ni = equation.n
            nt = div(ntime, nsave)
            nt = (nwrite == 0 ? nt : div(nwrite, nsave))
            nw = (nwrite == 0 ? ntime : nwrite)

            @assert nd > 0
            @assert nm > 0
            @assert ni > 0
            @assert nt ≥ 0
            @assert nw ≥ 0

            t = TimeSeries{TT}(nt, Δt, nsave)
            q = $TDataSeries(DT, nd, nt, ni)
            p = $TDataSeries(DT, nd, nt, ni)
            λ = $TDataSeries(DT, nm, nt, ni)
            s = $TSolution{DT,TT,N}(nd, nm, nt, ni, t, q, p, λ, ntime, nsave, nw)
            set_initial_conditions!(s, equation)

            if !isnothing(filename)
                isfile(filename) ? @warn("Overwriting existing HDF5 file.") : nothing
                create_hdf5!(s, filename)
            end

            return s
        end

        function $TSolution(t::TimeSeries{TT}, q::$TDataSeries{DT,N}, p::$TDataSeries{DT,N}, λ::$TDataSeries{DT,N}, ntime::Int) where {DT,TT,N}
            @assert q.nd == p.nd >= λ.nd
            @assert q.nt == p.nt == λ.nt
            @assert q.ni == p.ni == λ.ni

            # extract parameters
            nd = q.nd
            nm = λ.nd
            ni = q.ni
            nt = t.n
            ns = div(ntime, nt)

            @assert mod(ntime, nt) == 0

            # create solution
            $TSolution{DT,TT,N}(nd, nm, nt, ni, t, q, p, λ, ntime, ns, 0)
        end

        function $TSolution(file::String)
            # open HDF5 file
            get_config(:verbosity) > 1 ? @info("Reading HDF5 file ", file) : nothing
            h5 = h5open(file, "r")

            # read attributes
            ntime = read(attrs(h5)["ntime"])
            nsave = read(attrs(h5)["nsave"])

            # reading data arrays
            t = TimeSeries(read(h5["t"]), nsave)
            q = $TDataSeries(read(h5["q"]))
            p = $TDataSeries(read(h5["p"]))
            λ = $TDataSeries(read(h5["λ"]))

            # need to close the file
            close(h5)

            # create solution
            $TSolution(t, q, p, λ, ntime)
        end
    end
end

Base.:(==)(sol1::SolutionPDAE{DT1,TT1,N1}, sol2::SolutionPDAE{DT2,TT2,N2}) where {DT1,TT1,N1,DT2,TT2,N2} = (
                                DT1 == DT2
                             && TT1 == TT2
                             && N1  == N2
                             && sol1.nd == sol2.nd
                             && sol1.nm == sol2.nm
                             && sol1.nt == sol2.nt
                             && sol1.ni == sol2.ni
                             && sol1.t  == sol2.t
                             && sol1.q  == sol2.q
                             && sol1.p  == sol2.p
                             && sol1.λ  == sol2.λ
                             && sol1.ntime == sol2.ntime
                             && sol1.nsave == sol2.nsave
                             && sol1.nwrite == sol2.nwrite
                             && sol1.counter == sol2.counter
                             && sol1.woffset == sol2.woffset)

hdf5(sol::SolutionPDAE)  = sol.h5
timesteps(sol::SolutionPDAE)  = sol.t
ntime(sol::SolutionPDAE) = sol.ntime
nsave(sol::SolutionPDAE) = sol.nsave
offset(sol::SolutionPDAE) = sol.woffset

function set_initial_conditions!(sol::SolutionPDAE, equ::Union{IODE,VODE,PDAE,IDAE,VDAE})
    set_initial_conditions!(sol, equ.t₀, equ.q₀, equ.p₀, equ.λ₀)
end

function set_initial_conditions!(sol::SolutionPDAE{DT,TT}, t₀::TT, q₀::Union{Array{DT}, Array{TwicePrecision{DT}}}, p₀::Union{Array{DT}, Array{TwicePrecision{DT}}}, λ₀::Union{Array{DT}, Array{TwicePrecision{DT}}}) where {DT,TT}
    set_data!(sol.q, q₀, 0)
    set_data!(sol.p, p₀, 0)
    set_data!(sol.λ, λ₀, 0)
    compute_timeseries!(sol.t, t₀)
    sol.counter .= 1
end

function get_initial_conditions!(sol::SolutionPDAE{DT,TT}, asol::AtomicSolutionPDAE{DT,TT}, k, n=1) where {DT,TT}
    get_solution!(sol, asol.q, asol.p, asol.λ, n-1, k)
    asol.t  = sol.t[n-1]
    asol.q̃ .= 0
    asol.p̃ .= 0
end

function get_initial_conditions!(sol::SolutionPDAE{DT,TT}, q::SolutionVector{DT}, p::SolutionVector{DT}, λ::SolutionVector{DT}, k, n=1) where {DT,TT}
    get_solution!(sol, q, p, λ, n-1, k)
end

function get_initial_conditions!(sol::SolutionPDAE{DT,TT}, q::SolutionVector{DT}, p::SolutionVector{DT}, k, n=1) where {DT,TT}
    get_solution!(sol, q, p, n-1, k)
end

function get_initial_conditions(sol::SolutionPDAE, k, n=1)
    get_solution(sol, n-1, k)
end

function get_solution!(sol::SolutionPDAE{DT,TT}, q::SolutionVector{DT}, p::SolutionVector{DT}, λ::SolutionVector{DT}, n, k=1) where {DT,TT}
    q .= sol.q[:, n, k]
    p .= sol.p[:, n, k]
    λ .= sol.λ[:, n, k]
end

function get_solution!(sol::SolutionPDAE{DT,TT}, q::SolutionVector{DT}, p::SolutionVector{DT}, n, k=1) where {DT,TT}
    q .= sol.q[:, n, k]
    p .= sol.p[:, n, k]
end

function get_solution(sol::SolutionPDAE, n, k=1)
    (sol.t[n], sol.q[:, n, k], sol.p[:, n, k], sol.λ[:, n, k])
end

function set_solution!(sol::SolutionPDAE, t::Real, q::Vector, p::Vector, n::Int, k::Int=1)
    set_solution!(sol, q, p, n, k)
end

function set_solution!(sol::SolutionPDAE, t::Real, q::Vector, p::Vector, λ::Vector, n::Int, k::Int=1)
    set_solution!(sol, q, p, λ, n, k)
end

function set_solution!(sol::SolutionPDAE{DT,TT}, asol::AtomicSolutionPDAE{DT,TT}, n::Int, k::Int=1) where {DT,TT}
    set_solution!(sol, asol.t, asol.q, asol.p, asol.λ, n, k)
end

function set_solution!(sol::SolutionPDAE{DT,TT}, q::SolutionVector{DT}, p::SolutionVector{DT}, λ::SolutionVector{DT}, n::Int, k::Int=1) where {DT,TT}
    @assert n <= sol.ntime
    @assert k <= sol.ni
    if mod(n, sol.nsave) == 0
        if sol.counter[k] > sol.nt
            @error("Solution overflow. Call write_to_hdf5() and reset!() before continuing the simulation.")
        end
        set_data!(sol.q, q, sol.counter[k], k)
        set_data!(sol.p, p, sol.counter[k], k)
        set_data!(sol.λ, λ, sol.counter[k], k)
        sol.counter[k] += 1
    end
end

function set_solution!(sol::SolutionPDAE{DT,TT}, q::SolutionVector{DT}, p::SolutionVector{DT}, n::Int, k::Int=1) where {DT,TT}
    @assert n <= sol.ntime
    @assert k <= sol.ni
    if mod(n, sol.nsave) == 0
        if sol.counter[k] > sol.nt
            @error("Solution overflow. Call write_to_hdf5() and reset!() before continuing the simulation.")
        end
        set_data!(sol.q, q, sol.counter[k], k)
        set_data!(sol.p, p, sol.counter[k], k)
        sol.counter[k] += 1
    end
end

function CommonFunctions.reset!(sol::SolutionPDAE)
    reset!(sol.q)
    reset!(sol.p)
    reset!(sol.λ)
    compute_timeseries!(sol.t, sol.t[end])
    sol.counter .= 1
    sol.woffset += sol.nt
end

function Base.close(solution::SolutionPDAE)
    close(solution.h5)
end


text:
"""
    validate_grid(a::AbstractField, b::AbstractField)

Confirm that `a` and `b` are on the same grid if both are fields and return `a.grid`.
"""
function validate_grid(a::AbstractField, b::AbstractField)
    a.grid === b.grid || throw(ArgumentError("Fields in an AbstractOperation must be on the same grid."))
    return a.grid
end

"""Return `a.grid` when `b` is not an `AbstractField`."""
validate_grid(a::AbstractField, b) = a.grid

"""Return `b.grid` when `a` is not an `AbstractField`."""
validate_grid(a, b::AbstractField) = b.grid

"""Fallback when neither `a` nor `b` are `AbstractField`s."""
validate_grid(a, b) = nothing

"""
    validate_grid(a, b, c...)

Confirm that the grids associated with the 3+ long list `a, b, c...` are
consistent by checking each member against `a`.
This function is only correct when `a` is an `AbstractField`, though the
subsequent members `b, c...` may be anything.
"""
function validate_grid(a, b, c, d...)
    grids = []
    push!(grids, validate_grid(a, b))
    push!(grids, validate_grid(a, c))
    append!(grids, [validate_grid(a, di) for di in d])

    for g in grids
        if !(g === nothing)
            return g
        end
    end

    return nothing
end


text:
function check_derivative(f, g, h, θ)
       a = (f.(θ .+ h) .- f(θ)) ./ h
       b = g(θ)
       return a,b
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using ReverseDiff

# Set parameters:
# ngram_range = 2 will add bi-grams features
const NGRAM_RANGE = 1
const MAX_FEATURES = 20000
const MAXLEN = 400
const BATCH_SIZE = 32
const EMBEDDING_DIMS = 50
const NB_EPOCH = 5

model = Sequential()

# we start off with an efficient embedding layer which maps
# our vocab indices into embedding_dims dimensions
model.add(Embedding(MAX_FEATURES,
                    EMBEDDING_DIms,
                    input_length=MAXLEN))

# we add a GlobalAveragePooling1D, which will average the embeddings
# of all words in the document
model.add(GlobalAveragePooling1D())

# We project onto a single unit output layer, and squash it with a sigmoid:
model.add(Dense(1, activation='sigmoid'))

model.compile(loss='binary_crossentropy',
              optimizer='adam',
              metrics=['accuracy'])


text:
@testset "Transformer" begin
  gpt = Gpt(300, 10, 500, 3)
  Flux.testmode!(gpt)

  x = randn(Float32, 300, 5, 2)
  x1 = x[:,:,1]

  y = gpt(x)
  y1 = gpt(x1)

  @test size(y) == (300, 5, 2)
  @test size(y1) == (300, 5)
  @test y[:,:,1] ≈ y1

  _, yall = gpt(x; all=true)
  @test gpt.ts.models[1](x) ≈ yall[1]
  @test gpt.ts.models[2](yall[1]) ≈ yall[2]
  @test gpt.ts.models[3](yall[2]) ≈ yall[3]

  _, yall1 = gpt(x1; all=true)
  @test gpt.ts.models[1](x1) ≈ yall1[1]
  @test gpt.ts.models[2](yall1[1]) ≈ yall1[2]
  @test gpt.ts.models[3](yall1[2]) ≈ yall1[3]

  mask = cat(reshape([1,1,1,1,1], 1, 5, 1),
             reshape([1,1,1,0,0], 1, 5, 1); dims=3)

  ym, ymall = gpt(x, mask; all=true)
  topo = @nntopo ((x, m) => x:(x, m)) => 3
  @test topo(gpt.ts.models, x, mask) .* mask ≈ ym
end


text:
"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function kmer_path_to_sequence(kmer_path)
    sequence = BioSequences.LongDNASeq(first(kmer_path))
    for kmer in kmer_path[2:end]
        @assert kmer[1] == sequence[end-1]
        @assert kmer[2] == sequence[end]
        push!(sequence, kmer[end])
    end
    return sequence
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function random_fasta_record(;seed=rand(Int), L = rand(0:Int(typemax(UInt16))))
    id = Random.randstring(Int(ceil(log(L + 1))))
    seq = BioSequences.randdnaseq(Random.seed!(seed), L)
    return FASTX.FASTA.Record(id, seq)
end

# """
# $(DocStringExtensions.TYPEDSIGNATURES)

# A short description of the function

# ```jldoctest
# julia> 1 + 1
# 2
# ```
# """
# function observe(records::AbstractVector{R}; outfile = "", error_rate = 0.0) where {R <: Union{FASTX.FASTA.Record, FASTX.FASTQ.Record}}
#     if isempty(outfile)
#         error("no file name supplied")
#     end
#     io = open(outfile, 'w')
#     fastx_io = FASTX.FASTQ.Writer(io)
#     for record in records
#         new_seq = observe(FASTX.sequence(record), error_rate=error_rate)
#         new_seq_id = string(hash(new_seq)) * "-" * Random.randstring(32)
#     new_seq_description = FASTX.identifier(record)
#     quality = fill(UInt8(60), length(new_seq))
#     return FASTX.FASTQ.Record(new_seq_id, new_seq_description, new_seq, quality)
# end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function observe(record::R; error_rate = 0.0) where {R <: Union{FASTX.FASTA.Record, FASTX.FASTQ.Record}}

    new_seq = observe(FASTX.sequence(record), error_rate=error_rate)
    new_seq_id = string(hash(new_seq)) * "-" * Random.randstring(32)
    new_seq_description = FASTX.identifier(record)
    quality = fill(UInt8(60), length(new_seq))
    return FASTX.FASTQ.Record(new_seq_id, new_seq_description, new_seq, quality)
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function observe(record::R; error_rate = 0.0) where {R <: Union{FASTX.FASTA.Record, FASTX.FASTQ.Record}}

    new_seq = observe(FASTX.sequence(record), error_rate=error_rate)
    new_seq_id = string(hash(new_seq)) * "-" * Random.randstring(32)
    new_seq_description = FASTX.identifier(record)
    quality = fill(UInt8(60), length(new_seq))
    return FASTX.FASTQ.Record(new_seq_id, new_seq_description, new_seq, quality)
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function observe(sequence::BioSequences.LongSequence{T}; error_rate = 0.0) where T

    if T <: BioSequences.DNAAlphabet
        alphabet = DNA_ALPHABET
    elseif T <: BioSequences.RNAAlphabet
        alphabet = RNA_ALPHABET
    else
        @assert T <: BioSequences.AminoAcidAlphabet
        alphabet = AA_ALPHABET
    end

    new_seq = BioSequences.LongSequence{T}()
    for character in sequence
        if rand() > error_rate
            # match
            push!(new_seq, character)
        else
            error_type = rand(1:3)
            if error_type == 1
                # mismatch
                push!(new_seq, rand(setdiff(alphabet, character)))
            elseif error_type == 2
                # insertion
                total_insertions = 1 + rand(Distributions.Poisson(error_rate))
                for i in 1:total_insertions
                    push!(new_seq, rand(alphabet))
                end
                push!(new_seq, character)
            else
                # deletion
                continue
            end
        end
    end
    if (T <: BioSequences.DNAAlphabet || T <: BioSequences.RNAAlphabet) && rand(Bool)
        BioSequences.reverse_complement!(new_seq)
    end
    return new_seq
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

Return proportion of matched bases in alignment to total matches + edits.

0-1, not %

```jldoctest
julia> 1 + 1
2
```
"""
function assess_alignment_accuracy(alignment_result)
    return alignment_result.total_matches / (alignment_result.total_matches + alignment_result.total_edits)
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

Used to determine which orientation provides an optimal alignment for initiating path likelihood analyses in viterbi analysis

```jldoctest
julia> 1 + 1
2
```
"""
function assess_optimal_kmer_alignment(kmer, observed_kmer)

    forward_alignment_result = assess_alignment(kmer, observed_kmer)
    forward_alignment_accuracy = assess_alignment_accuracy(forward_alignment_result)

    reverse_alignment_result = assess_alignment(kmer, BioSequences.reverse_complement(observed_kmer))
    reverse_alignment_accuracy = assess_alignment_accuracy(reverse_alignment_result)

    if forward_alignment_accuracy > reverse_alignment_accuracy
        alignment_result = forward_alignment_result
        orientation = true
    elseif forward_alignment_accuracy < reverse_alignment_accuracy
        alignment_result = reverse_alignment_result
        orientation = false
    elseif forward_alignment_accuracy == reverse_alignment_accuracy
        alignment_result, orientation = rand(((forward_alignment_result, missing), (reverse_alignment_result, missing)))
    end

    return (alignment_result, orientation)
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function assess_alignment(a, b)
    pairwise_alignment = BioAlignments.pairalign(BioAlignments.LevenshteinDistance(), a, b)
    alignment_result = BioAlignments.alignment(pairwise_alignment)
    total_aligned_bases = BioAlignments.count_aligned(alignment_result)
    total_matches = Int(BioAlignments.count_matches(alignment_result))
    total_edits = Int(total_aligned_bases - total_matches)
    return (total_matches = total_matches, total_edits = total_edits)
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function count_canonical_kmers(::Type{KMER_TYPE}, sequence::BioSequences.LongSequence) where KMER_TYPE
    canonical_kmer_counts = DataStructures.OrderedDict{KMER_TYPE, Int}()
    canonical_kmer_iterator = (BioSequences.canonical(kmer.fw) for kmer in BioSequences.each(KMER_TYPE, sequence))
    for canonical_kmer in canonical_kmer_iterator
        canonical_kmer_counts[canonical_kmer] = get(canonical_kmer_counts, canonical_kmer, 0) + 1
    end
    return canonical_kmer_counts
end

function count_canonical_kmers(::Type{KMER_TYPE}, record::R) where {KMER_TYPE, R <: Union{FASTX.FASTA.Record, FASTX.FASTQ.Record}}
    return count_canonical_kmers(KMER_TYPE, FASTX.sequence(record))
end

function count_canonical_kmers(::Type{KMER_TYPE}, sequences::AbstractVector{T}) where {KMER_TYPE, T <: Union{FASTX.FASTA.Record, FASTX.FASTQ.Record}}
    joint_kmer_counts = DataStructures.OrderedDict{KMER_TYPE, Int}()
    for sequence in sequences
        sequence_kmer_counts = count_canonical_kmers(KMER_TYPE, sequence)
        merge!(+, joint_kmer_counts, sequence_kmer_counts)
    end
    sort!(joint_kmer_counts)
end

function count_canonical_kmers(::Type{KMER_TYPE}, sequences::R) where {KMER_TYPE, R <: Union{FASTX.FASTA.Reader, FASTX.FASTQ.Reader}}
    joint_kmer_counts = DataStructures.OrderedDict{KMER_TYPE, Int}()
    for sequence in sequences
        sequence_kmer_counts = count_canonical_kmers(KMER_TYPE, sequence)
        merge!(+, joint_kmer_counts, sequence_kmer_counts)
    end
    sort!(joint_kmer_counts)
end

function count_canonical_kmers(::Type{KMER_TYPE}, fastx_files::AbstractVector{S}) where {KMER_TYPE, S <: AbstractString}
    kmer_counts = count_canonical_kmers(KMER_TYPE, first(fastx_files))
    for fastx_file in fastx_files[2:end]
        _kmer_counts = count_canonical_kmers(KMER_TYPE, fastx_file)
        kmer_counts = merge!(+, kmer_counts, _kmer_counts)
    end
    return kmer_counts
end

function count_canonical_kmers(::Type{KMER_TYPE}, fastx_file::S) where {KMER_TYPE, S <: AbstractString}
    fastx_io = open_fastx(fastx_file)
    kmer_counts = Mycelia.count_canonical_kmers(KMER_TYPE, fastx_io)
    close(fastx_io)
    return kmer_counts
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function count_kmers(::Type{KMER_TYPE}, sequence::BioSequences.LongSequence) where KMER_TYPE
    kmer_counts = DataStructures.OrderedDict{KMER_TYPE, Int}()
    kmer_iterator = (kmer.fw for kmer in BioSequences.each(KMER_TYPE, sequence))
    for kmer in kmer_iterator
        kmer_counts[kmer] = get(kmer_counts, kmer, 0) + 1
    end
    return kmer_counts
end

function count_kmers(::Type{KMER_TYPE}, record::R) where {KMER_TYPE, R <: Union{FASTX.FASTA.Record, FASTX.FASTQ.Record}}
    return count_kmers(KMER_TYPE, FASTX.sequence(record))
end

function count_kmers(::Type{KMER_TYPE}, sequences::AbstractVector{T}) where {KMER_TYPE, T <: Union{FASTX.FASTA.Record, FASTX.FASTQ.Record}}
    joint_kmer_counts = DataStructures.OrderedDict{KMER_TYPE, Int}()
    for sequence in sequences
        sequence_kmer_counts = count_kmers(KMER_TYPE, sequence)
        merge!(+, joint_kmer_counts, sequence_kmer_counts)
    end
    sort!(joint_kmer_counts)
end

function count_kmers(::Type{KMER_TYPE}, sequences::R) where {KMER_TYPE, R <: Union{FASTX.FASTA.Reader, FASTX.FASTQ.Reader}}
    joint_kmer_counts = DataStructures.OrderedDict{KMER_TYPE, Int}()
    for sequence in sequences
        sequence_kmer_counts = count_kmers(KMER_TYPE, sequence)
        merge!(+, joint_kmer_counts, sequence_kmer_counts)
    end
    sort!(joint_kmer_counts)
end

function count_kmers(KMER_TYPE, fastx_files::AbstractVector{AbstractString})
    kmer_counts = count_kmers(KMER_TYPE, first(files))
    for file in files[2:end]
        _kmer_counts = count_kmers(KMER_TYPE, file)
        kmer_counts = merge!(+, kmer_counts, _kmer_counts)
    end
    return kmer_counts
end

function count_kmers(KMER_TYPE, fastx_file::AbstractString)
    fastx_io = open_fastx(fastx_file)
    kmer_counts = Mycelia.count_kmers(KMER_TYPE, fastx_io)
    close(fastx_io)
    return kmer_counts
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function assess_kmer_saturation(fastxs, kmer_type; kmers_to_assess=Inf, power=10)
    canonical_kmers = Set{kmer_type}()

    max_possible_kmers = max_canonical_kmers(kmer_type)

    if kmers_to_assess == Inf
        kmers_to_assess = max_possible_kmers
    end

    sampling_points = Int[0]
    i = 0
    while power^i <= kmers_to_assess
        push!(sampling_points, power^i)
        i += 1
    end

    unique_kmer_counts = zeros(Int, length(sampling_points))

    if length(sampling_points) < 3
        @info "increase the # of reads analyzed or decrease the power to acquire more data points"
        return (;sampling_points, unique_kmer_counts)
    end

    p = ProgressMeter.Progress(kmers_to_assess, 1)

    kmers_assessed = 0
    for fastx in fastxs
        for record in Mycelia.open_fastx(fastx)
            for kmer in BioSequences.each(kmer_type, FASTX.sequence(record))
                canonical_kmer = kmer.fw < kmer.bw ? kmer.fw : kmer.bw
                push!(canonical_kmers, canonical_kmer)
                kmers_assessed += 1
                if (length(canonical_kmers) == max_possible_kmers)
                    sampling_points = vcat(filter(s -> s < kmers_assessed, sampling_points), [kmers_assessed])
                    unique_kmer_counts = vcat(unique_kmer_counts[1:length(sampling_points)-1], length(canonical_kmers))
                    return (;sampling_points, unique_kmer_counts, eof = false)
                elseif kmers_assessed in sampling_points
                    i = findfirst(sampling_points .== kmers_assessed)
                    unique_kmer_counts[i] = length(canonical_kmers)
                    if i == length(sampling_points)
                        return (sampling_points = sampling_points, unique_kmer_counts = unique_kmer_counts, eof = false)
                    end
                end
                ProgressMeter.next!(p)
            end
        end
    end
    sampling_points = vcat(filter(s -> s < kmers_assessed, sampling_points), [kmers_assessed])
    unique_kmer_counts = vcat(unique_kmer_counts[1:length(sampling_points)-1], [length(canonical_kmers)])
    return (sampling_points = sampling_points, unique_kmer_counts = unique_kmer_counts, eof = true)
end

function assess_kmer_saturation(fastxs; outdir="", min_k=3, max_k=61)

    if isempty(outdir)
        outdir = joinpath(pwd(), "kmer-saturation")
    end
    mkpath(outdir)

    ks = Primes.primes(min_k, max_k)
    minimum_saturation = Inf
    midpoint = Inf
    for k in ks
        kmer_type = BioSequences.BigDNAMer{k}
        kmers_to_assess = 10_000_000
        sampling_points, kmer_counts, hit_eof = assess_kmer_saturation(fastxs, kmer_type, kmers_to_assess=kmers_to_assess)
        @show sampling_points, kmer_counts, hit_eof
        observed_midpoint_index = findfirst(i -> kmer_counts[i] > last(kmer_counts)/2, 1:length(sampling_points))
        observed_midpoint = sampling_points[observed_midpoint_index]
        initial_parameters = Float64[maximum(kmer_counts), observed_midpoint]
        @time fit = LsqFit.curve_fit(calculate_v, sampling_points, kmer_counts, initial_parameters)
        if hit_eof
            inferred_maximum = last(kmer_counts)
        else
            inferred_maximum = max(Int(ceil(fit.param[1])), last(kmer_counts))
        end

        inferred_midpoint = Int(ceil(fit.param[2]))
        predicted_saturation = inferred_maximum / max_canonical_kmers(kmer_type)
        @show k, predicted_saturation

        p = StatsPlots.scatter(
            sampling_points,
            kmer_counts,
            label="observed kmer counts",
            ylabel="# unique kmers",
            xlabel="# kmers assessed",
            title = "sequencing saturation @ k = $k",
            legend=:outertopright,
            size=(800, 400),
            margins=3Plots.PlotMeasures.mm
            )
        StatsPlots.hline!(p, [max_canonical_kmers(kmer_type)], label="absolute maximum")
        StatsPlots.hline!(p, [inferred_maximum], label="inferred maximum")
        StatsPlots.vline!(p, [inferred_midpoint], label="inferred midpoint")
        # xs = vcat(sampling_points, [last(sampling_points) * 2^i for i in 1:2])
        xs = sort([sampling_points..., inferred_midpoint])
        ys = calculate_v(xs, fit.param)
        StatsPlots.plot!(
            p,
            xs,
            ys,
            label="fit trendline")
        display(p)
        StatsPlots.savefig(p, joinpath(outdir, "$k.png"))
        StatsPlots.savefig(p, joinpath(outdir, "$k.svg"))

        if predicted_saturation < minimum_saturation
            minimum_saturation = predicted_saturation
            min_k = k
            midpoint = inferred_midpoint
        end
        if predicted_saturation < 0.1
            chosen_k_file = joinpath(outdir, "chosen_k.txt")
            println("chosen k = $k")
            open(chosen_k_file, "w") do io
                println(io, k)
            end
            return k
        end
    end
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function q_value_to_error_rate(q_value)
    error_rate = 10^(q_value/(-10))
    return error_rate
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function error_rate_to_q_value(error_rate)
    q_value = -10 * log10(error_rate)
    return q_value
end

# """
# $(DocStringExtensions.TYPEDSIGNATURES)

# A short description of the function

# ```jldoctest
# julia> 1 + 1
# 2
# ```
# """
# function path_to_sequence(kmer_graph, path)
#     sequence = BioSequences.LongDNASeq(oriented_kmer_to_sequence(kmer_graph, first(path)))
#     for oriented_kmer in path[2:end]
#         nucleotide = last(oriented_kmer_to_sequence(kmer_graph, oriented_kmer))
#         push!(sequence, nucleotide)
#     end
#     return sequence
# end

# """
# $(DocStringExtensions.TYPEDSIGNATURES)

# A short description of the function

# ```jldoctest
# julia> 1 + 1
# 2
# ```
# """
# function oriented_kmer_to_sequence(kmer_graph, oriented_kmer)
#     kmer_sequence = kmer_graph.kmers[oriented_kmer.index]
#     if !oriented_kmer.orientation
#         kmer_sequence = BioSequences.reverse_complement(kmer_sequence)
#     end
#     return kmer_sequence
# end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function is_equivalent(a, b)
    a == b || a == BioSequences.reverse_complement(b)
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function edge_path_to_sequence(kmer_graph, edge_path)
    edge = first(edge_path)
    sequence = BioSequences.LongDNASeq(kmers[edge.src])
    if !kmer_graph.eprops[edge][:orientations].source_orientation
        sequence = BioSequences.reverse_complement(sequence)
    end
    for edge in edge_path
        destination = BioSequences.LongDNASeq(kmers[edge.dst])
        if !kmer_graph.eprops[edge][:orientations].destination_orientation
            destination = BioSequences.reverse_complement(destination)
        end
        sequence_suffix = sequence[end-length(destination)+2:end]
        destination_prefix = destination[1:end-1]
        @assert sequence_suffix == destination_prefix
        push!(sequence, destination[end])
    end
    sequence
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

This turns a 4-line FASTQ entry into a single tab separated line,
adds a column with the length of each read, passes it to Unix sort,
removes the length column, and converts it back into a FASTQ file.

sorts longest to shortest!!

http://thegenomefactory.blogspot.com/2012/11/sorting-fastq-files-by-sequence-length.html
```jldoctest
julia> 1 + 1
2
```
"""
function sort_fastq(input_fastq, output_fastq="")

    if endswith(input_fastq, ".gz")
        p = pipeline(
                `gzip -dc $input_fastq`,
                `paste - - - -`,
                `perl -ne '@x=split m/\t/; unshift @x, length($x[1]); print join "\t",@x;'`,
                `sort -nr`,
                `cut -f2-`,
                `tr "\t" "\n"`,
                `gzip`
                )
    else
        p = pipeline(
                `cat $input_fastq`,
                `paste - - - -`,
                `perl -ne '@x=split m/\t/; unshift @x, length($x[1]); print join "\t",@x;'`,
                `sort -nr`,
                `cut -f2-`,
                `tr "\t" "\n"`
                )
    end
    run(pipeline(p, output_fastq))
    return output_fastq
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function count_records(fastx)
    n_records = 0
    for record in Mycelia.open_fastx(fastx)
        n_records += 1
    end
    return n_records
end

function determine_read_lengths(fastq_file; total_reads = Inf)
    if total_reads == Inf
        total_reads = count_reads(fastq_file)
    end
    read_lengths = zeros(Int, total_reads)
    @info "determining read lengths"
    p = ProgressMeter.Progress(total_reads, 1)
    for (i, record) in enumerate(Mycelia.open_fastx(fastq_file))
#         push!(read_lengths, length(FASTX.sequence(record)))
        read_lengths[i] = length(FASTX.sequence(record))
        ProgressMeter.next!(p)
    end
    return read_lengths
end

function max_canonical_kmers(kmer_type)
    k_size = last(kmer_type.parameters)
    # we only consider canonical kmers so cut in 1/2
    max_canonical_kmers = Int(4^k_size / 2)
    return max_canonical_kmers
end

# Michaelis–Menten
function calculate_v(s,p)
    vmax = p[1]
    km = p[2]
    v = (vmax .* s) ./ (km .+ s)
    return v
end

max_stars_repo_name: snystrom/Mycelia

max_issues_repo_name: snystrom/Mycelia

max_forks_repo_name: snystrom/Mycelia


text:
abstract type InformationMCTSPlanner{P} <: AbstractMCTSPlanner{P} end

struct InformationBelief{B <: AbstractParticleBelief}
    b::B
    i::Float64
end

Base.rand(rng::AbstractRNG, ib::InformationBelief) = rand(rng, ib.b)
Statistics.mean(ib::InformationBelief) = mean(ib.b)

ParticleFilters.pdf(ib::InformationBelief, s::S) where {S} = pdf(ib.b,s)

function ParticleFilters.resample(resampler, d::Sampleable, rng::AbstractRNG)
    n = n_init_samples(resampler)
    ParticleCollection([rand(rng, d) for i=1:n])
end


text:
using Base.Test
using QuantumOptics.sparsematrix

# SparseMatrix = quantumoptics.sparsematrix.SparseMatrix
const SparseMatrix = SparseMatrixCSC{Complex128, Int}

@testset "sparsematrix" begin

# Set up test matrices
A = rand(Complex128, 5, 5)
A_sp = sparse(A)

B = eye(Complex128, 5)
B_sp = speye(Complex128, 5)

C = rand(Complex128, 3, 3)
C[2,:] = 0
C_sp = sparse(C)

R_sp = A_sp + B_sp
R = A + B

# Test arithmetic
@test 0 ≈ norm(full(R_sp) - R)
@test 0 ≈ norm(full(Complex128(0.5,0)*A_sp) - 0.5*A)
@test 0 ≈ norm(full(A_sp/2) - A/2)
@test 0 ≈ norm(full(A_sp*B_sp) - A*B)

# Test kronecker product
@test 0 ≈ norm(full(kron(A_sp, C_sp)) - kron(A, C))
@test 0 ≈ norm(full(kron(A_sp, B_sp)) - kron(A, B))

end # testset


text:
# # includes -
# include("Include.jl")
#
# # setup some global consts -
# const path_to_cobra_mat_file = "$(pwd())/cobra/config/matlab_cobra_files/modelReg.mat"
# const model_file_name = "modelReg"
# const organism_id = :eco    # we use the KEGG organism symbols
# const path_to_measurements_file = "$(pwd())/experimental_data/test_data/Glucose.json"
#Added new function for regulated FBA

"""
TODO: Fill me in with some stuff ...
"""
function maximize_flux_at_index(flux_index_array::Array{Int64,1}, path_to_measurements_file::String; number_of_samples::Int64 = 100)

    # initalize -
    results_array = Array{VLOptimalFluxResult,1}()

    # Declare a progress meter for user feedback -
    p = Progress(number_of_samples,color=:yellow)

    # ok, so lets sample ...
    for sample_index = 1:number_of_samples

        # load the default data_dictionary -
        default_data_dictionary = generate_default_data_dictionary();

        # pass the default dictionary to a customization method -
        updated_data_dictionary = optimize_flux_at_index(flux_index_array, default_data_dictionary);

        # update dictionary with experimental data?
        updated_data_dictionary_fluxes = constrain_measured_fluxes(updated_data_dictionary, path_to_measurements_file);
        updated_data_dictionary_metabolites = constrain_measured_metabolites(updated_data_dictionary_fluxes, path_to_measurements_file);

        # estimate the optimal flux distrubution -
        (objective_value, calculated_flux_array, dual_value_array, uptake_array, exit_flag, status_flag) = calculate_optimal_flux_distribution(updated_data_dictionary_metabolites);

        # build a return type -
        fluxResult = VLOptimalFluxResult();
        fluxResult.objective_value = objective_value;
        fluxResult.flux_array = calculated_flux_array;
        fluxResult.dual_array = dual_value_array;
        fluxResult.uptake_array = uptake_array;
        fluxResult.exit_flag = exit_flag;
        fluxResult.status_flag = status_flag;

        # get some problem setup information -
        fluxResult.flux_bounds_array = updated_data_dictionary["flux_bounds_array"];

        # cache -
        push!(results_array, fluxResult);

        # user message -
        msg = "Completed $(sample_index) of $(number_of_samples) trials ...";

        # update the progress bar -
        ProgressMeter.next!(p; showvalues = [(:status,msg)]);
    end

    @info "Completed ...\r";

    # how many flux are there?
    number_of_fluxes = length(results_array[1].flux_array);

    # compute the flux values -
    flux_ensemble = zeros(number_of_fluxes,1);
    for flux_object in results_array

        # grab the flux -
        flux_array = flux_object.flux_array;

        # cache -
        flux_ensemble = [flux_ensemble flux_array];
    end

    # cut off the zeros -
    flux_ensemble = flux_ensemble[:,2:end];

    # compute the mean and std -
    µ = mean(flux_ensemble, dims=2);
    σ = std(flux_ensemble, dims=2);
    flux_distribution = [µ σ];

    # return -
    return (flux_distribution, results_array)
end

"""
TODO: Fill me in with some stuff ...
"""
function maximize_specific_growth_rate(path_to_measurements_file::String; number_of_samples::Int64 = 100)

    # initalize -
    results_array = Array{VLOptimalFluxResult,1}()

    # Declare a progress meter for user feedback -
    p = Progress(number_of_samples,color=:yellow)

    # ok, so lets sample ...
    for sample_index = 1:number_of_samples

        # load the default data_dictionary -
        default_data_dictionary = generate_default_data_dictionary();

        # pass the default dictionary to a customization method -
        updated_data_dictionary = optimize_specific_growth_rate(default_data_dictionary);

        # update dictionary with experimental data?
        updated_data_dictionary = constrain_measured_fluxes(updated_data_dictionary, path_to_measurements_file);

        # estimate the optimal flux distrubution -
        (objective_value, calculated_flux_array, dual_value_array, uptake_array, exit_flag, status_flag) = calculate_optimal_flux_distribution(updated_data_dictionary);

        # build a return type -
        fluxResult = VLOptimalFluxResult();
        fluxResult.objective_value = objective_value;
        fluxResult.flux_array = calculated_flux_array;
        fluxResult.dual_array = dual_value_array;
        fluxResult.uptake_array = uptake_array;
        fluxResult.exit_flag = exit_flag;
        fluxResult.status_flag = status_flag;

        # get some problem setup information -
        fluxResult.flux_bounds_array = updated_data_dictionary["flux_bounds_array"];

        # cache -
        push!(results_array, fluxResult);

        # user message -
        msg = "Completed $(sample_index) of $(number_of_samples) trials ...";

        # update the progress bar -
        ProgressMeter.next!(p; showvalues = [(:status,msg)]);
    end

    @info "Completed ...\r";

    # how many flux are there?
    number_of_fluxes = length(results_array[1].flux_array);

    # compute the flux values -
    flux_ensemble = zeros(number_of_fluxes,1);
    for flux_object in results_array

        # grab the flux -
        flux_array = flux_object.flux_array;

        # cache -
        flux_ensemble = [flux_ensemble flux_array];
    end

    # cut off the zeros -
    flux_ensemble = flux_ensemble[:,2:end];

    # compute the mean and std -
    µ = mean(flux_ensemble, dims=2);
    σ = (1/sqrt(number_of_samples))*std(flux_ensemble, dims=2);
    flux_distribution = [µ σ];

    # return -
    return (flux_distribution, results_array)
end

function sample_flux_space_with_experimental_constraints(solution_bounds_array::Array{Float64,2}, number_of_samples::Int64)

    # load the default data_dictionary -
    default_data_dictionary = generate_default_data_dictionary(path_to_cobra_mat_file, model_file_name, organism_id);

    # update dictionary with experimental data?
    updated_data_dictionary = constrain_measured_fluxes(default_data_dictionary, path_to_measurements_file)

    # call the sample method -
    return sample_flux_space(solution_bounds_array, updated_data_dictionary,number_of_samples)
end

function maximize_specific_growth_rate_regulated(path_to_measurements_file::String, initialReg)

    # initalize -
    results_array = Array{VLOptimalFluxResult,1}()

    # ok, so lets sample ...

        # load the default data_dictionary -
        default_data_dictionary = generate_default_data_dictionary();

        # pass the default dictionary to a customization method -
        updated_data_dictionary = optimize_specific_growth_rate(default_data_dictionary);

        # update dictionary with experimental data?
        updated_data_dictionary = constrain_measured_fluxes(updated_data_dictionary, path_to_measurements_file);

        # estimate the optimal flux distrubution -
        #(objective_value, calculated_flux_array, dual_value_array, uptake_array, exit_flag, status_flag) = calculate_optimal_flux_distribution(updated_data_dictionary);
        (regulated_data_dictionary, calculated_flux_array, dual_value_array, uptake_array, DRgenes, constrainedRxns, states, exit_flag, status_flag) = optimizeRegModel(updated_data_dictionary, initialReg)

        # build a return type -
        #fluxResult = VLOptimalFluxResult();
        fluxResult.objective_value = objective_value;
        fluxResult.flux_array = calculated_flux_array;
        fluxResult.dual_array = dual_value_array;
        fluxResult.uptake_array = uptake_array;
        fluxResult.exit_flag = exit_flag;
        fluxResult.status_flag = status_flag;
        fluxResults.DRgenes = DRgenes;
        fluxResults.constrainedRxns = constrainedRxns;
        fluxResults.states = states;

        # get some problem setup information -
        fluxResult.flux_bounds_array = regulated_data_dictionary["flux_bounds_array"];

        # cache -
        push!(results_array, fluxResult);

    # how many flux are there?
    #number_of_fluxes = length(results_array[1].flux_array);

    # compute the flux values -
    #flux_ensemble = zeros(number_of_fluxes,1);
    #for flux_object in results_array

        # grab the flux -
        #flux_array = flux_object.flux_array;

        # cache -
        #flux_ensemble = [flux_ensemble flux_array];
    #end

    # cut off the zeros -
    #flux_ensemble = flux_ensemble[:,2:end];

    # compute the mean and std -
    #µ = mean(flux_ensemble, dims=2);
    #σ = (1/sqrt(number_of_samples))*std(flux_ensemble, dims=2);
    #flux_distribution = [µ σ];

    # return -
    return (results_array)
end

#The 2 functions below are shortcuts to performing regulated and non-regulated FBA
function maximize_specific_growth_rate_regulated(path_to_measurements_file::String;) #using optimizeRegModel

    # initalize -
    results_array = Array{VLOptimalFluxResult,1}()

    # ok, so lets sample ...

        # load the default data_dictionary -
        default_data_dictionary = generate_default_data_dictionary();

        # pass the default dictionary to a customization method -
        updated_data_dictionary = optimize_specific_growth_rate(default_data_dictionary);

        # update dictionary with experimental data?
        updated_data_dictionary = constrain_measured_fluxes(updated_data_dictionary, path_to_measurements_file);
        #updated_data_dictionary["flux_bounds_array"][83,:] .= 0.0

        # estimate the optimal flux distrubution -
        #(objective_value, calculated_flux_array, dual_value_array, uptake_array, exit_flag, status_flag) = calculate_optimal_flux_distribution(updated_data_dictionary);
        results_array = optimizeRegModel(updated_data_dictionary)

        # build a return type -
        #fluxResult = VLOptimalFluxResult();
        #fluxResult.objective_value = objective_value;
        #fluxResult.flux_array = calculated_flux_array;
        #fluxResult.dual_array = dual_value_array;
        #fluxResult.uptake_array = uptake_array;
        #fluxResult.exit_flag = exit_flag;
        #fluxResult.status_flag = status_flag;
        #fluxResults.DRgenes = DRgenes;
        #fluxResults.constrainedRxns = constrainedRxns;
        #fluxResults.states = states;

        # get some problem setup information -
        #fluxResult.flux_bounds_array = regulated_data_dictionary["flux_bounds_array"];

        # cache -
        #push!(results_array, fluxResult);

    # how many flux are there?
    #number_of_fluxes = length(results_array[1].flux_array);

    # compute the flux values -
    #flux_ensemble = zeros(number_of_fluxes,1);
    #for flux_object in results_array

        # grab the flux -
        #flux_array = flux_object.flux_array;

        # cache -
        #flux_ensemble = [flux_ensemble flux_array];
    #end

    # cut off the zeros -
    #flux_ensemble = flux_ensemble[:,2:end];

    # compute the mean and std -
    #µ = mean(flux_ensemble, dims=2);
    #σ = (1/sqrt(number_of_samples))*std(flux_ensemble, dims=2);
    #flux_distribution = [µ σ];

    # return -
    return results_array
end

function maximize_specific_growth_rate_non_regulated(path_to_measurements_file::String;) #using calculate_optimal_flux_distribution

    # initalize -
    results_array = Array{VLOptimalFluxResult,1}()

    # ok, so lets sample ...

        # load the default data_dictionary -
        default_data_dictionary = generate_default_data_dictionary();

        # pass the default dictionary to a customization method -
        updated_data_dictionary = optimize_specific_growth_rate(default_data_dictionary);

        # update dictionary with experimental data?
        updated_data_dictionary = constrain_measured_fluxes(updated_data_dictionary, path_to_measurements_file);
        #updated_data_dictionary["flux_bounds_array"][83,:] .= 0.0

        # estimate the optimal flux distrubution -
        results_array = calculate_optimal_flux_distribution(updated_data_dictionary);

        flux_bounds_array = updated_data_dictionary["flux_bounds_array"]

        # build a return type -
        #fluxResult = VLOptimalFluxResult();
        #fluxResult.objective_value = objective_value;
        #fluxResult.flux_array = calculated_flux_array;
        #fluxResult.dual_array = dual_value_array;
        #fluxResult.uptake_array = uptake_array;
        #fluxResult.exit_flag = exit_flag;
        #fluxResult.status_flag = status_flag;
        #fluxResults.DRgenes = DRgenes;
        #fluxResults.constrainedRxns = constrainedRxns;
        #fluxResults.states = states;

        # get some problem setup information -
        #fluxResult.flux_bounds_array = regulated_data_dictionary["flux_bounds_array"];

        # cache -
        #push!(results_array, fluxResult);

    # how many flux are there?
    #number_of_fluxes = length(results_array[1].flux_array);

    # compute the flux values -
    #flux_ensemble = zeros(number_of_fluxes,1);
    #for flux_object in results_array

        # grab the flux -
        #flux_array = flux_object.flux_array;

        # cache -
        #flux_ensemble = [flux_ensemble flux_array];
    #end

    # cut off the zeros -
    #flux_ensemble = flux_ensemble[:,2:end];

    # compute the mean and std -
    #µ = mean(flux_ensemble, dims=2);
    #σ = (1/sqrt(number_of_samples))*std(flux_ensemble, dims=2);
    #flux_distribution = [µ σ];

    # return -
    return (results_array, flux_bounds_array)
end


text:
export Atom, vectorize

struct Atom{T} <: AbstractElement{T}
    e::Symbol; x::T; y::T; z::T
    Atom(e::Symbol, p::NTuple{3, T}) where T<:Real = new{T}(e, p[1], p[2], p[3])
end

# @code_warntype ✓
Base.getindex(tup::NTuple{N, Atom}, idx::Int, dir::Symbol) where N = getproperty(tup[idx], dir)
Base.getproperty(a::Atom, s::Symbol) = s ≡ :x ? getfield(a, :x) : s ≡ :y ? getfield(a, :y) : s ≡ :z ? getfield(a, :z) : error("Invalid :$s.")

# @code_warntype
vectorize(a::Atom) = (a.x,), (a.y,), (a.z,)
vectorize(a::Atom{T}, b::Atom{T}) where T<:Real = (a.x, b.x), (a.y, b.y), (a.z, b.z)

function vectorize(tup::NTuple{N, Atom}) where N
    if @generated
        tupX = Expr(:tuple)
        tupY = Expr(:tuple)
        tupZ = Expr(:tuple)
        argX = Vector{Expr}(undef, N)
        argY = Vector{Expr}(undef, N)
        argZ = Vector{Expr}(undef, N)

        @inbounds for i in 1:N
            argX[i] = :(tup[$i, :x])
            argY[i] = :(tup[$i, :y])
            argZ[i] = :(tup[$i, :z])
        end

        tupX.args = argX
        tupY.args = argY
        tupZ.args = argZ

        return quote
            $(Expr(:meta, :inline))
            @inbounds return ($tupX, $tupY, $tupZ)
        end
    else
        tupX = ntuple(i -> tup[i, :x], N)
        tupY = ntuple(j -> tup[j, :y], N)
        tupZ = ntuple(k -> tup[k, :z], N)
        return tupX, tupY, tupZ
    end
end


text:
"""
Positive definite diagonal matrix.
"""
struct PDiagMat{T<:Real,V<:AbstractVector} <: AbstractPDMat{T}
    dim::Int
    diag::V
    inv_diag::V

    PDiagMat{T,S}(d::Int,v::AbstractVector,inv_v::AbstractVector) where {T,S} =
        new{T,S}(d,v,inv_v)
end

function PDiagMat(v::AbstractVector,inv_v::AbstractVector)
    @check_argdims length(v) == length(inv_v)
    PDiagMat{eltype(v),typeof(v)}(length(v), v, inv_v)
end

PDiagMat(v::AbstractVector) = PDiagMat(v, inv.(v))

### Conversion
Base.convert(::Type{PDiagMat{T}},      a::PDiagMat) where {T<:Real} = PDiagMat(convert(AbstractArray{T}, a.diag))
Base.convert(::Type{AbstractArray{T}}, a::PDiagMat) where {T<:Real} = convert(PDiagMat{T}, a)

### Basics

dim(a::PDiagMat) = a.dim
Base.Matrix(a::PDiagMat) = Matrix(Diagonal(a.diag))
LinearAlgebra.diag(a::PDiagMat) = copy(a.diag)
LinearAlgebra.cholesky(a::PDiagMat) = cholesky(Diagonal(a.diag))

### Inheriting from AbstractMatrix

function Base.getindex(a::PDiagMat, i::Integer)
    ncol, nrow = fldmod1(i, a.dim)
    ncol == nrow ? a.diag[nrow] : zero(eltype(a))
end
Base.getindex(a::PDiagMat{T},i::Integer,j::Integer) where {T} = i == j ? a.diag[i] : zero(T)

### Arithmetics

function pdadd!(r::Matrix, a::Matrix, b::PDiagMat, c)
    @check_argdims size(r) == size(a) == size(b)
    if r === a
        _adddiag!(r, b.diag, c)
    else
        _adddiag!(copyto!(r, a), b.diag, c)
    end
    return r
end

*(a::PDiagMat, c::T) where {T<:Real} = PDiagMat(a.diag * c)
*(a::PDiagMat, x::AbstractVector) = a.diag .* x
*(a::PDiagMat, x::AbstractMatrix) = a.diag .* x
\(a::PDiagMat, x::AbstractVecOrMat) = a.inv_diag .* x
/(x::AbstractVector, a::PDiagMat) = a.inv_diag .* x
/(x::AbstractMatrix, a::PDiagMat) = transpose(a.inv_diag) .* x
Base.kron(A::PDiagMat, B::PDiagMat) = PDiagMat( vcat([A.diag[i] * B.diag for i in 1:dim(A)]...) )

### Algebra

Base.inv(a::PDiagMat) = PDiagMat(a.inv_diag, a.diag)
LinearAlgebra.logdet(a::PDiagMat) = sum(log, a.diag)
LinearAlgebra.eigmax(a::PDiagMat) = maximum(a.diag)
LinearAlgebra.eigmin(a::PDiagMat) = minimum(a.diag)

### whiten and unwhiten

function whiten!(r::StridedVector, a::PDiagMat, x::StridedVector)
    n = dim(a)
    @check_argdims length(r) == length(x) == n
    v = a.inv_diag
    for i = 1:n
        r[i] = x[i] * sqrt(v[i])
    end
    return r
end

function unwhiten!(r::StridedVector, a::PDiagMat, x::StridedVector)
    n = dim(a)
    @check_argdims length(r) == length(x) == n
    v = a.diag
    for i = 1:n
        r[i] = x[i] * sqrt(v[i])
    end
    return r
end

whiten!(r::StridedMatrix, a::PDiagMat, x::StridedMatrix) =
    broadcast!(*, r, x, sqrt.(a.inv_diag))

unwhiten!(r::StridedMatrix, a::PDiagMat, x::StridedMatrix) =
    broadcast!(*, r, x, sqrt.(a.diag))

### quadratic forms

quad(a::PDiagMat, x::AbstractVector) = wsumsq(a.diag, x)
invquad(a::PDiagMat, x::AbstractVector) = wsumsq(a.inv_diag, x)

function quad!(r::AbstractArray, a::PDiagMat, x::StridedMatrix)
    m, n = size(x)
    ad = a.diag
    @check_argdims m == length(ad) && length(r) == n
    @inbounds for j = 1:n
        s = zero(promote_type(eltype(ad), eltype(x)))
        for i in 1:m
            s += ad[i] * abs2(x[i,j])
        end
        r[j] = s
    end
    r
end

function invquad!(r::AbstractArray, a::PDiagMat, x::StridedMatrix)
    m, n = size(x)
    ainvd = a.inv_diag
    @check_argdims m == length(ainvd) && length(r) == n
    @inbounds for j = 1:n
        s = zero(promote_type(eltype(ainvd), eltype(x)))
        for i in 1:m
            s += ainvd[i] * abs2(x[i,j])
        end
        r[j] = s
    end
    r
end

### tri products

function X_A_Xt(a::PDiagMat, x::StridedMatrix)
    z = x .* reshape(sqrt.(a.diag), 1, dim(a))
    z * transpose(z)
end

function Xt_A_X(a::PDiagMat, x::StridedMatrix)
    z = x .* sqrt.(a.diag)
    transpose(z) * z
end

function X_invA_Xt(a::PDiagMat, x::StridedMatrix)
    z = x .* reshape(sqrt.(a.inv_diag), 1, dim(a))
    z * transpose(z)
end

function Xt_invA_X(a::PDiagMat, x::StridedMatrix)
    z = x .* sqrt.(a.inv_diag)
    transpose(z) * z
end


text:
# for GMM
function SVmoments(m, n, θ, η, ϵ)
    S = size(ϵ, 2)
    ms = zeros(S,size(m,1))
    Threads.@threads for s=1:S
        ms[s,:] = sqrt(n)*aux_stat(SVmodel(θ, n, η[:,s], ϵ[:,s])[1])
    end
    ms .- m'
end

max_stars_repo_name: mcreel/SV

max_issues_repo_name: mcreel/SV

max_forks_repo_name: mcreel/SV


text:
using Documenter, FrankWolfe
using SparseArrays
using LinearAlgebra

using Literate, Test

EXAMPLE_DIR = joinpath(dirname(@__DIR__), "examples")
DOCS_EXAMPLE_DIR = joinpath(@__DIR__, "src", "examples")
DOCS_REFERENCE_DIR = joinpath(@__DIR__, "src", "reference")

function file_list(dir, extension)
    return filter(file -> endswith(file, extension), sort(readdir(dir)))
end

function literate_directory(jl_dir, md_dir)
    for filename in file_list(md_dir, ".md")
        filepath = joinpath(md_dir, filename)
        rm(filepath)
    end
    for filename in file_list(jl_dir, ".jl")
        filepath = joinpath(jl_dir, filename)
        # `include` the file to test it before `#src` lines are removed. It is
        # in a testset to isolate local variables between files.
        if startswith(filename, "docs")
            Literate.markdown(
                filepath, md_dir; documenter=true, flavor=Literate.DocumenterFlavor()
            )
        end
    end
    return nothing
end

literate_directory(EXAMPLE_DIR, DOCS_EXAMPLE_DIR)

ENV["GKSwstype"] = "100"

generated_path = joinpath(@__DIR__, "src")
base_url = "https://github.com/ZIB-IOL/FrankWolfe.jl/blob/master/"
isdir(generated_path) || mkdir(generated_path)

open(joinpath(generated_path, "contributing.md"), "w") do io
    # Point to source license file
    println(
        io,
        """
        ```@meta
        EditURL = "$(base_url)CONTRIBUTING.md"
        ```
        """,
    )
    # Write the contents out below the meta block
    for line in eachline(joinpath(dirname(@__DIR__), "CONTRIBUTING.md"))
        println(io, line)
    end
end

open(joinpath(generated_path, "index.md"), "w") do io
    # Point to source license file
    println(
        io,
        """
        ```@meta
        EditURL = "$(base_url)README.md"
        ```
        """,
    )
    # Write the contents out below the meta block
    for line in eachline(joinpath(dirname(@__DIR__), "README.md"))
        println(io, line)
    end
end

makedocs(;
    modules=[FrankWolfe],
    sitename="FrankWolfe.jl",
    format=Documenter.HTML(; prettyurls=get(ENV, "CI", nothing) == "true", collapselevel=1),
    pages=[
        "Home" => "index.md",
        "How does it work?" => "basics.md",
        "Advanced features" => "advanced.md",
        "Examples" => [joinpath("examples", f) for f in file_list(DOCS_EXAMPLE_DIR, ".md")],
        "API reference" =>
            [joinpath("reference", f) for f in file_list(DOCS_REFERENCE_DIR, ".md")],
        "Contributing" => "contributing.md",
    ],
)

deploydocs(; repo="github.com/ZIB-IOL/FrankWolfe.jl.git", push_preview=true)


text:
module MatplotWrap

using Libdl
using Pkg.Artifacts

using CxxWrap

#const libmplxx = joinpath(@__DIR__, "..", "deps", "build", "lib", "libmplxx.$(dlext)")
const libmplxx = joinpath(artifact"libmplxx","libmplxx.$(Libdl.dlext)")
isfile(libmplxx) && @wrapmodule(libmplxx, :define_julia_module)

export plot, plot3, scatter, matplot_show

function plot(_x::AbstractVector{T1}, _y::AbstractVector{T2}, line_spec::AbstractString="") where {T1,T2 <: Real}
    x = convert.(Cdouble, _x) |> StdVector
    y = convert.(Cdouble, _y) |> StdVector
    plot(x, y, StdString(line_spec))
end

function plot(_y::AbstractVector{T}, line_spec::AbstractString="") where {T <: Real}
    _x = collect(1.0:length(_y))
    plot(_x, _y, line_spec)
end

function plot3(_x::AbstractVector{T1}, _y::AbstractVector{T2}, _z::AbstractVector{T3}) where {T1,T2,T3 <: Real}
    x = convert.(Cdouble, _x) |> StdVector
    y = convert.(Cdouble, _y) |> StdVector
    z = convert.(Cdouble, _z) |> StdVector
    plot3(x, y, z)
end

function scatter(_x::AbstractVector{T1}, _y::AbstractVector{T2}) where {T1,T2 <: Real}
    x = convert.(Cdouble, _x) |> StdVector
    y = convert.(Cdouble, _y) |> StdVector
    scatter(x, y)
end

    function scatter(_y::AbstractVector{T}) where {T <: Real}
    _x = collect(1.0:length(_y))
    scatter(_x, _y)
end

function __init__()
    @initcxx
end

end # module


text:
module TestGraphMLWiringDiagrams

using Test
using LightXML
using Catlab.WiringDiagrams

# Round trip wiring diagrams with dictionary box and wire data.

function roundtrip(f::WiringDiagram)
  xdoc = generate_graphml(f)
  parse_graphml(Dict, Nothing, Dict, xdoc)
end

ports(n::Int) = repeat([nothing], n)
diagram = WiringDiagram(ports(1), ports(1))
f = Box(Dict("name" => "f", "type" => "foo"), ports(1), ports(1))
fv = add_box!(diagram, f)
add_wires!(diagram, [
  Wire(Dict("name" => "w1", "type" => "woo"),
      (input_id(diagram), 1) => (fv, 1)),
  Wire(Dict("name" => "w2", "type" => "woo"),
       (fv, 1) => (output_id(diagram), 1)),
])
@test roundtrip(diagram) == diagram

# Round trip wiring diagrams with symbolic box and port values.

function roundtrip_symbolic(f::WiringDiagram)
  xdoc = generate_graphml(f)
  parse_graphml(Symbol, Symbol, Nothing, xdoc)
end

f = singleton_diagram(Box(:f, [:A], [:B]))
g = singleton_diagram(Box(:g, [:B], [:C]))

@test roundtrip_symbolic(f) == f
@test roundtrip_symbolic(compose(f,g)) == compose(f,g)
@test roundtrip_symbolic(otimes(f,g)) == otimes(f,g)

# Round trip nested, symbolic wiring diagrams.

inner = copy(f)
inner.value = :inner
outer = singleton_diagram(f)
outer.value = :outer
@test roundtrip_symbolic(outer) == outer

end


text:
using DataFrames, LowRankModels

# boolean example with only entries greater than threshold t observed
# ie, censored data
# example with only entries greater than threshold t observed
m,n,k,ktrue = 100,100,1,1
A = rand(m,ktrue)*rand(ktrue,n)
println("max value of A is ",maximum(maximum(A))," which is less than $ktrue")
B = round.(Int, ktrue*rand(m,n) .>= A) # Bernoulli samples with probability proportional to A
losses = fill(QuadLoss(),n)
r = QuadReg(.1)
obs = Array{Tuple{Int,Int}}(undef, 0)
for i=1:m
    for j=1:n
        if B[i,j] == 1
            push!(obs, (i,j))
        end
    end
end

(train_observed_features, train_observed_examples, test_observed_features,  test_observed_examples) =
    get_train_and_test(obs, m, n, .2)

train_glrm = GLRM(B,losses,r,r,k, observed_features=train_observed_features, observed_examples=train_observed_examples)

train_error, test_error, prec_at_k, train_time, reg_params, solution =
    precision_at_k(train_glrm, test_observed_features,
    					params=Params(1, max_iter=200, abs_tol=0.00001, min_stepsize=0.01),
                                 reg_params=exp10.(range(2,step=-2,length=9)))
df = DataFrame(train_error = train_error, prec_at_k = prec_at_k,
                   train_time = train_time, reg_param = reg_params, solution_1norm = [s[2] for s in solution]);
println(df)


text:
# HTML tables

export table

using DataFrames

# Model

immutable Table <: Tile
    head::AbstractVector
    body::AbstractDataFrame
end

table(body; head=Array(Tile, 0)) =
    Table(head, body)

# Behaviors

immutable SelectRows
    multi::Bool
    table::Table
end

immutable SelectCols
    multi::Bool
    table::Table
end

# Render logic

render_cell(x, state) =
    Elem(:td, render(x, state))

render(t::Table, state) = begin
    if length(t.head) > 0
        @assert length(t.head) == ncol(t.body)
        head = t.head
    else
        head = names(t.body)
    end

    Elem(:table,
        [Elem(:thead, Elem(:tr, map(x -> Elem(:th, render(x, state)), head))),
         Elem(:tbody,
            [Elem(:tr, [render_cell(v, state) for (k,v) in row])
                for row in eachrow(t.body)]
         )]
    )
end


text:
"""
Handles HttpServer related functionality, manages requests and responses and their logging.
"""
module AppServer

using HttpServer, Router, Genie, Millboard, Logger, Sessions, Genie.Configuration, MbedTLS, WebSockets, Channels, App, URIParser

struct WebServer
  task::Task
  server::HttpServer.Server
end

"""
    startup(port::Int = 8000) :: HttpServer.Server

Starts the web server on the configurated port.
Automatically invoked when Genie is started with the `s` or the `server:start` command line params.

# Examples
```julia
julia> AppServer.startup()
Listening on 0.0.0.0:8000...
```
"""
function startup(port::Int = 8000) :: WebServer
  http = setup_http_handler()
  App.config.websocket_server && (wsh = setup_ws_handler())
  App.config.lookup_ip && (http.events["connect"] = (http_client) -> handle_connect(http_client))

  server = App.config.websocket_server ? Server(http, wsh) : Server(http)
  server_task = @async run(server, port)

  if App.config.run_as_server
    while true
      sleep(1_000_000)
    end
  end

  WebServer(server_task, server)
end

"""
"""
function setup_http_handler() :: HttpHandler
  HttpHandler() do req::Request, res::Response
    try
      ip::IPv4 = App.config.lookup_ip ? task_local_storage(:ip) : ip"255.255.255.255"
      nworkers() == 1 ? handle_request(req, res, ip) : @fetch handle_request(req, res, ip)
    catch ex
      Logger.log(string(ex), :critical)
      Logger.log(sprint(io->Base.show_backtrace(io, catch_backtrace() )), :critical)
      Logger.log("$(@__FILE__):$(@__LINE__)", :critical)

      message = Configuration.is_prod() ?
                  "The error has been logged and we'll look into it ASAP." :
                  string(ex, " in $(@__FILE__):$(@__LINE__)", "\n\n", sprint(io->Base.show_backtrace(io, catch_backtrace())))

      Router.serve_error_file(500, message, Router.@params)
    end
  end
end

"""
"""
function setup_ws_handler() :: WebSocketHandler
  WebSocketHandler() do req::Request, ws_client::WebSockets.WebSocket
    while true
      response =  try
                    msg = read(ws_client)
                    ip::IPv4 = App.config.lookup_ip ? task_local_storage(:ip) : ip"255.255.255.255"

                    nworkers() == 1 ? handle_ws_request(req, String(msg), ws_client, ip) : @fetch handle_ws_request(req, String(msg), ws_client, ip)
                  catch ex
                    if typeof(ex) == WebSockets.WebSocketClosedError
                      Channels.unsubscribe_client(ws_client)

                      break
                    end

                    try
                      Logger.log(string(ex), :critical)
                      Logger.log("$(@__FILE__):$(@__LINE__)", :critical)

                      Configuration.is_prod() ? "The error has been logged and we'll look into it ASAP." : string(ex)
                    catch exx
                      print_with_color(:red, "One can not simply log an error")
                    end

                    break
                  end

      try
        write(ws_client, response)
      catch socket_exception
        Channels.unsubscribe_client(ws_client)

        break
      end
    end
  end
end

"""
    handle_connect(client::HttpServer.Client) :: Void

Connection callback for HttpServer. Stores the Request IP in the current task's local storage.
"""
function handle_connect(client::HttpServer.Client) :: Void
  try
    ip, port = getsockname(isa(client.sock, MbedTLS.SSLContext) ? client.sock.bio : client.sock)
    task_local_storage(:ip, ip)
  catch ex
    Logger.log("Failed getting IP address of request", :err)
    Logger.log(string(ex), :err)
    Logger.log("$(@__FILE__):$(@__LINE__)", :err)

    task_local_storage(:ip, ip"255.255.255.255")
  end

  nothing
 end

"""
    handle_request(req::Request, res::Response, ip::IPv4 = ip"0.0.0.0") :: Response

HttpServer handler function - invoked when the server gets a request.
"""
function handle_request(req::Request, res::Response, ip::IPv4 = ip"0.0.0.0") :: Response
  App.config.log_requests && log_request(req)
  App.config.server_signature != "" && sign_response!(res)

  app_response::Response = set_headers!(req, res, Router.route_request(req, res, ip))

  App.config.log_responses && log_response(req, app_response)

  app_response
end

function set_headers!(req::Request, res::Response, app_response::Response) :: Response
  if req.method == OPTIONS
    App.config.cors_headers["Access-Control-Allow-Origin"] = strip(App.config.cors_headers["Access-Control-Allow-Origin"])

    ! isempty(App.config.cors_allowed_origins) &&
      in(req.headers["Origin"], App.config.cors_allowed_origins) &&
      (App.config.cors_headers["Access-Control-Allow-Origin"] == "" ||
        App.config.cors_headers["Access-Control-Allow-Origin"] == "*") &&
      (App.config.cors_headers["Access-Control-Allow-Origin"] = req.headers["Origin"])

    app_response.headers = merge(res.headers, App.config.cors_headers)
  end

  app_response.headers = merge(res.headers, app_response.headers)

  app_response.cookies = merge(res.cookies, app_response.cookies)

  app_response
end

"""
    handle_ws_request(req::Request, client::Client, ip::IPv4 = ip"0.0.0.0") :: String

HttpServer handler function - invoked when the server gets a request.
"""
function handle_ws_request(req::Request, msg::String, ws_client::WebSockets.WebSocket, ip::IPv4 = ip"0.0.0.0") :: String
  App.config.log_requests && log_request(req)

  Router.route_ws_request(req, msg, ws_client, ip)
end

"""
    sign_response!(res::Response) :: Response

Adds a signature header to the response using the value in `App.config.server_signature`.
If `App.config.server_signature` is empty, the header is not added.
"""
function sign_response!(res::Response) :: Response
  if ! isempty(App.config.server_signature)
    res.headers["Server"] = App.config.server_signature
  end

  res
end

"""
    log_request(req::Request) :: Void

Logs information about the request.
"""
function log_request(req::Request) :: Void
  if Router.is_static_file(req.resource)
    App.config.log_resources && log_request_response(req)
  elseif App.config.log_responses
    log_request_response(req)
  end

  nothing
end

"""
    log_response(req::Request, res::Response) :: Void

Logs information about the response.
"""
function log_response(req::Request, res::Response) :: Void
  if Router.is_static_file(req.resource)
    App.config.log_resources && log_request_response(res)
  elseif App.config.log_responses
    log_request_response(res)
  end

  nothing
end

"""
    log_request_response(req_res::Union{Request,Response}) :: Void

Helper function that logs `Request` or `Response` objects.
"""
function log_request_response(req_res::Union{Request,Response}) :: Void
  req_data = Dict{String,String}()
  response_is_error = false

  for f in fieldnames(req_res)
    f = string(f)
    v = getfield(req_res, Symbol(f))

    f == "status" && (req_res.status == 404 || req_res.status == 500) && (response_is_error = true)

    req_data[f] = if f == "data" && ! isempty(v)
                    mapreduce(x -> string(Char(Int(x))), *, v) |> Logger.truncate_logged_output
                  elseif isa(v, Dict) && App.config.log_formatted
                    Millboard.table(parse_inner_dict(v)) |> string
                  else
                    string(v) |> Logger.truncate_logged_output
                  end
  end

  Logger.log(string(req_res) * "\n" * string(App.config.log_formatted ? Millboard.table(req_data) : req_data), response_is_error ? :err : :debug, showst = false)

  nothing
end

"""
    parse_inner_dict{K,V}(d::Dict{K,V}) :: Dict{String,String}

Helper function that knows how to parse a `Dict` containing `Request` or `Response` data and prepare it for being logged.
"""
function parse_inner_dict{K,V}(d::Dict{K,V}) :: Dict{String,String}
  r = Dict{String,String}()
  for (k, v) in d
    k = string(k)
    if k == "Cookie" && App.config.log_verbosity == Configuration.LOG_LEVEL_VERBOSITY_VERBOSE
      cookie = Dict{String,String}()
      cookies = split(v, ";")
      for c in cookies
        cookie_part = split(c, "=")
        cookie[cookie_part[1]] = cookie_part[2] |> Logger.truncate_logged_output
      end

      r[k] = (App.config.log_formatted ? Millboard.table(cookie) : cookie) |> string
    else
      r[k] = Logger.truncate_logged_output(string(v))
    end
  end

  r
end

"""
"""
function shutdown(w::WebServer) :: Void
  close(w.server)
end

end


text:
using Test
using OSMMakie
using OSMMakie.LightOSM
using CairoMakie
using StatsBase

## Get data

# download OpenStreetMap data and load as OSMGraph
osm = graph_from_download(:bbox; # rectangular area
    minlat = 51.5015, # bottom left corner
    minlon = -0.0921,
    maxlat = 51.5154, # top right corner
    maxlon = -0.0662,
    network_type = :drive, # download motorways
    graph_type = :light, # SimpleDiGraph
    weight_type = :distance
)

# download OpenStreetMap buildings and load as buildings dict
buildings = buildings_from_download(:bbox;
    minlat = 51.5015, # bottom left corner
    minlon = -0.0921,
    maxlat = 51.5154, # top right corner
    maxlon = -0.0662,
    metadata = true,
    download_format = :osm,
);

## Tests

@testset "Basic" begin
    @test begin # OSMPlot recipe
        osmplot(osm)
        true
    end

    @test begin # DataInspector
        fig, ax, p = osmplot(osm)
        DataInspector(fig)
        true
    end

    # wrong recipe usage, for whatever reasons this is different between CI and local
    if get(ENV, "CI", nothing) == "true"
        @test_throws ArgumentError osmplot(; osm)
    else
        @test_throws MethodError osmplot(; osm)
    end

    @test begin # recipe with buildings and limits/autolimitaspect keywords
        fig = Figure()
        ax = fig[1,1] = Axis(fig;
            limits = ((-0.0921, -0.0662), (51.5015, 51.5154)),
            autolimitaspect = 1 + abs(mean([51.5015, 51.5154])) / 90
        )
        p = osmplot!(osm; buildings)
        true
    end
end


text:
# Nystrom Approximation ====================================================================

for orientation in (:row, :col)
    (dim, S, fulldim) = orientation == :row ? (1, :S, :(:)) : (2, :(:), :S)

    @eval begin
        function samplematrix(
                σ::Val{$(Meta.quot(orientation))},
                X::Matrix,
                r::T
            ) where {T<:AbstractFloat}
            0 < r <= 1 || error("Sample rate must be in range (0,1]")
            n = size(X, $dim)
            s = max(Int64(trunc(n*r)),1)
            S = [rand(1:n) for i = 1:s]

        end

        function nystrom_sample(
                σ::Val{$(Meta.quot(orientation))},
                κ::Kernel{T},
                X::Matrix{T},
                S::Vector{U}
            ) where {T<:AbstractFloat,U<:Integer}
            Xs = getindex(X, $S, $fulldim)
            C = kernelmatrix(σ, κ, Xs, X)   # kernel matrix of X and sampled X
            Cs = getindex(C, :, S)  # purely sampled component of C
            return (C, Cs)
        end

    end
end

function nystrom_pinv!(
        Cs::Matrix{T},
        tol::T = eps(T)*size(Cs,1)
    ) where {T<:LinearAlgebra.BlasReal}
    # Compute eigendecomposition of sampled component of C
    QΛQᵀ = LinearAlgebra.eigen!(LinearAlgebra.Symmetric(Cs))

    # Solve for D = Λ^(-1/2) (pseudo inverse - use tolerance from before factorization)
    D = QΛQᵀ.values
    λ_tol = maximum(D)*tol

    for i in eachindex(D)
        @inbounds D[i] = abs(D[i]) <= λ_tol ? zero(T) : one(T)/sqrt(D[i])
    end

    # Scale eigenvectors by D
    Q = QΛQᵀ.vectors
    QD = LinearAlgebra.rmul!(Q, LinearAlgebra.Diagonal(D))  # Scales column i of Q by D[i]

    # W := (QD)(QD)ᵀ = (QΛQᵀ)^(-1)  (pseudo inverse)
    W = LinearAlgebra.syrk_wrapper!(similar(QD), 'N', QD)

    return LinearAlgebra.copytri!(W, 'U')
end

"""
    NystromFact

Type for storing a Nystrom factorization. The factorization contains two fields: `W` and
`C` as described in the `nystrom` documentation.
"""
struct NystromFact{T<:LinearAlgebra.BlasReal}
    W::Matrix{T}
    C::Matrix{T}
end

@doc raw"""
    nystrom([σ::Orientation,] κ::Kernel, X::Matrix, [S::Vector])

Computes a factorization of Nystrom approximation of the square kernel matrix of data
matrix `X` with respect to kernel `κ`. Returns a `NystromFact` struct which stores a
Nystrom factorization satisfying:

```math
\mathbf{K} \approx \mathbf{C}^{\intercal}\mathbf{WC}
```
"""
function nystrom(
        σ::Orientation,
        κ::Kernel{T},
        X::Matrix{T},
        S::Vector{U} = samplematrix(σ, X, convert(T,0.15))
    ) where {T<:LinearAlgebra.BlasReal,U<:Integer}
    C, Cs = nystrom_sample(σ, κ, X, S)
    W = nystrom_pinv!(Cs)
    NystromFact{T}(W, C)
end

function nystrom(
        κ::Kernel{T},
        X::Matrix{T},
        S::Vector{U} = samplematrix(:row(), X, convert(T,0.15))
    ) where {T<:LinearAlgebra.BlasReal,U<:Integer}
    nystrom(:row(), κ, X, S)
end

function kernelmatrix(CᵀWC::NystromFact{T}) where {T<:LinearAlgebra.BlasReal}
    W = CᵀWC.W
    C = CᵀWC.C
    At_mul_B(C,W)*C
end


text:
module MonteCarloX
# dependencies
using Random
using StatsBase
using LinearAlgebra
using StaticArrays
# exports that are relevant to run simulations in MonteCarloX
# using Reexports
# @reexport using Random
# @reexport using StatsBase

include("utils.jl")
include("event_handler.jl")
include("reweighting.jl")
include("rng.jl")

# Equilibrium
include("importance_sampling.jl")

# Non-equilibrium
include("kinetic_monte_carlo.jl")
include("poisson_process.jl")
include("gillespie.jl")

# TODO: move to external SpinSystems.jl package
include("cluster_wolff.jl")

# algorithms
export  Metropolis,
        Gillespie,
        KineticMonteCarlo,
        InhomogeneousPoisson,
        InhomogeneousPoissonPiecewiseDecreasing

# functions
export  # base
        # equilibrium
        accept,
        sweep,
        # update -> will be moved to test/utils.jl for now and later to SpinSystems.jl
        # non-equilibrium
        next_event,
        next_time,
        next,
        advance!

# helper
export  log_sum,
        binary_search,
        random_element,
        # event handler
        AbstractEventHandlerRate,
        ListEventRateSimple,
        ListEventRateActiveMask,
        # rng
        MutableRandomNumbers

# Check georges 2nd workshop notebook on github
#    i::Int = StatsBase.binindex(d.h, x) -> remember this for custom things similar to EmpiricalDistributions.jl (not well documented though)

end # module

# Maybe embedd this into StatisticalPhysics.jl the including SpinSystems.jl PolymerSystems.jl etc ;)
# TODO: external modules obviously in external modules
# ("DirectedPercolation.jl") [inlcuding ContactProcess, CellularAutomatoa, etc but not as modules but as models]
# ("NeuralNetworks.jl) -> maybe NeuralSystems.jl?
# ("ComplexNetworks.jl)


text:
module PowerSimulations

#################################################################################
# Exports

# Base Models
export Stage
export Simulation
export OperationModel
export ModelReference
export InitialCondition

#Network Relevant Exports
export StandardPTDFModel
export CopperPlatePowerModel

######## Device Models ########
export DeviceModel
######## Service Models ########
export ServiceModel
######## Branch Models ########
export StaticLine
export StaticTransformer
export TapControl
export StaticLineUnbounded
#export PhaseControl
export HVDCLossless
export HVDCDispatch
#export VoltageSourceDC
######## Load Models ########
export StaticPowerLoad
export InterruptiblePowerLoad
export DispatchablePowerLoad
######## Renewable Formulations ########
export RenewableFixed
export RenewableFullDispatch
export RenewableConstantPowerFactor
######## Hydro Formulations ########
export HydroFixed
######## Renewable Formulations ########
export BookKeeping
export BookKeepingwReservation
######## Thermal Formulations ########
export ThermalStandardUnitCommitment
export ThermalBasicUnitCommitment
export ThermalDispatch
export ThermalRampLimited
export ThermalDispatchNoMin

# feedforward sequences
export RecedingHorizon
export Synchronize
export Sequential

# feedforward models
export UpperBoundFF
export SemiContinuousFF
export RangeFF

# Initial Conditions Quantities
export DevicePower
export DeviceStatus
export TimeDurationON
export TimeDurationOFF
export DeviceEnergy

# cache_models
export TimeStatusChange

#operation_models
#export UnitCommitment
#export EconomicDispatch
#export OptimalPowerFlow

# Functions
## Construction Exports
export construct_device!
export construct_network!
## Op Model Exports
export solve_op_model!
export get_initial_conditions
export set_transmission_ref!
export set_devices_ref!
export set_branches_ref!
export set_services_ref!
export set_device_model!
export set_branch_model!
export set_device_model!
## Sim Model Exports
export run_sim_model!
## Utils Exports
export write_op_model
export write_model_results
export load_operation_results
export load_simulation_results
export get_all_constraint_index
export get_all_var_index
export get_con_index
export get_var_index
# Plotting Utils
export sort_data
export get_stacked_plot_data
export get_bar_plot_data
export get_stacked_generation_data
export bar_plot
export stack_plot
export report

#################################################################################
# Imports
#Modeling Imports
import JuMP
# so that users do not need to import JuMP to use a solver with PowerModels
import JuMP: with_optimizer
export with_optimizer
import MathOptInterface
import ParameterJuMP
import LinearAlgebra
import PowerSystems
import InfrastructureSystems
import PowerModels
import RecipesBase
import Requires

#TimeStamp Management Imports
import Dates
import TimeSeries

#I/O Imports
import MathOptFormat
import DataFrames
import Feather

include("core/definitions.jl")

#################################################################################
##### JuMP methods overloading
JuMP.Model(optimizer::Nothing; kwargs...) = JuMP.Model(kwargs...)

################################################################################
# Includes

#Abstract Models
include("network_models/networks.jl")
include("service_models/services.jl")

#Core Models and constructors
include("core/core_structs/aux_structs.jl")
include("core/core_structs/cache_models.jl")
include("core/core_structs/feedforward_model.jl")
include("core/core_structs/device_model.jl")
include("core/core_structs/initial_conditions.jl")
include("core/core_structs/canonical_model.jl")
include("core/core_structs/service_model.jl")
include("core/core_structs/operation_model.jl")
include("core/core_structs/chronology.jl")
include("core/core_structs/simulations_stages.jl")
include("core/core_structs/simulation_model.jl")
include("core/core_structs/results_model.jl")
include("core/build_cache.jl")
include("core/build_operations.jl")
include("core/build_simulations.jl")

#FeedForward Model Files
include("simulation_models/feedforward_models.jl")

#Device Modeling components
include("devices/device_models/common.jl")
include("devices/device_models/renewable_generation.jl")
include("devices/device_models/thermal_generation.jl")
include("devices/device_models/electric_loads.jl")
include("devices/device_models/AC_branches.jl")
include("devices/device_models/DC_branches.jl")
include("devices/device_models/storage.jl")
include("devices/device_models/hydro_generation.jl")

#Network models
include("network_models/copperplate_model.jl")
include("network_models/powermodels_interface.jl")
include("network_models/ptdf_model.jl")

#Device constructors
include("devices/device_constructors/device_constructors.jl")

#Network constructors
include("network_models/network_constructor.jl")

#Services Models
#include("service_models/reserves.jl")

#Services constructors
include("service_models/services_constructor.jl")

#Operational Model Constructors
include("operation_models/operation_models.jl")

#Simulations Model Files
include("simulation_models/stage_update.jl")

#Routines
include("routines/make_initial_conditions.jl")
include("routines/get_results.jl")
include("routines/solve_routines.jl")
include("routines/write_model.jl")

#Utils
include("utils/optimization_debugging.jl")
include("utils/printing.jl")
include("utils/plot_results.jl")
include("utils/plot_recipes.jl")

#Initialization

function __init__()
   Requires.@require Weave = "44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9" include("utils/make_report.jl")
end

end


text:
#LogitMarginLoss

using LaTeXStrings, LossFunctions, Measures, Latexify

Plots.plot(LogitMarginLoss(), c=:olive, leg=:topright, linewidth=2, annotation=(1, 1.5, L"L\left( a \right) = ln(1+e^{-a})", 18))

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations