text:
using DocStringExtensions
"""
Computes the evaporative and runoff flux as well as the storage of the interception component of the model.

$(SIGNATURES)

The function returns the amount of water leaving the component as runoff and evaporation and the amount of water stored in the component.
Function needs inputs to be in area independent units (e.g. mm)
"""
function interception(Potential_Evaporation::Float64, Precipitation::Float64, Temp::Float64, Interceptionstorage::Float64, Interceptionstoragecapacity::Float64, Temp_Thresh::Float64)
    #print(Interceptionstoragecapacity - Interceptionstorage, "\n")
    #print(Interceptionstorage <= Interceptionstoragecapacity, "\n")
    @assert Potential_Evaporation >= 0
    @assert Precipitation >= 0
    @assert Interceptionstorage >= 0
    @assert Interceptionstorage <= Interceptionstoragecapacity
    @assert Interceptionstoragecapacity >= 0
    @assert - 2.5 <= Temp_Thresh  <= 2.5
    @assert -60 <= Temp <= 60

    # if interception capacity is 0

    if Temp > Temp_Thresh
        if Interceptionstoragecapacity == 0.0
            Effective_Precipitation = Precipitation
            Interception_Evaporation = 0.0
            Interceptionstorage = 0.0
        #if the temperature is higher than freezing temp, precipitation falls as rain
        elseif Precipitation > 0
            #if it rains
            # the amount stored in Interception Reservoir will increase by amount of precipitation
            Interceptionstorage = Interceptionstorage + Precipitation
            # the excess precipitation will leave the reservoir directly
            Effective_Precipitation = max(0.0, Interceptionstorage - Interceptionstoragecapacity)
            #change in storage if Effective_Precipitation > 0
            if Effective_Precipitation > 0.0
                Interceptionstorage = Interceptionstoragecapacity
            end
            # after excess water has left interception storage evaporation occurs which can be 50% of total potential evaporation
            Interception_Evaporation = min(Interceptionstorage, Potential_Evaporation / 2.0)
            if Interception_Evaporation > 0.0
                #change in storage
                Interceptionstorage = Interceptionstorage - Interception_Evaporation
            end
        else
            # if it does not rain, no excess water leaves storage
            Effective_Precipitation = 0.0
            # the Interception Evporation will be either the amount stored or 50% the potential evaporation
            Interception_Evaporation = min(Interceptionstorage, Potential_Evaporation /2.0)
            # the amount stored in the Interception Reservoir will decrease because of evaporation
            if Interception_Evaporation > 0.0
                #change in storage
                Interceptionstorage = Interceptionstorage - Interception_Evaporation
            end

        end
    else
        # snow accumulates
        # evporation is 0 and no effective precipitation is released
        Interception_Evaporation = 0.0
        Effective_Precipitation = 0.0
        #Interceptionstorage = Interceptionstorage #amount stored does not change??!
    end
    #print("interception", Potential_Evaporation / 2.0, " ", Interception_Evaporation, "\n")
    @assert Interception_Evaporation <= Potential_Evaporation / 2.0
    @assert Effective_Precipitation >= 0
    @assert Interception_Evaporation >= 0
    @assert Interceptionstorage >= 0
    @assert Interceptionstorage <= Interceptionstoragecapacity
    return Effective_Precipitation::Float64, Interception_Evaporation::Float64, Interceptionstorage::Float64
end

"""
Computes the accumulation and melt of snow in the model.

$(SIGNATURES)

The function returns the amount of water leaving the component as runoff amount of water stored as snow in the model.
Function needs inputs to be in area independent units (e.g. mm)
"""
function snow(Area_Glacier::Float64, Precipitation::Float64, Temp::Float64, Snowstorage::Float64, Meltfactor::Float64, Mm::Float64, Temp_Thresh::Float64)

    @assert 0.0 <= Area_Glacier <= 1.0
    @assert Precipitation >= 0
    @assert Snowstorage >= 0
    @assert Meltfactor >= 0 #within the parameter range
    @assert Mm > 0 #within the parameter range
    @assert -2.5 <= Temp_Thresh <= 2.5 # it should be within the parameter range
    @assert -60 <= Temp <= 60

    if Temp > Temp_Thresh
        # if temperature higher than freezing temperature, melting takes place
        Melt = Meltfactor * Mm * ( (Temp - Temp_Thresh) / Mm + log(1 + exp(- (Temp - Temp_Thresh) / Mm) ) )
        # it can only melt as much as it is stored in snow storage
        Melt_Snow = min(Melt, Snowstorage)
        Melt_Glacier = Melt
        # the total melt is the combination of snow and glacier melt and the areal extent
        Melt_Total = Melt_Snow * (1.0 - Area_Glacier) + Melt_Glacier * Area_Glacier
        # the amount of snow stored decreases by amount melted
        Snowstorage = Snowstorage - Melt_Snow
        Melt_Glacier = Melt * Area_Glacier - Melt_Snow * Area_Glacier
    else
        # the amount of snow stored increases by Precipitation
        Snowstorage = Snowstorage + Precipitation
        # no snow melts
        Melt_Total = 0.0
        Melt_Glacier = 0.0
    end
    @assert Melt_Total >= 0
    @assert Snowstorage >= 0

    return Melt_Glacier::Float64, Melt_Total::Float64, Snowstorage::Float64
end

"""
Computes the processes in the soil component the hillslope hydrological response units of the model.

$(SIGNATURES)

The function returns the amount of water leaving the component as evaporation, overlandflow or preferential flow and the amount of water stored in soil component.
Function needs inputs to be in area independent units (e.g. mm)
beta: factor accounting for nonlinearity
Ce: Evapotranspiration control factor
"""
function soilstorage(Effective_Precipitation::Float64, Interception_Evaporation::Float64, Potential_Evaporation::Float64, Soilstorage::Float64, beta::Float64, Ce::Float64, Ratio_Pref::Float64, Soilstoragecapacity::Float64)
    @assert Effective_Precipitation >= 0
    @assert Interception_Evaporation >= 0
    @assert Potential_Evaporation >= 0
    #print("soil", Potential_Evaporation / 2.0, " ", Interception_Evaporation, "\n")
    @assert Potential_Evaporation / 2.0 - Interception_Evaporation >= -eps(Float64) * 100000
    @assert Soilstorage >= 0
    @assert Soilstoragecapacity - Soilstorage >= 0
    @assert Soilstoragecapacity >= 0 #within the parameter range
    @assert beta > 0 #within the parameter range
    @assert Ce > 0 #within the parameter range
    @assert Ratio_Pref >= 0

    if Effective_Precipitation > 0
        # rho represents the non linear process that only part of precipitation enters soil
        Ratio_Soil = 1 - (1 - (Soilstorage/Soilstoragecapacity))^beta
        #Ratio_Soil = min(Ratio_Soil, 1)
        @assert 0 <= Ratio_Soil <= 1
        # part of the water enters the soil, it cannot exceed the soil storage capacity
        Unused_Capacity = Soilstoragecapacity - Soilstorage
        Q_Soil = (1 - Ratio_Soil) * Effective_Precipitation
        if Unused_Capacity > Q_Soil
            Soilstorage = Soilstorage + Q_Soil
        else
            Q_Soil = Unused_Capacity
            Soilstorage = Soilstoragecapacity
        end
        Overlandflow = (Effective_Precipitation - Q_Soil) * Ratio_Pref
        # or flows into the groundwater
        Preferentialflow = (Effective_Precipitation - Q_Soil) * (1 - Ratio_Pref)
    else
        # if it does not rain no overland flow occurs
        Overlandflow = 0.0
        Preferentialflow = 0.0
    end
    # Transpiration in soil, only the part that not evaporated in interception reservoir can evaporate
    Potential_Soilevaporation = max(Potential_Evaporation - Interception_Evaporation, 0)
    # transpiration can maximum be the amount stored in soil, or a percentage of potential evaporation
    Soil_Evaporation = Potential_Soilevaporation * min(Soilstorage / (Soilstoragecapacity * Ce), 1.0)
    Soil_Evaporation = min(Soilstorage, Soil_Evaporation)
    Soilstorage = Soilstorage - Soil_Evaporation
    # if Soilstorage > Soilstoragecapacity
    #     print("evap", Soilstorage, " ", Soilstoragecapacity, "\n")
    # end
    @assert Overlandflow >= 0
    @assert Preferentialflow >= 0
    @assert Soil_Evaporation <= Potential_Evaporation - Interception_Evaporation
    @assert Soil_Evaporation >= 0
    @assert Soilstorage >= 0
    @assert Soilstoragecapacity - Soilstorage >= 0
    return Overlandflow::Float64, Preferentialflow::Float64, Soil_Evaporation::Float64, Soilstorage::Float64
end
"""
Computes the processes in the soil component the riparian hydrological response units of the model.

$(SIGNATURES)

The function returns the amount of water leaving the component as evaporation, overlandflow flow and the amount of water stored in soil component.
Function needs inputs to be in area independent units (e.g. mm)
beta: factor accounting for nonlinearity
Ce: Evapotranspiration control factor
Riparian Discharge: Input into Soil from ground water
"""
function ripariansoilstorage(Effective_Precipitation, Interception_Evaporation, Potential_Evaporation, Riparian_Discharge, Soilstorage, beta, Ce, Drainagecapacity, Soilstoragecapacity)
    @assert Effective_Precipitation >= 0
    @assert Interception_Evaporation >= 0
    @assert Potential_Evaporation >= 0
    @assert - eps(Float64) * 10 <= Potential_Evaporation * 0.5 - Interception_Evaporation
    @assert Potential_Evaporation / 2.0 - Interception_Evaporation >= -eps(Float64) * 100000
    @assert Riparian_Discharge >= 0
    #@assert Soil_Evaporation >= 0 #or should it be zero?
    @assert Soilstorage >= 0
    @assert Soilstoragecapacity - Soilstorage >= 0
    @assert Soilstoragecapacity > 0 #within the parameter range
    @assert beta > 0 #within the parameter range
    @assert Ce > 0 #within the parameter range
    @assert Drainagecapacity >= 0

    Ratio_Soil = 1 - (1 - (Soilstorage/Soilstoragecapacity))^beta

    Unused_Capacity = Soilstoragecapacity - Soilstorage
    Q_Soil = (1 - Ratio_Soil) * (Effective_Precipitation + Riparian_Discharge)
    if Unused_Capacity > Q_Soil
        Soilstorage = Soilstorage + Q_Soil
    else
        Q_Soil = Unused_Capacity
        Soilstorage = Soilstoragecapacity
    end
    # the other part does not enter the soil but flows into the fast reservoir
    Overlandflow = (Effective_Precipitation + Riparian_Discharge - Q_Soil)
    # Transpiration in soil, only the part that not evaporated in interception reservoir can evaporate
    Potential_Soilevaporation = max(Potential_Evaporation - Interception_Evaporation,0)
    # transpiration can maximum be the amount stored in soil, or a percentage of potential evaporation
    Soil_Evaporation = Potential_Soilevaporation * min(Soilstorage / (Soilstoragecapacity * Ce), 1)
    Soil_Evaporation = min(Soilstorage, Soil_Evaporation)
    Soilstorage = Soilstorage - Soil_Evaporation
    # Part of the water stored in soil will drain to a maximum capacity, which is routed into the fast response reservoir
    if Drainagecapacity > 0
        Fastdrainage = (Soilstorage / Soilstoragecapacity) * Drainagecapacity
        Soilstorage = Soilstorage - Fastdrainage
        Overlandflow = Overlandflow + Fastdrainage
    end

    @assert Overlandflow >= 0
    @assert Soil_Evaporation <= max(Potential_Evaporation - Interception_Evaporation,0)
    @assert Soil_Evaporation >= 0
    @assert Soilstorage >= 0
    @assert Soilstoragecapacity - Soilstorage >= 0
    return Overlandflow, Soil_Evaporation, Soilstorage
end

