text:
"""
```
sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) + sum(pconvac[c] for c in bus_convs) - pd - gs*v^2
sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) + sum(qconvac[c] for c in bus_convs) - qd + bs*v^2
```
"""
function constraint_power_balance_ac(pm::_PM.AbstractACPModel, n::Int,  i::Int, bus_arcs, bus_arcs_dc, bus_gens, bus_convs_ac, bus_loads, bus_shunts, pd, qd, gs, bs)
    vm = _PM.var(pm, n,  :vm, i)
    p = _PM.var(pm, n,  :p)
    q = _PM.var(pm, n,  :q)
    pg = _PM.var(pm, n,  :pg)
    qg = _PM.var(pm, n,  :qg)
    pconv_grid_ac = _PM.var(pm, n,  :pconv_tf_fr)
    qconv_grid_ac = _PM.var(pm, n,  :qconv_tf_fr)

    JuMP.@NLconstraint(pm.model, sum(p[a] for a in bus_arcs) + sum(pconv_grid_ac[c] for c in bus_convs_ac)  == sum(pg[g] for g in bus_gens)   - sum(pd[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*vm^2)
    JuMP.@NLconstraint(pm.model, sum(q[a] for a in bus_arcs) + sum(qconv_grid_ac[c] for c in bus_convs_ac)  == sum(qg[g] for g in bus_gens)  - sum(qd[d] for d in bus_loads) + sum(bs[s] for s in bus_shunts)*vm^2)
end
"""
Creates Ohms constraints for DC branches

```
p[f_idx] == p * g[l] * vmdc[f_bus] * (vmdc[f_bus] - vmdc[t_bus])
```
"""
function constraint_ohms_dc_branch(pm::_PM.AbstractACPModel, n::Int,  f_bus, t_bus, f_idx, t_idx, r, p)
    p_dc_fr = _PM.var(pm, n,  :p_dcgrid, f_idx)
    p_dc_to = _PM.var(pm, n,  :p_dcgrid, t_idx)
    vmdc_fr = _PM.var(pm, n,  :vdcm, f_bus)
    vmdc_to = _PM.var(pm, n,  :vdcm, t_bus)

    if r == 0
        JuMP.@constraint(pm.model, p_dc_fr + p_dc_to == 0)
    else
        g = 1 / r
        JuMP.@NLconstraint(pm.model, p_dc_fr == p * g * vmdc_fr * (vmdc_fr - vmdc_to))
        JuMP.@NLconstraint(pm.model, p_dc_to == p * g * vmdc_to * (vmdc_to - vmdc_fr))
    end
end
"`vdc[i] == vdcm`"
function constraint_dc_voltage_magnitude_setpoint(pm::_PM.AbstractACPModel, n::Int,  i, vdcm)
    v = _PM.var(pm, n,  :vdcm, i)
    JuMP.@constraint(pm.model, v == vdcm)
end

function constraint_dc_branch_current(pm::_PM.AbstractACPModel, n::Int,  f_bus, f_idx, ccm_max, p)
# do nothing
end

############## TNEP constraints ###############
"""
```
sum(p[a] for a in bus_arcs) + sum(pconv_grid_ac[c] for c in bus_convs_ac) + sum(pconv_grid_ac_ne[c] for c in bus_convs_ac_ne) == sum(pg[g] for g in bus_gens)  - pd - gs*v^2
sum(q[a] for a in bus_arcs) + sum(qconv_grid_ac[c] for c in bus_convs_ac) + sum(qconv_grid_ac_ne[c] for c in bus_convs_ac_ne) == sum(qg[g] for g in bus_gens)  - qd + bs*v^2
```
"""
function constraint_power_balance_ac_dcne(pm::_PM.AbstractACPModel, n::Int, i::Int, bus_arcs, bus_arcs_dc, bus_gens, bus_convs_ac, bus_convs_ac_ne, bus_loads, bus_shunts, pd, qd, gs, bs)
    vm = _PM.var(pm, n, :vm, i)
    p = _PM.var(pm, n, :p)
    q = _PM.var(pm, n, :q)
    pg = _PM.var(pm, n, :pg)
    qg = _PM.var(pm, n, :qg)
    pconv_grid_ac = _PM.var(pm, n, :pconv_tf_fr)
    qconv_grid_ac = _PM.var(pm, n, :qconv_tf_fr)
    pconv_grid_ac_ne = _PM.var(pm, n, :pconv_tf_fr_ne)
    qconv_grid_ac_ne = _PM.var(pm, n, :qconv_tf_fr_ne)
    JuMP.@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(pconv_grid_ac[c] for c in bus_convs_ac) + sum(pconv_grid_ac_ne[c] for c in bus_convs_ac_ne)  == sum(pg[g] for g in bus_gens)  - sum(pd[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*vm^2)
    JuMP.@constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(qconv_grid_ac[c] for c in bus_convs_ac) + sum(qconv_grid_ac_ne[c] for c in bus_convs_ac_ne)  == sum(qg[g] for g in bus_gens)  - sum(qd[d] for d in bus_loads) + sum(bs[s] for s in bus_shunts)*vm^2)
end
"""
```
sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) + sum(pconv_grid_ac[c] for c in bus_convs_ac) + sum(pconv_grid_ac_ne[c] for c in bus_convs_ac_ne) == sum(pg[g] for g in bus_gens)  - pd - gs*v^2
sum(q[a] for a in bus_arcs) + sum(q_ne[a] for a in bus_arcs_ne) + sum(qconv_grid_ac[c] for c in bus_convs_ac) + sum(qconv_grid_ac_ne[c] for c in bus_convs_ac_ne) == sum(qg[g] for g in bus_gens)  - qd + bs*v^2
```
"""
function constraint_power_balance_acne_dcne(pm::_PM.AbstractACPModel, n::Int, i::Int, bus_arcs, bus_arcs_ne, bus_arcs_dc, bus_gens, bus_convs_ac, bus_convs_ac_ne, bus_loads, bus_shunts, pd, qd, gs, bs)
    vm = _PM.var(pm, n, :vm, i)
    p = _PM.var(pm, n, :p)
    p_ne = _PM.var(pm,n, :p_ne)
    q = _PM.var(pm, n, :q)
    q_ne = _PM.var(pm,n, :q_ne)
    pg = _PM.var(pm, n, :pg)
    qg = _PM.var(pm, n, :qg)
    pconv_grid_ac = _PM.var(pm, n, :pconv_tf_fr)
    qconv_grid_ac = _PM.var(pm, n, :qconv_tf_fr)
    pconv_grid_ac_ne = _PM.var(pm, n, :pconv_tf_fr_ne)
    qconv_grid_ac_ne = _PM.var(pm, n, :qconv_tf_fr_ne)
    JuMP.@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) + sum(pconv_grid_ac[c] for c in bus_convs_ac) + sum(pconv_grid_ac_ne[c] for c in bus_convs_ac_ne)  == sum(pg[g] for g in bus_gens)  - sum(pd[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*vm^2)
    JuMP.@constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_ne[a] for a in bus_arcs_ne) + sum(qconv_grid_ac[c] for c in bus_convs_ac) + sum(qconv_grid_ac_ne[c] for c in bus_convs_ac_ne)  == sum(qg[g] for g in bus_gens)  - sum(qd[d] for d in bus_loads) + sum(bs[s] for s in bus_shunts)*vm^2)
end
"""
Creates Ohms constraints for DC branches

```
p[f_idx] == p * g[l] * vmdc[f_bus] * (vmdc[f_bus] - vmdc[t_bus])
```
"""
function constraint_ohms_dc_branch_ne(pm::_PM.AbstractACPModel, n::Int,  f_bus, t_bus, f_idx, t_idx, r, p)
    p_dc_fr = _PM.var(pm, n, :p_dcgrid_ne, f_idx)
    p_dc_to = _PM.var(pm, n, :p_dcgrid_ne, t_idx)
    vmdc_fr = []
    vmdc_to = []
    z = _PM.var(pm, n, :branchdc_ne, f_idx[1])
    vmdc_to, vmdc_fr = contraint_ohms_dc_branch_busvoltage_structure(pm, n, f_bus, t_bus, vmdc_to, vmdc_fr)
    if r == 0
        JuMP.@constraint(pm.model, p_dc_fr + p_dc_to == 0)
    else
        g = 1 / r

        JuMP.@NLconstraint(pm.model, p_dc_fr == z * p * g * vmdc_fr * (vmdc_fr - vmdc_to))
        JuMP.@NLconstraint(pm.model, p_dc_to == z * p * g * vmdc_to * (vmdc_to - vmdc_fr))
    end
end

# V or W as input, the function returns V or W with existing dc bus and dc_ne bus voltages
function contraint_ohms_dc_branch_busvoltage_structure(pm::_PM.AbstractACPModel, n::Int, f_bus, t_bus, vmdc_to, vmdc_fr)
    for i in _PM.ids(pm, n, :busdc_ne)
        if t_bus == i
            vmdc_to = _PM.var(pm, n, :vdcm_ne, t_bus)
        end
        if f_bus == i
            vmdc_fr = _PM.var(pm, n, :vdcm_ne, f_bus)
        end
    end
    for i in _PM.ids(pm, n, :busdc)
        if t_bus == i
            vmdc_to = _PM.var(pm, n, :vdcm, t_bus)
        end
        if f_bus == i
            vmdc_fr = _PM.var(pm, n, :vdcm, f_bus)
        end
    end
    return vmdc_to, vmdc_fr
end


text:
abstract type BasicMonteCarloAlgorithm <: DiffEqBase.MonteCarloAlgorithm end
struct MonteThreads <: BasicMonteCarloAlgorithm end
struct MonteDistributed <: BasicMonteCarloAlgorithm end
struct MonteSplitThreads <: BasicMonteCarloAlgorithm end
struct MonteSerial <: BasicMonteCarloAlgorithm end

#=
if (kwargs[:parallel_type] == != :none && kwargs[:parallel_type] == != :threads)
  error("Distributed arrays cannot be generated via none or threads")
end
(batch_size != num_monte) && warn("batch_size and reductions are ignored when !collect_result")

elapsed_time = @elapsed u = DArray((num_monte,)) do I
    solve_batch(prob,alg,kwargs[:parallel_type] ==,I[1],pmap_batch_size,kwargs...)
end
return MonteCarloSolution(u,elapsed_time,false)
=#

function DiffEqBase.__solve(prob::DiffEqBase.AbstractMonteCarloProblem,
                            alg::Union{DiffEqBase.DEAlgorithm,Nothing};
                            kwargs...)
    if alg isa DiffEqBase.MonteCarloAlgorithm
      @error "You forgot to pass a DE solver algorithm! Only a MonteCarloAlgorithm has been supplied. Exiting"
    end
    if :parallel_type ∈ keys(kwargs)
      if kwargs[:parallel_type] == :none
        montealg = MonteSerial()
      elseif kwargs[:parallel_type] == :pmap || kwargs[:parallel_type] == :parfor
        montealg = MonteDistributed()
      elseif kwargs[:parallel_type] == :threads
        montealg = MonteThreads()
      elseif kwargs[:parallel_type] == :split_threads
        montealg = MonteSplitThreads()
      else
        @error "parallel_type value not recognized"
      end
    else
      montealg = MonteSerial()
    end
    if :parallel_type ∈ keys(kwargs)
      @warn "parallel_type has been deprecated. Please refer to the docs for the new dispatch-based system."
    end
    DiffEqBase.__solve(prob,alg,montealg;kwargs...)
end

function DiffEqBase.__solve(prob::DiffEqBase.AbstractMonteCarloProblem,
                 alg::Union{DiffEqBase.DEAlgorithm,Nothing},
                 montealg::BasicMonteCarloAlgorithm;
                 num_monte, batch_size = num_monte,
                 pmap_batch_size = batch_size÷100 > 0 ? batch_size÷100 : 1, kwargs...)

  num_batches = num_monte ÷ batch_size
  num_batches * batch_size != num_monte && (num_batches += 1)

  u = deepcopy(prob.u_init)
  converged = false
  elapsed_time = @elapsed for i in 1:num_batches
    if i == num_batches
      I = (batch_size*(i-1)+1):num_monte
    else
      I = (batch_size*(i-1)+1):batch_size*i
    end
    batch_data = solve_batch(prob,alg,montealg,I,pmap_batch_size,kwargs...)
    u,converged = prob.reduction(u,batch_data,I)
    converged && break
  end
  if typeof(u) <: Vector{Any}
    _u = convert(Array{typeof(u[1])},u)
  else
    _u = u
  end
  return MonteCarloSolution(_u,elapsed_time,converged)
end

function batch_func(i,prob,alg,I,kwargs...)
  iter = 1
  new_prob = prob.prob_func(deepcopy(prob.prob),i,iter)
  rerun = true
  x = prob.output_func(solve(new_prob,alg;kwargs...),i)
  if !(typeof(x) <: Tuple)
      @warn("output_func should return (out,rerun). See docs for updated details")
      _x = (x,false)
  else
    _x = x
  end
  rerun = _x[2]
  while rerun
      iter += 1
      new_prob = prob.prob_func(deepcopy(prob.prob),i,iter)
      x = prob.output_func(solve(new_prob,alg;kwargs...),i)
      if !(typeof(x) <: Tuple)
          @warn("output_func should return (out,rerun). See docs for updated details")
          _x = (x,false)
      else
        _x = x
      end
      rerun = _x[2]
  end
  _x[1]
end

function solve_batch(prob,alg,::MonteDistributed,I,pmap_batch_size,kwargs...)
  wp=CachingPool(workers())
  batch_data = let
    pmap(wp,I,batch_size=pmap_batch_size) do i
      batch_func(i,prob,alg,I,kwargs...)
    end
  end
  _batch_data = convert(Array{typeof(batch_data[1])},batch_data)
end

function solve_batch(prob,alg,::MonteSerial,I,pmap_batch_size,kwargs...)
  batch_data = let
    map(I) do i
      batch_func(i,prob,alg,I,kwargs...)
    end
  end
  _batch_data = convert(Array{typeof(batch_data[1])},batch_data)
end

function solve_batch(prob,alg,::MonteThreads,I,pmap_batch_size,kwargs...)
  batch_data = Vector{Any}(undef,length(I))
  let
    Threads.@threads for batch_idx in axes(batch_data, 1)
        i = I[batch_idx]
        iter = 1
        new_prob = prob.prob_func(deepcopy(prob.prob),i,iter)
        x = prob.output_func(solve(new_prob,alg;kwargs...),i)
        if !(typeof(x) <: Tuple)
            warn("output_func should return (out,rerun). See docs for updated details")
            _x = (x,false)
        else
          _x = x
        end
        rerun = _x[2]

        while rerun
            iter += 1
            new_prob = prob.prob_func(deepcopy(prob.prob),i,iter)
            x = prob.output_func(solve(new_prob,alg;kwargs...),i)
            if !(typeof(x) <: Tuple)
                warn("output_func should return (out,rerun). See docs for updated details")
                _x = (x,false)
            else
              _x = x
            end
            rerun = _x[2]
        end
        batch_data[batch_idx] = _x[1]
    end
  end
  _batch_data = convert(Array{typeof(batch_data[1])},batch_data)