"""
Computes the processes in the fast component of the model. The fast component is represented by a linear response reservoir.

$(SIGNATURES)

The function returns the amount of water leaving the component as runoff towards the stream and the amount of water stored in fast component.
Function needs inputs to be in area independent units (e.g. mm)
OVerlandflow: amount of water entering the component
Kf: Fast storage Coefficient
Riparian Discharge: Input into Soil from ground water
"""
function faststorage(Overlandflow, Faststorage, Kf)
    @assert Overlandflow >= 0
    @assert Faststorage >= 0
    @assert Kf >=0 and <= 1
    # the fast storage increases with the overland flow
    Faststorage = Faststorage + Overlandflow
    # a part of the fast storage gets redirected into discharge depending on the reservoir constant (linear response)
    if Kf * Faststorage < Faststorage
        Fast_Discharge = Kf * Faststorage
        Faststorage = Faststorage - Fast_Discharge
    else
        Fast_Discharge = Faststorage
        Faststorage = 0
    end
    @assert Fast_Discharge >= 0
    @assert Faststorage >= 0
    return Fast_Discharge, Faststorage
end
"""
Computes the processes in the slow component of the model. The slow component is represented by a linear response reservoir.

$(SIGNATURES)

The function returns the amount of water leaving the component as runoff towards the stream and the amount of water stored in slow component (ground water)
Function needs inputs to be in area independent units (e.g. mm)
GWflow: flow towards slow component
Area_Riparian: areal percentage of riparian zone
Ks: Slow storage coefficient
Ratio_Riparian: Share of riparian flow of total outflow of slow reservoir
"""
function slowstorage(GWflow, Slowstorage, Area_Riparian::Float64, Ks::Float64, Ratio_Riparian::Float64)
    @assert GWflow >= 0
    @assert Slowstorage >= 0
    @assert Ks >=0 and <= 1
    @assert Ratio_Riparian >=0 and <= 1

    Slowstorage = Slowstorage + GWflow
    Slow_Discharge = Ks * Slowstorage * (1 - Ratio_Riparian)
    # the riparian discharge is the areal percentage of the total possible riparian discharge
    Riparian_Discharge = Ks * Slowstorage * Ratio_Riparian * Area_Riparian
    Slowstorage = Slowstorage - Slow_Discharge - Riparian_Discharge

    @assert Riparian_Discharge >= 0
    @assert Slow_Discharge >= 0
    @assert Slowstorage >= 0
    return Riparian_Discharge, Slow_Discharge, Slowstorage
end

export interception
export snow
export ripariansoilstorage
export soilstorage
export faststorage
export slowstorage

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
function variablecost(m::JuMP.AbstractModel, devices::Array{T,1}, device_formulation::Type{D}, system_formulation::Type{S}) where {T <: PowerSystems.ThermalGen, D <: AbstractThermalDispatchForm, S <: PM.AbstractPowerFormulation}

    p_th = m[:p_th]
    time_index = m[:p_th].axes[2]
    name_index = m[:p_th].axes[1]

    var_cost = JuMP.AffExpr()

    for  (ix, name) in enumerate(name_index)
        if name == devices[ix].name
            c = gencost(m, p_th[name,:], devices[ix].econ.variablecost)
        else
            @error "Bus name in Array and variable do not match"
        end

        # TODO: Add no load costs?

        # TODO: Move this to a new function
        if isa(var_cost,JuMP.AffExpr) && isa(c,JuMP.AffExpr)
            JuMP.add_to_expression!(var_cost,c)
        elseif isa(var_cost,JuMP.GenericQuadExpr) && isa(c,JuMP.GenericQuadExpr)
            JuMP.add_to_expression!(var_cost,c)
        else
            var_cost += c
        end
    end

    return var_cost

end

ext: jl


text:
using SearchLight, SearchLight.Migrations, SearchLight.Relationships

cd(@__DIR__)

connection_file = joinpath(@__DIR__,"mysql_connection.yml")
conn_info_postgres = SearchLight.Configuration.load(connection_file)
const conn = SearchLight.connect(conn_info_postgres)

try
  SearchLight.Migrations.status()
catch _
  SearchLight.Migrations.create_migrations_table()
end

isempty(SearchLight.Migrations.downed_migrations()) || SearchLight.Migrations.all_up!!()

Base.@kwdef mutable struct User <: AbstractModel
  id::DbId = DbId()
  username::String = ""
  password::String = ""
  name::String = ""
  email::String = ""
end

Base.@kwdef mutable struct Role <: AbstractModel
  id::DbId = DbId()
  name::String = ""
end

Base.@kwdef mutable struct Ability <: AbstractModel
  id::DbId = DbId()
  name::String = ""
end

u1 = findone_or_create(User, username = "a") |> save!
r1 = findone_or_create(Role, name = "abcd") |> save!

for x in 'a':'d'
  findone_or_create(Ability, name = "$x") |> save!
end

Relationships.Relationship!(u1, r1)

for a in all(Ability)
  Relationships.Relationship!(r1, a)
end

Relationships.related(u1, Role)
Relationships.related(findone(Role, id = 1), Ability)
Relationships.related(u1, Ability, through = [Role])

max_stars_repo_name: GenieFramework/SearchLightMySQL

max_issues_repo_name: GenieFramework/SearchLightMySQL

max_forks_repo_name: GenieFramework/SearchLightMySQL


text:
scalaradd(a, b) = a + b

function sum(v)
  s = 0.0
  N = length(v)
  for i = 1:N
      s += v[i]
  end
  return s
end

# y = dot(a,b)
function dot(a, b)
  s = 0.0
  N = length(a)
  for i = 1:N
      s += a[i] * b[i]
  end
  return s
end

function tuplereturn(x)
  a = x * 1.2
  b = x * 1.3
  return a, b
end

function tuplecall(x)
  ab = tuplereturn(x)
  r = ab[1]
  #r = getfield(ab,1)
  return r
end

# y = W*X
mul(W, X) = mul_(zeros(size(W,1), size(X,2)), W, X)
function mul_(y, W, X)
    szW = size(W)
    szX = size(X)
    s = 0.0
    for n = 1:szX[2]
        for i = 1:szW[1]
            s = 0.0
            for j = 1:szW[2]
                s += W[i,j] * X[j,n]
            end
            y[i,n] = s
        end
    end
    return y
end

# y = W*X .+ b
muladd(W, X, b) = muladd_(zeros(size(W,1), size(X,2)), W, X, b)
function muladd_(y, W, X, b)
    szW = size(W)
    szX = size(X)
    s = 0.0
    for n = 1:szX[2]
        for i = 1:szW[1]
            s = 0.0
            for j = 1:szW[2]
                s += W[i,j] * X[j,n]
            end
            y[i,n] = s + b[i]
        end
    end
    return y
end

# y = a .+ b
add(a, b) = add_(zeros(size(a)), a, b)
function add_(r, a, b)
  N = length(r)
  for i = 1:N
      r[i] = a[i] + b[i]
  end
  return r
end

function exports()
  scalaradd(1.0, 2.0)

  a=randn(3,1)
  b=randn(3,1)
  dot(a,b)
  res = add(a,b)

  tuplereturn(1.2)
  tuplecall(2.1)

  M1=randn(4,2)
  M2=randn(2,5)
  r = mul(M1,M2)

  b = randn(4,1)
  r2 = muladd(M1,M2, b)
  sum(a)
end

exports()


text:
# We define several Returns to Scale (RS) types. The RS constraint used in DDF depends on these RTS type
abstract type RS end

type CRS <: RS
end

type VRS <: RS
end

type NIRS <: RS
end

type NDRS <: RS
end

# The default RTS constraint returns an all zero vector.
# Constant Returns to Scale (CRS) only imposes \lambda_i >= 0
function RSconstraint{T<:RS}(K,RStype::T)
	return zeros(Float64,1,K+1),0,'='
end

# Variable Returns to Scale (VRS) imposes sum_i{\lambda_i} = 1
function RSconstraint(K,RStype::VRS)
	return [0 ones(1,K)],1,'='
end

# Non-Increasing Returns to Scale (NIRS) imposes sum_i{\lambda_i} <= 1
function RSconstraint(K,RStype::NIRS)
	return [0 ones(1,K)],1,'<'
end

# Non-Decreasing Returns to Scale (NDRS) imposes sum_i{\lambda_i} >= 1
function RSconstraint(K,RStype::NDRS)
	return [0 ones(1,K)],1,'>'
end

max_stars_repo_name: kepiej/NEPA


text:
# ===================
# Inhomogeneous case
# ===================

#using LinearAlgebra.BLAS

# inhomogeneous case without invariant, non-recursive implementation
function reach_inhomog_BOX!(F::Vector{ReachSet{N, Hyperrectangle{N, VNC, VNR}}},
                            Ω0::Hyperrectangle{N, VNC, VNR},
                            Φ::AbstractMatrix,
                            NSTEPS::Integer,
                            δ::N,
                            X::Universe,
                            U::Hyperrectangle,
                            recursive::Val{false},
                            Δt0::TN) where {N, TN, VNC, VNR}

    # initial reach set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = ReachSet(Ω0, Δt)

    # preallocations
    c = Vector{VNC}(undef, NSTEPS)
    r = Vector{VNR}(undef, NSTEPS)

    c[1] = Ω0.center
    r[1] = Ω0.radius

    # input sequence
    Wk₊ = copy(U)

    # cache for powers of Φ
    Φ_power_k = similar(Φ)
    copyto!(Φ_power_k, Φ)
    Φ_power_k_cache = similar(Φ_power_k)

    k = 2
    @inbounds while k <= NSTEPS
        # Hₖ = overapproximate(Φ_power_k * Ω0, Hyperrectangle)
        abs_Φ_power_k = abs.(Φ_power_k)
        c[k] = Φ_power_k * c[1] + Wk₊.center
        r[k] = abs_Φ_power_k * r[1] + Wk₊.radius
        Hₖ = Hyperrectangle(c[k], r[k], check_bounds=false)

        Δt += δ
        F[k] = ReachSet(Hₖ, Δt)

        # in-place computation of Wk₊ <- Wk₊ + Φ_power_k * U,
        # by overapproximating the second term with a hyperrectangle
        Wk₊.center .+= Φ_power_k * U.center
        Wk₊.radius .+= abs_Φ_power_k * U.radius
        #gemv!('N', 1.0, Φ_power_k, U.center, 1.0, Wk₊.center)
        #gemv!('N', 1.0, abs_Φ_power_k, U.radius, 1.0, Wk₊.radius)

        mul!(Φ_power_k_cache, Φ_power_k, Φ)
        copyto!(Φ_power_k, Φ_power_k_cache)
        k += 1
    end
    return F
end

# homogeneous case with an invariant, non-recursive implementation
function reach_inhomog_BOX!(F::Vector{ReachSet{N, Hyperrectangle{N, VNC, VNR}}},
                            Ω0::Hyperrectangle{N, VNC, VNR},
                            Φ::AbstractMatrix,
                            NSTEPS::Integer,
                            δ::N,
                            X::LazySet,
                            U::Hyperrectangle,
                            recursive::Val{false},
                            Δt0::TN) where {N, TN, VNC, VNR}

    # initial reach set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = ReachSet(Ω0, Δt)

    # preallocations
    c = Vector{VNC}(undef, NSTEPS)
    r = Vector{VNR}(undef, NSTEPS)

    c[1] = Ω0.center
    r[1] = Ω0.radius

    # input sequence
    Wk₊ = copy(U)

    # cache for powers of Φ
    Φ_power_k = similar(Φ)
    copyto!(Φ_power_k, Φ)
    Φ_power_k_cache = similar(Φ_power_k)

    k = 2
    @inbounds while k <= NSTEPS
      # Hₖ = overapproximate(Φ_power_k * Ω0, Hyperrectangle)
      abs_Φ_power_k = abs.(Φ_power_k)
      c[k] = Φ_power_k * c[1] + Wk₊.center
      r[k] = abs.(Φ_power_k) * r[1] + Wk₊.radius
      Hₖ = Hyperrectangle(c[k], r[k], check_bounds=false)

      _is_intersection_empty(X, Hₖ) && break

      Δt += δ
      F[k] = ReachSet(Hₖ, Δt)

      # in-place computation of Wk₊ <- Wk₊ + Φ_power_k * U,
      # by overapproximating the second term with a hyperrectangle
      Wk₊.center .+= Φ_power_k * U.center
      Wk₊.radius .+= abs_Φ_power_k * U.radius

      mul!(Φ_power_k_cache, Φ_power_k, Φ)
      copyto!(Φ_power_k, Φ_power_k_cache)
      k += 1
    end
    if k < NSTEPS + 1
        resize!(F, k-1)
    end
    return F
end