end

function solve_batch(prob,alg,::MonteSplitThreads,I,pmap_batch_size,kwargs...)
  wp=CachingPool(workers())
  batch_data = let
    pmap(wp,1:nprocs(),batch_size=pmap_batch_size) do i
      thread_monte(prob,I,alg,i,kwargs...)
    end
  end
  _batch_data = vector_batch_data_to_arr(batch_data)
end

function thread_monte(prob,I,alg,procid,kwargs...)
  start = I[1]+(procid-1)*length(I)
  stop = I[1]+procid*length(I)-1
  portion = start:stop
  batch_data = Vector{Any}(undef,length(portion))
  let
    Threads.@threads for i in portion
      iter = 1
      new_prob = prob.prob_func(deepcopy(prob.prob),i,iter)
      rerun = true
      x = prob.output_func(solve(new_prob,alg;kwargs...),i)
      if !(typeof(x) <: Tuple)
          warn("output_func should return (out,rerun). See docs for updated details")
          _x = (x,false)
      else
        _x = x
      end
      rerun = _x[2]
      while rerun
          iter += 1
          new_prob = prob.prob_func(deepcopy(prob.prob),i,iter)
          x = prob.output_func(solve(new_prob,alg;kwargs...),i)
          if !(typeof(x) <: Tuple)
              warn("output_func should return (out,rerun). See docs for updated details")
              _x = (x,false)
          else
            _x = x
          end
          rerun = _x[2]
      end
      batch_data[i - start + 1] = _x[1]
    end
  end
  batch_data
end

function vector_batch_data_to_arr(batch_data)
  _batch_data = Vector{typeof(batch_data[1][1])}(undef,sum((length(x) for x in batch_data)))
  idx = 0
  @inbounds for a in batch_data
    for x in a
      idx += 1
      _batch_data[idx] = x
    end
  end
  _batch_data
end

ext: jl


text:
module PSlibrary

using DataFrames: DataFrame, groupby
using AcuteML: UN, parsehtml, root, nextelement, nodecontent
using JLD2: load
using LazyArtifacts
using MLStyle: @match
using Pseudopotentials:
    CoreHole, ExchangeCorrelationFunctional, ValenceCoreState, Pseudization
using Pkg.Artifacts: @artifact_str
using REPL.TerminalMenus: RadioMenu, request

using ..UnifiedPseudopotentialFormat: UPFFile, analyzename

import Pseudopotentials: download_potential

export list_elements,
    list_potentials, download_potentials, download_potential, potential_exists

const ELEMENTS = (
    "h",
    "he",
    "li",
    "be",
    "b",
    "c",
    "n",
    "o",
    "f",
    "ne",
    "na",
    "mg",
    "al",
    "si",
    "p",
    "s",
    "cl",
    "ar",
    "k",
    "ca",
    "sc",
    "ti",
    "v",
    "cr",
    "mn",
    "fe",
    "co",
    "ni",
    "cu",
    "zn",
    "ga",
    "ge",
    "as",
    "se",
    "br",
    "kr",
    "rb",
    "sr",
    "y",
    "zr",
    "nb",
    "mo",
    "tc",
    "ru",
    "rh",
    "pd",
    "ag",
    "cd",
    "in",
    "sn",
    "sb",
    "te",
    "i",
    "xe",
    "cs",
    "ba",
    "la",
    "ce",
    "pr",
    "nd",
    "pm",
    "sm",
    "eu",
    "gd",
    "tb",
    "dy",
    "ho",
    "er",
    "tm",
    "yb",
    "lu",
    "hf",
    "ta",
    "w",
    "re",
    "os",
    "ir",
    "pt",
    "au",
    "hg",
    "tl",
    "pb",
    "bi",
    "po",
    "at",
    "rn",
    "fr",
    "ra",
    "ac",
    "th",
    "pa",
    "u",
    "np",
    "pu",
)
const PERIODIC_TABLE = raw"""
H                                                  He
Li Be                               B  C  N  O  F  Ne
Na Mg                               Al Si P  S  Cl Ar
K  Ca Sc Ti V  Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr
Rb Sr Y  Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I  Xe
Cs Ba    Hf Ta W  Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn
Fr Ra
      La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu
      Ac Th Pa U  Np Pu
"""

"""
    list_elements(pt=true)

List all elements that has pseudopotentials available in `PSlibrary`. Print the periodic table if `pt` is `true`.
"""
function list_elements(pt = true)
    if pt
        println(PERIODIC_TABLE)
    end
    return loaddb()
end

"""
    list_potentials(element::Union{AbstractString,AbstractChar,Integer})

List all pseudopotentials in `PSlibrary` for a specific element (abbreviation or index).
"""
function list_potentials(element::Union{AbstractString,AbstractChar})
    element = lowercase(string(element))
    @assert element in ELEMENTS "element $element is not recognized!"
    return loaddb(element)
end
function list_potentials(atomic_number::Integer)
    @assert 1 <= atomic_number <= 94
    element = ELEMENTS[atomic_number]
    return list_potentials(element)
end

"""
    download_potential(element::Union{AbstractString,AbstractChar,Integer})

Download one or multiple pseudopotentials from `PSlibrary` for a specific element.
"""
function download_potentials(element)
    df = list_potentials(element)
    display(df)
    paths, finished = String[], false
    while !finished
        printstyled("Enter its index (integer) to download a potential: "; color = :green)
        i = parse(Int, readline())
        printstyled(
            "Enter the file path to save the potential (press enter to skip): ";
            color = :green,
        )
        str = readline()
        path = abspath(expanduser(isempty(str) ? tempname() : strip(str)))  # `abspath` is necessary since the path will depend on where you run it
        download(df.src[i], path)
        push!(paths, path)
        finished = request("Finished?", RadioMenu(["yes", "no"])) == 1
    end
    return paths
end

function potential_exists(name::AbstractString)
    x = analyzename(UPFFile(name))
    db = list_potentials(x.element)
    return name in db.name
end

function download_potential(name::AbstractString, path)
    x = analyzename(UPFFile(name))
    df = list_potentials(x.element)
    if name in df.name
        i = findfirst(==(name), df.name)
        return download(df.src[i], path)
    else
        throw("potential '$name' is not in the database!")
    end
end

function loaddb(element::AbstractString)
    artifact_path = joinpath(@artifact_str(element), "$element.jld2")
    return load(artifact_path)["database"]
end
loaddb(i::Integer) = loaddb(ELEMENTS[i])
function loaddb()
    artifact_path = joinpath(artifact"all", "all.jld2")
    return load(artifact_path)["database"]
end

end


text:
using Crux, Flux, POMDPs, POMDPGym, Distributions, POMDPPolicies, HDF5
init_mujoco_render() # Required for visualization

include("d4rl_utils.jl")

# Construct the Mujoco environment
mdp = GymPOMDP(Symbol("hopper-medium"), version=:v0)
S = state_space(mdp)
adim = length(POMDPs.actions(mdp)[1])
amin = -1*ones(Float32, adim)
amax = 1*ones(Float32, adim)

filename = "/home/anthonycorso/.d4rl/datasets/hopper_medium.hdf5"

## (Load from bson) get expert trajectories
𝒟_train = h5toEB(filename)
# μ_s = mean(𝒟_train[:s], dims=2)[:]
# σ_s = std(𝒟_train[:s], dims=2)[:] .+ 1f-3
# S = state_space(mdp, μ=μ_s, σ=σ_s)
S = state_space(mdp)

## Behavioral Cloning
A() = ContinuousNetwork(Chain(Dense(S.dims[1], 256, relu), Dense(256, 256, relu), Dense(256, 6, tanh)))
function SAC_A()
    base = Chain(Dense(S.dims[1], 256, relu), Dense(256, 256, relu))
    mu = ContinuousNetwork(Chain(base..., Dense(256, adim)))
    logΣ = ContinuousNetwork(Chain(base..., Dense(256, adim)))
    SquashedGaussianPolicy(mu, logΣ)
end

Q() = ContinuousNetwork(Chain(Dense(S.dims[1] + adim, 256, relu), Dense(256, 256, relu), Dense(256, 1)))

## Behavioral Cloning
𝒮_bc = BC(π=SAC_A(),
          𝒟_demo=𝒟_train,
          S=S,
          opt=(epochs=100000, batch_size=1024),
          log=(period=100, fns=[log_undiscounted_return(1)]),
          max_steps=1000)
solve(𝒮_bc, mdp)

## Batch SAC
𝒮_sac = BatchSAC(;π=ActorCritic(SAC_A(), DoubleNetwork(Q(), Q())),
                     𝒟_train=𝒟_train,
                     S=S,
                     log = (;period=100, fns=[log_undiscounted_return(1)]))
solve(𝒮_sac, mdp)

## CQL
𝒮_cql = CQL(;π=ActorCritic(SAC_A(), DoubleNetwork(Q(), Q())),
             𝒟_train=𝒟_train,
             S=S,
             max_steps=1000,
             CQL_α_thresh = 10f0,
             a_opt=(epochs=100, batch_size=256, optimizer=ADAM(1e-4)),
             c_opt=(batch_size=256, optimizer=ADAM(3f-4),),
             SAC_α_opt=(batch_size=256, optimizer=ADAM(1f-4),),
             CQL_α_opt=(batch_size=256, optimizer=ADAM(1f-4),),
             log=(period=500, fns=[log_undiscounted_return(1)]),)
solve(𝒮_cql, mdp)

max_stars_repo_path: examples/offline rl/hopper_medium.jl

max_issues_repo_path: examples/offline rl/hopper_medium.jl

max_forks_repo_path: examples/offline rl/hopper_medium.jl


text:
using Unitful: Energy, Frequency, Length, Mass, Wavenumber, c0, h, ħ, unit, ustrip

"""
    MassEnergy()

Equivalence to convert between mass and energy according to the relation ``E = mc^2``, where
* ``E`` is the energy,
* ``m`` is the mass and
* ``c`` is the speed of light in vacuum.

# Example

```jldoctest
julia> uconvert(u"keV", 1u"me", MassEnergy()) # electron rest mass is equivalent to ≈511 keV
510.9989499961642 keV
```
"""
struct MassEnergy <: Equivalence end

const c² = Quantity(convert(Int64, ustrip(c0)), unit(c0))^2

@eqrelation MassEnergy Energy/Mass = c²

"""
    PhotonEnergy(; frequency=:linear, wavelength=:linear, wavenumber=:linear)

Equivalence that relates the energy of a photon to its (linear or angular) frequency,
wavelength, and wavenumber. Whether to convert to linear or angular quantities is
determined by optional keyword arguments, `:linear` is the default for all quantities.

Equivalent quantities are converted according to the relations
``E = hf = ħω = hc/λ = ħc/ƛ = hcν̃ = ħck``, where
* ``E`` is the photon energy,
* ``f`` is the (temporal) frequency (`frequency=:linear`),
* ``ω`` is the angular frequency (`frequency=:angular`),
* ``λ`` is the wavelength (`wavelength=:linear`),
* ``ƛ`` is the angular (also called reduced) wavelength (`wavelength=:angular`),
* ``ν̃`` is the spectroscopic wavenumber (`wavenumber=:linear`),
* ``k`` is the angular wavenumber (`wavelength=:angular`),
* ``h`` is the Planck constant,
* ``ħ`` is the reduced Planck constant and
* ``c`` is the speed of light in vacuum.

# Examples

```jldoctest
julia> uconvert(u"nm", 13.6u"eV", PhotonEnergy()) # photon wavelength needed to ionize hydrogen
91.16485178911785 nm

julia> uconvert(u"Hz", 589u"nm", PhotonEnergy(frequency=:angular)) # angular frequency of sodium D line
3.1980501991661345e15 Hz
```
"""
struct PhotonEnergy{freq, len, num} <: Equivalence
    function PhotonEnergy{freq, len, num}() where {freq, len, num}
        check_photonarg(freq)
        check_photonarg(len)
        check_photonarg(num)
        new()
    end
end

@inline check_photonarg(s::Symbol) =
    s === :linear || s === :angular || throw(ArgumentError("PhotonEnergy parameter must be :linear or :angular"))

PhotonEnergy(; frequency=:linear, wavelength=:linear, wavenumber=:linear) =
    PhotonEnergy{frequency, wavelength, wavenumber}()

Base.show(io::IO, e::PhotonEnergy{freq, len, num}) where {freq, len, num} =
    print(io, PhotonEnergy, "(frequency=", repr(freq), ", wavelength=", repr(len), ", wavenumber=", repr(num), ")")

@eqrelation PhotonEnergy{:linear}  Energy/Frequency = h
@eqrelation PhotonEnergy{:angular} Energy/Frequency = ħ

@eqrelation PhotonEnergy{F,:linear}  where F Energy*Length = h*c0
@eqrelation PhotonEnergy{F,:angular} where F Energy*Length = ħ*c0

@eqrelation PhotonEnergy{F,L,:linear}  where {F,L} Energy/Wavenumber = h*c0
@eqrelation PhotonEnergy{F,L,:angular} where {F,L} Energy/Wavenumber = ħ*c0

function edconvert(::dimtype(Frequency), x::Length, ::PhotonEnergy{F,L}) where {F,L}
    F === L       ? c0/x      :
    F === :linear ? c0/(π*2x) : 2c0*(π/x)
end
function edconvert(::dimtype(Length), x::Frequency, ::PhotonEnergy{F,L}) where {F,L}
    F === L       ? c0/x      :
    F === :linear ? c0/(π*2x) : 2c0*(π/x)
end

function edconvert(::dimtype(Frequency), x::Wavenumber, ::PhotonEnergy{F,L,N}) where {F,L,N}
    F === N       ? c0*x       :
    F === :linear ? c0/2*(x/π) : 2c0*(π*x)
end
function edconvert(::dimtype(Wavenumber), x::Frequency, ::PhotonEnergy{F,L,N}) where {F,L,N}
    F === N       ? x/c0       :
    F === :linear ? 2*(π*x)/c0 : (x/π)/2c0
end

function edconvert(::dimtype(Length), x::Wavenumber, ::PhotonEnergy{F,L,N}) where {F,L,N}
    L === N       ? inv(x)  :
    L === :linear ? 2*(π/x) : inv(2*(π*x))
end
function edconvert(::dimtype(Wavenumber), x::Length, ::PhotonEnergy{F,L,N}) where {F,L,N}
    L === N       ? inv(x)  :
    L === :linear ? 2*(π/x) : inv(2*(π*x))
end


text:
@testset "Struct accessors" begin
    @test PH.value(PH.ScenarioID(1)) == 1
    @test PH.value(PH.StageID(2)) == 2
    @test PH.value(PH.Index(3)) == 3
    @test PH.value(PH.NodeID(4)) == 4