text:
"""
We update constraints doing:

1. ...
"""
function update_constraints!(constraints::Constraints, trial::Trial)
    lower_bound = LetterCounts()

    for i in 1:5
        char = trial.candidate.chars[i]

        # If we see a PRESENT_NOWHERE code, we update the lower bound and forbid the character with some caveats.
        if trial.result.codes[i] === PRESENT_NOWHERE
            # Don't set this if there's a PRESENT_ELSEWHERE code set or an PRESENT_HERE set for this character.
            not_repeated = true
            for j in 1:5
                if (trial.result.codes[j] === PRESENT_ELSEWHERE && trial.candidate.chars[j] === char) ||
                   (trial.result.codes[j] === PRESENT_HERE && trial.candidate.chars[j] === char)
                   not_repeated = false
                end
            end
            if not_repeated
                constraints.global_ban[char] = true
                lower_bound[char] = 0
            end
        # If we see an PRESENT_ELSEWHERE code, we store a Boolean flag forbidding that character.
        elseif trial.result.codes[i] === PRESENT_ELSEWHERE
            constraints.local_ban[i][char] = true
            lower_bound[char] += 1
        # If we see an PRESENT_HERE code, we store the matched character and a Boolean flag.
        elseif trial.result.codes[i] == PRESENT_HERE
            constraints.is_exact_match[i] = true
            constraints.match[i] = char
            lower_bound[char] += 1
        else
            error("Invariant failed")
        end
    end

    for c in 'a':'z'
        constraints.lower_bound[c] = max(constraints.lower_bound[c], lower_bound[c])
    end
end


text:
module GAIO

using Parameters: @with_kw
using Memoize: @memoize

import Combinatorics: combinations

include("Trees.jl")
import .Trees: Node

export BoxNode, generate_sites!, eval!, subdivide!, init_box_tree,
    results, sites, contains;

export children, parent, Leaves;
export ≺ ;

include("interfaces.jl");

# implementation of the `SomeDataBase` interface
include("DataBase.jl");

include("Boxes.jl");
include("mo_utils.jl")

end


text:
using ModelParameters
using Test
import Base.@kwdef

@testset "Get Properties" begin
    x1 = Parameter(1.0, "x1")
    x2 = BoundedParameter(1.0, 0.0, 10.0, "x2")

    x3 = Parameter(1, "x3")
    x4 = BoundedParameter(1, 0, 10, "x4")

    @kwdef struct Params <: ParameterGroup
        y1::Parameter{Float64,String} = x1
        y2::BoundedParameter{Float64,String} = x2

        y3::Parameter{Int64,String} = x3
        y4::BoundedParameter{Int64,String} = x4
    end

    z1 = x1.value
    @test z1 ≈ 1.0
    z2 = x2.value
    @test z2 ≈ 1.0
    z3 = x3.value
    @test z3 == 1
    z4 = x2.value
    @test z4 == 1

    T = Params()

    z1 = T.y1
    @test z1 ≈ 1.0
    z2 = T.y2
    @test z2 ≈ 1.0
    z3 = T.y3
    @test z3 == 1
    z4 = T.y4
    @test z4 == 1
end

@testset "Set Properties" begin

    x1 = Parameter(1.0, "x1")
    x2 = BoundedParameter(1.0, 0.0, 10.0, "x2")

    x3 = Parameter(1, "x3")
    x4 = BoundedParameter(1, 0, 10, "x4")

    @kwdef struct Params <: ParameterGroup
        y1::Parameter{Float64,String} = x1
        y2::BoundedParameter{Float64,String} = x2

        y3::Parameter{Int64,String} = x3
        y4::BoundedParameter{Int64,String} = x4
    end
    T = Params()
    T = set(T, :y1, 0.3)
    @test T.y1 ≈ 0.3
    T = T ← ((:y1, 0.2), (:y3, 0))
    @test T.y1 ≈ 0.2
    @test T.y3 == 0

end

@testset "Get Bounds" begin
    x1 = Parameter(1.0, "x1")
    x2 = BoundedParameter(1.0, 0.0, 10.0, "x2")

    x3 = Parameter(1, "x3")
    x4 = BoundedParameter(1, 0, 10, "x4")

    @kwdef struct Params <: ParameterGroup
        y1::Parameter{Float64,String} = x1
        y2::BoundedParameter{Float64,String} = x2

        y3::Parameter{Int64,String} = x3
        y4::BoundedParameter{Int64,String} = x4
    end
    T = Params()

    initial_guess = get_initial_guess([x2, x4])
    @test initial_guess == [1.0, 1]

    bounds = get_bounds(T,:y2)
    @test bounds == (0.0, 10.0)

    lower_bounds = get_lower_bounds([x2, x4])
    @test lower_bounds == [0.0, 0]

    upper_bounds = get_upper_bounds([x2, x4])
    @test upper_bounds == [10.0, 10]

    descr = description(T,:y1)
    @test descr == "x1"

end


text:
using CSV
using DataFrames
using Pipe
using Dates

function read_in_TP(date, state)

    # read in the reff file
    TP_file_name = "results/soc_mob_R"*date*".csv"
    # drop the first column
    df = CSV.read(TP_file_name, DataFrame, drop=[1])
    # extract the unique states
    unique_states = unique(df[!,"state"])

    TP_dict_local = Dict{String, Array}()
    # get the dates
    TP_dict_local["date"] = @pipe df |>
        filter(:state => ==(state), _) |>
        filter(:type => ==("R_L"), _) |>
        select(_, :date) |>
        Matrix(_)[:]
    # create a vector of the indices
    TP_dict_local["t_indices"] = collect(1:length(TP_dict_local["date"]))

    TP_dict_import = Dict{String, Array}()
    # get the dates
    TP_dict_import["date"] = @pipe df |>
        filter(:state => ==(state), _) |>
        filter(:type => ==("R_L"), _) |>
        select(_, :date) |>
        Matrix(_)[:]
    # create a vector of the indices
    TP_dict_import["t_indices"] = collect(1:length(TP_dict_import["date"]))

    for state in unique_states
        # filter by state and RL and keep only TP trajectories
        df_state_matrix = @pipe df |>
            filter(:state => ==(state), _) |>
            filter(:type => ==("R_L"), _) |>
            select(_, Not(1:8)) |>
            Matrix(_)
        TP_dict_local[state] = df_state_matrix
        # filter by state and RL and keep only TP trajectories
        df_state_matrix = @pipe df |>
            filter(:state => ==(state), _) |>
            filter(:type => ==("R_I"), _) |>
            select(_, Not(1:8)) |>
            Matrix(_)
        TP_dict_import[state] = df_state_matrix
    end

    return (TP_dict_local, TP_dict_import)

end

function create_state_TP_matrices(forecast_start_date, date, state)
    """
    Read in the TP for a given state from csv with filedate, date. Adjust the indices
    according to the forecast_start_date and then create matrices of local and import TP's.
    The indices vector is used to allow for pre-forecast window infection dates.
    """

    # read in the TP
    (TP_dict_local, TP_dict_import) = read_in_TP(date, state)
    # adjust the indices so that the forecast start date is t = 0
    TP_indices = convert.(Int, Dates.value.(TP_dict_local["date"] - forecast_start_date))
    # make a matrix for a given state
    TP_local = TP_dict_local[state]
    TP_import = TP_dict_import[state]

    return (TP_indices, TP_local, TP_import)
end


text:
using DashBootstrapComponents, DashHtmlComponents

row = html_div([
    dbc_row([
        dbc_col(html_div("One of three columns"), md=4),
        dbc_col(html_div("One of three columns"), md=4),
        dbc_col(html_div("One of three columns"), md=4),
    ]),
    dbc_row([
        dbc_col(html_div("One of four columns"), width=6, lg=3),
        dbc_col(html_div("One of four columns"), width=6, lg=3),
        dbc_col(html_div("One of four columns"), width=6, lg=3),
        dbc_col(html_div("One of four columns"), width=6, lg=3),
    ]),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
module MPGenbb

using Random
using Distributions
using Documenter
using DataFrames

include("MomentumFuncs.jl")
export  get_pMag,            # returns the magnitude of the momentum vector calculated from kinetic energy T
        get_px,              # x-component of momentum vector
        get_py,              # y-component of momentum vector
        get_pz,              # z-component of momentum vector
        get_first_vector,    # returns the momentum vector of the first electron (independent, direction generated uniformaly)
        get_second_vector    # returns the momentum vector of the second electron (dependent on first, direction generated based on θdif distribuition)

include("StringFuncs.jl")
export get_particle_string,
       get_header_string,
       get_event_string

include("SamplingFuncs.jl")
export sample_phi,           # returns ϕ azimuthal angle uniformly
       sample_theta,         # returns θ polar angle uniformly
       sample_energies,      # returns !!tuple!! of electron kintetic energies T1, T2
       sample_theta_dif,     # returns θdif as the angle between the 2 electrons (based on 1 - cos(θ) angular distribution)
       solve_quadratic       # returns x1, x2 of a quadratic equation

end # module MPgenbb


text:
using LinearAlgebra
using Distributions
using Zygote
using Random
using RandomNumbers.PCG
using Distributed
@everywhere using SharedArrays

struct PSPoint
    q::Vector{Float64}
    p::Vector{Float64}
end

function maybeupdate!(c::Dict, s::Symbol, v)
    if !haskey(c, s)
        merge!(c, Dict(s => v))
    end
    return nothing
end

function checkcontrol(c::Dict)

    c = convert(Dict{Any, Any}, c)

    # TODO (ear) add checks on all values
    maybeupdate!(c, :iterations, 2000)
    maybeupdate!(c, :iterations_warmup, c[:iterations] ÷ 2)
    maybeupdate!(c, :chains, 4)
    maybeupdate!(c, :chainid, 1)
    maybeupdate!(c, :rng, PCGStateOneseq(UInt64, PCG_RXS_M_XS))
    maybeupdate!(c, :maxtreedepth, 10)
    maybeupdate!(c, :μ, log(10))
    maybeupdate!(c, :adaptinitpercent, 0.15)
    maybeupdate!(c, :adapttermpercent, 0.1)
    maybeupdate!(c, :adaptwindowsize, 25)
    maybeupdate!(c, :initradius, 2)
    maybeupdate!(c, :initattempts, 100)
    maybeupdate!(c, :δ, 0.8)
    maybeupdate!(c, :γ, 0.05)
    maybeupdate!(c, :t0, 10)
    maybeupdate!(c, :κ, 0.75)
    maybeupdate!(c, :regularize, true)
    maybeupdate!(c, :skewsymmetric, false)

    return c
end

function stan(U, ndim;
              M::AbstractArray{Float64} = ones(ndim), control::Dict = Dict())

    control = checkcontrol(control)
    I = control[:iterations] + control[:iterations_warmup]
    chains = control[:chains]

    samples = SharedArray{Float64}(control[:iterations], chains, ndim)
    leaps = SharedArray{Float64}(I, chains)
    acceptstats = SharedArray{Float64}(I, chains)
    stepsizes = SharedArray{Float64}(I, chains)
    treedepths = SharedArray{Float64}(I, chains)
    energies = SharedArray{Float64}(I, chains)
    divergences = SharedArray{Bool}(I, chains)
    massmatrices = SharedArray{Float64}(size(M)..., chains)

    @sync @distributed for chain in 1:control[:chains]
        control = merge(control, Dict(:chainid => chain))
        samples[:, chain, :], d = hmc(U, ndim; M = M, control = control)

        leaps[:, chain] = d[:leapfrog]
        acceptstats[:, chain] = d[:acceptstat]
        stepsizes[:, chain] = d[:stepsize]
        treedepths[:, chain] = d[:treedepth]
        energies[:, chain] = d[:energy]
        divergences[:, chain] = d[:divergent]
        # TODO (ear) there's got to be a better way
        massmatrices[fill((:), length(size(M)))..., chain] = d[:massmatrix]
    end

    return samples, Dict(:leapfrog => leaps,
                         :acceptstat => acceptstats,
                         :stepsize => stepsizes,
                         :treedepth => treedepths,
                         :energy => energies,
                         :divergent => divergences,
                         :massmatrix => massmatrices,
                         :control => control)
end

function hmc(U, ndim;
             M::AbstractArray{Float64} = ones(ndim), control::Dict = Dict())

    control = checkcontrol(control)
    I = control[:iterations] + control[:iterations_warmup]

    # advance RandomNumbers RNG
    # TODO (ear) figure out their type structure
    # and document this requirement.
    advance!(control[:rng], convert(UInt64, 1 << 50 * control[:chainid]))

    ∇U = q -> first(Zygote.gradient(U, q))
    q = initializesample(ndim, U, ∇U, control)

    samples = zeros(I, ndim)
    samples[1, :] = q
    zsample = PSPoint(q, generatemomentum(control[:rng], ndim, M))

    ε = findepsilon(1.0, zsample, U, ∇U, M, control)
    μ = control[:μ]
    εbar = 0.0
    sbar = 0.0
    xbar = 0.0

    W = WelfordState(zeros(ndim), zero(M), 0)
    openwindow, closewindow, lastwindow, windowstep =
        adaptionwindows(control)

    leaps = zeros(Int, I)
    treedepths = zeros(Int, I)
    energies = zeros(I)
    acceptstats = zeros(I)
    εs = zeros(I)
    divergences = falses(I)
    stepsizecounter = 0

    for i = 2:I
        z = PSPoint(samples[i - 1, :], generatemomentum(control[:rng], ndim, M))
        H0 = hamiltonian(U, z, M)
        H = H0

        zf = z
        zb = z
        zsample = z
        zpr = z

        # Momentum and sharp momentum at forward end of forward subtree
        pff = z.p
        psharpff = control[:skewsymmetric] ? z.p : rhosharp(z.p, M) # dtau_dp

        # Momentum and sharp momentum at backward end of forward subtree
        pfb = z.p
        psharpfb = psharpff

        # Momentum and sharp momentum at forward end of backward subtree
        pbf = z.p
        psharpbf = psharpff

        # Momentum and sharp momentum at backward end of backward subtree
        pbb = z.p
        psharpbb = psharpff

        # Integrated momenta along trajectory
        rho = z.p

        α = 0.0
        depth = 0
        nleapfrog = 0
        lsw = 0.0

        while depth < control[:maxtreedepth]

            rhof = zeros(length(rho))
            rhob = zeros(length(rho))
            lswsubtree = -Inf

            if rand(control[:rng]) > 0.5
                rhob = rho
                pbf = pff
                psharpbf = psharpff

                zf, zpr, validsubtree,
                psharpfb, psharpff, rhof, pfb, pff,
                nleapfrog, lswsubtree, α  =
                    buildtree(depth, zf,
                              psharpfb, psharpff, rhof, pfb, pff,
                              H0, 1 * ε, U, ∇U, M,
                              nleapfrog, lswsubtree, α, control)
            else
                rhof = rho
                pfb = pbb
                psharpfb = psharpbb

                zb, zpr, validsubtree,
                psharpbf, psharpbb, rhob, pbf, pbb,
                nleapfrog, lswsubtree, α =
                    buildtree(depth, zb,
                              psharpbf, psharpbb, rhob, pbf, pbb,
                              H0, -1 * ε, U, ∇U, M,
                              nleapfrog, lswsubtree, α, control)
            end

            if !validsubtree
                divergences[i] = true
                break
            end
            depth += 1

            if lswsubtree > lsw
                zsample = zpr
            else
                if rand(control[:rng]) < exp(lswsubtree - lsw)
                    zsample = zpr
                end
            end

            lsw = logsumexp(lsw, lswsubtree)

            # Demand satisfication around merged subtrees
            rho = rhob + rhof
            persist = stancriterion(psharpbb, psharpff, rho)

            # Demand satisfaction between subtrees
            rhoextended = rhob + pfb
            persist &= stancriterion(psharpbb, psharpfb, rhoextended)

            rhoextended = rhof + pbf
            persist &= stancriterion(psharpbf, psharpff, rhoextended)

            if !persist
                break
            end
        end # end while

        samples[i, :] = zsample.q
        energies[i] = hamiltonian(U, zsample, M)
        treedepths[i] = depth
        leaps[i] = nleapfrog
        εs[i] = ε
        acceptstats[i] = α / nleapfrog

        if i < control[:iterations_warmup]
            ε, εbar, stepsizecounter, sbar, xbar =
                adaptstepsize(acceptstats[i], stepsizecounter, sbar, xbar, μ, control)

            if openwindow <= i <= lastwindow
                W = accmoments(W, samples[i, :])
            end

            if i == closewindow
                # reset var
                M = samplevariance(W, control)
                W = WelfordState(zeros(ndim), zero(M), 0)

                # reset stepsize
                ε = findepsilon(ε, zsample, U, ∇U, M, control)
                μ = log(10 * ε)
                sbar = 0.0
                xbar = 0.0
                stepsizecounter = 0

                # update windows
                windowstep *= 2
                closewindow += windowstep
                nextclosewindow = closewindow + windowstep * 2
                if nextclosewindow > lastwindow
                    closewindow = lastwindow
                end
            end
        else
            ε = εbar
        end

    end # end for

    postwarmup = control[:iterations_warmup]+1:I
    return samples[postwarmup, :, :], Dict(:stepsize => εs,
                                           :leapfrog => leaps,
                                           :acceptstat => acceptstats,
                                           :treedepth => treedepths,
                                           :energy => energies,
                                           :divergent => divergences,
                                           :massmatrix => M)
end

function buildtree(depth::Int, z::PSPoint,
                    psharpbeg, psharpend, rho, pbeg, pend,
                    H0::Float64, ε::Float64, U, ∇U, M::AbstractArray{Float64},
                    nleapfrog::Int, logsumweight::Float64, α::Float64, control)
    if depth == 0
        zpr = leapfrog(z, ∇U, ε, M)
        z = zpr
        nleapfrog += 1

        H = hamiltonian(U, zpr, M)
        if isnan(H)
            H = Inf
        end

        divergent = false
        if H - H0 > 1000.0
            divergent = true
        end

        Δ = H0 - H
        logsumweight = logsumexp(logsumweight, Δ)
        α += Δ > 0.0 ? 1.0 : exp(Δ)

        psharpbeg = control[:skewsymmetric] ? z.p : rhosharp(z.p, M) # dtau_dp
        psharpend = psharpbeg

        rho += zpr.p
        pbeg = zpr.p
        pend = pbeg

        return z, zpr, !divergent, psharpbeg, psharpend, rho, pbeg, pend, nleapfrog, logsumweight, α
    end

    lswinit = -Inf

    psharpinitend = similar(z.p)
    rhoinit = zeros(length(rho))
    pinitend = similar(z.p)

    z, zpr, validinit,
    psharpbeg, psharpinitend, rhoinit, pbeg, pinitend,
    nleapfrog, lswinit, α =
        buildtree(depth - 1, z,
                   psharpbeg, psharpinitend, rhoinit, pbeg, pinitend,
                   H0, ε, U, ∇U, M, nleapfrog, lswinit, α, control)

    if !validinit
        return zpr, zpr, false,
        psharpbeg, psharpend, rho, pbeg, pend,
        nleapfrog, logsumweight, α
    end

    lswfinal = -Inf

    psharpfinalbeg = similar(z.p)
    rhofinal = zeros(length(rho))
    pfinalbeg = similar(z.p)

    z, zfinalpr, validfinal,
    psharpfinalbeg, psharpend, rhofinal, pfinalbeg, pend,
    nleapfrog, lswfinal, α =
        buildtree(depth - 1, z,
                   psharpfinalbeg, psharpend, rhofinal, pfinalbeg, pend,
                   H0, ε, U, ∇U, M, nleapfrog, lswfinal, α, control)

    if !validfinal
        return zfinalpr, zfinalpr, false,
        psharpbeg, psharpend, rho, pbeg, pend,
        nleapfrog, logsumweight, α
    end

    lswsubtree = logsumexp(lswinit, lswfinal)

    if lswfinal > lswsubtree
        zpr = zfinalpr
    else
        if rand(control[:rng]) < exp(lswfinal - lswsubtree)
            zpr = zfinalpr
        end
    end

    logsumweight = logsumexp(logsumweight, lswsubtree)

    rhosubtree = rhoinit + rhofinal
    rho += rhosubtree

    # Demand satisfaction around merged subtrees
    persist = stancriterion(psharpbeg, psharpend, rhosubtree)

    # Demand satisfaction between subtrees
    rhosubtree = rhoinit + pfinalbeg
    persist &= stancriterion(psharpbeg, psharpfinalbeg, rhosubtree)

    rhosubtree = rhofinal + pinitend
    persist &= stancriterion(psharpinitend, psharpend, rhosubtree)

    return z, zpr, persist,
    psharpbeg, psharpend, rho, pbeg, pend,
    nleapfrog, logsumweight, α
end

function leapfrog(z::PSPoint, ∇U, ε::Float64, M::Array{Float64})::PSPoint
    p_ = z.p - 0.5 * ε * ∇U(z.q)
    q = z.q + ε * rhosharp(p_, M)
    p = p_ - 0.5 * ε * ∇U(q)
    return PSPoint(q, p)
end

function leapfrog(z::PSPoint, ∇U, ε::Float64,
                  L::LowerTriangular{Float64,Matrix{Float64}})::PSPoint
    p_ = z.p - 0.5 * ε * (L' * ∇U(z.q))
    q = z.q + ε * rhosharp(p_, L)
    p = p_ - 0.5 * ε * (L' * ∇U(q))
    return PSPoint(q, p)
end

function hamiltonian(U, z::PSPoint, M::Array{Float64})::Float64
    return U(z.q) + 0.5 * (z.p' * rhosharp(z.p, M))
end

function hamiltonian(U, z::PSPoint, L::LowerTriangular{Float64,Matrix{Float64}})::Float64
    return U(z.q) + 0.5 * (z.p' * z.p)
end

function rhosharp(p::Vector{Float64}, M::Vector{Float64})::Vector{Float64}
    return M .* p
end

function rhosharp(p::Vector{Float64}, M::Matrix{Float64})::Vector{Float64}
    return M * p
end

function rhosharp(p::Vector{Float64}, L::LowerTriangular{Float64,Matrix{Float64}})::Vector{Float64}
    return L * p
end

function stancriterion(psharp_m::Vector{Float64}, psharp_p::Vector{Float64},
                       rho::Vector{Float64})::Bool
    return psharp_p' * rho > 0 && psharp_m' * rho > 0
end

function generatemomentum(rng, ndim::Int, M::Vector{Float64})::Vector{Float64}
    return randn(rng, ndim) ./ sqrt.(M)
end

function generatemomentum(rng, ndim::Int, M::Matrix{Float64})::Vector{Float64}
    return cholesky(Symmetric(M)).U \ randn(rng, ndim)
end

function generatemomentum(rng, ndim::Int,
                          L::LowerTriangular{Float64,Matrix{Float64}})::Vector{Float64}
    return randn(rng, ndim)
end

function findepsilon(ε::Float64, z::PSPoint, U, ∇U,
                     M::AbstractArray{Float64}, control)::Float64
    ndim = length(z.q)
    H0 = hamiltonian(U, z, M)

    zp = leapfrog(z, ∇U, ε, M)
    H = hamiltonian(U, zp, M)
    isnan(H) && (H = Inf)

    ΔH = H0 - H
    direction = ΔH > log(0.8) ? 1 : -1

    while true
        rp = generatemomentum(control[:rng], ndim, M)
        H0 = hamiltonian(U, z, M)

        zp = leapfrog(z, ∇U, ε, M)
        H = hamiltonian(U, zp, M)
        isnan(H) && (H = Inf)

        ΔH = H0 - H
        if direction == 1 && !(ΔH > log(0.8))
            break
        elseif direction == -1 && !(ΔH < log(0.8))
            break
        else
            ε = direction == 1 ? 2.0 * ε : 0.5 * ε
        end

        @assert ε <= 1.0e7 "Posterior is impropoer.  Please check your model."
        @assert ε >= 0.0 "No acceptable small step size could be found.  Perhaps the posterior is not continuous."
    end

    return ε
end

function adaptionwindows(control)

    # TODO(ear) check, Madapt * termpercent > 20
    openwindow::Int = floor(control[:iterations_warmup] *
        control[:adaptinitpercent])

    lastwindow::Int = ceil(control[:iterations_warmup] *
        (1 - control[:adapttermpercent]))

    return openwindow, openwindow + control[:adaptwindowsize],
    lastwindow, control[:adaptwindowsize]
end

function initializesample(ndim::Int, U, ∇U, control)::Vector{Float64}
    q = zeros(ndim)
    initialized = false
    a = 0

    while a < control[:initattempts] && !initialized
        q = generateuniform(control[:rng], length(q), control[:initradius])

        lq = U(q)
        if isfinite(lq) && !isnan(lq)
            initialized = true
        end

        gq = sum(∇U(q))
        if isfinite(gq) && !isnan(gq)
            initialized &= true
        end

        a += 1
    end

    @assert a <= control[:initattempts] && initialized "Failed to find inital values in $(initattempts) attempts."

    return q
end

function generateuniform(rng, N, radius)
    return rand(rng, N) .* (radius * 2) .- radius
end

function adaptstepsize(adaptstat, counter, sbar, xbar, μ, control)
    counter += 1
    adaptstat = adaptstat > 1 ? 1 : adaptstat
    eta = 1.0 / (counter + control[:t0])
    sbar = (1.0 - eta) * sbar + eta * (control[:δ] - adaptstat)
    x = μ - sbar * sqrt(counter) / control[:γ]
    xeta = counter ^ -control[:κ]
    xbar = (1.0 -  xeta) * xbar + xeta * x
    ε = exp(x)
    εbar = exp(xbar)
    return ε, εbar, counter, sbar, xbar
end

function log1pexp(a)
    if a > 0.0
        return a + log1p(exp(-a))
    end
    return log1p(exp(a))
end

function logsumexp(a, b)
    if a == -Inf
        return b
    end

    if a == Inf && b == Inf
        return Inf
    end

    if a > b
        return a + log1pexp(b - a)
    end

    return b + log1pexp(a - b)
end

abstract type AbstractWelfordState{S, T} end

struct WelfordState{S<:Vector{Float64}, T<:AbstractArray{Float64}} <: AbstractWelfordState{S, T}
    m::S
    s::T
    n::Int
end

function accmoments(ws::AbstractWelfordState, x::Vector{Float64})::WelfordState
    n = ws.n + 1
    d = similar(ws.m)
    M = similar(ws.m)
    @. d = x - ws.m
    @. M = ws.m + d / n
    S = updateCov(ws.s, d, x .- M)
    return WelfordState(M, S, n)
end

function updateCov(s::Vector{Float64}, d::Vector{Float64}, dp::Vector{Float64})
    return @. s + d * dp
end

function updateCov(s::AbstractArray{Float64}, d::Vector{Float64}, dp::Vector{Float64})
    return s .+ dp * d'
end

function samplevariance(ws::AbstractWelfordState{Vector{Float64}, Vector{Float64}},
                        control)
    if ws.n > 1
        σ = similar(ws.s)
        @. σ = ws.s / (ws.n - 1)
        if control[:regularize]
            w = ws.n / (ws.n + 5.0)
            σ = w * σ + 1e-3 * (1 - w) * ones(length(σ))
        end
        return σ
    end
    return ones(length(ws.m))
end

function samplevariance(ws::AbstractWelfordState{Vector{Float64}, Matrix{Float64}},
                        control)
    if ws.n > 1
        σ = similar(ws.s)
        @. σ = ws.s / (ws.n - 1)
        if control[:regularize]
            w = ws.n / (ws.n + 5.0)
            σ = w * σ + 1e-3 * (1 - w) * one(σ)
        end
        if control[:skewsymmetric]
            return cholesky(Symmetric(σ)).L
        end
        return σ
    end
    return one(ws.s)
end

max_stars_repo_name: roualdes/yastan

max_issues_repo_name: roualdes/yastan

max_forks_repo_name: roualdes/yastan


text:
module Simulation

include("Gates.jl")

end  # module Simulatio


text:
export l_bfgs_rcst
export Steighaug
export BhaskaraTop
export LimitedMemory

"""
l_bfgs_rcst(nlp; options...) - Ainda não chegamos nesse!

Este método é chamado L-BFGS com região de confiança por Steihaug-Toint.

Tenta-se resolver B_k d = - ∇f(xₖ) usando Gradientes Conjugados.
Se em algum momento a direção ficar maior que a região de confiança,
ela é truncada. Bₖ é a aproximação de memória limitada de BFGS.
Porém, deve-se guardar apenas os últimos p vetores sₖ e yₖ,
e não criar a matriz Bₖ explicitamente.

Options:
- atol: absolute tolerance for the first order condition (default: 1e-6)
- rtol: relative tolerance for the first order condition (default: 1e-6)
- max_eval: maximum number of [functions] evaluations, use ≤ 0 for unlimited (default: 1000)
- max_iter: maximum number of iterations, use ≤ 0 for unlimited (default: 0)
- max_time: maximum elapsed time in seconds, use ≤ 0 for unlimited (default: 10)

Disclaimers for the developer:
  - nlp should be the only mandatory argument
  - these five options are the current default for other JSO-compliant solvers
  - always return a GenericExecutionStats
  """

  function Steighaug(gx, B, Δ)
    # ϵ = 1e-4
    ϵ = min(0.5,(norm(gx))^0.5)*norm(gx) # tolerância sugerida Algor 7.1
    m = length(gx)
    zx = zeros(m)
    r = gx
    d = -r

    rx = r
    z = zx
    normr = norm(r)
    if normr < ϵ
      return z
    end

    # enquanto estiver dentro da região de confiança
    k=0
    while normr > ϵ || k < m # CG vai em no máximo m-variáveis direções
      dotdBd = dot(d,B*d)
      if dotdBd ≤ 0
        # pára o método se a direção dⱼ é direção de curvatura não positiva
        # encontrar τ >=0 tal que pk e satisfaz ||pk|| = Δk
        # ou seja, a interseção da direção com a região de confiança
        m1, m2 = BhaskaraTop(z, d, Δ)
        return z + m1*d
      end

      dotrr = dot(r,r)
      α = dotrr/dotdBd
      zx = z + α*d
      if norm(zx) ≥ Δ
        #para se zⱼ₊₁ viola os limites da região de confiança
        # encontrar τ >=0 tal que pk e satisfaz ||pk|| = Δk
        # ou seja, a interseção da direção com a região de confiança
        m1, m2 = BhaskaraTop(z, d, Δ)
        return z + m1*d
      end

      rx = r + α*B*d #Conjugate Gradient
      #se α = 0 => rx = r , β = 1 (logo abaixo)

      if norm(rx) < ϵ
        return zx
      end
      β = dot(rx,rx)/dotrr
      d = -rx+ β*d
      r = rx
      normr = norm(r)
      z = zx
      k+=1
    end
    if k >= m
      #@error ("não encontrou a borda, direção multipla da própria direção")
      return zeros(m)
    end

  end

  function BhaskaraTop(z, d, Δ)
    a=dot(d,d)
    b=2dot(z,d)
    c=dot(z,z) - Δ^2
    Delta = b^2-4*a*c
    if Delta < 0
      @warn("Δ<0")
    else
      t1 = (-b+sqrt(Delta))/2a
      t2 = (-b-sqrt(Delta))/2a
      m1 = max(t1,t2)
      m2 = min(t1,t2)
      return m1, m2
    end
  end

  # L-BGFS-RC-Steighaug-Toint
  function l_bfgs_rcst(
    nlp::AbstractNLPModel;
    atol::Real = 1e-6,
    rtol::Real = 1e-6,
    max_eval::Int = 1000,
    max_iter::Int = 0,
    max_time::Float64 = 10.0
    )

    if !unconstrained(nlp)
      error("Problem is not unconstrained")
    end

    #Given the starting point
    x = copy(nlp.meta.x0)
    f(x) = obj(nlp, x)
    ∇f(x) = grad(nlp, x)

    fx = f(x)
    ∇fx = ∇f(x)
    # initial Hessian approximation
    n = length(x)

    #trust region radius
    Δ = 1.0

    B = Matrix(1.0*I, n, n)
    s = Steighaug(∇fx, B, Δ)
    y = ∇f(x.+s) - ∇fx
    γ = dot(y,y)/dot(s,y) #ver Nocedal p. 178
    B = Matrix(γ*I, n, n)

    # Criar matrizes para armazenar s e y
    if n < 5
      m_vetores = n
    else
      m_vetores = 5
    end

    S = zeros(n, m_vetores)
    Y = zeros(n, m_vetores)

    #convergence tolerance
    ϵ = atol + rtol * norm(∇fx)

    #parameters η e r
    η = 0.01 #∈ (0,1e-3)
    r = 0.001 #∈ (0,1)

    t₀ = time()

    iter = 0
    Δt = time() - t₀
    solved = norm(∇fx) < ϵ # First order stationary
    tired = neval_obj(nlp) ≥ max_eval > 0 ||
                      iter ≥ max_iter > 0 ||
                      Δt ≥ max_time > 0

    status = :unknown
    # log_header is up for some rewrite in the future. For now, it simply prints the column names with some spacing
    @info log_header(
    [:iter, :fx, :ngx, :nf, :Δt],
    [Int, Float64, Float64, Int, Float64],
    hdr_override=Dict(:fx => "f(x)", :ngx => "‖∇f(x)‖", :nf => "#f")
    )
    # log_row uses the type information of each value, thus we use `Any` here.
    @info log_row(
    Any[iter, fx, norm(∇fx), neval_obj(nlp), Δt]
    )

    # Aqui começa o show

    while !(solved || tired)
      #compute sₖ by solving the subproblem
      s = Steighaug(∇fx, B, Δ)
      y = ∇f(x.+s) - ∇fx
      if norm(y) < 10e-18
        status=:small_step
      end
      ared = fx - f(x.+s)
      pred = -(dot(∇fx, s) + 1/2 * dot(s, B*s))

      ρ = ared/pred
      if ρ < η
        Δ = Δ/2
        if Δ < 10e-50
          status =:user
        end
      else
        x = x + s
        fx = f(x)
        ∇fx = ∇f(x)
        if ρ > 0.75 && norm(s) > 0.8 * Δ
          Δ = 2*Δ
          if Δ > 10e50 #evita que o raio aumente muito
            status =:not_desc
          end
        end
      end

      if status != :unknown #small_step
        break
      end
      if iter < m_vetores
        yBs = y .-B*s
        if abs(dot(s,yBs)) > r*norm(s,2)*norm(yBs,2)
          #6.26 - evita atualizações se o denominador é pequeno
          B = B + (yBs*yBs')/dot(yBs,s)
        end
        S[:,iter+1] = s
        Y[:,iter+1] = y
      else
        yBs = y .-B*s
        if abs(dot(s,yBs)) > r*norm(s,2)*norm(yBs,2) #evita situações em que são ambas zero
          #6.26 - evita atualizações se o denominador é pequeno
          #aqui tem que ter mágica - LBFGS
          δₖ = dot(y,y)/dot(s, y)
          B = LimitedMemory(δₖ, S, Y, iter, n, m_vetores)
        end

        cpatu = mod(iter, m_vetores) + 1
        S[:,cpatu] = s
        Y[:,cpatu] = y
      end

      iter+= 1

      Δt = time() - t₀
      solved = norm(∇fx) < ϵ # First order stationary
      tired = neval_obj(nlp) ≥ max_eval > 0|| iter ≥ max_iter > 0 ||
                          Δt ≥ max_time > 0 # Excess time, iteration, evaluations

      @info log_row(
      Any[iter, fx, norm(∇fx), neval_obj(nlp), Δt]
      )
    end

    if solved
      status = :first_order
    elseif tired
      if neval_obj(nlp) ≥ max_eval > 0
        status = :max_eval
      elseif iter ≥ max_iter > 0
        status = :max_iter
      elseif Δt ≥ max_time > 0
        status = :max_time
      end
    end

    return GenericExecutionStats(
    status,
    nlp,
    solution=x,
    objective=f(x),
    dual_feas=norm(∇fx),
    elapsed_time=Δt,
    iter=iter
    )
  end

  function LimitedMemory(δₖ, S, Y, iter, n, m)
    Dₖ = zeros(m,m)
    Lₖ = zeros(m,m)

    B1 = δₖ.*I(n)
    inicio = mod(iter, m) + 1
    Sₖ= [S[:, inicio:m] S[:, 1:inicio-1]]
    Yₖ = [Y[:, inicio:m] Y[:, 1:inicio-1]]
    # Sₖ= [S[:, inicio-1:-1:1] S[:, m:-1:inicio] ] #sey_inverso
    # Yₖ = [Y[:, inicio-1:-1:1] Y[:, m:-1:inicio] ]

    δS = δₖ.*Sₖ
    B2 = [δS Yₖ]
    for i=1:m, j=1:m
      if i>j
        Lₖ[i,j] = dot(S[:,i], Y[:,j])
      end
      if i==j
        Dₖ[i,j] = dot(S[:,i], Y[:,j])
      end
    end
    # B3 = inv([δS'*Sₖ Lₖ; Lₖ' -Dₖ])

    B4 = [δS';Yₖ']

    return B1 - B2*(([δS'*Sₖ Lₖ; Lₖ' -Dₖ])\B4)

  end

  function Unrolling(B, S, Y, iter, n, m)
    #Unrolling the BFGS formula - Nocedal p.184
    inicio = mod(iter, m) + 1
    Sₖ= [S[:, inicio:m] S[:, 1:inicio-1]]
    Yₖ = [Y[:, inicio:m] Y[:, 1:inicio-1]]
    # Sₖ= [S[:, inicio-1:-1:1] S[:, m:-1:inicio] ] #sey_inverso
    # Yₖ = [Y[:, inicio-1:-1:1] Y[:, m:-1:inicio] ]
    a = zeros(n,m)
    b = zeros(n,m)
    for i = 1:m
      b[:,i] = Y[:,i]/sqrt(dot(Y[:,i], S[:,i]))
      a[:,i] = B*S[:,i] + sum(dot(b[:,j],S[:,i])*b[:,j] -
                                    dot(a[:,j],S[:,i])*a[:,j]
                                    for j=1:m)
    a[:,i] = a[:,i]/sqrt(dot(S[:,i], a[:,i]))
    end
    return B + sum(b[:,i]*b[:,i]' - a[:,i]*a[:,i]' for i=1:m)
  end


text:
@recipe function dummy(curve::ROC)
  x, y = get_curve(curve)
	ticks = [0.2, 0.4, 0.6, 0.8, 1]
	xlims := (0,1)
	ylims := (0,1)
	xguide := "False Positive Rate (%)"
	yguide := "True Positive Rate (%)"
	title --> "Receiver Operator Characteristic (ROC)"
	legend --> :outerleft
	ticks --> (ticks, string.(ticks.*100))
	@series begin
		seriescolor --> :black
		linestyle --> :dash
		label := ""
		[0, 1], [0, 1]
	end
  EER, idx = eer(curve)
	@series begin
		seriescolor --> :red
    seriestype := :scatter
		label := ""
    [x[idx+1]], [y[idx+1]]
	end
	@series begin
    x, y
	end
end

@recipe function dummy(curve::PR)
  x, y = get_curve(curve)
	ticks = [0.2, 0.4, 0.6, 0.8, 1]
	xlims := (0,1)
	ylims := (0,1)
	xguide := "Recall (%)"
	yguide := "Precision (%)"
	title --> "Precision Recall (PR)"
	legend --> :outerleft
	ticks --> (ticks, string.(ticks.*100))
	@series begin # this is just to get the same colors with ROC
		seriescolor --> :black
		linestyle --> :dash
		label := ""
		[0, 0], [0, 0]
	end
  EER, idx = eer(curve)
	@series begin
		seriescolor --> :red
    seriestype := :scatter
		label := ""
    [x[idx+1]], [y[idx+1]]
	end
	@series begin
    x, y
	end
end

@recipe function dummy(curve::DET)
  x, y = get_curve(curve)
	ticks = [0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 40, 60]
	xguide := "False Positive Rate (%)"
	yguide := "False Negative Rate (%)"
	title --> "Detection Error Tradeoff (DET)"
	legend --> :outerleft
	ticks --> (det_axis.( ticks ./ 100), string.(ticks))
	xlims --> (det_axis(0.001), det_axis(0.6))
	ylims --> (det_axis(0.001), det_axis(0.6))
	@series begin # this is just to get the same colors with ROC
		seriescolor --> :black
		linestyle --> :dash
		label := ""
		[0, 0], [0, 0]
	end
  EER, idx = eer(curve)
	@series begin
    seriestype := :scatter
		seriescolor --> :red
		label := ""
    [x[idx]], [y[idx]]
	end
	@series begin
		x, y
	end
end


text:
using ColorSchemes, CairoMakie, FileIO, JuMP
using ProgressMeter

include(joinpath("src", "EMP-model", "model.jl"))
import .GlnModel

imgpath = joinpath("docs", "imgs", "EMP-model")

#=
Warm up problem to get all the symbols in the model.
Create arrays using these symbols as the names.
This is the fast way of doing:
mu = Float64[]
atp = Float64[]
etc.
=#
test_model = GlnModel.gln_model();
syms = Symbol.(all_variables(test_model))
for sym in syms
    eval(:($sym=Float64[]))
end

prev_sol=Dict{Symbol, Float64}() # a variable to store the previous solution

ratio_lb = 0.05
ratio_ub = 10
@showprogress for var in reverse(range(log(ratio_lb); stop=log(ratio_ub), length=20))
    nadh_nad_ratio = exp(var)

    gln_model = GlnModel.gln_model(
        glc_ext = 50e-3,
        lac_ext = 1e-4,
        nh4_ext = 10e-3,
        ac_ext = 1e-16,
        etoh_ext = 1e-4,
        # atp_adp_ratio = -1, # unconstrained
        nadh_nad_ratio = nadh_nad_ratio,
        num_ms = 10,
    )

    solve_1 = GlnModel.max_mu!(gln_model, prev_sol)

    if !solve_1
        continue # if mu not found, skip
    end

    solve_2 = GlnModel.find_nearest!(gln_model, prev_sol) # find a solution close to the previous one, smoothen things

    if solve_1 && solve_2
        for sym in syms
            push!(eval(sym), value(gln_model[sym]))
        end
        # update old previous solution
        GlnModel.update!(gln_model, prev_sol)
    end
end

function plot_proteome(x, fracs, fraclabels; legendlabel="",xlabel="", ylabel="", xscale=identity)
    # NB: can't display more than 10 and not less than 3
    z = cumsum(fracs'; dims = 1) # no zeros
    z = [zeros(size(z, 2))'; z ./ z[end, :]']
    N = length(fraclabels)
    cs = eval(Meta.parse("ColorSchemes.Spectral_$N"))

    f = Figure()
    ax = Axis(f[1, 1], xscale=xscale, xlabel=xlabel, ylabel=ylabel)

    for i = 1:N
        band!(ax, x, z[i, :], z[i+1, :], label = fraclabels[i], color=cs[i])
    end

    f[1, 2] = Legend(f, ax, legendlabel, unique = true, framevisible = false)

    return f
end

#=
Plot cofactor metabolites
=#
fracs = exp.(hcat([adp, atp, nad, nadh]...))
fraclabels = ["adp", "atp", "nad", "nadh"]
x = exp.(nadh)./exp.(nad)

f = plot_proteome(x, fracs, fraclabels;
    legendlabel="Metabolites",
    xlabel="NADH/NAD ratio",
    ylabel="Relative abundance",
    xscale=log10,
)

FileIO.save(joinpath(imgpath, "Cofactor_Metabolites_NADH_NAD.pdf"), f)

#=
Plot upper glycolysis metabolites
=#
fracs = exp.(hcat([g6p, f6p, fdp, dhap, g3p, dpg13]...))
fraclabels = ["g6p", "f6p", "fdp", "dhap", "g3p", "dpg13"]
x = exp.(nadh)./exp.(nad)

f = plot_proteome(x, fracs, fraclabels;
    legendlabel="Metabolites",
    xlabel="NADH/NAD ratio",
    ylabel="Relative abundance",
    xscale=log10,
)

FileIO.save(joinpath(imgpath, "Upper_Glycolysis_Metabolites_NADH_NAD.pdf"), f)

#=
Plot lower glycolysis metabolites
=#
fracs = exp.(hcat([pg3, pg2, pep]...))
fraclabels = ["pg3", "pg2", "pep"]
x = exp.(nadh)./exp.(nad)

f = plot_proteome(x, fracs, fraclabels;
    legendlabel="Metabolites",
    xlabel="NADH/NAD ratio",
    ylabel="Relative abundance",
    xscale=log10,
)

FileIO.save(joinpath(imgpath, "Lower_Glycolysis_Metabolites_NADH_NAD.pdf"), f)

#=
Plot Gibbs dissipation of reactions in upper glycolysis
=#
fracs = hcat([dg_pts.*v_pts, dg_pgi .* v_pgi, dg_pfk .* v_pfk, dg_fba .*v_fba, dg_tpi .* v_tpi, dg_gapd .* v_gapd]...)
# fracs = hcat([dg_pts, dg_pgi, dg_pfk, dg_fba, dg_tpi, dg_gapd]...)
fraclabels = ["pts", "pgi", "pfk", "fba", "tpi", "gapd"]
x = exp.(nadh)./exp.(nad)

f = plot_proteome(x, fracs, fraclabels;
    legendlabel="Reaction",
    xlabel="NADH/NAD ratio",
    ylabel="Relative Gibbs dissipation",
    xscale=log10,
)

FileIO.save(joinpath(imgpath, "Upper_Glycolysis_Gibbs_Dissapation_NADH_NAD.pdf"), f)

#=
Plot upper glycolysis proteome
=#
fracs = hcat([pts, pgi, pfk, fba, tpi, gapd]...)
fraclabels = ["pts", "pgi", "pfk", "fba", "tpi", "gapd"]
x = exp.(nadh)./exp.(nad)

f = plot_proteome(x, fracs, fraclabels;
    legendlabel="Enzyme",
    xlabel="NADH/NAD ratio",
    ylabel="Relative Enzyme concentration",
    xscale=log10,
)

FileIO.save(joinpath(imgpath, "Upper_Glycolysis_Proteome_NADH_NAD.pdf"), f)

#=
Plot lower glycolysis proteome
=#
fracs = hcat([pgk, pgm, pyk, eno, ldh]...)
fraclabels = ["pgk", "pgm", "pyk", "eno", "ldh"]
x = exp.(nadh)./exp.(nad)

f = plot_proteome(x, fracs, fraclabels;
    legendlabel="Enzyme",
    xlabel="NADH/NAD ratio",
    ylabel="Relative Enzyme concentration",
    xscale=log10,
)

FileIO.save(joinpath(imgpath, "Lower_Glycolysis_Proteome_NADH_NAD.pdf"), f)

#=
Plot glutamine synthesis proteome
=#
fracs = hcat([ppc, pfl, pdh, cs, aconta, acontb, icdh, gludy, glns]...)
fraclabels = ["ppc", "pfl", "pdh", "cs", "aconta", "acontb", "icdh", "gludy", "glns"]
x = exp.(nadh)./exp.(nad)

f = plot_proteome(x, fracs, fraclabels;
    legendlabel="Enzyme",
    xlabel="NADH/NAD ratio",
    ylabel="Relative Enzyme concentration",
    xscale=log10,
)

FileIO.save(joinpath(imgpath, "Glutamine_synthesis_Proteome_NADH_NAD.pdf"), f)

#=
Plot remaining proteome (transporters; ethanol -> accoa; acetate -> accoa)
=#
fracs = hcat([lact, nh4t, fort, etoht, alcd, acald, act, ackr, ptar]...)
fraclabels = ["lact", "nh4t", "fort", "etoht", "alcd", "acald", "act", "ackr", "ptar"]
x = exp.(nadh)./exp.(nad)

f = plot_proteome(x, fracs, fraclabels;
    legendlabel="Enzyme",
    xlabel="NADH/NAD ratio",
    ylabel="Relative Enzyme concentration",
    xscale=log10,
)

FileIO.save(joinpath(imgpath, "Remaining_Proteome_NADH_NAD.pdf"), f)

#=
Plot Fluxes of reactions in upper glycolysis
=#
fracs = hcat([v_pts,v_pgi, v_pfk, v_fba, v_tpi, v_gapd]...)
fraclabels = ["v_pts", "v_pgi", "v_pfk", "v_fba", "v_tpi", "v_gapd"]
x = exp.(nadh)./exp.(nad)

f = plot_proteome(x, fracs, fraclabels;
    legendlabel="Reaction",
    xlabel="NADH/NAD ratio",
    ylabel="Relative Flux",
    xscale=log10,
)

FileIO.save(joinpath(imgpath, "Upper_Glycolysis_Fluxes_NADH_NAD.pdf"), f)

#=
Plot Fluxes of reactions in glutamine synthesis
=#
fracs = hcat([v_ppc, v_pfl, v_pdh, v_cs, v_aconta, v_acontb, v_icdh, v_gludy, v_glns]...)
fraclabels = ["v_ppc", "v_pfl", "v_pdh", "v_cs", "v_aconta", "v_acontb", "v_icdh", "v_gludy", "v_glns"]
x = exp.(nadh)./exp.(nad)

f = plot_proteome(x, fracs, fraclabels;
    legendlabel="Reaction",
    xlabel="NADH/NAD ratio",
    ylabel="Relative Flux",
    xscale=log10,
)

FileIO.save(joinpath(imgpath, "Glutamine_synthesis_Fluxes_NADH_NAD.pdf"), f)

#=
Plot dG of reactions in upper glycolysis
=#

fracs = hcat([dg_pts, dg_pgi, dg_pfk, dg_fba, dg_tpi, dg_gapd]...)
fraclabels = ["dG_pts", "dG_pgi", "dG_pfk", "dG_fba", "dG_tpi", "dG_gapd"]
x = exp.(nadh)./exp.(nad)

f = plot_proteome(x, fracs, fraclabels;
    legendlabel="Reaction",
    xlabel="NADH/NAD ratio",
    ylabel="Relative dG of reaction",
    xscale=log10,
)

FileIO.save(joinpath(imgpath, "Upper_Glycolysis_dG_NADH_NAD.pdf"), f)

#=
Plot dG of reactions in glycolysis synthesis
=#
fracs = hcat([dg_ppc, dg_pfl, dg_pdh, dg_cs, dg_aconta, dg_acontb, dg_icdh, dg_gludy, dg_glns]...)
fraclabels = ["dg_ppc", "dg_pfl", "dg_pdh", "dg_cs", "dg_aconta", "dg_aconta", "dg_icdh", "dg_gludy", "dg_glns"]
x = exp.(nadh)./exp.(nad)

f = plot_proteome(x, fracs, fraclabels;
    legendlabel="Reaction",
    xlabel="NADH/NAD ratio",
    ylabel="Relative dG of reaction",
    xscale=log10,
)

FileIO.save(joinpath(imgpath, "Glutamine_synthesis_dG_NADH_NAD.pdf"), f)

### Absolte plots ###

#=
Plot cofactor metabolites
=#
x = exp.(nadh)./exp.(nad)
a = exp.([adp]...)
b = exp.([atp]...)
c = exp.([nad]...)
d = exp.([nadh]...)

f = Figure()
ax = Axis(f[1, 1])
scatter1 = scatter!(ax, x, a)
line1 = lines!(ax, x, a)
scatter2 = scatter!(ax, x, b)
line2 = lines!(ax, x, b)
scatter3 = scatter!(ax, x, c)
line3 = lines!(ax, x, c)
scatter4 = scatter!(ax, x, d)
line4 = lines!(ax, x, d)

ax.xlabel = "NADH/NAD ratio"
ax.ylabel = "Metabolite concentration [M]"
xscale = log10
f[1, 2] = Legend(
    f,
    legendlabel="Metabolite",
    [[scatter1, line1], [scatter2, line2], [scatter3, line3], [scatter4, line4]],
    ["ADP", "ATP", "NAD", "NADH"],
)
f

FileIO.save(joinpath(imgpath, "Cofactor_Metabolites_NADH_NAD_abs.pdf"), f)

#=
Plot upper glycolysis metabolites
=#
x = exp.(nadh)./exp.(nad)
a = exp.([g6p]...)
b = exp.([f6p]...)
c = exp.([fdp]...)
d = exp.([dhap]...)
e = exp.([g3p]...)
g = exp.([dpg13]...)

f = Figure()
ax = Axis(f[1, 1])
scatter1 = scatter!(ax, x, a)
line1 = lines!(ax, x, a)
scatter2 = scatter!(ax, x, b)
line2 = lines!(ax, x, b)
scatter3 = scatter!(ax, x, c)
line3 = lines!(ax, x, c)
scatter4 = scatter!(ax, x, d)
line4 = lines!(ax, x, d)
scatter5 = scatter!(ax, x, e)
line5 = lines!(ax, x, e)
scatter6 = scatter!(ax, x, g)
line6 = lines!(ax, x, g)

ax.xlabel = "NADH/NAD ratio"
ax.ylabel = "Metabolite concentration [M]"
xscale = log10
f[1, 2] = Legend(
    f,
    legendlabel="Metabolite",
    [[scatter1, line1], [scatter2, line2], [scatter3, line3], [scatter4, line4], [scatter5, line5], [scatter6, line6]],
    ["g6p", "f6p", "fdp", "dhap", "g3p", "dpg13"],
)
f

FileIO.save(joinpath(imgpath, "Upper_Glycolysis_Metabolites_NADH_NAD_abs.pdf"), f)

#=
Plot lower glycolysis metabolites
=#
x = exp.(nadh)./exp.(nad)
a = exp.([pg3]...)
b = exp.([pg2]...)
c = exp.([pep]...)

f = Figure()
ax = Axis(f[1, 1])
scatter1 = scatter!(ax, x, a)
line1 = lines!(ax, x, a)
scatter2 = scatter!(ax, x, b)
line2 = lines!(ax, x, b)
scatter3 = scatter!(ax, x, c)
line3 = lines!(ax, x, c)

ax.xlabel = "NADH/NAD ratio"
ax.ylabel = "Metabolite concentration [M]"
xscale = log10
f[1, 2] = Legend(
    f,
    legendlabel="Metabolite",
    [[scatter1, line1], [scatter2, line2], [scatter3, line3]],
    ["p3g", "pg2", "pep"],
)
f

FileIO.save(joinpath(imgpath, "Lower_Glycolysis_Metabolites_NADH_NAD_abs.pdf"), f)

#=
Plot upper glycolysis proteome
=#
x = exp.(nadh)./exp.(nad)
a = (pts)
b = (pgi)
c = (pfk)
d = (fba)
e = (tpi)
g = (gapd)

f = Figure()
ax = Axis(f[1, 1])
scatter1 = scatter!(ax, x, a)
line1 = lines!(ax, x, a)
scatter2 = scatter!(ax, x, b)
line2 = lines!(ax, x, b)
scatter3 = scatter!(ax, x, c)
line3 = lines!(ax, x, c)
scatter4 = scatter!(ax, x, d)
line4 = lines!(ax, x, d)
scatter5 = scatter!(ax, x, e)
line5 = lines!(ax, x, e)
scatter6 = scatter!(ax, x, g)
line6 = lines!(ax, x, g)

ax.xlabel = "NADH/NAD ratio"
ax.ylabel = "Enzyme concentration [g enz / g DW]"
xscale = log10
f[1, 2] = Legend(
    f,
    legendlabel="Enzyme",
    [[scatter1, line1], [scatter2, line2], [scatter3, line3], [scatter4, line4], [scatter5, line5], [scatter6, line6]],
    ["pts", "pgi", "pfk", "fba", "tpi", "gapd"],
)
f

FileIO.save(joinpath(imgpath, "Upper_Glycolysis_Proteome_NADH_NAD_abs.pdf"), f)

#=
Plot lower glycolysis proteome
=#
x = exp.(nadh)./exp.(nad)
a = (pgk)
b = (pgm)
c = (pyk)
d = (eno)
e = (ldh)

f = Figure()
ax = Axis(f[1, 1])
scatter1 = scatter!(ax, x, a)
line1 = lines!(ax, x, a)
scatter2 = scatter!(ax, x, b)
line2 = lines!(ax, x, b)
scatter3 = scatter!(ax, x, c)
line3 = lines!(ax, x, c)
scatter4 = scatter!(ax, x, d)
line4 = lines!(ax, x, d)
scatter5 = scatter!(ax, x, e)
line5 = lines!(ax, x, e)

ax.xlabel = "NADH/NAD ratio"
ax.ylabel = "Enzyme concentration [g enz / g DW]"
xscale = log10
f[1, 2] = Legend(
    f,
    legendlabel="Enzyme",
    [[scatter1, line1], [scatter2, line2], [scatter3, line3], [scatter4, line4], [scatter5, line5]],
    ["pgk", "pgm", "pyk", "eno", "ldh"],
)
f

FileIO.save(joinpath(imgpath, "Lower_Glycolysis_Proteome_NADH_NAD_abs.pdf"), f)

#=
Plot Fluxes of reactions in upper glycolysis
=#
x = exp.(nadh)./exp.(nad)
a = (v_pts)
b = (v_pgi)
c = (v_pfk)
d = (v_fba)
e = (v_tpi)
g = (v_gapd)

f = Figure()
ax = Axis(f[1, 1])
scatter1 = scatter!(ax, x, a)
line1 = lines!(ax, x, a)
scatter2 = scatter!(ax, x, b)
line2 = lines!(ax, x, b)
scatter3 = scatter!(ax, x, c)
line3 = lines!(ax, x, c)
scatter4 = scatter!(ax, x, d)
line4 = lines!(ax, x, d)
scatter5 = scatter!(ax, x, e)
line5 = lines!(ax, x, e)
scatter6 = scatter!(ax, x, g)
line6 = lines!(ax, x, g)

ax.xlabel = "NADH/NAD ratio"
ax.ylabel = "Fluxes [mmol/gDW/h]"
xscale = log10
f[1, 2] = Legend(
    f,
    legendlabel="Flux",
    [[scatter1, line1], [scatter2, line2], [scatter3, line3], [scatter4, line4], [scatter5, line5], [scatter6, line6]],
    ["v_pts", "v_pgi", "v_pfk", "v_fba", "v_tpi", "v_gapd"],
)
f

FileIO.save(joinpath(imgpath, "Upper_Glycolysis_Fluxes_NADH_NAD_abs.pdf"), f)

#=
Plot dG of reactions in upper glycolysis
=#
x = exp.(nadh)./exp.(nad)
a = (dg_pts)
b = (dg_pgi)
c = (dg_pfk)
d = (dg_fba)
e = (dg_tpi)
g = (dg_gapd)

f = Figure()
ax = Axis(f[1, 1])
scatter1 = scatter!(ax, x, a)
line1 = lines!(ax, x, a)
scatter2 = scatter!(ax, x, b)
line2 = lines!(ax, x, b)
scatter3 = scatter!(ax, x, c)
line3 = lines!(ax, x, c)
scatter4 = scatter!(ax, x, d)
line4 = lines!(ax, x, d)
scatter15 = scatter!(ax, x, e)
line5 = lines!(ax, x, e)
scatter6 = scatter!(ax, x, g)
line6 = lines!(ax, x, g)

ax.xlabel = "NADH/NAD ratio"
ax.ylabel = "dG of reaction"
xscale = log10
f[1, 2] = Legend(
    f,
    legendlabel="dG of reaction",
    [[scatter1, line1], [scatter2, line2], [scatter3, line3], [scatter4, line4],[scatter5, line5], [scatter6, line6]],
    ["dg_pts", "dg_pgi", "dg_pfk", "dg_fba", "dg_tpi", "dg_gapd"],
)
f

FileIO.save(joinpath(imgpath, "Upper_Glycolysis_dG_NADH_NAD_abs.pdf"), f)


text:
# visualization
function visualize!(vis, model::DoublePendulum, x;
        color=RGBA(0.0, 0.0, 0.0, 1.0),
        r = 0.1, Δt = 0.1)

    i = 1
    l1 = Cylinder(Point3f0(0.0, 0.0, 0.0), Point3f0(0.0, 0.0, model.l1),
        convert(Float32, 0.025))
    setobject!(vis["l1$i"], l1, MeshPhongMaterial(color = color))
    l2 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l2),
        convert(Float32, 0.025))
    setobject!(vis["l2$i"], l2, MeshPhongMaterial(color = color))

    setobject!(vis["elbow$i"], Sphere(Point3f0(0.0),
        convert(Float32, 0.05)),
        MeshPhongMaterial(color = color))
    setobject!(vis["ee$i"], Sphere(Point3f0(0.0),
        convert(Float32, 0.05)),
        MeshPhongMaterial(color = color))

    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

    T = length(x)
    for t = 1:T

        MeshCat.atframe(anim,t) do
            p_mid = [kinematics_elbow(model, x[t])[1], 0.0, kinematics_elbow(model, x[t])[2]]
            p_ee = [kinematics(model, x[t])[1], 0.0, kinematics(model, x[t])[2]]

            settransform!(vis["l1$i"], cable_transform(zeros(3), p_mid))
            settransform!(vis["l2$i"], cable_transform(p_mid, p_ee))

            settransform!(vis["elbow$i"], Translation(p_mid))
            settransform!(vis["ee$i"], Translation(p_ee))
        end
    end

    # settransform!(vis["/Cameras/default"],
    #    compose(Translation(0.0 , 0.0 , 0.0), LinearMap(RotZ(pi / 2.0))))

    MeshCat.setanimation!(vis, anim)
end

# # visualization
# function visualize_elbow!(vis, model::DoublePendulum, x;
#         color=RGBA(0.0, 0.0, 0.0, 1.0),
#         r = 0.1, Δt = 0.1)
#
#     i = 1
#     l1 = Cylinder(Point3f0(0.0, 0.0, 0.0), Point3f0(0.0, 0.0, model.l1),
#         convert(Float32, 0.025))
#     setobject!(vis["l1"], l1, MeshPhongMaterial(color = color))
#     l2 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l2),
#         convert(Float32, 0.025))
#     setobject!(vis["l2"], l2, MeshPhongMaterial(color = color))
#
#     setobject!(vis["elbow_nominal"], Sphere(Point3f0(0.0),
#         convert(Float32, 0.05)),
#         MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))
# 	setobject!(vis["elbow_limit"], Sphere(Point3f0(0.0),
#         convert(Float32, 0.05)),
#         MeshPhongMaterial(color = RGBA(0.0, 1.0, 0.0, 1.0)))
#     setobject!(vis["ee"], Sphere(Point3f0(0.0),
#         convert(Float32, 0.05)),
#         MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))
#
#     anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))
#
# 	ϵ = 1.0e-1
#     T = length(x)
#     for t = 1:T
#
#         MeshCat.atframe(anim,t) do
#             p_mid = [kinematics_elbow(model, x[t])[1], 0.0, kinematics_elbow(model, x[t])[2]]
#             p_ee = [kinematics(model, x[t])[1], 0.0, kinematics(model, x[t])[2]]
#
#             settransform!(vis["l1"], cable_transform(zeros(3), p_mid))
#             settransform!(vis["l2"], cable_transform(p_mid, p_ee))
#
#             settransform!(vis["elbow_nominal"], Translation(p_mid))
# 			settransform!(vis["elbow_limit"], Translation(p_mid))
#             settransform!(vis["ee"], Translation(p_ee))
#
# 			if x[t][2] <= -0.5 * π + ϵ || x[t][2] >= 0.5 * π - ϵ
# 				setvisible!(vis["elbow_nominal"], false)
# 				setvisible!(vis["elbow_limit"], true)
# 			else
# 				setvisible!(vis["elbow_nominal"], true)
# 				setvisible!(vis["elbow_limit"], false)
# 			end
#         end
#     end
#
#     # settransform!(vis["/Cameras/default"],
#     #    compose(Translation(0.0 , 0.0 , 0.0), LinearMap(RotZ(pi / 2.0))))
#
#     MeshCat.setanimation!(vis, anim)
# end

# visualization
function _create_acrobot!(vis, model::DoublePendulum;
        tl = 1.0,
        limit_color = RGBA(0.0, 1.0, 0.0, tl),
        color = RGBA(0.0, 0.0, 0.0, tl),
        i = 0,
        r = 0.1)

    l1 = Cylinder(Point3f0(0.0, 0.0, 0.0), Point3f0(0.0, 0.0, model.l1),
        convert(Float32, 0.025))
    setobject!(vis["l1_$i"], l1, MeshPhongMaterial(color = color))
    l2 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l2),
        convert(Float32, 0.025))
    setobject!(vis["l2_$i"], l2, MeshPhongMaterial(color = color))

    setobject!(vis["elbow_nominal_$i"], Sphere(Point3f0(0.0),
        convert(Float32, 0.05)),
        MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, tl)))
	setobject!(vis["elbow_limit_$i"], Sphere(Point3f0(0.0),
        convert(Float32, 0.05)),
        MeshPhongMaterial(color = limit_color))
    setobject!(vis["ee_$i"], Sphere(Point3f0(0.0),
        convert(Float32, 0.05)),
        MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, tl)))