end

@testset "Convenience constructors" begin
    @test PH.scid(1) == PH.ScenarioID(1)
    @test PH.stid(1) == PH.StageID(1)
    @test PH.index(1) == PH.Index(1)
end

@testset "Struct less than operators" begin
    @test PH.NodeID(3) < PH.NodeID(7)
    @test PH.scid(0) < PH.scid(2)
    @test PH.stid(1) < PH.stid(4)
    @test PH._increment(PH.stid(1)) == PH.stid(2)
    @test PH.index(10) < PH.index(16)
    @test PH._increment(PH.index(15)) == PH.index(16)

    uvid1 = PH.VariableID(PH.scid(4), PH.stid(1), PH.index(347))
    uvid2 = PH.VariableID(PH.scid(0), PH.stid(2), PH.index(2))
    uvid3 = PH.VariableID(PH.scid(1), PH.stid(2), PH.index(2))
    uvid4 = PH.VariableID(PH.scid(1), PH.stid(2), PH.index(3))
    @test uvid1 < uvid2
    @test uvid2 < uvid3
    @test uvid3 < uvid4

    xhid1 = PH.XhatID(PH.NodeID(0), PH.Index(4))
    xhid2 = PH.XhatID(PH.NodeID(1), PH.Index(2))
    xhid3 = PH.XhatID(PH.NodeID(1), PH.Index(3))
    @test xhid1 < xhid2
    @test xhid2 < xhid3
end

@testset "Consensus variable functions" begin
    tv = PH.HatVariable(true)
    @test PH.is_integer(tv)
    vid = PH.VariableID(PH.scid(3), PH.stid(2), PH.index(392))
    PH.add_variable(tv, vid)
    @test PH.variables(tv) == Set([vid])
    rv = rand()
    PH.set_value(tv, rv)
    @test tv.value == rv
    @test PH.value(tv) == rv
end

@testset "Convenience functions" begin
    st = two_stage_tree(3)
    @test length(st.tree_map) == 4
    @test length(st.prob_map) == 3
    @test isapprox(sum(values(st.prob_map)), 1.0)

    p = [0.8, 0.2]
    st = two_stage_tree(p)
    @test st.prob_map[PH.scid(0)] == p[1]
    @test st.prob_map[PH.scid(1)] == p[2]
end

nv1 = 3
nv2 = 4
phd = fake_phdata(nv1, nv2)
nscen = length(PH.scenarios(phd))

@testset "Accessor Utilities" begin
    scid = PH.scid(0)
    stid = PH.stid(1)
    index = PH.index(2)
    vid = PH.VariableID(scid, stid, index)
    @test PH.name(phd, vid) == "a2"
    val = convert(Float64, PH.value(index) - PH.value(scid))
    @test PH.value(phd, vid) == val
    @test PH.value(phd, scid, stid, index) == val
    @test PH.branch_value(phd, vid) == val
    @test PH.branch_value(phd, scid, stid, index) == val

    scid = PH.scid(1)
    stid = PH.stid(2)
    index = PH.index(3)
    vid = PH.VariableID(scid, stid, index)
    @test PH.name(phd, vid) == "b4"
    val = convert(Float64, nscen + (PH.value(index) + 1)*(PH.value(scid) + 1))
    @test PH.value(phd, vid) == val
    @test PH.value(phd, scid, stid, index) == val
    @test PH.leaf_value(phd, vid) == val
    @test PH.leaf_value(phd, scid, stid, index) == val

    @test_throws ErrorException PH.name(phd, PH.VariableID(PH.scid(8303), stid, index))

    @test scenario_tree(phd) === phd.scenario_tree

    time_str_1 = @capture_out(print_timing(phd))
    @test length(time_str_1) > 0
    io = IOBuffer()
    print_timing(io, phd)
    time_str_2 = String(take!(io))
    @test length(time_str_2) > 0
end

@testset "Conversion utilities" begin

    # Branch variables
    for i in 1:nv1
        vid1 = PH.VariableID(PH.scid(0), PH.stid(1), PH.index(i))
        vid2 = PH.VariableID(PH.scid(1), PH.stid(1), PH.index(i+1))
        xhid = PH.convert_to_xhat_id(phd, vid1)
        @test PH.convert_to_xhat_id(phd, vid2) == xhid
        @test PH.convert_to_variable_ids(phd, xhid) == Set{PH.VariableID}([vid1, vid2])
    end

    # Leaf variables
    for j in 1:nv2
        vid1 = PH.VariableID(PH.scid(0), PH.stid(2), PH.index(j-1))
        vid2 = PH.VariableID(PH.scid(1), PH.stid(2), PH.index(j-1))
        @test PH.convert_to_variable_ids(phd, PH.convert_to_xhat_id(phd, vid1)) == Set{PH.VariableID}([vid1])
        @test PH.convert_to_variable_ids(phd, PH.convert_to_xhat_id(phd, vid2)) == Set{PH.VariableID}([vid2])
        @test PH.convert_to_xhat_id(phd, vid1) != PH.convert_to_xhat_id(phd, vid2)
    end

end

@testset "PH variable accessors" begin
    scid = PH.scid(0)
    stid = PH.stid(1)
    index = PH.index(1)
    vid = PH.VariableID(scid, stid, index)
    rval = rand()

    phd.scenario_map[scid].w_vars[vid] = rval
    @test PH.w_value(phd, vid) == rval
    @test PH.w_value(phd, scid, stid, index) == rval

    rval = rand()
    xhid = PH.convert_to_xhat_id(phd, vid)

    phd.xhat[xhid].value = rval
    @test PH.xhat_value(phd, xhid) == rval
    @test PH.xhat_value(phd, vid) == rval
    @test PH.xhat_value(phd, scid, stid, index) == rval

    @test PH.is_leaf(phd, xhid) == false
    vid = PH.VariableID(PH.scid(0), PH.stid(2), PH.index(0))
    xhid = PH.convert_to_xhat_id(phd, vid)
    @test PH.is_leaf(phd, xhid) == true
end


text:
######################################################
### This test implements example 13.5 in Greene (2012)
######################################################
using MomentBasedEstimators
using LinearAlgebra
using SpecialFunctions

# this is the data from Example C.1 in Greene 2012
income_data = [20.5, 31.5, 47.7, 26.2, 44, 8.28, 30.8, 17.2, 19.9, 9.96,
               55.8, 25.2, 29, 85.5, 15.1, 28.5, 21.4, 17.7, 6.42, 84.9]

function g_gamma(θ)
    P, λ = θ[1], θ[2]
    hcat(income_data .- P/λ,
         income_data.^2 .- P*(P+1)/(λ^2),
         log.(income_data) .+ log(λ) .- SpecialFunctions.digamma(P),
         1 ./income_data .- λ/(P-1))
end

two_step = GMMEstimator(g_gamma, [2.40, 0.08], initialW = Matrix(I, 4, 4),  mgr=TwoStepGMM())

MomentBasedEstimators.setparLB!(two_step, [0.0, 0.0])
MomentBasedEstimators.setparUB!(two_step, [+Inf, +Inf])

estimate!(two_step)


text:
mutable struct xf_indiv <: AbstractSolution # Single Objective
    x::Vector{Float64}
    f::Float64
end

mutable struct xfgh_indiv <: AbstractSolution # Single Objective Constraied
    x::Vector{Float64}
    f::Float64
    g::Vector{Float64}
    h::Vector{Float64}
    sum_violations::Float64 # ∑ max(0,g) + ∑|h|
    is_feasible::Bool

end

function xfgh_indiv(
    x::Vector{Float64},
    f::Float64,
    g::Vector{Float64},
    h::Vector{Float64};
    sum_violations = 0.0,
    ε = 0.0
)
    if sum_violations <= 0.0
        sum_violations = violationsSum(g, h; ε=ε)
    end

    xfgh_indiv(x, f, g, h, sum_violations, sum_violations == 0.0)
end

mutable struct xFgh_indiv <: AbstractSolution# Single Objective Constraied
    x::Vector{Float64}
    f::Vector{Float64}
    g::Vector{Float64}
    h::Vector{Float64}
    rank::Int
    crowding::Float64
    sum_violations::Float64 # ∑ max(0,g) + ∑|h|
    is_feasible::Bool
end

function xFgh_indiv(
    x::Vector{Float64},
    f::Vector{Float64},
    g::Vector{Float64},
    h::Vector{Float64};
    rank = 0,
    crowding = 0.0,
    sum_violations = 0.0,
    ε = 0.0
)

    if sum_violations <= 0
        sum_violations = violationsSum(g, h;ε=ε)
    end
    xFgh_indiv(x, f, g, h, Int(rank), crowding, sum_violations, sum_violations == 0.0)
end

const Solution = Union{xf_indiv, xfgh_indiv, xFgh_indiv}
const Population = Array{Solution, 1}

############################################################
# Generate solutions depending on the objective function
# output
#############################################################

function generateChild(x::Vector{Float64}, fResult::Float64;ε=0.0)
    return xf_indiv(x, fResult)
end

function generateChild(x::Vector{Float64},
        fResult::Tuple{Float64,Array{Float64,1},Array{Float64,1}};
        ε=0.0
    )
    f, g, h = fResult
    return xfgh_indiv(x, f, g, h; ε=ε)
end

function generateChild(x::Vector{Float64},
        fResult::Tuple{Array{Float64,1},Array{Float64,1},Array{Float64,1}};
        ε = 0.0
    )
    f, g, h = fResult
    return xFgh_indiv(x, f, g, h;ε=ε)
end

function generate_population(N::Int, problem;ε=0.0)
    a = view(problem.bounds, 1, :)'
    b = view(problem.bounds, 2, :)'
    D = length(a)

    X = a .+ (b - a) .* rand(N, D)

    population = [ create_solution(X[i,:], problem; ε=ε) for i in 1:N]

    return population
end

"""
    Metaheuristics.create_child(x, fx)

Constructor for a solution depending on the result of `fx`.

## Example

```
julia> import Metaheuristics

julia> Metaheuristics.create_child(rand(3), 1.0)
| f(x) = 1
| solution.x = [0.2700437125780806, 0.5233263210622989, 0.12871108215859772]

julia> Metaheuristics.create_child(rand(3), (1.0, [2.0, 0.2], [3.0, 0.3]))
| f(x) = 1
| g(x) = [2.0, 0.2]
| h(x) = [3.0, 0.3]
| x = [0.9881102595664819, 0.4816273348099591, 0.7742585077942159]

julia> Metaheuristics.create_child(rand(3), ([-1, -2.0], [2.0, 0.2], [3.0, 0.3]))
| f(x) = [-1.0, -2.0]
| g(x) = [2.0, 0.2]
| h(x) = [3.0, 0.3]
| x = [0.23983577719146854, 0.3611544510766811, 0.7998754930109109]

julia> population = [ Metaheuristics.create_child(rand(2), (randn(2),  randn(2), rand(2))) for i = 1:100  ]
                           F space
          ┌────────────────────────────────────────┐
        2 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠄⠀⠀⠂⠀⠀⠀⠀⠀⡇⠈⡀⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠘⠀⡇⠀⠀⠘⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠂⠀⠀⢀⠠⠐⠀⡇⠄⠁⠀⠀⠀⡀⠀⢁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠂⢈⠀⠈⡇⠀⡐⠃⠀⠄⠄⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠄⢐⠠⠀⡄⠀⠀⡇⠀⠂⠈⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠉⠉⠉⠉⠋⠉⠉⠉⠉⠉⠉⠙⢉⠉⠙⠉⠉⡏⠉⠉⠩⠋⠉⠩⠉⠉⠉⡉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉│
   f_2    │⠀⠀⠀⠀⠀⡀⠀⠀⠀⠄⠀⠀⡀⠀⠀⠂⠀⡇⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠄⠀⠀⠐⡇⠠⠀⠀⠀⠈⢀⠄⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠄⠀⡀⠀⠂⡇⠐⠘⠈⠂⠀⠈⡀⠀⠀⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠄⠀⠀⠀⠀⠀⠂⠀⠂⠀⠀⡇⠀⠈⢀⠐⠀⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⠀⠀⠀⠀⢁⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       -3 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          └────────────────────────────────────────┘
          -3                                       4
                             f_1
```

"""
create_child(x, fx) = generateChild(x, fx)

function create_solution(x, problem::Problem; ε=0.0)
    problem.f_calls += 1
    return generateChild(x, problem.f(x), ε=ε)
end

# getters for the above structures

"""
    get_position(solution)

Get the position vector.
"""
get_position(solution::Solution) = solution.x

function positions(population::AbstractArray)
    if isempty(population)
        zeros(0,0)
    end

    D = length(get_position(population[1]))
    [ get_position(sol)[i] for sol in population, i in 1:D]
end
"""
    fval(solution)

Get the objective function value (fitness) of a solution.
"""
fval(solution::Solution) = solution.f
sum_violations(solution::xfgh_indiv) = solution.sum_violations
sum_violations(solution::xFgh_indiv) = solution.sum_violations

function fvals(population::AbstractArray{xFgh_indiv})
    if isempty(population)
        return zeros(0, 0)
    end

    M = length(fval(population[1]))
    [ fval(sol)[i] for sol in population, i in 1:M]
end

fvals(population::AbstractArray) = fval.(population)

"""
    pareto_front(state::State)
Returns the non-dominated solutions in state.population.
"""
pareto_front(st::State) = pareto_front(st.population)

"""
    pareto_front(population::Array)
Returns non-dominated solutions.
"""
pareto_front(population::Array) = fvals(get_non_dominated_solutions(population))


text:
#floyd(n) creates an n-row floyd's triangle counting from 1 to (n/2+.5)*n
function floyd(n)
	x = 1
	dig(x,line,n) = (while line < n; x+=line; line+= 1 end; return ndigits(x)+1)
	for line = 1:n, i = 1:line; print(lpad(x,dig(x,line,n)," ")); x+=1; i==line && print("\n") end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Nemo
using Hecke
using ToyFHE
using Primes
using Test
using GaloisFields

const ℤ = ZZ
const ℤx, x = ZZ["x"]
const ℤpx, xp = ResidueRing(ZZ, 2)["xp"]
ℛ = ResidueRing(ℤpx, cyclotomic(7, x))

plain = PolyCRTEncoding(zero(ℛ))
β = GaloisFields.gen(eltype(plain))

plain[1] = β+1
plain[2] = β^2+1

p = convert(ℛ, plain)
enc = PolyCRTEncoding(p)
@test enc[1] == β+1
@test enc[2] == β^2+1


text:
function create_atomic_number_mapping()
  mapping::Dict{String,Int64} = Dict(
    "H" => 1,
    "He" => 2,
    "Li" => 3,
    "Be" => 4,
    "B" => 5,
    "C" => 6,
    "N" => 7,
    "O" => 8,
    "F" => 9,
    "Ne" => 10,
    "Na" => 11,
    "Mg" => 12,
    "Al" => 13,
    "Si" => 14,
    "P" => 15,
    "S" => 16,
    "Cl" => 17,
    "Ar" => 18,
    "K" => 19,
    "Ca" => 20,
    "Sc" => 21,
    "Ti" => 22,
    "V" => 23,
    "Cr" => 24,
    "Mn" => 25,
    "Fe" => 26,
    "Co" => 27,
    "Ni" => 28,
    "Cu" => 29,
    "Zn" => 30,
    "Ga" => 31,
    "Ge" => 32,
    "As" => 33,
    "Se" => 34,
    "Br" => 35,
    "Kr" => 36,
    "Rb" => 37,
    "Sr" => 38,
    "Y" => 39,
    "Zr" => 40,
    "Nb" => 41,
    "Mo" => 42,
    "Tc" => 43,
    "Ru" => 44,
    "Rh" => 45,
    "Pd" => 46,
    "Ag" => 47,
    "Cd" => 48,
    "In" => 49,
    "Sn" => 50,
    "Sb" => 51,
    "Te" => 52,
    "I" => 53,
    "Xe" => 54
  )
  return mapping
end

function create_atomic_mass_mapping()
  mapping::Dict{String,Float64} = Dict(
    "H" => 1.00794,
    "He" => 4.002602,
    "Li" => 6.941,
    "Be" => 9.012182,
    "B" => 10.811,
    "C" => 12.0107,
    "N" => 14.0067,
    "O" => 15.9994,
    "F" => 18.9984032,
    "Ne" => 20.1797,
    "Na" => 22.98976928,
    "Mg" => 24.3050,
    "Al" => 26.9815386,
    "Si" => 28.0855,
    "P" => 30.973762,
    "S" => 32.0655,
    "Cl" => 35.453,
    "Ar" => 39.948,
    "K" => 39.0983,
    "Ca" => 40.078,
    "Sc" => 44.955912,
    "Ti" => 47.867,
    "V" => 50.9415,
    "Cr" => 51.9961,
    "Mn" => 54.938045,
    "Fe" => 55.845,
    "Co" => 58.933195,
    "Ni" => 58.6934,
    "Cu" => 63.546,
    "Zn" => 65.409,
    "Ga" => 69.723,
    "Ge" => 72.64,
    "As" => 74.92160,
    "Se" => 78.96,
    "Br" => 79.804,
    "Kr" => 83.798,
    "Rb" => 85.4678,
    "Sr" => 87.62,
    "Y" => 88.90585,
    "Zr" => 91.224,
    "Nb" => 92.90638,
    "Mo" => 95.94,
    "Tc" => 98.0,
    "Ru" => 101.07,
    "Rh" => 102.90550,
    "Pd" => 106.42,
    "Ag" => 107.8682,
    "Cd" => 112.411,
    "In" => 114.818,
    "Sn" => 118.710,
    "Sb" => 121.760,
    "Te" => 127.60,
    "I" => 126.90447,
    "Xe" => 131.293
  )
  return mapping
end

function create_shell_am_mapping()
  mapping::Dict{String,Int64} = Dict(
    "L" => -1,
    "S" => 1,
    "P" => 2,
    "D" => 3,
    "F" => 4,
    "G" => 5,
    "H" => 6,
    "I" => 7,
  )
  return mapping
end


text:
"""
    βs(X::T, Y::T; dims::Union{Nothing,Integer}=nothing) where {T<:BinaryNetwork}

Components of β-diversity as measured on species.

#### References

Koleff, P., Gaston, K.J., Lennon, J.J., 2003. Measuring beta diversity for
presence–absence data. Journal of Animal Ecology 72, 367–382.
https://doi.org/10.1046/j.1365-2656.2003.00710.x
"""
function βs(X::T, Y::T; dims::Union{Nothing,Integer}=nothing) where {T<:BinaryNetwork}
    a = richness(intersect(X,Y); dims=dims)
    b = richness(Y; dims=dims)-a
    c = richness(X; dims=dims)-a
    return (a=a, b=b, c=c)
end

"""
    βos(X::T, Y::T) where {T<:BipartiteNetwork}

Overlapping species (bipartite)

#### References

Canard, E.F., Mouquet, N., Mouillot, D., Stanko, M., Miklisova, D., Gravel, D.,
2014. Empirical evaluation of neutral interactions in host-parasite networks.
The American Naturalist 183, 468–479. https://doi.org/10.1086/675363

Poisot, T., Canard, E., Mouillot, D., Mouquet, N., Gravel, D., 2012. The
dissimilarity of species interaction networks. Ecol. Lett. 15, 1353–1361.
https://doi.org/10.1111/ele.12002

Poisot, T., Cirtwill, A.R., Cazelles, K., Gravel, D., Fortin, M.-J., Stouffer,
D.B., 2016. The structure of probabilistic networks. Methods in Ecology and
Evolution 7, 303–312. https://doi.org/10.1111/2041-210X.12468
"""
function βos(X::T, Y::T) where {T<:BipartiteNetwork}
    core = intersect(X,Y)
    Xi = X[species(core; dims=1),species(core; dims=2)]
    Yi = Y[species(core; dims=1),species(core; dims=2)]
    a = links(intersect(Xi,Yi))
    b = links(Xi)-a
    c = links(Yi)-a
    return (a=a, b=b, c=c)
end

"""
    βos(X::T, Y::T) where {T<:UnipartiteNetwork}

Overlapping species (unipartite)

#### References

Canard, E.F., Mouquet, N., Mouillot, D., Stanko, M., Miklisova, D., Gravel, D.,
2014. Empirical evaluation of neutral interactions in host-parasite networks.
The American Naturalist 183, 468–479. https://doi.org/10.1086/675363

Poisot, T., Canard, E., Mouillot, D., Mouquet, N., Gravel, D., 2012. The
dissimilarity of species interaction networks. Ecol. Lett. 15, 1353–1361.
https://doi.org/10.1111/ele.12002

Poisot, T., Cirtwill, A.R., Cazelles, K., Gravel, D., Fortin, M.-J., Stouffer,
D.B., 2016. The structure of probabilistic networks. Methods in Ecology and
Evolution 7, 303–312. https://doi.org/10.1111/2041-210X.12468
"""
function βos(X::T, Y::T) where {T<:UnipartiteNetwork}
    core = intersect(X,Y)
    Xi = X[species(core)]
    Yi = Y[species(core)]
    a = links(intersect(Xi,Yi))
    b = links(Xi)-a
    c = links(Yi)-a
    return (a=a, b=b, c=c)
end

"""
    βwn(X::T, Y::T) where {T<:BinaryNetwork}

Whole network

#### References

Canard, E.F., Mouquet, N., Mouillot, D., Stanko, M., Miklisova, D., Gravel, D.,
2014. Empirical evaluation of neutral interactions in host-parasite networks.
The American Naturalist 183, 468–479. https://doi.org/10.1086/675363

Poisot, T., Canard, E., Mouillot, D., Mouquet, N., Gravel, D., 2012. The
dissimilarity of species interaction networks. Ecol. Lett. 15, 1353–1361.
https://doi.org/10.1111/ele.12002

Poisot, T., Cirtwill, A.R., Cazelles, K., Gravel, D., Fortin, M.-J., Stouffer,
D.B., 2016. The structure of probabilistic networks. Methods in Ecology and
Evolution 7, 303–312. https://doi.org/10.1111/2041-210X.12468
"""
function βwn(X::T, Y::T) where {T<:BinaryNetwork}
    a = links(intersect(X,Y))
    c = links(X)-a
    b = links(Y)-a
    return (a=a, b=b, c=c)
end


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7, x8, x9)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7, x8, x9]);F = [25448*x1^4-19984*x1^3*x3-14784*x1^3*x6-1980*x1^2*x2*x3-4048*x1^2*x2*x5+7456*x1^2*x3^2+11968*x1^2*x3*x5+3360*x1^2*x3*x6-4860*x1^2*x3*x7+6664*x1^2*x3*x8-7020*x1^2*x4*x5-17248*x1^2*x4*x7-9016*x1^2*x5*x8+14980*x1^2*x6^2-9604*x1^2*x6*x9-13720*x1^2*x8*x9+920*x1*x2*x3*x5+1964*x1*x2*x5*x6-2720*x1*x3^2*x5-5712*x1*x3^2*x8+3920*x1*x3*x4*x7-11424*x1*x3*x5*x6+7728*x1*x3*x5*x8+8232*x1*x3*x6*x9+11760*x1*x3*x8*x9+16464*x1*x4*x6*x7-5776*x1*x5*x6*x9+6084*x2^4+121*x2^2*x3^2-11388*x2^2*x3*x4-12324*x2^2*x3*x5+2652*x2^2*x4*x5+1154*x2^2*x5^2-2808*x2^2*x6*x9+5776*x2^2*x9^2+594*x2*x3^2*x7-12768*x2*x3^2*x9+858*x2*x3*x4*x5-3128*x2*x3*x5^2-792*x2*x3*x6^2-3496*x2*x3*x9^2+13224*x2*x4^2*x9+4508*x2*x4*x5*x7+3800*x2*x5^2*x9-3800*x2*x7*x9^2+7056*x3^4+3864*x3^3*x9-9287*x3^2*x4^2+11534*x3^2*x4*x5+10865*x3^2*x5^2+729*x3^2*x7^2+4200*x3^2*x7*x9+1156*x3^2*x8^2+529*x3^2*x9^2-2482*x3*x4^2*x5-4002*x3*x4^2*x9-2686*x3*x4*x5^2-11222*x3*x4*x5*x7+2628*x3*x4*x6*x9+2844*x3*x5*x6*x9-3128*x3*x5*x8^2-1944*x3*x6^2*x7-3332*x3*x6*x8*x9+1150*x3*x7*x9^2-4760*x3*x8^2*x9+7569*x4^4+1810*x4^2*x5^2+9604*x4^2*x7^2-4350*x4^2*x7*x9-2808*x4*x5*x6^2-612*x4*x5*x6*x9+2116*x5^2*x8^2+5776*x5^2*x9^2+4508*x5*x6*x8*x9+6440*x5*x8^2*x9+1296*x6^4+2725*x6^2*x9^2+6860*x6*x8*x9^2+625*x7^2*x9^2+4900*x8^2*x9^2-4680*x1^2*x5-6384*x1^2*x6+4200*x1*x2*x5-6992*x1*x5*x6+12768*x1*x5*x9+6764*x1*x6^2+572*x2*x3*x5+4600*x2*x5^2-4450*x2*x5*x6-10944*x2*x6*x9+12096*x3^2*x6+1404*x3*x5*x7+3312*x3*x6*x9-12528*x4^2*x6+2028*x4*x5^2+13984*x5^2*x9-1872*x5*x6^2-13528*x5*x6*x9+3600*x6*x7*x9+7056*x1^2+15456*x1*x5-14952*x1*x6-2028*x2^2+1898*x3*x4+2054*x3*x5-442*x4*x5+9140*x5^2-16376*x5*x6+13105*x6^2+468*x6*x9+169, -1980*x1^2*x3-4048*x1^2*x5+920*x1*x3*x5+1964*x1*x5*x6+24336*x2^3+242*x2*x3^2-22776*x2*x3*x4-24648*x2*x3*x5+5304*x2*x4*x5+2308*x2*x5^2-5616*x2*x6*x9+11552*x2*x9^2+594*x3^2*x7-12768*x3^2*x9+858*x3*x4*x5-3128*x3*x5^2-792*x3*x6^2-3496*x3*x9^2+13224*x4^2*x9+4508*x4*x5*x7+3800*x5^2*x9-3800*x7*x9^2+4200*x1*x5+572*x3*x5+4600*x5^2-4450*x5*x6-10944*x6*x9-4056*x2, -19984*x1^3-1980*x1^2*x2+14912*x1^2*x3+11968*x1^2*x5+3360*x1^2*x6-4860*x1^2*x7+6664*x1^2*x8+920*x1*x2*x5-5440*x1*x3*x5-11424*x1*x3*x8+3920*x1*x4*x7-11424*x1*x5*x6+7728*x1*x5*x8+8232*x1*x6*x9+11760*x1*x8*x9+242*x2^2*x3-11388*x2^2*x4-12324*x2^2*x5+1188*x2*x3*x7-25536*x2*x3*x9+858*x2*x4*x5-3128*x2*x5^2-792*x2*x6^2-3496*x2*x9^2+28224*x3^3+11592*x3^2*x9-18574*x3*x4^2+23068*x3*x4*x5+21730*x3*x5^2+1458*x3*x7^2+8400*x3*x7*x9+2312*x3*x8^2+1058*x3*x9^2-2482*x4^2*x5-4002*x4^2*x9-2686*x4*x5^2-11222*x4*x5*x7+2628*x4*x6*x9+2844*x5*x6*x9-3128*x5*x8^2-1944*x6^2*x7-3332*x6*x8*x9+1150*x7*x9^2-4760*x8^2*x9+572*x2*x5+24192*x3*x6+1404*x5*x7+3312*x6*x9+1898*x4+2054*x5, -7020*x1^2*x5-17248*x1^2*x7+3920*x1*x3*x7+16464*x1*x6*x7-11388*x2^2*x3+2652*x2^2*x5+858*x2*x3*x5+26448*x2*x4*x9+4508*x2*x5*x7-18574*x3^2*x4+11534*x3^2*x5-4964*x3*x4*x5-8004*x3*x4*x9-2686*x3*x5^2-11222*x3*x5*x7+2628*x3*x6*x9+30276*x4^3+3620*x4*x5^2+19208*x4*x7^2-8700*x4*x7*x9-2808*x5*x6^2-612*x5*x6*x9-25056*x4*x6+2028*x5^2+1898*x3-442*x5, -4048*x1^2*x2+11968*x1^2*x3-7020*x1^2*x4-9016*x1^2*x8+920*x1*x2*x3+1964*x1*x2*x6-2720*x1*x3^2-11424*x1*x3*x6+7728*x1*x3*x8-5776*x1*x6*x9-12324*x2^2*x3+2652*x2^2*x4+2308*x2^2*x5+858*x2*x3*x4-6256*x2*x3*x5+4508*x2*x4*x7+7600*x2*x5*x9+11534*x3^2*x4+21730*x3^2*x5-2482*x3*x4^2-5372*x3*x4*x5-11222*x3*x4*x7+2844*x3*x6*x9-3128*x3*x8^2+3620*x4^2*x5-2808*x4*x6^2-612*x4*x6*x9+4232*x5*x8^2+11552*x5*x9^2+4508*x6*x8*x9+6440*x8^2*x9-4680*x1^2+4200*x1*x2-6992*x1*x6+12768*x1*x9+572*x2*x3+9200*x2*x5-4450*x2*x6+1404*x3*x7+4056*x4*x5+27968*x5*x9-1872*x6^2-13528*x6*x9+15456*x1+2054*x3-442*x4+18280*x5-16376*x6, -14784*x1^3+3360*x1^2*x3+29960*x1^2*x6-9604*x1^2*x9+1964*x1*x2*x5-11424*x1*x3*x5+8232*x1*x3*x9+16464*x1*x4*x7-5776*x1*x5*x9-2808*x2^2*x9-1584*x2*x3*x6+2628*x3*x4*x9+2844*x3*x5*x9-3888*x3*x6*x7-3332*x3*x8*x9-5616*x4*x5*x6-612*x4*x5*x9+4508*x5*x8*x9+5184*x6^3+5450*x6*x9^2+6860*x8*x9^2-6384*x1^2-6992*x1*x5+13528*x1*x6-4450*x2*x5-10944*x2*x9+12096*x3^2+3312*x3*x9-12528*x4^2-3744*x5*x6-13528*x5*x9+3600*x7*x9-14952*x1-16376*x5+26210*x6+468*x9, -4860*x1^2*x3-17248*x1^2*x4+3920*x1*x3*x4+16464*x1*x4*x6+594*x2*x3^2+4508*x2*x4*x5-3800*x2*x9^2+1458*x3^2*x7+4200*x3^2*x9-11222*x3*x4*x5-1944*x3*x6^2+1150*x3*x9^2+19208*x4^2*x7-4350*x4^2*x9+1250*x7*x9^2+1404*x3*x5+3600*x6*x9, 6664*x1^2*x3-9016*x1^2*x5-13720*x1^2*x9-5712*x1*x3^2+7728*x1*x3*x5+11760*x1*x3*x9+2312*x3^2*x8-6256*x3*x5*x8-3332*x3*x6*x9-9520*x3*x8*x9+4232*x5^2*x8+4508*x5*x6*x9+12880*x5*x8*x9+6860*x6*x9^2+9800*x8*x9^2, -9604*x1^2*x6-13720*x1^2*x8+8232*x1*x3*x6+11760*x1*x3*x8-5776*x1*x5*x6-2808*x2^2*x6+11552*x2^2*x9-12768*x2*x3^2-6992*x2*x3*x9+13224*x2*x4^2+3800*x2*x5^2-7600*x2*x7*x9+3864*x3^3+4200*x3^2*x7+1058*x3^2*x9-4002*x3*x4^2+2628*x3*x4*x6+2844*x3*x5*x6-3332*x3*x6*x8+2300*x3*x7*x9-4760*x3*x8^2-4350*x4^2*x7-612*x4*x5*x6+11552*x5^2*x9+4508*x5*x6*x8+6440*x5*x8^2+5450*x6^2*x9+13720*x6*x8*x9+1250*x7^2*x9+9800*x8^2*x9+12768*x1*x5-10944*x2*x6+3312*x3*x6+13984*x5^2-13528*x5*x6+3600*x6*x7+468*x6];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
using Distributions: Dirichlet