end

function _set_acrobot!(vis, model::DoublePendulum, x;
        i = 0, ϵ = 1.0e-1)

    p_mid = [kinematics_elbow(model, x)[1], 0.0, kinematics_elbow(model, x)[2]]
    p_ee = [kinematics(model, x)[1], 0.0, kinematics(model, x)[2]]

    settransform!(vis["l1_$i"], cable_transform(zeros(3), p_mid))
    settransform!(vis["l2_$i"], cable_transform(p_mid, p_ee))

    settransform!(vis["elbow_nominal_$i"], Translation(p_mid))
	settransform!(vis["elbow_limit_$i"], Translation(p_mid))
    settransform!(vis["ee_$i"], Translation(p_ee))

	if x[2] <= -0.5 * π + ϵ || x[2] >= 0.5 * π - ϵ
		setvisible!(vis["elbow_nominal_$i"], false)
		setvisible!(vis["elbow_limit_$i"], true)
	else
		setvisible!(vis["elbow_nominal_$i"], true)
		setvisible!(vis["elbow_limit_$i"], false)
	end
end

# visualization
function visualize_elbow!(vis, model::DoublePendulum, x;
        tl = 1.0,
        i = 0,
        color = RGBA(0.0, 0.0, 0.0, 1.0),
        limit_color = RGBA(0.0, 1.0, 0.0, tl),
        r = 0.1, Δt = 0.1,
        ϵ = 1.0e-1)

    _create_acrobot!(vis, model,
        tl = tl,
        color = color,
        limit_color = limit_color,
        i = i,
        r = r)

    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

    T = length(x)
    for t = 1:T
        MeshCat.atframe(anim,t) do
            _set_acrobot!(vis, model, x[t], i = i, ϵ = ϵ)
        end
    end

    # settransform!(vis["/Cameras/default"],
    #    compose(Translation(0.0 , 0.0 , 0.0), LinearMap(RotZ(pi / 2.0))))

    MeshCat.setanimation!(vis, anim)
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
if VERSION < v"0.7.0"
    cd(Pkg.dir("Coverage"))