export crossover,
    crossover!,
    CrossoverOperator,
    CrossoverResult,
    setup!

export ArithmeticCrossover,
    LiftedCrossover,
    NoCrossover,
    UniformCrossover

abstract type CrossoverOperator{G, K, P} end

setup!(operator::CrossoverOperator{G}, model::AbstractPopulationModel{G}) where {G} = operator

"""
    LiftedCrossover{G, C}(args)

Lifts a `CrossoverOperator{I, K, P}` on `I`, to a new `CrossoverOperator{G, K, P}` on `G`,
by storing `C(args)`.  The actual application/lifting needs to be defined manually.
"""
struct LiftedCrossover{G, C, I, K, P} <: CrossoverOperator{G, K, P}
    inner::C

    LiftedCrossover{G}(operator::C) where {G, I, K, P, C<:CrossoverOperator{I, K, P}} =
        new{G, C, I, K, P}(operator)
    LiftedCrossover{G, C}(args...) where {G, I, K, P, C<:CrossoverOperator{I, K, P}} =
        new{G, C, I, K, P}(C(args...))
end

"""
    crossover!(parents::Union{Family, NTuple}, operator, generation) -> children

Perform crossover between `parents`.  You only need to define the method for `parents` being an
`NTuple{K}`, returning an `NTuple{P}`.  The method will be automatically lifted to `Family`s (i.e.,
tuples of `Individual`s) of the same sizes.
"""
crossover!(parents::Family{G, K},
           operator::CrossoverOperator{G, K, P},
           generation::Int) where {G, K, P} =
    Individual.(crossover!(genome.(parents), operator, generation))

struct NoCrossover{N} <: CrossoverOperator{Any, N, N} end

crossover!(parents::NTuple{N}, operator::NoCrossover{N}, generation::Integer) where {N} =
    parents

struct ArithmeticCrossover{G<:AbstractVector, K, P} <: CrossoverOperator{G, K, P}
    rate::Rate
end

function crossover!(parents::NTuple{2, G},
                    operator::ArithmeticCrossover{G, 2, 2},
                    generation::Integer) where {G}
    if rand() < operator.rate(generation)
        mixing = rand()
        return ((1 - mixing) .* parents[1] .+ mixing .* parents[2],
                (1 - mixing) .* parents[2] .+ mixing .* parents[1])
    else
         return parents
    end
end

struct UniformCrossover{G<:AbstractVector, K, P} <: CrossoverOperator{G, K, P}
    p::Rate

    UniformCrossover{G, K, P}(p::Rate = ConstantRate(0.5)) where {G, K, P} = new{G, K, P}(p)
end

function crossover!(parents::NTuple{2, G},
                    operator::UniformCrossover{G, 2, 1},
                    generation::Integer) where {G}
    crossover_points = rand(length(parents[1])) .≤ operator.p(generation)
    (map(ifelse, crossover_points, parents...),)
end

function crossover!(parents::NTuple{2, G},
                    operator::UniformCrossover{G, 2, 2},
                    generation::Integer) where {G}
    crossover_points = rand(length(parents[1])) .≤ operator.p(generation)
    map(ifelse, crossover_points, parents...), map(ifelse, .~crossover_points, parents...)
end

# function crossover!(parents::Family{G, 1}, operator::ArithmeticCrossover{G, 1}) where {G}
#     if rand() < operator.rate
#         mixing = rand()
#         return ((1 - mixing) .* parents[1] .+ mixing .* parents[2],
#                 (1 - mixing) .* parents[2] .+ mixing .* parents[1])
#     else
#          return parents
#     end
# end

## TODO: move D to operator struct, use something different than "N random parent permutations"?
# function crossover(parents::NTuple{N, <:AbstractArray{G}}, operator::ArithmeticCrossover{G, N}) where {G}
#     D = Dirichlet(N, 1)

#     if rand() < operator.rate
#         mixing = rand(D)
#         return ntuple(sum(parents[randperm(N)] .* mixing), N)

#         return ((1 - mixing) .* parents[1] .+ mixing .* parents[2],
#                 (1 - mixing) .* parents[2] .+ mixing .* parents[1])
#     else
#          return parents
#     end
# end

# struct UniformCrossover{G} <: CrossoverOperator{NTuple{2, <:AbstractArray{G}}}
#     p::Float64
#     _dist::Bernoulli
# end

# function crossover!((p₁, p₂)::NTuple{2, <:AbstractArray{G}}, operator::UniformCrossover{G}) where {G}
#     @assert length(p₁) == length(p₂)
#     l = length(p₁)

#     crossover_points = rand(Uniform(operator.p), l)
#     p₁[crossover_points], p₂[crossover_points] = p₂[crossover_points], p₁[crossover_points]

#     p₁, p₂
# end


text:
@testset "Power series matrix inverse" begin

    T, t = Nemo.PowerSeriesRing(Nemo.GF(2^31 - 1), 50, "t"; model=:capped_absolute)

    for d in 1:5
        S = Nemo.MatrixSpace(T, d, d)
        for case in 1:20
            M = S([random_ps(T) for i in 1:d, j in 1:d])
            while isequal(StructuralIdentifiability.LinearAlgebra.det(StructuralIdentifiability.ps_matrix_const_term(M)), 0)
                M = S([random_ps(T) for i in 1:d, j in 1:d])
            end
            invM = ps_matrix_inv(M)
            prod = invM * M
            @test prod == one(S)
        end
    end
end


text:
################################################################################
# Primal-dual algorithms based on Asymmetric Forward-Backward-Adjoint
#
# [1] Latafat, Patrinos. "Asymmetric forward–backward–adjoint splitting for solving monotone inclusions involving three operators"  Computational Optimization and Applications, pages 1–37, 2017.
# [2] Condat. "A primal–dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms" Journal of Optimization Theory and Applications 158.2 (2013): 460-479.
# [3] Vũ. "A splitting algorithm for dual monotone inclusions involving cocoercive operators"" Advances in Computational Mathematics, 38(3), pp.667-681.

struct AFBAIterator{I <: Integer, R <: Real, T1 <: AbstractArray, T2 <: AbstractArray} <: ProximalAlgorithm{I, Tuple{T1, T2}}
    x::T1
    y::T2
    g
    hconj
    f
    lconj
    L
    theta::R
    mu::R
    lam::R
    betaQ::R # Lipschitz constant of gradient of f
    betaR::R # Lipschitz constant of gradient of l conjugate
    gamma1::R
    gamma2::R
    maxit::I
    tol::R
    verbose::I
    verbose_freq::I
    xbar::T1
    ybar::T2
    gradf::T1
    gradl::T2
    FPR_x::T1
    FPR_y::T2
    temp_x::T1
    temp_y::T2
end

################################################################################
# Constructor(s)

function AFBAIterator(x0::T1, y0::T2; g=IndFree(), h=IndFree(), f=IndFree(), l=IndZero(), L=Identity(size(x0)), theta=1, mu=1, lam=1, betaQ=0, betaR=0, gamma1=-1, gamma2=-1, maxit::I=10000, tol::R=1e-5, verbose=1, verbose_freq = 100) where {I, R, T1, T2}
    x = copy(x0)
    xbar = copy(x0)
    y = copy(y0)
    ybar = copy(y0)
    gradf = copy(x0)
    gradl = copy(y0)
    FPR_x = copy(x0)
    FPR_x .= Inf
    FPR_y = copy(y0)
    FPR_y .= Inf
    temp_x = copy(x0)
    temp_y = copy(y0)
    hconj = Conjugate(h)
    lconj = Conjugate(l)

    # default stepsizes
    par = 4 #  scale parameter for comparing Lipschitz constant and opnorm(L)
    nmL = opnorm(L)
    alpha = 1
    if isa(h, ProximalOperators.IndFree) && (gamma1<0 || gamma2<0)
        # mu=0 is the only case where stepsizes matter
        alpha = 1000/(betaQ+1e-5) # the speed is determined by gamma1 since bary=0
        temp = theta^2-3*theta+3
        gamma1 = 0.99/(betaQ/2+ temp*nmL/alpha) # in this case R=0
        gamma2 = 1/(nmL*alpha)
    elseif gamma1<0 || gamma2<0
        if theta==2 # default stepsize for Vu-Condat
            if betaQ > par*nmL && betaR > par*nmL
                alpha = 1
            elseif betaQ > par*nmL
                alpha = 2*nmL/betaQ
            elseif betaR > par*nmL
                alpha = betaR/(2*nmL)
            end
            gamma1 = 1/(betaQ/2+nmL/alpha)
            gamma2 = 0.99/(betaR/2+nmL*alpha)
        elseif theta==1 && mu==1 # default stepsize for theta=1, mu=1
            if betaQ > par*nmL && betaR > par*nmL
                alpha = 1
            elseif betaQ > par*nmL
                alpha = 2*nmL/betaQ
            elseif  betaR > par*nmL
                alpha = betaR/(2*nmL)
            end
            gamma1 = 1/(betaQ/2+nmL/alpha)
            gamma2 = 0.99/(betaR/2+gamma1*nmL^2)
        elseif theta==0 && mu==1 # default stepsize for theta=0, mu=1
            temp=3
            if betaQ > par*nmL && betaR > temp*par*nmL # denominator for Sigma involves 3α opnorm(L) in this case
                alpha = 1
            elseif betaQ > par*nmL
                alpha = 2*nmL/betaQ
            elseif betaR > temp*par*nmL # denominator for Sigma involves 3α opnorm(L) in this case
                alpha = betaR/(temp*2*nmL)
            end
            gamma1 = 1/(betaQ/2+nmL/alpha)
            gamma2 = 0.99/(betaR/2+2*gamma1*nmL^2+ alpha*nmL)
        elseif mu==0 # default stepsize for  mu=0
            temp = theta^2-3*theta+3
            if betaQ > temp*par*nmL && betaR > par*nmL
                alpha = 1
            elseif betaR > par*nmL
                alpha = betaR/(2*nmL)
            elseif betaQ > temp*par*nmL # denominator for Sigma involves 3α opnorm(L) in this case
                alpha = 2*nmL*temp/betaQ
            end
            gamma2 = 1/(betaR/2+ alpha*nmL)
            gamma1 = 0.99/(betaQ/2+ (temp-1)*gamma2*nmL^2+ nmL/alpha)
        else
            error("this choice of theta and mu is not supported!")
        end
    end
    return AFBAIterator{I, R, T1, T2}(x0, y0, g, hconj, f, lconj, L, theta, mu, lam, betaQ, betaR, gamma1, gamma2, maxit, tol, verbose, verbose_freq, xbar, ybar, gradf, gradl, FPR_x,FPR_y, temp_x, temp_y)
end

################################################################################
# Utility methods
maxit(sol::AFBAIterator) = sol.maxit

converged(sol::AFBAIterator, it) = it > 0 && norm(sol.FPR_x)+norm(sol.FPR_y) <= sol.tol

verbose(sol::AFBAIterator) = sol.verbose > 0
verbose(sol::AFBAIterator, it) = sol.verbose > 0 && (sol.verbose == 2 ? true : (it == 1 || it%sol.verbose_freq == 0))

function display(sol::AFBAIterator)
    @printf("%6s | %22s | %10s |\n ", "it", "gamma", "fpr")
    @printf("------|------------------------|------------|\n")
end

function display(sol::AFBAIterator, it)
    @printf("%6d | %7.4e, %7.4e | %7.4e |\n", it, sol.gamma1, sol.gamma2, norm(sol.FPR_x)+norm(sol.FPR_y))
end

function Base.show(io::IO, sol::AFBAIterator)
    println(io, "Asymmetric Forward-Backward-Adjoint Splitting" )
    if sol.theta==2
        println(io, "theta               : $(sol.theta)" )
    else
        println(io, "theta, mu           : $(sol.theta), $(sol.mu)" )
    end
    println(io, "fpr                 : $(norm(sol.FPR_x)+norm(sol.FPR_y))")
    print(  io, "gamma1, gamma2      : $(sol.gamma1), $(sol.gamma2)")
end

################################################################################
# Initialization

function initialize!(sol::AFBAIterator)
    return sol.x, sol.y
end

################################################################################
# Iteration

function iterate!(sol::AFBAIterator{I, R, T1, T2}, it::I) where {I, R, T1, T2}
    # perform xbar-update step
    gradient!(sol.gradf, sol.f, sol.x)
    mul!(sol.temp_x, sol.L', sol.y)
    sol.temp_x .+= sol.gradf
    sol.temp_x .*= -sol.gamma1
    sol.temp_x .+= sol.x
    prox!(sol.xbar, sol.g, sol.temp_x, sol.gamma1)
    # perform ybar-update step
    gradient!(sol.gradl, sol.lconj, sol.y)
    sol.temp_x .=  (sol.theta * sol.xbar) .+ ((1-sol.theta) * sol.x)
    mul!(sol.temp_y, sol.L, sol.temp_x)
    sol.temp_y .-= sol.gradl
    sol.temp_y .*= sol.gamma2
    sol.temp_y .+= sol.y
    prox!(sol.ybar, sol.hconj, sol.temp_y, sol.gamma2)
    # the residues
    sol.FPR_x .= sol.xbar .- sol.x
    sol.FPR_y .= sol.ybar .- sol.y
    # perform x-update step
    sol.temp_y .= sol.mu*(2-sol.theta)*sol.gamma1*sol.FPR_y
    mul!(sol.temp_x, sol.L', sol.temp_y)
    sol.x .+= sol.lam *(sol.FPR_x .- sol.temp_x)
    # perform y-update step
    sol.temp_x .= (1-sol.mu)*(2-sol.theta)*sol.gamma2*sol.FPR_x
    mul!(sol.temp_y, sol.L, sol.temp_x)
    sol.y .+= sol.lam *(sol.FPR_y .+ sol.temp_y)
    return sol.x, sol.y
end

################################################################################
# Solver interface(s)

"""
**Asymmetric forward-backward-adjoint algorithm**

    AFBA(x0, y0; kwargs)

Solves convex optimization problems of the form

    minimize f(x) + g(x) + (h □ l)(L x),

where `f` is smooth, `g` and `h` are possibly nonsmooth and `l` is strongly convex.
Symbol `□` denotes the infimal convolution, and `L` is a linear mapping.
Points `x0` and `y0` are the initial primal and dual iterates, respectively.

Keyword arguments are as follows:
* `f`: smooth, convex function (default: zero)
* `g`: convex function (possibly nonsmooth, default: zero)
* `h`: convex function (possibly nonsmooth, default: zero)
* `l`: strongly convex function (possibly nonsmooth, default: indicator of {0})
* `L`: linear operator (default: identity)
* `betaQ`: Lipschitz constant of gradient of f (default: zero)
* `betaR`: Lipschitz constant of gradient of l conjugate (default: zero)
* `gamma1`: stepsize corresponding to the primal updates (default: see [1] for each case)
* `gamma2`: stepsize corresponding to the dual updates (default: see [1] for each case)
* `theta`: nonnegative algorithm parameter (default: '1')
* `mu`: algorithm parameter in the range [0,1] (default: '1')
* `tol`: primal-dual termination tolerance (default: `1e-5`)
* `maxit`: maximum number of iterations (default: `10000`)
* `verbose`, verbosity level (default: `1`)
* `verbose_freq`, verbosity frequency for `verbose = 1` (default: `100`)

The iterator implements Algorithm 3 of [1] with constant stepsize (α_n=λ) for several prominant special cases:
1) θ = 2          ==>   Corresponds to the Vu-Condat Algorithm [2,3].
2) θ = 1, μ=1
3) θ = 0, μ=1
4) θ ∈ [0,∞), μ=0

See [1, Figure 1] for other special cases and relation to other algorithms.

[1] Latafat, Patrinos. "Asymmetric forward–backward–adjoint splitting for solving monotone inclusions involving three operators"  Computational Optimization and Applications, pages 1–37, 2017.
[2] Condat. "A primal–dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms" Journal of Optimization Theory and Applications 158.2 (2013): 460-479.
[3] Vũ. "A splitting algorithm for dual monotone inclusions involving cocoercive operators"" Advances in Computational Mathematics, 38(3), pp.667-681.
"""
function AFBA(x0, y0; kwargs...)
    # Create iterable
    sol = AFBAIterator(x0, y0; kwargs...)
    # Run iterations
    (it, (primal, dual)) = run!(sol)
    return (it, primal, dual, sol)
end

"""
**Vũ-Condat primal-dual algorithm**

    VuCondat(x0, y0; kwargs)

Solves convex optimization problems of the form

    minimize f(x) + g(x) + (h □ l)(L x).

where `f` is smooth, `g` and `h` are possibly nonsmooth and `l` is strongly convex.

Symbol `□` denotes the infimal convolution, and `L` is a linear mapping.
Points `x0` and `y0` are the initial primal and dual iterates, respectively.

See documentation of `AFBA` for the list of keyword arguments.
"""
function VuCondat(x0, y0; kwargs...)
    # Create iterable
    sol = AFBAIterator(x0, y0; kwargs..., theta=2)
    # Run iterations
    (it, (primal, dual)) = run!(sol)
    return (it, primal, dual, sol)
end

"""
**Chambolle-Pock primal-dual algorithm**

    ChambollePock(x0, y0; kwargs)

Solves convex optimization problems of the form

    minimize g(x) + (h □ l)(L x).

where `g` and `h` are possibly nonsmooth and `l` is strongly convex.
Symbol `□` denotes the infimal convolution, and `L` is a linear mapping.
Points `x0` and `y0` are the initial primal and dual iterates, respectively.

See documentation of `AFBA` for the list of keyword arguments.
"""
function ChambollePock(x0, y0; kwargs...)
    # Create iterable
    sol = AFBAIterator(x0, y0; kwargs..., f=IndFree(), theta=2)
    # Run iterations
    (it, (primal, dual)) = run!(sol)
    return (it, primal, dual, sol)
end


text:
##############################
## by Qin Yu, Nov 2018
## using Julia 1.0.1
##############################

############################## Starting here:
# If this gives error, please call:
# using Pkg # Julia 1.0.1 only
# Pkg.add("Plots")
# Pkg.add("GR")

using LinearAlgebra
using Plots
using Printf

############################## Linear Regression with Polynomial Basis on Simple Dataset:
# Prepare Data
# Let's try to see the effect of polynomial basis,
# on this extremely simple dataset, S:
S = [(1,3),(2,2),(3,0),(4,5)]
plot(S, line=:scatter, lab="data set")

x = [x[1] for x in S]
y = [y[2] for y in S]

# Formulae - Find 𝝎:
phik(x, k) = [xi^k for xi in x]  # ϕₖ(𝒙) = ..., for vector 𝒙, all inputs
phi1tok(x, k) = [phik(x, i) for i in 0:k-1]  # Φ(𝒙) = 𝒙 ⋅ 𝝋 = 𝒙 ⋅ (ϕ₁(), ϕ₂(), ..., ϕₖ())
transformed_x_kk(x, k) = hcat(phi1tok(x, k)...)  # Φ(𝒙) as a matrix
w_k(x, y, k) = transformed_x_kk(x, k) \ y  # 𝝎 = Φ\𝒚

# Formulae - Equation of Fitted Regression Line:
w_phi_k(x_test, x_train, y, k) = dot(transformed_x_kk(x_test, k), w_k(x_train, y, k))  # ̂𝑦 = Φ(𝑥) ⋅ 𝝎

plot(S, line=:scatter, lab="data set", legend=:bottomright)
W = zeros(4, 4)
for i = 1:4
    𝒘Φ(x_test) = w_phi_k(x_test, x, y, i)
    display(plot!(𝒘Φ, 0, 4, lab="k = $i"))
    W[i,1:i] = w_k(x, y, i)

    # print the formulae for these curves
    print("k = $i, f(x) = ")
    for j = 1:i
        j > 1 && print(" + ")
        @printf "%.2f" W[i,j]
        j > 1 && print("x^$(j-1)")
    end
    print("\n")
end
savefig("./graph/1.1.pdf")

############################## Training Error
# SSE = 𝚺ᵢ(𝑦ᵢ - ̂𝑦ᵢ)² = 𝚺ᵢ(𝑦ᵢ - Φ(𝑥ᵢ) ⋅ 𝝎)²
# MSE = SSE/N, where N = number_of_rows(input_data_set), here is S
# (here I use N as the book ESLII uses, standing for m in qestions)
sse_k(x, y, k) = sum((y - transformed_x_kk(x, k) * w_k(x, y, k)).^2)
mse_k(x, y, k) = sse_k(x, y, k) / first(size(y))

MSE = [mse_k(x, y, i) for i = 1:4]
plot(MSE, xlabel="k", lab="MSE")
savefig("./graph/1.2.pdf")

max_stars_repo_name: qin-yu/julia-regression-boston-housing

max_issues_repo_name: qin-yu/ml-julia-boston-housing

max_forks_repo_name: qin-yu/ml-julia-boston-housing


text:
import Base: sin, cos, tanh, exp, log

"""
    Sin <: Func
"""
@func mutable struct Sin end

foward(f::Sin, x) = sin.(x)

backward(f::Sin, gy) = cos.(gy)

sin(x::Var) = Sin()(x)

"""
    Cos <: Func
"""
@func mutable struct Cos end

foward(f::Cos, x) = cos.(x)

backward(f::Cos, gy) = -sin.(gy)

cos(x::Var) = Cos()(x)

"""
    Tanh <: Func
"""
@func mutable struct Tanh end

forward(f::Tanh, x) = tanh.(x)

function backward(f::Tanh, gy)
    y = f.outputs[1].value
    return gy .* (1 - y .* y)
end

tanh(x::Var) = Tanh()(x)

"""
    Exp <: Func
"""
@func mutable struct Exp end

forward(f::Exp, x) = exp.(x)

backward(f::Exp, gy) = gy .* f.outputs[1].value

exp(x::Var) = Exp()(x)

"""
    Log <: Func
"""
@func mutable struct Log end

forward(f::Log, x) = log.(x)

function backward(f::Log, gy)
    x, = f.args
    return gy / x
end

log(x::Var) = Log()(x)


text:
@testset "test_freqresp" begin
## EVALFR ##
H = [tf(0) tf([3, 0],[1, 1, 10]) ; tf([1, 1],[1, 5]) tf([2],[1, 6])]
G = ss([-5 0 0 0; 0 -1 -2.5 0; 0 4 0 0; 0 0 0 -6], [2 0; 0 1; 0 0; 0 2],
       [0 3 0 0; -2 0 0 1], [0 0; 1 0])

@test evalfr(H, -6) == [0.0 -0.45; 5.0 Inf]
@test evalfr(H, -5) == [0.0 -0.5; Inf 2.0]
@test evalfr(H, -1) == [0.0 -0.3; 0.0 0.4]
@test evalfr(H, 0) ≈ [0.0 0.0; 0.2 1/3]

@test evalfr(G, -6) == [Inf Inf; Inf Inf]
@test evalfr(G, -5) == [Inf Inf; Inf Inf]
@test evalfr(G, -1) ≈ [0.0 -0.3; 0.0 0.4]
@test evalfr(G, 0) ≈ [0.0 0.0; 0.2 1/3]

## Constant system
w = exp10.(range(-1, stop=1, length=50))

sys1 = ss(1)
G1 = tf(1)
H1 = zpk(1)
resp1 = ones(ComplexF64, length(w), 1, 1)

@test evalfr(sys1, im*w[1]) == fill(resp1[1], 1, 1)
@test evalfr(G1, im*w[1]) == fill(resp1[1], 1, 1)
@test evalfr(H1, im*w[1]) == fill(resp1[1], 1, 1)

@test freqresp(sys1, w) == resp1
@test freqresp(G1, w) == resp1
@test freqresp(H1, w) == resp1

## First order system

sys2 = ss(-1, 1, 1, 1)
G2 = tf([1, 2], [1,1])
H2 = zpk([-2], [-1.0], 1.0)
resp2 = reshape((im*w .+ 2)./(im*w  .+ 1), length(w), 1, 1)

@test evalfr(sys2, im*w[1]) ≈ fill(resp2[1], 1, 1)
@test evalfr(G2, im*w[1]) == fill(resp2[1], 1, 1)
@test evalfr(H2, im*w[1]) == fill(resp2[1], 1, 1)

@test freqresp(sys2, w) ≈ resp2 rtol=1e-15
@test freqresp(G2, w) == resp2
@test freqresp(H2, w) == resp2

## Complex-coefficient system
sys3 = ss(-1+im, 1, (1-im), (1-im))
G3 = (1-im)*tf([1,2-im], [1,1-im])
H3 = zpk([-2+im], [-1+im], (1-im))
resp3 = reshape((1 .- im)*(2 .+ im*(w .- 1))./(1 .+ im*(w .- 1)), length(w), 1, 1)

@test evalfr(sys3, im*w[1]) ≈ fill(resp3[1], 1, 1) rtol=1e-16
@test evalfr(G3, im*w[1]) == fill(resp3[1], 1, 1)
@test evalfr(H3, im*w[1]) == fill(resp3[1], 1, 1)

@test freqresp(sys3, w) ≈ resp3 rtol=1e-15
@test freqresp(G3, w) ≈ resp3 rtol=1e-16
@test freqresp(H3, w) == resp3

## Shortcut notation for evalfr ##
F = tf([1],[1,0.5],-1)
omega = 2
z = 0.5(1+im)
# This test is not correct if we dont have a sample time
# @test_throws F(omega,true)[1] == 1/(exp(-im*2)+0.5)
@test F(z,false)[1] == 1/(z+0.5)
@test_throws ErrorException F(z,true)

F = [tf([1],[1,0.5],2.0) 3*tf([1],[1,0.5],2.0)]
omegas = [1,2]
z = 0.5(1+im)
@test F(omegas[1],true) ≈ [1 3].*1/(exp(im*2)+0.5) atol=1e-14
@test F(omegas[2],true) == [1 3].*1/(exp(2*im*2)+0.5)
@test F(omegas,true) ≈ [k/(exp(omega*im*2)+0.5) for omega=omegas, o=1:1, k=[1,3]] atol=1e-14
@test F(omegas,false) ≈ [k/(omega+0.5) for omega=omegas, o=1:1, k=[1,3]] atol=1e-14
@test F(z,false)[1] == 1/(z+0.5)
@test_throws ErrorException F(z,true)

## Test bode, nyquist and sigma
sys = [tf([1,-1], [1,1,1]) 0; 0 tf([1],[1,1])]
f(s) = [(s-1)./(s.^2+s+1) 0; 0 1 ./(1+s)]
ws = exp10.(range(-2, stop=2, length=50))
resp = Array{ComplexF64}(undef, 50,2,2)
for (i,w) in enumerate(ws)
    resp[i,:,:] = f(im*w)
end

@test bode(sys, ws)[1:2] == (abs.(resp), rad2deg.(angle.(resp)))
@test nyquist(sys, ws)[1:2] == (real(resp), imag(resp))
sigs = Array{Float64}(undef, 50,2)
for i in eachindex(ws)
    sigs[i,:] =  svdvals(resp[i,:,:])
end
@test sigma(sys, ws)[1] == sigs

#Test default freq vector contains at least half a decade more than all features
p, z = 100, 1/1000
sys2 = [tf([1],[1/p,1]) tf([1/z, 2/z, 1],[1])]
mag, phase, ws2 = bode(sys2)
@test maximum(ws2) >= 2max(p,z)
@test minimum(ws2) <= 0.5min(p,z)

#Test default freq vector not too small
p, z = 1, 1.2
sys2 = [tf([1],[1/p,1]) tf([1/z, 2/z, 1],[1])]
mag, mag, ws2 = bode(sys2)
@test maximum(ws2) >= 5max(p,z)
@test minimum(ws2) <= 0.2min(p,z)
@test length(ws2) > 100
end


text:
#=
# Type definition for a FeatureVector.
# Wrapper around a Dict type.
# Restricted to Any type [key] => [value] of Number type.
=#

type FeatureVector{K,V<:Number} <: FeatureSpace
  map::Dict{K,V}
  total::Number
  length::Number
  mdata::Any
  FeatureVector() = new(Dict{K,V}(),0,0)
  FeatureVector{K,V}(map::Dict{K,V}) = new(map,get_total(map),get_length(map))
  function get_total(fv::Dict{K,V})
    total = 0
    for value in values(fv)
      total += value
    end
    return total
  end
  function get_length(fv::Dict{K,V})
    length = 0
    for value in values(fv)
      length += value^2
    end
    return sqrt(length)
  end
end
FeatureVector() = FeatureVector{Any,Number}()
FeatureVector{K,V}(map::Dict{K,V}) = FeatureVector{K,V}(sanitize!(Base.copy(map)))

# gets value of [key] in a FeatureVector
function getindex(fv::FeatureVector, key)
  if haskey(fv,key)
    return fv.map[key]
  end
  return 0
end

# sets value of [key] in a FeatureVector. Must be subtype of number/dict type
function setindex!(fv::FeatureVector, value, key)
  fv.total -= fv[key]
  fv.total += value

  fv.length = sqrt((fv.length^2-fv[key]^2)+value^2)

  if value == 0
    Base.delete!(fv.map, key)
  else
    fv.map[key] = value
  end
end

# gets all keys of a FeatureVector
function keys(fv::FeatureVector)
  return Base.keys(fv.map)
end

# gets all values of a FeaturVector
function values(fv::FeatureVector)
  return Base.values(fv.map)
end

# check to see if the FeatureVector is empty.
function isempty(fv::FeatureVector)
  return Base.isempty(fv.map)
end

# check to see if the FeatureVector has key.
function haskey(fv::FeatureVector, key)
  return Base.haskey(fv.map,key)
end

# returns length of FeatureVector
function length(fv::FeatureVector)
  return Base.length(fv.map)
end

# copies selected fv, and makes a new one.
function copy{K,V}(fv::FeatureVector{K,V})
  new_fv = FeatureVector{K,V}()

  for key in keys(fv)
    new_fv[key] = fv[key]
  end

  return new_fv
end

# deletes 0 value keys in a dictionary
function sanitize!(dict::Dict)
  for key in keys(dict)
    if dict[key] == 0
      Base.delete!(dict,key)
    end
  end
  return dict
end

# returns an Array of the (key, value) tuples
# in fv using the provided boolean expression
# by default it sorts by largest value
function freq_list(fv::FeatureVector, expression::Function = (a,b) -> a[2]>b[2])
  words = Array(find_common_type(fv,fv),length(fv))
  i = 1
  for word in fv.map
    words[i] = word
    i+=1
  end
  return Base.sort!(words,lt=expression)
end

# finds common type of two FeatureVectors
function find_common_type(fv1::FeatureVector,fv2::FeatureVector)
  if isempty(fv1) && isempty(fv2)
    commonType = (Any,Number)
  elseif isempty(fv1)
    commonType = typeof(first(fv2.map))
  elseif isempty(fv2)
    commonType = typeof(first(fv1.map))
  else
    fv1_type = typeof(first(fv1.map))
    fv2_type = typeof(first(fv2.map))
    commonType = (promote_type(fv1_type[1],fv2_type[1]),
    promote_type(fv1_type[2],fv2_type[2]))
  end

  return commonType
end

# adds two FeatureVectors together
function +(fv1::FeatureVector, fv2::FeatureVector)
  dict_type = find_common_type(fv1,fv2)
  dict = Dict{dict_type[1],dict_type[2]}(fv1.map)
  fv1_keys = keys(fv1)
  fv2_keys = keys(fv2)

  for key in fv2_keys
    if key in fv1_keys
      dict[key] += fv2[key]
    else
      dict[key] = fv2[key]
    end
  end

  return FeatureVector(dict)
end

# adds two FeatureVectors together in place
function add!(fv1::FeatureVector, fv2::FeatureVector)
  for key in keys(fv2)
      fv1[key] += fv2[key]
  end
end

# subtracts two FeatureVectors
function -(fv1::FeatureVector, fv2::FeatureVector)
  dict_type = find_common_type(fv1,fv2)
  dict = Dict{dict_type[1],dict_type[2]}(fv1.map)
  fv1_keys = keys(fv1)
  fv2_keys = keys(fv2)

  for key in fv2_keys
    if key in fv1_keys
      dict[key] -= fv2[key]
    else
      dict[key] = -fv2[key]
    end
  end

  return FeatureVector(dict)
end

# subtracts two FeatureVectors in place
function subtract!(fv1::FeatureVector, fv2::FeatureVector)
  for key in keys(fv2)
      fv1[key] -= fv2[key]
  end
end

# multiplies a FeatureVector by a scalar
function *(fv::FeatureVector, value::Number)
  if isempty(fv)
    return fv
  end

  fv_keys = keys(fv)
  fv_type = typeof(first(fv.map)[2])
  dict = Dict{typeof(first(fv_keys)), promote_type(fv_type,typeof(value))}()

  for key in fv_keys
    dict[key] = fv[key]*value
  end

  return FeatureVector(dict)
end

# handles commutativity of multiplication
function *(value::Number, fv::FeatureVector)
  return fv*value
end

# in place multiplication by a scalar
function multiply!(fv::FeatureVector, value::Number)
  for key in keys(fv)
    fv[key] = fv[key]*value
  end
end

# divides a FeatureVector by a scalar
function /(fv::FeatureVector, value::Number)
  if isempty(fv)
    return fv
  end

  fv_keys = keys(fv)
  dict = Dict{Any, typeof(fv[first(fv_keys)]/value)}()

  for key in fv_keys
    dict[key] = fv[key]/value
  end

  return FeatureVector(dict)
end

# in place division by a scalar
function divide!(fv::FeatureVector, value::Number)
  for key in keys(fv)
    fv[key] = fv[key]/value
  end
end

# rationalizes a FeatureVectors values
function //(fv::FeatureVector, value::Integer)
  if isempty(fv)
    return fv
  end

  fv_keys = keys(fv)
  dict = Dict{typeof(first(fv_keys)), Rational}()
  for key in fv_keys
    dict[key] = fv[key]//value
  end

  return FeatureVector(dict)
end

# in place rationalization by a scalar
function rationalize!(fv::FeatureVector, value::Integer)
  for key in keys(fv)
    fv[key] = fv[key]//value
  end
end

# finds the cosine between two vectors and returns 1-cos
function dist_cos(fv1::FeatureVector, fv2::FeatureVector, norm::Bool=false)
  fv1_magnitude = 0
  fv2_magnitude = 0
  dot_product = 0

  for key in keys(fv1)
    fv1_value = fv1[key]
    fv1_magnitude += fv1_value*fv1_value
    dot_product += fv1_value*fv2[key]
  end

  for key in keys(fv2)
    fv2_value = fv2[key]
    fv2_magnitude += fv2_value*fv2_value
  end

  cosine = dot_product/(sqrt(fv1_magnitude)*sqrt(fv2_magnitude))
  if cosine > 1 - 1e-15
    cosine = 1
  end

  return 1 - cosine
end

# number of disjoint nonzero dimensions between vectors
function dist_zero(fv1::FeatureVector, fv2::FeatureVector, norm::Bool=false)
  distance = 0

  for key in keys(fv1)
    if fv1[key] != 0 && fv2[key] == 0
      distance += 1
    end
  end

  for key in keys(fv2)
    if fv2[key] != 0 && fv1[key] == 0
      distance += 1
    end
  end

  return distance
end

# zero distance divided by number of dimensions
function dist_zero_weighted(fv1::FeatureVector, fv2::FeatureVector, norm::Bool=false)
  d1 = 0
  d2 = 0
  l1 = 0
  l2 = 0

  for key in keys(fv1)
    if fv1[key] != 0
      l1 += 1
      if fv2[key] == 0
	d1 += 1
      end
    end
  end

  for key in keys(fv2)
    if fv2[key] != 0
      l2 += 1
      if fv1[key] == 0
	d2 += 1
      end
    end
  end

  return d1/l1 + d2/l2
end

# sum of absolute distance between dimensions
function dist_taxicab(fv1::FeatureVector, fv2::FeatureVector, norm::Bool = true)
  distance = 0

  for key in keys(fv1)
    if norm
      distance += abs((fv1[key]/fv1.total)-(fv2[key]/fv2.total))
    else
      distance += abs(fv1[key]-fv2[key])
    end
  end

  for key in keys(fv2)
    if !haskey(fv1, key)
      if norm
        distance += abs(fv2[key]/fv2.total)
      else
        distance += abs(fv2[key])
      end
    end
  end

  return distance
end

# ordinary distance between vectors
function dist_euclidean(fv1::FeatureVector, fv2::FeatureVector, norm::Bool = true)
  distance = 0

  for key in keys(fv1)
    if norm
      distance += ((fv1[key]/fv1.length)-(fv2[key]/fv2.length))^2
    else
      distance += (fv1[key]-fv2[key])^2
    end
  end

  for key in keys(fv2)
    if !haskey(fv1, key)
      if norm
        distance += (fv2[key]/fv2.length)^2
      else
        distance += fv2[key]^2
      end
    end
  end

  return sqrt(distance)
end

# maximum absolute difference between dimensions
function dist_infinite(fv1::FeatureVector, fv2::FeatureVector, norm::Bool = true)
  max_dist = 0
  max_fv1 = 0
  max_fv2
  if norm
    for value in values(fv1)
      if value > max_fv1
        max_fv1 = value
      end
    end
    for value in values(fv2)
      if value > max_fv2
        max_fv2 = value
      end
    end
  end

  for key in keys(fv1)
    if norm
      current_dist = abs((fv1[key]/max_fv1)-(fv2[key]/max_fv2))
    else
      current_dist = abs(fv1[key]-fv2[key])
    end

    if current_dist > max_dist
      max_dist = current_dist
    end
  end

  for key in keys(fv2)
    current_dist = 0
    if !haskey(fv1, key)
      if norm
        current_dist = abs(fv2[key]/max_fv2)
      else
        current_dist = abs(fv2[key])
      end
    end
    if current_dist > max_dist
      max_dist = current_dist
    end
  end

  return max_dist
end

# prints out the fv cleanly.
function display(fv::FeatureVector)
  print_with_color(:white,string(typeof(fv)))
  fv_list = freq_list(fv)
  k = length(fv_list)
  print_with_color(:white," with $k features")
  if k > 0
    print_with_color(:white,":\n")
    if k > 5
      k = 5
    end

    print_with_color(:white,"  Top $k Features:\n")
    for key in fv_list[1:k]
      print_with_color(:white,"    $key\n")
    end
  end
end

function show(io::IO, fv::FeatureVector)
  print(io,string(typeof(fv)))
  k = length(fv.map)
  print(io," with $k features")
end


text:
"""
$(SIGNATURES)

Given `{{ ... }}` blocks, identify what kind of blocks they are and return a
vector of qualified blocks of type `AbstractBlock`.
"""
function qualify_html_hblocks(blocks::Vector{OCBlock})::Vector{AbstractBlock}

    qb = Vector{AbstractBlock}(undef, length(blocks))

    # TODO improve this (if there are many blocks this would be slow)
    for (i, β) ∈ enumerate(blocks)
        # if block {{ if v }}
        m = match(HBLOCK_IF_PAT, β.ss)
        isnothing(m) || (qb[i] = HIf(β.ss, m.captures[1]); continue)
        # else block {{ else }}
        m = match(HBLOCK_ELSE_PAT, β.ss)
        isnothing(m) || (qb[i] = HElse(β.ss); continue)
        # else if block {{ elseif v }}
        m = match(HBLOCK_ELSEIF_PAT, β.ss)
        isnothing(m) || (qb[i] = HElseIf(β.ss, m.captures[1]); continue)
        # end block {{ end }}
        m = match(HBLOCK_END_PAT, β.ss)
        isnothing(m) || (qb[i] = HEnd(β.ss); continue)
        # ---
        # isdef block
        m = match(HBLOCK_ISDEF_PAT, β.ss)
        isnothing(m) || (qb[i] = HIsDef(β.ss, m.captures[1]); continue)
        # isnotdef block
        m = match(HBLOCK_ISNOTDEF_PAT, β.ss)
        isnothing(m) || (qb[i] = HIsNotDef(β.ss, m.captures[1]); continue)
        # ---
        # isempty block
        m = match(HBLOCK_ISEMPTY_PAT, β.ss)
        isnothing(m) || (qb[i] = HIsEmpty(β.ss, m.captures[1]); continue)
        # isnotempty block
        m = match(HBLOCK_ISNOTEMPTY_PAT, β.ss)
        isnothing(m) || (qb[i] = HIsNotEmpty(β.ss, m.captures[1]); continue)
        # ---
        # ispage block
        m = match(HBLOCK_ISPAGE_PAT, β.ss)
        isnothing(m) || (qb[i] = HIsPage(β.ss, split(m.captures[1])); continue)
        # isnotpage block
        m = match(HBLOCK_ISNOTPAGE_PAT, β.ss)
        isnothing(m) || (qb[i] = HIsNotPage(β.ss, split(m.captures[1])); continue)
        # ---
        # for block
        m = match(HBLOCK_FOR_PAT, β.ss)
        if !isnothing(m)
            v, iter = m.captures
            check_for_pat(v)
            qb[i] = HFor(β.ss, v, iter);
            continue
        end
        # ---
        # function block {{ fname v1 v2 ... }}
        m = match(HBLOCK_FUN_PAT, β.ss)
        if !isnothing(m)
            if isnothing(m.captures[2]) || isempty(strip(m.captures[2]))
                ps = String[]
            else
                # merge \"abc def\" (issue #731)
                ps = split_hfun_parameters(m.captures[2])
            end
            qb[i] = HFun(β.ss, m.captures[1], ps)
            continue
        end
        # ---
        throw(HTMLBlockError("I found a HBlock that did not match anything, " *
                             "verify '$(β.ss)'"))
    end
    return qb
end

"""Blocks that can open a conditional block which are special."""
const HTML_OPEN_COND_SP =
    Union{HIsDef, HIsNotDef, HIsEmpty, HIsNotEmpty, HIsPage, HIsNotPage}

"""Blocks that can open a conditional block."""
const HTML_OPEN_COND = Union{HIf, HTML_OPEN_COND_SP}

"""
$SIGNATURES

Internal function to balance conditional blocks. See [`process_html_qblocks`](@ref).
"""
hbalance(::HTML_OPEN_COND) = 1
hbalance(::HEnd) = -1
hbalance(::AbstractBlock) = 0

function split_hfun_parameters(s)
    matches = [e.match for e in eachmatch(HFUN_PARAMS_PAT, s; overlap=false)]
    return strip.(matches, '\"')
end


text:
#This is the place tu put all the variables needed for the project
#All the variable will be sent from here to required places in the programme

##Parameters which are needed for setting up the plant and Differential Equations
function Model_Par(m1,m2,m3,m4,m5)
    myDict = Dict(
        :x2_f   => m1,
        :km     => m2,
        :k1     => m3,
        :Y      => m4,
        :μmax   => m5
    )
    return myDict
end

## Parameters which are needed to set up the simulation of the plant.

function Simulation_Data(Tfsim, TfMPC, Ncp, states, inputs)
    T0_sim  = 0.0;
    Tf_sim  = Tfsim;
    dt_sim  = 1;
    Tf_MPC  = TfMPC;
    dt_MPC  = 1;
    NCP     = Ncp;
    xk      = states;
    uk      = inputs;

    NFE_sim = convert(Int, (Tf_sim - T0_sim) / dt_sim);
    NFE_MPC = convert(Int, (Tf_MPC) / dt_MPC);
    xₛₚ     = hcat(1.5032 * ones(1, convert(Int, NFE_sim / 3)),     0.9951 * ones(1,convert(Int, NFE_sim / 3)),    0 * ones(1,convert(Int, NFE_sim / 3)), 0 * ones(1,convert(Int, NFE_MPC)));
    myDict = Dict(
        :T0_sim  => T0_sim,
        :Tf_sim  => Tf_sim,
        :dt_sim  => dt_sim,
        :Tf_MPC  => Tf_MPC,
        :dt_MPC  => dt_MPC,
        :NCP     => NCP,
        :NFE_sim => NFE_sim,
        :NFE_MPC => NFE_MPC,
        :xₛₚ     => xₛₚ,
        :xk      => xk,
        :uk      => uk
    )
    return myDict
end


text:
using Weave
using Test
import Plots

function publish_test(outfile, format)
  outfile = joinpath("documents/publish", outfile)
  infile = "documents/publish_test.jmd"
  weave(infile, doctype = format, out_path = outfile, template = "templates/mini.tpl")
  result =  read(outfile, String)
  ref =  read(outfile * ".ref", String)
  @test result == ref
  rm(outfile)
end

#Test formatters
publish_test("publish_tex.tex", "md2tex")
#!is_windows() && publish_test("publish_test.html", "md2html")


text:
export HostCall, hostcall!

const SENTINEL_COUNTER = Ref{Int}(2)
const SENTINEL_LOCK = ReentrantLock()
const DEFAULT_HOSTCALL_LATENCY = 0.01

"""
    HostCall{S,RT,AT}

GPU-compatible struct for making hostcalls.
"""
struct HostCall{S,RT,AT}
    signal::S
    host_sentinel::Int
    device_sentinel::Int
    buf_ptr::LLVMPtr{UInt8,AS.Global}
    buf_len::UInt
end
function HostCall(RT::Type, AT::Type{<:Tuple}, signal::S;
                  agent=get_default_agent(), buf_len=nothing) where {S<:UInt64}
    # FIXME: atomic_add!
    host_sentinel, device_sentinel = lock(SENTINEL_LOCK) do
        host_sentinel = SENTINEL_COUNTER[]
        SENTINEL_COUNTER[] += 1
        device_sentinel = SENTINEL_COUNTER[]
        SENTINEL_COUNTER[] += 1
        if SENTINEL_COUNTER[] > typemax(Int64)-2
            SENTINEL_COUNTER[] = 2
            @debug "Sentinel counter wrapped around!"
        end
        return host_sentinel, device_sentinel
    end
    if buf_len === nothing
        buf_len = 0
        for T in AT.parameters
            @assert isbitstype(T) "Hostcall arguments must be bits-type"
            buf_len += sizeof(T)
        end
    end
    buf_len = max(sizeof(UInt64), buf_len) # make room for return buffer pointer
    buf = Mem.alloc(agent, buf_len; coherent=true)
    buf_ptr = LLVMPtr{UInt8,AS.Global}(Base.unsafe_convert(Ptr{UInt8}, buf))
    HostCall{S,RT,AT}(signal, host_sentinel, device_sentinel, buf_ptr, buf_len)
end

## device signal functions
# TODO: device_signal_load, device_signal_add!, etc.
@inline device_signal_store!(signal::HSA.Signal, value::Int64) =
    device_signal_store!(signal.handle, value)
@inline @generated function device_signal_store!(signal::UInt64, value::Int64)
    Context() do ctx
        T_nothing = convert(LLVMType, Nothing; ctx)
        T_i32 = LLVM.Int32Type(ctx)
        T_i64 = LLVM.Int64Type(ctx)

        # create a function
        llvm_f, _ = create_function(T_nothing, [T_i64, T_i64])
        mod = LLVM.parent(llvm_f)

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            T_signal_store = LLVM.FunctionType(T_nothing, [T_i64, T_i64, T_i32])
            signal_store = LLVM.Function(mod, "__ockl_hsa_signal_store", T_signal_store)
            call!(builder, signal_store, [parameters(llvm_f)[1],
                                          parameters(llvm_f)[2],
                                          # __ATOMIC_RELEASE == 3
                                          ConstantInt(Int32(3); ctx)])

            ret!(builder)
        end

        call_function(llvm_f, Nothing, Tuple{UInt64,Int64}, :signal, :value)
    end
end
@inline device_signal_wait(signal::HSA.Signal, value::Int64) =
    device_signal_wait(signal.handle, value)
@inline @generated function device_signal_wait(signal::UInt64, value::Int64)
    Context() do ctx
        T_nothing = convert(LLVMType, Nothing; ctx)
        T_i32 = LLVM.Int32Type(ctx)
        T_i64 = LLVM.Int64Type(ctx)

        # create a function
        llvm_f, _ = create_function(T_nothing, [T_i64, T_i64])
        mod = LLVM.parent(llvm_f)

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            signal_match = BasicBlock(llvm_f, "signal_match"; ctx)
            signal_miss = BasicBlock(llvm_f, "signal_miss"; ctx)

            position!(builder, entry)
            br!(builder, signal_miss)

            position!(builder, signal_miss)
            T_sleep = LLVM.FunctionType(T_nothing, [T_i32])
            sleep_f = LLVM.Function(mod, "llvm.amdgcn.s.sleep", T_sleep)
            call!(builder, sleep_f, [ConstantInt(Int32(1); ctx)])
            T_signal_load = LLVM.FunctionType(T_i64, [T_i64, T_i32])
            signal_load = LLVM.Function(mod, "__ockl_hsa_signal_load", T_signal_load)
            loaded_value = call!(builder, signal_load, [parameters(llvm_f)[1],
                                                        # __ATOMIC_ACQUIRE == 2
                                                        ConstantInt(Int32(2); ctx)])
            cond = icmp!(builder, LLVM.API.LLVMIntEQ, loaded_value, parameters(llvm_f)[2])
            br!(builder, cond, signal_match, signal_miss)

            position!(builder, signal_match)
            ret!(builder)
        end

        call_function(llvm_f, Nothing, Tuple{UInt64,Int64}, :signal, :value)
    end
end

## Device

"Calls the host function stored in `hc` with arguments `args`."
@inline function hostcall!(hc::HostCall{UInt64,RT,AT}, args...) where {RT,AT}
    _hostcall_write_args!(hc, args...)
    _hostcall!(hc)
end
@inline @generated function _hostcall_write_args!(hc::HostCall{UInt64,RT,AT}, args...) where {RT,AT}
    ex = Expr(:block)

    # Copy arguments into buffer
    # Modified from CUDAnative src/device/cuda/dynamic_parallelism.jl
    off = 1
    for i in 1:length(args)
        T = args[i]
        sz = sizeof(T)
        # TODO: Should we do what CUDAnative does instead?
        ptr = :(reinterpret(LLVMPtr{$T,AS.Global}, hc.buf_ptr+$off-1))
        push!(ex.args, :(Base.unsafe_store!($ptr, args[$i])))
        off += sz
    end

    return ex
end
@inline @generated function _hostcall!(hc::HostCall{UInt64,RT,AT}) where {RT,AT}
    ex = Expr(:block)

    # Ring the doorbell
    push!(ex.args, :($device_signal_store!(hc.signal, hc.device_sentinel)))

    # Wait on doorbell
    push!(ex.args, :($device_signal_wait(hc.signal, hc.host_sentinel)))

    # Get return buffer and load first value
    ptr = :(reinterpret(LLVMPtr{LLVMPtr{$RT,AS.Global},AS.Global}, hc.buf_ptr))
    push!(ex.args, :(unsafe_load(unsafe_load($ptr))))

    return ex
end

## Host

@generated function _hostcall_args(hc::HostCall{UInt64,RT,AT}) where {RT,AT}
    ex = Expr(:tuple)

    # Copy arguments into buffer
    off = 1
    for i in 1:length(AT.parameters)
        T = AT.parameters[i]
        sz = sizeof(T)
        # TODO: Should we do what CUDAnative does instead?
        ptr = :(reinterpret(LLVMPtr{$T,AS.Global}, hc.buf_ptr+$off-1))
        # FIXME: We should not be using a device intrinsic here, even though it works...
        push!(ex.args, :(Base.unsafe_load($ptr)))
        off += sz
    end

    return ex
end

struct HostCallException <: Exception
    reason::String
    err::Union{Exception,Nothing}
    bt::Union{Vector,Nothing}
end
HostCallException(reason) = HostCallException(reason, nothing, backtrace())
HostCallException(reason, err) = HostCallException(reason, err, catch_backtrace())
function Base.showerror(io::IO, err::HostCallException)
    print(io, "HostCallException: $(err.reason)")
    if err.err !== nothing || err.bt !== nothing
        print(io, ":\n")
    end
    err.err !== nothing && Base.showerror(io, err.err)
    err.bt !== nothing && Base.show_backtrace(io, err.bt)
end

"""
    HostCall(func, return_type::Type, arg_types::Type{Tuple}) -> HostCall

Construct a `HostCall` that executes `func` with the arguments passed from the
calling kernel. `func` must be passed arguments of types contained in
`arg_types`, and must return a value of type `return_type`, or else the
hostcall will fail with undefined behavior.

Note: This API is currently experimental and is subject to change at any time.
"""
function HostCall(func, rettype, argtypes; return_task=false,
                  agent=get_default_agent(), maxlat=DEFAULT_HOSTCALL_LATENCY,
                  timeout=nothing, continuous=false, buf_len=nothing)
    signal = HSASignal()
    hc = HostCall(rettype, argtypes, signal.signal[].handle; agent=agent, buf_len=buf_len)

    tsk = @async begin
        ret_buf = Ref{Mem.Buffer}()
        while true
            try
                try
                    if !_hostwait(signal, hc.device_sentinel; maxlat=maxlat, timeout=timeout)
                        error("Hostwait timeout")
                    end
                catch err
                    throw(HostCallException("Error during hostwait", err))
                end
                if length(argtypes.parameters) > 0
                    args = try
                        _hostcall_args(hc)
                    catch err
                        throw(HostCallException("Error getting arguments", err))
                    end
                    @debug "Hostcall: Got arguments of length $(length(args))"
                else
                    args = ()
                end
                ret = try
                    func(args...,)
                catch err
                    throw(HostCallException("Error executing host function", err))
                end
                if typeof(ret) != rettype
                    throw(HostCallException("Host function result of wrong type: $(typeof(ret)), expected $rettype"))
                end
                if !isbits(ret)
                    throw(HostCallException("Host function result not isbits: $(typeof(ret))"))
                end
                @debug "Hostcall: Host function returning value of type $(typeof(ret))"
                try
                    ret_len = sizeof(ret)
                    ret_buf = Mem.alloc(agent, ret_len; coherent=true) # FIXME: Don't be coherent
                    ret_buf_ptr = Base.unsafe_convert(Ptr{typeof(ret)}, ret_buf)
                    Base.unsafe_store!(ret_buf_ptr, ret)
                    ret_buf_ptr = Base.unsafe_convert(Ptr{UInt64}, ret_buf)
                    args_buf_ptr = reinterpret(Ptr{UInt64}, hc.buf_ptr)
                    Base.unsafe_store!(args_buf_ptr, ret_buf_ptr)
                    HSA.signal_store_release(signal.signal[], hc.host_sentinel)
                catch err
                    throw(HostCallException("Error returning hostcall result", err))
                end
                @debug "Hostcall: Host function return completed"
                continuous || break
            catch err
                @error "Hostcall error" exception=(err,catch_backtrace())
                # TODO: Gracefully terminate just the immediate waiters
                # FIXME: Get the right queue
                kill_queue!(get_default_queue(agent))
                continuous || rethrow(err)
            end
        end
        # TODO: This is probably a bad idea to free while the kernel might
        # still hold a reference. We should probably have some way to
        # reference count these buffers from kernels (other than just
        # using a ROCArray, which can't currently be serialized to the
        # GPU).
        if isassigned(ret_buf)
            Mem.free(ret_buf)
        end
        Mem.free(hc.buf_ptr)
    end

    if return_task
        return hc, tsk
    else
        return hc
    end
end

# CPU functions
get_value(hc::HostCall{UInt64,RT,AT} where {RT,AT}) =
    HSA.signal_load_scacquire(HSA.Signal(hc.signal))
function _hostwait(signal, sentinel; maxlat=DEFAULT_HOSTCALL_LATENCY, timeout=nothing)
    @debug "Hostcall: Waiting on signal for sentinel: $sentinel"
    start_time = time_ns()
    while true
        value = HSA.signal_load_scacquire(signal.signal[])
        if value == sentinel
            @debug "Hostcall: Signal triggered with sentinel: $sentinel"
            return true
        end
        if timeout !== nothing
            now_time = time_ns()
            diff_time = (now_time - start_time) / 10^9
            if diff_time > timeout
                @debug "Hostcall: Signal wait timeout with sentinel: $sentinel"
                return false
            end
        end
        sleep(maxlat)
    end
end