end

using Coverage
cov_res = process_folder()
Coveralls.submit(cov_res)
Codecov.submit(cov_res)


text: abstract type AbstractDiscount end


text:
module kalman;
import StateSpace;
import Distributions;

export newModel,initialGuess,predict,update,extractMeanFromState,extractVarianceFromState;

sampleTimeDeltaTable = Dict{UTF8String,Float64}();

function initialGuess(observations::Array{Float64,1},varianceEstimate::Array{Float64,1})
   state = [observations; [0.0,0.0,0.0]]
   variance  = [varianceEstimate; [0.0001,0.0001,0.0001]]
   return Distributions.MvNormal(state,variance)
end

function fThunk(id)
   function f(t)
      dt = sampleTimeDeltaTable[id];
      F = eye(6)
      F[1,4] = dt
      F[2,5] = dt
      F[3,6] = dt
      return F;
   end
   return f
end

function qThunk(id,accVar)
   function q(t)
      dt = sampleTimeDeltaTable[id];
      G = [ dt^2/2; dt^2/2; dt^2/2; dt; dt; dt ]
      Q = eye(6).*(G*G'*accVar);
      return Q;
   end
   return q
end

function recordSampleTimeDelta(id,sampleTimeDelta)
   sampleTimeDeltaTable[id] = sampleTimeDelta;
end

function newModel(id,accVar::Float64,obVar::Float64,dt0)
   H = [ eye(3) zeros(3,3) ];
   R = eye(3)*obVar;
   recordSampleTimeDelta(id,dt0);
   return StateSpace.LinearGaussianSSM(fThunk(id),(t)->zeros(Float64,6,1),qThunk(id,accVar),(t)->H,(t)->R);
end

function update(id,model,predictedState,observations,dt)
   recordSampleTimeDelta(id,dt);
   return StateSpace.update(model,predictedState,copy(observations));
end

function predict(id,model,state)
   return StateSpace.predict(model,state)
end

function extractMeanFromState(state)
   return Distributions.mean(state)
end

function extractVarianceFromState(state)
   return Distributions.var(state)
end

end

max_stars_repo_name: waTeim/gpsFlow

max_issues_repo_name: waTeim/flying-monkey

max_forks_repo_name: waTeim/flying-monkey


text:
using LinearAlgebra

export loss

"""
    prior_KL(gp::SparseGaussianProcess{<:Any,<:Any,<:Any,<:Any,
                                       <:MarginalInducingPoints})

Computes the prior Kullback-Leibler divergence for a Gaussian process with
marginal inducing points, given by the expression

``KL(q(u) \\mathbin{||} p(u)) = \\frac{1}{2} \\left(
                                              \\ln\\frac{|K_{zz}|}{|\\Sigma|} +
                                              tr(K_{zz}^{-1}\\Sigma) +
                                              \\mu^T K_{zz} \\mu \\right)``

where the mean is re-parameterized according to
``\\mu = \\mathbb{E}( (K_{zz} + \\xi I)^{-1} u )``.
"""
function prior_KL(gp::SparseGaussianProcess{<:Any,<:Any,<:Any,<:Any,<:MarginalInducingPoints})
  (k,B) = (gp.kernel, gp.inter_domain_operator)
  (z,mu,U,V,QC) = gp.inducing_points()
  K = (B*k*B)(z,z)
  Q = K + V'*V
  UQ = QC.U
  logdet_term = 2 .* sum(log.(diag(UQ)) .- log.(diag(U)); dims=1)
  invUQ = inv(UQ)
  trace_term = sum((invUQ * invUQ') .* (U' * U); dims=(1,2))
  reparameterized_quadratic_form_term = sum(mu .* (Q * mu); dims=1)
  (logdet_term .- length(mu) .+ trace_term .+ reparameterized_quadratic_form_term) ./ 2
end

"""
    prior_KL(gp::SparseGaussianProcess{<:Any,<:Any,<:Any,<:Any,
                                       <:PseudoDataInducingPoints})

Computes the prior Kullback-Leibler divergence for a Gaussian process with
pseudo-data inducing points, given by the expression

``KL(q(u) \\mathbin{||} p(u)) = \\frac{1}{2} \\left(
                                              \\ln\\frac{|K_{zz} + \\Sigma|}{|\\Sigma|} +
                                              tr((K_{zz}+\\Sigma)^{-1}K_{zz}) +
                                              \\mu^T K_{zz} \\mu \\right)``

where the mean is re-parameterized according to
``\\mu = \\mathbb{E}( (K_{zz} + \\Sigma)^{-1} u )``.
"""
function prior_KL(gp::SparseGaussianProcess{<:Any,<:Any,<:Any,<:Any,<:PseudoDataInducingPoints})
  (k,B) = (gp.kernel, gp.inter_domain_operator)
  (z,mu,U,V,QC) = gp.inducing_points()
  K = (B*k*B)(z,z)
  Q = K + V'*V
  UQ = QC.U
  logdet_term = 2 .* sum(log.(diag(UQ)) .- log.(diag(V)); dims=1)
  invUQ = inv(UQ)
  trace_term = sum((invUQ * invUQ') .* K; dims=(1,2))
  reparameterized_quadratic_form_term = sum(mu .* (Q * mu); dims=1)
  (logdet_term .- length(mu) .+ trace_term .+ reparameterized_quadratic_form_term) ./ 2
end

"""
    loss(gp::GaussianProcess, x::AbstractMatrix, y::AbstractMatrix;
                                                 n_data::Int = size(x,ndims(x)))

Computes the Kullback-Leibler divergence of the variational family from the
posterior Gaussian process, up to an additive constant. Minimizing this
function trains the Gaussian process.
"""
function loss(gp::GaussianProcess, x::AbstractMatrix, y::AbstractMatrix; n_data::Int = size(x,ndims(x)))
  # sample the GP
  # rand!(gp.prior_basis)
  rand!(gp)

  # KL term
  kl = prior_KL(gp)

  # likelihood term
  n_samples = num_samples(gp)
  n_batch = size(x,ndims(x))
  f = gp(x)
  c = n_data ./ (n_batch .* n_samples .* 2)
  l = n_data .* gp.log_error .+ c .* sum(((y .- f) ./ exp.(gp.log_error)).^2; dims=(1,2))

  # regularizer term for hyperparameters
  r = gp.hyperprior.log_error(gp.log_error) .+ hyperprior_logpdf(gp.kernel)

  kl .+ l .+ r
end