text:
# Unit tests for Region type.
# This file is supposed to be included from runtests.jl.

@testset "Range" begin
    @test (0:-1).start == 0
    @test (0:-1).stop == -1
    @test length(0:-1) == 0
    @test isempty(0:-1)

    @test (0:0).start == 0
    @test (0:0).stop == 0
    @test length(0:0) == 1
    @test !isempty(0:0)

    @test (0:1).start == 0
    @test (0:1).stop == 1
    @test length(0:1) == 2
    @test !isempty(0:1)

    @test 0:1 == 0:1
    @test 0:1 ≠ 1:2

    @test 0:1 ≤ 0:1

    @test 0:1 < 1:2
    @test 0:1 ≤ 1:2
    @test 0:1 ≤ 0:1

    @test 0:1 ≥ 0:1

    @test 1:2 > 0:1
    @test 1:2 ≥ 0:1
    @test 0:1 ≥ 0:1

    @test translate(1:2, -5) == -4:-3
    @test translate(1:2, 5) == 6:7
    @test (1:2) - 5 == -4:-3
    @test (1:2) + 5 == 6:7

    @test invert(1:2) == -2:-1
    @test invert(-1:0) == 0:1
    @test invert(invert(-1:1)) == -1:1
    @test invert(invert(1:2)) == 1:2

    @test !contains(0:-1, 0)
    @test !contains(0:0, -1)
    @test contains(0:0, 0)
    @test !contains(0:0, 1)
    @test !contains(0:1, -1)
    @test contains(0:1, 0)
    @test contains(0:1, 1)
    @test !contains(0:1, 2)

    @test 0 ∉ 0:-1
    @test -1 ∉ 0:-1
    @test 0 ∈ 0:0
    @test 1 ∉ 0:0
    @test -1 ∉ 0:1
    @test 0 ∈ 0:1
    @test 1 ∈ 0:1
    @test 2 ∉ 0:1

    @test !isoverlapping(0:1, 6:7)
    @test !isoverlapping(1:2, 5:6)
    @test !isoverlapping(2:3, 4:5)
    @test isoverlapping(2:3, 3:4)
    @test isoverlapping(3:4, 3:4)
    @test !isoverlapping(4:5, 2:3)
    @test !isoverlapping(5:6, 1:2)
    @test !isoverlapping(6:7, 0:1)

    @test !istouching(0:1, 6:7)
    @test !istouching(1:2, 5:6)
    @test istouching(2:3, 4:5)
    @test istouching(2:3, 3:4)
    @test istouching(3:4, 3:4)
    @test istouching(4:5, 2:3)
    @test !istouching(5:6, 1:2)
    @test !istouching(6:7, 0:1)

    @test !isclose(0:1, 6:7, 0)
    @test !isclose(1:2, 5:6, 0)
    @test !isclose(2:3, 4:5, 0)
    @test isclose(2:3, 3:4, 0)
    @test isclose(3:4, 3:4, 0)
    @test !isclose(4:5, 2:3, 0)
    @test !isclose(5:6, 1:2, 0)
    @test !isclose(6:7, 0:1, 0)

    @test !isclose(0:1, 6:7, 1)
    @test !isclose(1:2, 5:6, 1)
    @test isclose(2:3, 4:5, 1)
    @test isclose(2:3, 3:4, 1)
    @test isclose(3:4, 3:4, 1)
    @test isclose(4:5, 2:3, 1)
    @test !isclose(5:6, 1:2, 1)
    @test !isclose(6:7, 0:1, 1)

    @test !isclose(0:1, 6:7, 3)
    @test isclose(1:2, 5:6, 3)
    @test isclose(2:3, 4:5, 3)
    @test isclose(2:3, 3:4, 3)
    @test isclose(3:4, 3:4, 3)
    @test isclose(4:5, 2:3, 3)
    @test isclose(5:6, 1:2, 3)
    @test !isclose(6:7, 0:1, 3)
end # "Range"


text:
#
# Scenario
#

@generated function Base.setindex!(s::Scenario, val::Number, sc::CashRiskFactor)
    if val <: Currencies.Cash
        return quote
            error("setindex!($s, $val, $sc) not implemented.")
        end
    else
        return quote
            error("The value for a `Core.CashRiskFactor` must be a `Currencies.Cash`.")
        end
    end
end

function Base.haskey(scenario::Scenario, c::Currencies.Currency) :: Bool
    return haskey(scenario, SpotCurrency(c))
end

function Base.getindex(scenario::Scenario, c::Currencies.Currency) :: Currencies.Cash
    return scenario[SpotCurrency(c)]
end

function Base.setindex!(scenario::Scenario, val::Currencies.Cash, c::Currencies.Currency)
    Base.setindex!(scenario, val, SpotCurrency(c))
    nothing
end

for T in (:FixedCashRiskFactor, :FixedNonCashRiskFactor), S in (:FixedScenario, :ActualScenario)
    @eval begin
        Base.haskey(scenario::$S, rf::$T) = true
        Base.getindex(scenario::$S, rf::$T) = rf.val
    end
end

#
# FixedScenario
#

function Base.setindex!(s::FixedScenario, val::Currencies.Cash, rf::CashRiskFactor)
    s.cash_risk_factors[rf] = val
    nothing
end

function Base.haskey(s::FixedScenario, rf::CashRiskFactor) :: Bool
    return haskey(s.cash_risk_factors, rf)
end

function Base.getindex(s::FixedScenario, rf::CashRiskFactor) :: Currencies.Cash
    return s.cash_risk_factors[rf]
end

function Base.setindex!(s::FixedScenario, val::Float64, rf::NonCashRiskFactor)
    s.non_cash_risk_factors[rf] = val
    nothing
end

# risk factor values are currently stored as Float64
Base.setindex!(s::FixedScenario, val::Number, rf::NonCashRiskFactor) = s[rf] = Float64(val)

function Base.haskey(s::FixedScenario, rf::NonCashRiskFactor) :: Bool
    return haskey(s.non_cash_risk_factors, rf)
end

function Base.getindex(s::FixedScenario, rf::NonCashRiskFactor) :: Float64
    return s.non_cash_risk_factors[rf]
end

Base.length(scenario::FixedScenario) = length(scenario.cash_risk_factors) + length(scenario.non_cash_risk_factors)
Base.isempty(scenario::FixedScenario) = length(scenario) == 0

#
# ActualScenario
#

@inline function get_locf_option(scenario::ActualScenario, sym::Symbol) :: Bool
    return scenario.locf_option(scenario, sym)
end

@inline function Base.haskey(act::ActualScenario, rf::RiskFactor)
    result = MarketData.has_serie(act.provider, market_data_symbol(rf))
    if !result
        @warn("Provider doesn't have time series for $(market_data_symbol(rf))")
    end
    return result
end

@inline function Base.getindex(scenario::ActualScenario, rf::SpotCurrency{C}) :: Currencies.Cash where {C<:Currencies.Currency}
    sym = market_data_symbol(rf)
    result = MarketData.get_cash(scenario.provider, sym, scenario.date, scenario.date; locf=get_locf_option(scenario, sym))
    if ismissing(result)
        error("Error getting value for risk factor $rf: value is missing for date $(scenario.date)")
    else
        return result
    end
end

@inline function  Base.getindex(scenario::ActualScenario, rf::Stock) :: Currencies.Cash
    sym = market_data_symbol(rf)
    result = MarketData.get_cash(scenario.provider, sym, scenario.date, scenario.date; locf=get_locf_option(scenario, sym))
    if ismissing(result)
        error("Error getting value for risk factor $rf: value is missing for date $(scenario.date)")
    else
        return result
    end
end

@inline function Base.getindex(scenario::ActualScenario, rf::DiscountFactor) :: Float64
    @assert scenario.date <= rf.maturity "Why get value a DiscountFactor value for a past date? Scenario date $(scenario.date); DiscountFactor maturity: $(rf.maturity)."
    sym = market_data_symbol(rf)
    curve = MarketData.get_value(scenario.provider, sym, scenario.date, scenario.date; locf=get_locf_option(scenario, sym))

    if ismissing(curve)
        error("Error getting value for risk factor $rf: curve $sym is missing for date $(scenario.date).")
    else
        return InterestRates.discountfactor(curve, rf.maturity)
    end
end

#
# CompositeScenario
#

@inline function Base.haskey(composite::CompositeScenario, rf::RiskFactor)
    for scenario in composite.scenarios
        if haskey(scenario, rf)
            return true
        end
    end

    return false
end

@inline function Base.getindex(composite::CompositeScenario, rf::RiskFactor) :: Float64
    for scenario in composite.scenarios
        if haskey(scenario, rf)
            return scenario[rf]
        end
    end
    error("Scenario has no value for risk factor $rf.")
end

@inline function Base.getindex(composite::CompositeScenario, sc::SpotCurrency) :: Currencies.Cash
    for scenario in composite.scenarios
        if haskey(scenario, sc)
            return scenario[sc]
        end
    end
    error("Scenario has no value for risk factor $rf.")
end

function Base.length(composite::CompositeScenario)
    l = 0

    for scenario in composite.scenarios
        l += length(scenario)
    end

    return l
end

function Base.isempty(composite::CompositeScenario)
    for scenario in composite.scenarios
        if !isempty(scenario)
            return false
        end
    end

    return true
end

#
# DebugScenario
#

Base.haskey(s::DebugScenario, key::RiskFactor) = haskey(s.scenario, key)

function Base.getindex(s::DebugScenario, key::RiskFactor)
    # tracks which risk factors are used
    item = getindex(s.scenario, key)
    s.record[key] = item
    return item
end

#
# ScenarioMap
#

Base.haskey(s::ScenarioMap, rf::RiskFactor) = haskey(s.scenario, rf)
Base.getindex(s::ScenarioMap, rf::RiskFactor) = s.f(rf, s.scenario[rf])

ext: jl


text:
module EntityManager

export Entity

using ..Management
using CoreExtended
using MatrixMath
using JLGEngine

""" TODO """
type Entity <: JLGEngine.IComponent
	id::Symbol
	enabled::Bool

  position::Vec3f
  rotation::Vec3f
	scaling::Vec3f

	Entity(id=:_) = new(id,true,zeros(Vec3f),zeros(Vec3f),ones(Vec3f))
end

# preset for manager: includes code here
#presetManager(Entity)
m=Management.presetManager(Entity)

end #EntityManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
using Random, HCubature, QuadGK, Distributions, Statistics, LinearAlgebra

function test_t_integral_1D_gaussian()

    rg = MersenneTwister(1)

    # easy case: standard univeriate Gaussian

    σ = 3.1

    μ = 0.1

    g = Normal(μ, σ)

    x = randn(rg)

    f(t) = pdf(g, t*x)

    # numerical integral

    num1 = hquadrature(f, -50, 50)[1]

    num2 = quadgk(f, -50, 50)[1]

    # analytical solution

    ana = 1/sqrt(2π*σ*σ) * sqrt(π) / sqrt(x*x/(2*σ*σ))

    num1, num2, ana, abs(num1-ana)

end

function test_t_integral_2D_diagonal_gaussian(seed)

    rg = MersenneTwister(seed)

    D = 3

    Σ = Diagonal(rand(rg, D)*2)

    U = svd(Σ).U

    μ = randn(rg, D)

    x = randn(rg, D)

    num1 = numericalintegraloft(x, μ, Σ)

    ana = analyticalintegraloft(x, μ, Σ)

    num1, ana, abs(num1-ana)

end

function test_t_integral_2D_full_gaussian(seed=1)

    rg = MersenneTwister(seed)

    D = 3

    Σ = randn(rg,D, D)*0.5; Σ = Σ*Σ' + 0.5*I

    U = svd(Σ).U

    μ = randn(rg, D)

    x = randn(rg, D)

    num1 = numericalintegraloft(x, μ, Σ)

    ana = analyticalintegraloft(x, μ, Σ)

    num1, ana, abs(num1-ana)

end

function numericalintegraloft(x, μ, Σ)

    g = MvNormal(μ, Σ)

    hquadrature(t -> pdf(g, t*x), -50, 50)[1]

end

function analyticalintegraloft(x, μ, Σ)

    D = length(μ) ; @assert(size(Σ, 1) == size(Σ, 2) == D == length(x))

    # see https://en.wikipedia.org/wiki/Gaussian_integral#Generalizations
    a =  0.5*x'*(Σ\x)
    b =  x'*(Σ\μ)
    c = -0.5*μ'*(Σ\μ)

    1/sqrt((2π)^D*det(Σ)) * sqrt(π/a) * exp(b^2/(4*a)+c)

end


text:
# some type alias
const RGBAf0 = RGBA{Float32}

"""
A simple iterator that returns a new, unique color when `next(::UniqueColorIter)` is called.
"""
mutable struct UniqueColorIter{T}
    colors::T
    state::Int
end

UniqueColorIter(x::Union{Symbol, String}) = UniqueColorIter(to_colormap(x), 1)

function Base.getindex(iter::UniqueColorIter, idx::Int)
    # TODO make out of bounds more graceful? But hooow
    iter.colors[mod1(idx, length(iter.colors))]
end

Base.start(iter::UniqueColorIter) = (iter.state = 1; (iter, iter))

function Base.next(iter::UniqueColorIter)
    result = iter[iter.state]
    iter.state += 1
    result
end
# make the base interface happy
Base.next(iter::UniqueColorIter, iter2) = (next(iter), iter2)


text:
function table{T}(d::AbstractDataArray{T})
    counts = Dict{Union(T, NAtype), Int}()
    for i = 1:length(d)
        if haskey(counts, d[i])
            counts[d[i]] += 1
        else
            counts[d[i]] = 1
        end
    end
    return counts
end

const letters = convert(Vector{ASCIIString}, split("abcdefghijklmnopqrstuvwxyz", ""))
const LETTERS = convert(Vector{ASCIIString}, split("ABCDEFGHIJKLMNOPQRSTUVWXYZ", ""))

# Like string(s), but preserves Vector{String} and converts
# Vector{Any} to Vector{String}.
_vstring{T <: String}(s::T) = s
_vstring{T <: String}(s::AbstractVector{T}) = s
_vstring(s::AbstractVector) = String[_vstring(x) for x in s]
_vstring(s::Any) = string(s)
vcatstring(x) = vcat(_vstring(x))

function paste(s...)
    s = map(vcatstring, {s...})
    sa = {s...}
    N = maximum(length, sa)
    res = fill("", N)
    for i in 1:length(sa)
        Ni = length(sa[i])
        k = 1
        for j = 1:N
            res[j] = string(res[j], sa[i][k])
            if k == Ni   # This recycles array elements.
                k = 1
            else
                k += 1
            end
        end
    end
    res
end

function paste_columns(d::AbstractDataFrame, sep)
    res = fill("", nrow(d))
    for j in 1:ncol(d)
        for i in 1:nrow(d)
            res[i] *= string(d[i,j])
            if j != ncol(d)
                res[i] *= sep
            end
        end
    end
    res
end
paste_columns(d::AbstractDataFrame) = paste_columns(d, "_")

function cut{S, T}(x::Vector{S}, breaks::Vector{T})
    if !issorted(breaks)
        sort!(breaks)
    end
    min_x = minimum(x)
    max_x = maximum(x)
    if breaks[1] > min_x
        unshift!(breaks, min_x)
    end
    if breaks[end] < max_x
        push!(breaks, max_x)
    end
    refs = fill(zero(DEFAULT_POOLED_REF_TYPE), length(x))
    for i in 1:length(x)
        if x[i] == min_x
            refs[i] = 1
        else
            refs[i] = searchsortedfirst(breaks, x[i]) - 1
        end
    end
    n = length(breaks)
    from = map(x -> sprint(showcompact, x), breaks[1:(n - 1)])
    to = map(x -> sprint(showcompact, x), breaks[2:n])
    pool = Array(ASCIIString, n - 1)
    if breaks[1] == min_x
        pool[1] = string("[", from[1], ",", to[1], "]")
    else
        pool[1] = string("(", from[1], ",", to[1], "]")
    end
    for i in 2:(n - 1)
        pool[i] = string("(", from[i], ",", to[i], "]")
    end
    PooledDataArray(RefArray(refs), pool)
end
cut(x::Vector, ngroups::Int) = cut(x, quantile(x, [1 : ngroups - 1] / ngroups))

##############################################################################
##
## rep()
##
##############################################################################

function rep(x::AbstractVector, lengths::Vector{Int})
    if length(x) != length(lengths)
        error("vector lengths must match")
    end
    res = similar(x, sum(lengths))
    i = 1
    for idx in 1:length(x)
        tmp = x[idx]
        for kdx in 1:lengths[idx]
            res[i] = tmp
            i += 1
        end
    end
    res
end
rep(x::AbstractVector, times::AbstractVector{Int}, each::Int) = rep(x, times)

function rep(x::AbstractVector, times::Int = 1, each::Int = 1)
    res = similar(x, each * sum(times) * length(x))
    i = 1
    for jdx in 1:times
        for idx in 1:length(x)
            for kdx in 1:each
                res[i] = x[idx]
                i += 1
            end
        end
    end
    res
end

rep(x, times) = fill(x, times)

function rep(x; times = 1, each::Int = 1)
    rep(x, times, each)
end


text:
MOI.get(model::Optimizer, ::MOI.NumberOfVariables) = length(model.inner.variable_info)
MOI.get(model::Optimizer, ::MOI.ListOfVariableIndices) = VI.(1 : length(model.inner.variable_info))

function MOI.add_variable(model::Optimizer)
    push!(model.inner.variable_info, VariableInfo())
    return VI(length(model.inner.variable_info))
end

function MOI.add_variables(model::Optimizer, nvars::Integer)
    return [MOI.add_variable(model) for i in 1:nvars]
end

function MOI.add_constraint(model::Optimizer, vi::MOI.VariableIndex, lt::MOI.LessThan{Float64})
    check_variable_indices(model, vi)
    set_upper_bound(model.inner.variable_info[vi.value], lt.upper)
    return MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}}(vi.value)
end

function MOI.add_constraint(model::Optimizer, vi::MOI.VariableIndex, gt::MOI.GreaterThan{Float64})
    check_variable_indices(model, vi)
    set_lower_bound(model.inner.variable_info[vi.value], gt.lower)
    return MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{Float64}}(vi.value)
end

function MOI.add_constraint(model::Optimizer, vi::MOI.VariableIndex, eq::MOI.EqualTo{Float64})
    check_variable_indices(model, vi)
    set_lower_bound(model.inner.variable_info[vi.value], eq.value)
    set_upper_bound(model.inner.variable_info[vi.value], eq.value)
    return MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}}(vi.value)
end

# see comment in: write_bar_file
# MOI.supports(::Optimizer, ::MOI.VariableName, ::Type{VI}) = true
# function MOI.set(model::Optimizer, attr::MOI.VariableName, vi::VI, value)
#     check_variable_indices(model, vi)
#     model.inner.variable_info[vi.value].name = value
# end
# function MOI.get(model::Optimizer, ::MOI.VariableName, vi::VI)
#      return model.inner.variable_info[vi.value].name
# end
# function MOI.get(model::Optimizer, ::Type{MathOptInterface.VariableIndex}, name::String)
#     for (i,var) in enumerate(model.inner.variable_info)
#         if name == var.name
#             return VI(i)
#         end
#     end
#     error("Unrecognized variable name $name.")
# end

MOI.supports(::Optimizer, ::MOI.VariablePrimalStart, ::Type{VI}) = true

function MOI.set(model::Optimizer, ::MOI.VariablePrimalStart, vi::VI, value::Union{Real, Nothing})
    check_variable_indices(model, vi)
    model.inner.variable_info[vi.value].start = value
    return
end


text:
# %% imports

include("../../../Modules/save_param.jl")
using Statistics
using PyCall
gym = pyimport("gym")
include("../../../Modules/Networks/net_arch.jl")
include("../../../Modules/Poisson/PoissonStateSpike.jl")
include("../../../Modules/simviz.jl")

# %% Read the results

results = load_param("Julia/CartPole/LIF_base/results/results_noalpharandn.jld")

# %% sim functions

function count_spikes(arch::Array{Int64}, spiked::Array{Int64}, layer::Int64)
    count_layer = zeros(arch[layer])
    interval = [sum(arch[1:layer-1])+1, sum(arch[1:layer])]
    spikes_in = filter(x-> interval[1] <= x <= interval[2], spiked)
    spikes_in .-= sum(arch[1:layer-1])
    count_layer[spikes_in] .+= 1
    return count_layer
end

function left_right(count_spikes::Array{Float64}, n_left::Int64, n_right::Int64)
    sl = sum(count_spikes[1:n_left])
    sr = sum(count_spikes[n_left+1:end])
    (sl >= sr) ? action = 0 : action = 1
    return action, sl, sr
end

function borne(x,a,b)
    y = a + 0.5*(b-a)*(1+cos(pi*x/10))
    return y
end

function sig(z::Float64, alpha::Float64)
    return 1.0/(1.0+exp(-z/alpha))
end

function norm_obs(obs::Array{Float64})
    pos = obs[1]
    theta = obs[3]%pi
    vel = [obs[2], obs[4]]
    pos_scaled = (pos + 2.4)/4.8
    cone = 41.8*pi/180
    theta_scaled = (theta + cone)/(2*cone)
    vel_scaled = (vel .+ 7.5)/15
    obs_scaled = [pos_scaled, vel_scaled[1], theta_scaled, vel_scaled[2]]
    return obs_scaled
end

function norm_sig(obs::Array{Float64})
    pos = obs[1] # in [-2.4,2.4]
    theta = obs[3]%pi # in [-pi, pi]
    vel = [obs[2], obs[4]] # in [-15,15]
    alpha_pos = 2.4/5
    alpha_theta = pi/5
    alpha_vel = 15/5
    pos_scaled = sig(pos,alpha_pos)
    theta_scaled = sig(theta,alpha_theta)
    vel_scaled = sig.(vel,alpha_vel)
    obs_scaled = [pos_scaled, vel_scaled[1], theta_scaled, vel_scaled[2]]
    return obs_scaled
end

function obs_input(obs::Array{Float64}, gain::Float64)
    input = Float64[]
    obs_s = norm_obs(obs)
    for idx in eachindex(obs_s)
        anti_obs = gain*(1.0-obs_s[idx])
        push!(input, obs_s[idx])
        push!(input, anti_obs)
    end
    return input
end

# %% main functions

function play_episode(net::Network, n_steps::Int64 = 500)
    env = gym.make("CartPole-v1")
    env.seed(0)
    #obs = swingup_reset!(env)
    obs = env.reset()
    rng_action = Random.MersenneTwister(0)
    win_size = 100 # ms
    tot_reward = 0
    N_rand = 0
    N_tot = 0
    matalpha =  randn((4,4))
    @inbounds for j in 1:n_steps
        count_tot = zeros(net.arch[end])
        obs_alpha = matalpha*obs
        input = obs_input(obs_alpha, 1.0)
        spikes = int_spikes(input) # state2spike function
        #reset voltage
        net.v .= net.param["c"]
        influence = 2
        @inbounds for msec in 1:win_size
            spiked = step_LIF!(net, spikes, false)
            if msec >= length(net.arch)-1
                count_tot .+= influence*count_spikes(net.arch, spiked, length(net.arch))
                influence *= 0.995
                #push!(spike_hist,spiked)
            end
        end
        n_left = div(net.arch[end],2)
        n_right = net.arch[end] - n_left
        action, sl, sr = left_right(count_tot, n_left, n_right)
        if sl == sr
            action = rand(rng_action, 0:1)
            tot_reward -= 0.9
            N_rand += 1
        end
        N_tot += 1
        #push!(score_hist[1], sl)
        #push!(score_hist[2], sr)
        #obs_new, reward, done, _ = swingup_step(env,action)
        obs_new, reward, done, _ = env.step(action)
        tot_reward += reward
        if done
            break
        end
        obs = obs_new
    end
    println("Rand = ", N_rand, "/", N_tot)
    env.close()
    env = nothing
    Base.GC.gc()
    return tot_reward
end

function set_weight!(net::Network, genes_w::Array{Float64})
    w = net.s.e
    for mat_idx in 1:length(w)
        for idx in eachindex(w[mat_idx])
            offset = sum([length(mat) for mat in w[1:mat_idx-1]])
            net.s.e[mat_idx][idx] = borne(genes_w[idx+offset], 0.0, 1.1)
        end
    end
end

function sim(genes, render::Bool = false, n_steps::Int64 = 500)
    #load the optimisation parameters
    optim_param = YAML.load(open("Julia/CartPole/LIF_base/cfg_cmaes.yml"))
    arch = optim_param["arch"]
    genes_w = genes

    params = YAML.load(open("Julia/CartPole/LIF_base/noalpha/cfglifnoalpha.yml"))
    n = Network(arch, params)
    set_weight!(n,genes_w)
    reward = play_episode(n, n_steps)
    return reward, n.arch
end

# %% test robustness

function test_random(results::Dict)
    mean_reward = Float64[]
    for i in 1:length(results["fit"])
        genes = results["weights"][i]
        tmp_reward = []
        for ep in 1:30
            reward, _, _ = sim(genes, false, 200)
            push!(tmp_reward, reward)
        end
        push!(mean_reward, mean(tmp_reward))
    end
    return mean_reward
end

function test_random(results::Dict, indiv::Int64)
    genes = results["weights"][indiv]
    tmp_reward = []
    for ep in 1:30
        reward, _ = sim(genes, false, 200)
        push!(tmp_reward, reward)
    end
    mean_r = mean(tmp_reward)
    return mean_r
end

means = test_random(results,20)

# %% visualize the activity in the network

function activiz(results::Dict, indiv::Int64)
    genes = results["weights"][indiv]
    r, spike_hist, arch = sim(genes)
    spike_histogram(arch, spike_hist)
end

for idx in [20]
    activiz(results, idx)
end

# %% run one episode

genes = results["weights"][20]
r, spike_hist, arch = sim(genes, true, 500)

# %%
using Plots

for indiv in 20
    display(histogram(map(x->borne(x,0.0, 1.1), results["weights"][indiv]), xlabel = "indiv $(indiv)",  bins = collect(0.0:0.1:1.1)))
end

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
#------------------------------------------------------------------------------
# Author: Edward Tate <edward.tate@erlang-solutions.com>
#------------------------------------------------------------------------------
module FsFootball
include("simulation.jl")
end

max_stars_repo_name: esl/FsFootball


text:
struct Projection{S, I}
    set::S
    indices::I
end
Polyhedra.project(set, indices) = Projection(set, indices)
function space_variables(p::Projection)
    vars = space_variables(p.set)
    if vars === nothing
        return nothing
    else
        return vars[p.indices]
    end
end
dimension(p::Projection) = length(p.indices)

# A^{-1} * S
struct LinearImage{S, T, MT <: AbstractMatrix{T}} <: AbstractSet{T}
    set::S
    A::MT
end
perspective_variable(li::LinearImage) = perspective_variable(li.set)
space_variables(li::LinearImage) = nothing
dimension(li::LinearImage) = size(li.A, 1)

# A^{-1} * S
struct LinearPreImage{S, T, MT <: AbstractMatrix{T}} <: AbstractSet{T}
    set::S
    A::MT
end

# S + c
struct Translation{S, T, VT <: AbstractVector{T}} <: AbstractSet{T}
    set::S
    c::VT
end
dimension(t::Translation) = length(t.c)
space_variables(t::Translation) = space_variables(t.set)
function Polyhedra.project(t::Translation, I)
    return Translation(Polyhedra.project(t.set, I), t.c[I])
end

"""
    householder(x)

Householder reflection
```math
I - 2 v v^T / (v^T v)
```
It is symmetric and orthogonal.
"""
function householder(x)
    y = copy(x)
    t = LinearAlgebra.reflector!(y)
    v = [1; y[2:end]]
    I - t * v * v'
end
_householder(h) = householder([1.0; h]) # We add 1, for perspective variable z

struct Householder{T, S <: AbstractSet{T}, U} <: AbstractSet{T}
    set::S
    p::DynamicPolynomials.Polynomial{true, U}
    h::Vector{Float64}
    z::SpaceVariable
    x::Vector{SpaceVariable}
end
perspective_gauge0(set::Householder) = set.p
perspective_variable(set::Householder) = set.z
space_variables(set::Householder) = set.x
convexity_proof(set::Householder) = convexity_proof(set.set)

const HouseDualOf{S, T, U} = PerspectiveDualOf{Householder{T, S, U}}

function Polyhedra.project(set::Polar{T}, I) where T
    return Polyhedra.polar(zero_eliminate(Polyhedra.polar(set), setdiff(1:dimension(set), I)))
end

function zero_eliminate(set::Householder{T}, I) where T
    J = setdiff(1:dimension(set), I)
    p = subs(set.p, set.x[I] => zeros(T, length(I)))
    return Householder(UnknownSet{T}(), p, set.h[J], set.z, set.x[J])
end
function Polyhedra.project(set::PerspectiveDual, I)
    return perspective_dual(zero_eliminate(perspective_dual(set),
                                           setdiff(1:dimension(set), I)))
end

function _HPH(set::Householder)
    H = _householder(set.h)
    return H * _HPH(set.set) * H
end


text:
facts("Error-free add") do
  x, y = err_add(1.0, 2.0)
  @fact x => 1.0 + 2.0
  @fact y => 0.0

  x, y = err_add(1.0, 2e-50)
  @fact x => 1.0
  @fact y => 2e-50

  x, y = err_add(2e-50, 1.0)
  @fact x => 1.0
  @fact y => 2e-50
end

facts("Fast error-free add") do
  x, y = err_fast_add(2.0, 1.0)
  @fact x => 2.0 + 1.0
  @fact y => 0.0

  x, y = err_fast_add(1.0, 2e-50)
  @fact x => 1.0
  @fact y => 2e-50
end


text:
# TODO: Add support for uploading a CuSparseMatrixBSR

Base.@kwdef mutable struct AMGXMatrix <: AMGXObject
    handle::API.AMGX_matrix_handle = API.AMGX_matrix_handle(C_NULL)
    mode::Union{Mode, Nothing} = nothing
    resources::Union{Resources, Nothing} = nothing
    function AMGXMatrix(handle::API.AMGX_matrix_handle, mode::Union{Mode, Nothing},
                       resources::Union{Resources, Nothing})
        m = new(handle, mode, resources)
        finalizer(warn_not_destroyed_on_finalize, m)
        return m
    end
end
function Base.show(io::IO, mime::MIME"text/plain", m::AMGXMatrix)
    invoke(show, Tuple{IO, MIME"text/plain", AMGXObject}, io, mime, m)
    m.handle == C_NULL && return
    n, block_dims = matrix_get_size(m)
    if n !== 0
        nz = amgx_nnz(m)
        if block_dims == (1, 1)
            print(io, " of size $n×$n with $nz stored entries")
        else
            block_dim_x, block_dim_y = block_dims
            print(io, " of size $n⋅$block_dim_x×$n⋅$block_dim_y with $nz stored block entries")
        end
    end
end
get_api_destroy_call(::Type{AMGXMatrix}) = API.AMGX_matrix_destroy
function dec_refcount_parents(m::AMGXMatrix)
    dec_refcount!(m.resources)
    m.resources = nothing
    nothing
end

function create!(m::AMGXMatrix, res::Resources, mode::Mode)
    mat_handle_ptr = Ref{API.AMGX_matrix_handle}()
    amgx_mode = API.AMGX_Mode(Int(mode))
    @checked API.AMGX_matrix_create(mat_handle_ptr, res.handle, amgx_mode)
    m.handle = mat_handle_ptr[]
    m.mode = mode
    m.resources = res
    inc_refcount!(m.resources)
    return m
end
AMGXMatrix(res::Resources, mode::Mode) = create!(AMGXMatrix(), res, mode)

function upload!(m::AMGXMatrix, n::Int, nnz::Int, row_ptrs_ptr::PtrOrCuPtrUnion{Cint}, col_indices_ptr::PtrOrCuPtrUnion{Cint},
                 data_ptr::PtrOrCuPtrUnion{T}; block_dims::Tuple{Int, Int}=(1,1),
                 diag_data_ptr::PtrOrCuPtrUnion{T}=Ptr{T}(C_NULL)) where {T <: Union{Float64, Float32}}
    mT = matrix_type(m.mode)
    if mT != T
        throw(ArgumentError("inconsistent AMGX matrix mode ($mT) with element type of upload ($T)"))
    end
    # TODO the CuPtr conversions should be handled in the ccall
    diag_data_ptr = reinterpret(Ptr{T}, diag_data_ptr)
    col_indices_ptr = reinterpret(Ptr{Cint}, col_indices_ptr)
    data_ptr = reinterpret(Ptr{T}, data_ptr)
    row_ptrs_ptr = reinterpret(Ptr{Cint}, row_ptrs_ptr)
    block_dim_x, block_dim_y = block_dims
    @checked API. AMGX_matrix_upload_all(m.handle, n, nnz, block_dim_x, block_dim_y, row_ptrs_ptr,
                                            col_indices_ptr, data_ptr, diag_data_ptr)
    return m
end

# Note: Data within the block is assumed to be arranged in row-major scanline order.
function upload!(m::AMGXMatrix, row_ptrs::VectorOrCuVector{Cint}, col_indices::VectorOrCuVector{Cint},
                 data::VectorOrCuVector{T}; block_dims::Tuple{Int, Int}=(1,1),
                 diag_data::Union{VectorOrCuVector{T}, Nothing}=nothing) where {T <: Union{Float64, Float32}}
    nnz, rem = divrem(length(data), prod(block_dims))
    if rem != 0
        throw(ArgumentError("length of `data` ($length(data)) is not an integer multiple of block dimension size $block_dims"))
    end
    n = length(row_ptrs) - 1
    if nnz != length(col_indices)
        throw(ArgumentError("length of `data` ($(length(data))) is inconsistent with length of `row_ptrs` ($(length(row_ptrs))"))
    end
    # Disable this check if we have a CUArray since it does a scalar index (perhaps harmless here though)
    if row_ptrs isa Vector
        if row_ptrs[end] != nnz
            throw(ArgumentError("length of data ($(length(data)) is inconsitent with last element of `row_ptrs` ($(row_ptrs[end]))"))
        end
    end
    if diag_data !== nothing
        if length(diag_dat) != n * prod(block_dims)
            throw(ArgumentError("length of `diag_data` ($(length(diag_data))) is not equal to number of elements on diagonal"))
        end
    end
    GC.@preserve row_ptrs col_indices data diag_data begin
        diag_data_ptr = diag_data === nothing ? Ptr{T}(C_NULL) : pointer(diag_data)
        upload!(m, n, nnz, pointer(row_ptrs), pointer(col_indices), pointer(data); block_dims, diag_data_ptr)
    end
    return m
end

function upload!(matrix::AMGXMatrix, cu_matrix::CUDA.CUSPARSE.CuSparseMatrixCSR)
    upload!(matrix, cu_matrix.rowPtr, cu_matrix.colVal, cu_matrix.nzVal)
end

function matrix_get_size(matrix::AMGXMatrix)
    n_ptr, block_dim_x_ptr, block_dim_y_ptr = Ref{Cint}(), Ref{Cint}(), Ref{Cint}()
    @checked API.AMGX_matrix_get_size(matrix.handle, n_ptr, block_dim_x_ptr, block_dim_y_ptr)
    return Int(n_ptr[]), (Int(block_dim_x_ptr[]), Int(block_dim_y_ptr[]))
end
function Base.size(matrix::AMGXMatrix)
    n, block_dims = matrix_get_size(matrix)
    return n * block_dims[1], n * block_dims[2]
end

function replace_coefficients!(m::AMGXMatrix, data::VectorOrCuVector{T}, diag_data::Union{VectorOrCuVector{T}, Nothing}=nothing) where {T <: Union{Float64, Float32}}
    n, block_dims = matrix_get_size(m)
    _amgx_nnz = amgx_nnz(m)
    if length(data) != _amgx_nnz * prod(block_dims)
        throw(ArgumentError("can not change the number of nnz entries in `replace_coefficients!`"))
    end
    mT = matrix_type(m.mode)
    if mT != T
        throw(ArgumentError("inconsistent AMGX matrix mode ($mT) with element type of upload ($T)"))
    end
        GC.@preserve data diag_data begin
            diag_data_ptr = diag_data === nothing ? Ptr{T}(C_NULL) : pointer(diag_data)
            replace_coefficients!(m, n, _amgx_nnz, pointer(data), diag_data_ptr)
        end
    return m
end

function replace_coefficients!(m::AMGXMatrix, n::Int, nnz::Int, data_ptr::PtrOrCuPtrUnion{T},
                               diag_data_ptr::PtrOrCuPtrUnion{T}=Ptr{T}(C_NULL)) where {T <: Union{Float64, Float32}}
    # TODO the CuPtr conversions should be handled in the ccall
    data_ptr = reinterpret(Ptr{T}, data_ptr)
    diag_data_ptr = reinterpret(Ptr{T}, diag_data_ptr)
    @checked API.AMGX_matrix_replace_coefficients(m.handle, n, nnz, data_ptr, diag_data_ptr)
    return m
end

function amgx_nnz(matrix)
    nnz_ptr = Ref{Cint}()
    @checked API.AMGX_matrix_get_nnz(matrix.handle, nnz_ptr)
    return Int(nnz_ptr[])
end

function SparseArrays.nnz(matrix::AMGXMatrix)
    _, block_dims = matrix_get_size(matrix)
    return amgx_nnz(matrix) * prod(block_dims)
end


text:
using Zygote
using HMCUtilities:
    VariableConstraint,
    free_dimension,
    constrain_dimension,
    free,
    constrain,
    constrain_with_pushlogpdf,
    constrain_with_pushlogpdf_grad,
    _logabsdet

_size(x) = size(x)
_size(x::Number) = (length(x),)

function logabsdetjac(c, y)
    J = HMCUtilities.constrain_jacobian(c, y)
    return _logabsdet(J'J)[1] / 2
end

_format_grad(z, x) = Zygote.accum(z, x)

function test_constraint(
        c,
        x,
        y,
        logπx,
        ∇x_logπx,
        y_exp,
        logπy_exp,
        ∇y_logπy_exp;
        cvtypes=[Vector{Float64}, Vector{Float32}],
        fvtypes=[Vector{Float64}, Vector{Float32}],
        test_type_stability = true,
        test_logdetJ_consistency = true
    )
    @test isa(c, VariableConstraint)

    @testset "$CVType, $FVType" for (CVType, FVType) in zip(cvtypes, fvtypes)
        _test_constraint(
            c,
            convert(CVType, x),
            convert(FVType, y),
            convert(eltype(CVType), logπx),
            convert(CVType, ∇x_logπx),
            convert(FVType, y_exp),
            convert(eltype(FVType), logπy_exp),
            convert(FVType, ∇y_logπy_exp);
            test_type_stability = test_type_stability,
            test_logdetJ_consistency = test_logdetJ_consistency
        )
    end
end

function _test_constraint(
    c,
    x,
    y,
    logπx,
    ∇x_logπx,
    y_exp,
    logπy_exp,
    ∇y_logπy_exp;
    test_type_stability = true,
    test_logdetJ_consistency = true
)

    @testset "free" begin
        y2 = free(c, x)
        @test _size(y2) == (free_dimension(c),)
        @test y2 ≈ y_exp
        x3 = constrain(c, y2)
        @test _size(x3) == (constrain_dimension(c),)
        @test x3 ≈ x
    end

    @testset "constrain" begin
        x2 = constrain(c, y)
        @test _size(x2) == (constrain_dimension(c),)
        @test x2 ≈ x
    end

    @testset "constrain_with_pushlogpdf" begin
        x2, pushlogpdf = constrain_with_pushlogpdf(c, y)
        @test x2 ≈ constrain(c, y)
        logπy = pushlogpdf(logπx)
        @test isreal(logπy)
        @test logπy ≈ logπy_exp
    end

    @testset "constrain_with_pushlogpdf_grad" begin
        x2, pushlogpdf_grad = constrain_with_pushlogpdf_grad(c, y)
        @test x2 ≈ constrain(c, y)
        logπy, ∇y_logπy = pushlogpdf_grad(logπx, ∇x_logπx)
        @test isreal(logπy)
        @test logπy ≈ logπy_exp
        @test _size(∇y_logπy) == (free_dimension(c),)
        @test ∇y_logπy ≈ ∇y_logπy_exp
    end

    @testset "jacobian consistency" begin
        J = HMCUtilities.constrain_jacobian(c, y)
        J2, back = Zygote.pullback(HMCUtilities.constrain_jacobian, c, y)
        @test J ≈ J2
    end

    if test_logdetJ_consistency
        @testset "logabsdet jacobian consistency" begin
            logdetJ = logabsdetjac(c, y)
            _, pushlogpdf = constrain_with_pushlogpdf(c, y)
            @test logdetJ ≈ pushlogpdf(zero(eltype(y)))

            ∇y_logdetJ = _format_grad(zero(∇y_logπy_exp), Zygote.gradient(logabsdetjac, c, y)[2])
            _, pushlogpdf_grad = constrain_with_pushlogpdf_grad(c, y)
            logdetJ2, ∇y_logdetJ2 = pushlogpdf_grad(zero(logπx), zero(∇x_logπx))
            @test logdetJ ≈ logdetJ2
            @test isapprox(∇y_logdetJ, ∇y_logdetJ2; atol=1e-6)
        end
    end

    if test_type_stability
        @testset "type stability" begin
            @inferred free(c, x)
            @inferred constrain(c, y)

            @inferred (y -> constrain_with_pushlogpdf(c, y)[1])(y)
            _, pushlogpdf = constrain_with_pushlogpdf(c, y)
            @inferred pushlogpdf(logπx)

            @inferred (y -> constrain_with_pushlogpdf_grad(c, y)[1])(y)
            _, pushlogpdf_grad = constrain_with_pushlogpdf_grad(c, y)
            @inferred pushlogpdf_grad(logπx, ∇x_logπx)
        end
    end
end


text:
# TODO: Document this.
macro where(src, exprs...)
    for expr in exprs
        # TODO: Add proper error messages here.
        @assert !has_alias(expr)
    end

    e = fuse_conjunction(exprs)

    Expr(
        :call,
        :Selection,
        esc(src),
        expression_expr(e),
    )
end


text:
#module Tmp
using Random, BenchmarkTools, Test, InteractiveUtils, JLD, PyPlot
using Base.Iterators, StatsPlots, StatsBase
#
# include("src/params.jl")
# include("src/functions.jl")
# import .Params
# import .Func

data = load("examplereturns.jld")
stocks = data["stockreturns"][:,100:end]
stdevs = StatsBase.std(stockreturns, dims=2)
stocks = stocks ./ stdevs

stocks = abs.(stocks)
stocks = collect(flatten(stocks))
stockshist = PyPlot.hist(stocks, 100)
returns = stockshist[1]
returnsprob = returns ./ sum(returns)
bins = stockshist[2][2:end]
logbins = log10.(bins)
logreturns = log10.(returnsprob)
StatsPlots.plot(logbins, logreturns, label=)
#end # module


text:
using JuMP, MathOptInterface
using PowerModels
#
function getBranchIds(pm_data)
    branch_ids=Int64[]
    for ll in keys(pm_data["branch"])
        push!(branch_ids,parse(Int64,ll))
    end
    sort!(branch_ids)
    return branch_ids
end

function ConvertModelToDualizableForm(model::JuMP.Model)
    soc_model = MOI.Utilities.Model{Float64}()
    soc_bridged_model = MOI.Bridges.Constraint.QuadtoSOC{Float64}(soc_model)
    MOI.copy_to(soc_bridged_model,backend(model))

    dualizable_model = JuMP.Model()
    bridged_model = MOI.Bridges.Constraint.Square{Float64}(backend(dualizable_model))
    MOI.copy_to(bridged_model,soc_model)
	return dualizable_model
end

function removeRedundantConstraints(pm::AbstractPowerModel)
    # "So far, do nothing"
end
function removeRedundantConstraints(pm::AbstractConicModels)
    con_types=list_of_constraint_types(pm.model)
    n_con_types=length(con_types)
    for cc=1:n_con_types
        if con_types[cc][2]==MathOptInterface.RotatedSecondOrderCone
	        ### These constraints are presumably associated with defining the quadratic cost function for the OPF and are usually not needed here.
	        rsoc_con = all_constraints(pm.model, con_types[cc][1], con_types[cc][2])
	        n_rsoc_con = length(rsoc_con)
	        for nn=1:n_rsoc_con
		        rsoc_expr = constraint_object(rsoc_con[nn]).func
		        n_vars=length( rsoc_expr )
		        if n_vars == 3  ### "This condition is necessary for the RSOC to have the nonproductive constraint in question"
                    io2=IOBuffer()
                    print(io2,rsoc_expr[2])
                    io3=IOBuffer()
                    print(io3,rsoc_expr[3])
                    if occursin("pg", String(take!(io2))) && occursin("pg", String(take!(io3)))
                        #println("Removing constraint: ",rsoc_con[nn])
		                JuMP.delete(pm.model,rsoc_con[nn])
                    end
                end
	        end
	    end
    end
end

function replaceThermalLineLimits(pm::AbstractPowerModel)
    con_types=list_of_constraint_types(pm.model)
    n_con_types=length(con_types)
    for cc=1:n_con_types
        correct_con_type = (con_types[cc][1]==GenericQuadExpr{Float64,VariableRef} && con_types[cc][2]==MathOptInterface.LessThan{Float64})
        if correct_con_type
	        for con in all_constraints(pm.model, con_types[cc][1], con_types[cc][2])
		        expr = constraint_object(con).func
                io_con=IOBuffer()
                print(io_con,expr)
                con_str=String(take!(io_con))
                if occursin("p[(", con_str) && occursin("q[(", con_str)
                    #println("Removing constraint: ",con)
		            JuMP.delete(pm.model,con)
                end
	        end
	    end
    end
    pm_type = typeof(pm)
    for l in ids(pm,pm.cnw,:branch)
        branch = ref(pm, pm.cnw, :branch, l)
        if haskey(branch, "rate_a")
            f_bus,t_bus = branch["f_bus"],branch["t_bus"]
            f_idx,t_idx = (l, f_bus, t_bus),(l, t_bus, f_bus)
            #re-add constraints using auxiliary variable proxies for power flow variables
            pf_m = var(pm, pm.cnw, :pt_br)[f_idx,0]
            pf_p = var(pm, pm.cnw, :pt_br)[f_idx,1]
            pt_m = var(pm, pm.cnw, :pt_br)[t_idx,0]
            pt_p = var(pm, pm.cnw, :pt_br)[t_idx,1]
            qf_m = var(pm, pm.cnw, :qt_br)[f_idx,0]
            qf_p = var(pm, pm.cnw, :qt_br)[f_idx,1]
            qt_m = var(pm, pm.cnw, :qt_br)[t_idx,0]
            qt_p = var(pm, pm.cnw, :qt_br)[t_idx,1]
            cref_f=JuMP.@constraint(pm.model, (pf_p-pf_m)^2 + (qf_p-qf_m)^2 <= branch["rate_a"]^2)
            JuMP.set_name(cref_f,string("th_l_lim_f[",f_idx,"]"))
            #println("Replacing the 'from' thermal line limit with: ",cref_f)
            cref_t=JuMP.@constraint(pm.model, (pt_p-pt_m)^2 + (qt_p-qt_m)^2 <= branch["rate_a"]^2)
            JuMP.set_name(cref_t,string("th_l_lim_t[",t_idx,"]"))
            #println("Replacing the 'to' thermal line limit with: ",cref_t)
        end
    end
end
function replaceThermalLineLimits(pm::AbstractConicModels)
    con_types=list_of_constraint_types(pm.model)
    n_con_types=length(con_types)
    for cc=1:n_con_types
        correct_con_type = con_types[cc][2]==MathOptInterface.SecondOrderCone
        if correct_con_type
	        for con in all_constraints(pm.model, con_types[cc][1], con_types[cc][2])
		        expr = constraint_object(con).func
                io_con=IOBuffer()
                print(io_con,expr)
                con_str=String(take!(io_con))
                if occursin("p[(", con_str) && occursin("q[(", con_str)
                    #println("Removing constraint: ",con)
		            JuMP.delete(pm.model,con)
                end
	        end
	    end
    end
    for l in ids(pm,pm.cnw,:branch)
        branch = ref(pm, pm.cnw, :branch, l)
        if haskey(branch, "rate_a")
            f_bus,t_bus = branch["f_bus"],branch["t_bus"]
            f_idx,t_idx = (l, f_bus, t_bus),(l, t_bus, f_bus)
            #re-add constraints using auxiliary variable proxies for power flow variables
            pf_m = var(pm, pm.cnw, :pt_br)[f_idx,0]
            pf_p = var(pm, pm.cnw, :pt_br)[f_idx,1]
            pt_m = var(pm, pm.cnw, :pt_br)[t_idx,0]
            pt_p = var(pm, pm.cnw, :pt_br)[t_idx,1]
            qf_m = var(pm, pm.cnw, :qt_br)[f_idx,0]
            qf_p = var(pm, pm.cnw, :qt_br)[f_idx,1]
            qt_m = var(pm, pm.cnw, :qt_br)[t_idx,0]
            qt_p = var(pm, pm.cnw, :qt_br)[t_idx,1]
            cref_f=JuMP.@constraint(pm.model, [branch["rate_a"], pf_p-pf_m, qf_p-qf_m] in JuMP.SecondOrderCone())
            JuMP.set_name(cref_f,string("th_l_lim_f[",f_idx,"]"))
            #println("Replacing the 'from' thermal line limit with: ",cref_f)
            cref_t=JuMP.@constraint(pm.model, [branch["rate_a"], pt_p-pt_m, qt_p-qt_m] in JuMP.SecondOrderCone())
            JuMP.set_name(cref_t,string("th_l_lim_t[",t_idx,"]"))
            #println("Replacing the 'to' thermal line limit with: ",cref_t)
        end
    end
end

# "Precondition: replaceThermalLineLimits() has been called"
function removeThermalLineLimits(pm::AbstractPowerModel)
    for l in ids(pm,pm.cnw,:branch)
        branch = ref(pm, pm.cnw, :branch, l)
        if haskey(branch, "rate_a")
            f_bus,t_bus = branch["f_bus"],branch["t_bus"]
            f_idx,t_idx = (l, f_bus, t_bus),(l, t_bus, f_bus)
            JuMP.delete(pm.model,JuMP.constraint_by_name(pm.model,string("th_l_lim_f[",f_idx,"]")))
            JuMP.delete(pm.model,JuMP.constraint_by_name(pm.model,string("th_l_lim_t[",t_idx,"]")))
        end
    end
end

function getMaxShedding(pm_data)
    maxShed=0.0
    for ll in keys(pm_data["load"])
        maxShed += abs(pm_data["load"][ll]["pd"]) + abs(pm_data["load"][ll]["qd"])
    end
    for gg in keys(pm_data["gen"])
        maxShed += max(pm_data["gen"][gg]["pmin"],0) + max(pm_data["gen"][gg]["qmin"],0)
    end
    return maxShed
end

function convertSOCtoPSD(model::JuMP.Model)
    #BRIDGE SOC CONSTRAINTS
    model_rsoc_moi = MOI.Utilities.Model{Float64}()
    rsoc_bridged_model = MOI.Bridges.Constraint.SOCtoPSD{Float64}(model_rsoc_moi)
    MOI.copy_to(rsoc_bridged_model,backend(model))

    model_psd_moi = MOI.Utilities.Model{Float64}()
    psd_bridged_model = MOI.Bridges.Constraint.RSOCtoPSD{Float64}(model_psd_moi)
    MOI.copy_to(psd_bridged_model,model_rsoc_moi)
    model_psd = JuMP.Model()
    MOI.copy_to(backend(model_psd),model_psd_moi)
    return model_psd
end

function unfix_vars(model,branch_ids)
    for l in branch_ids
        if is_fixed(variable_by_name(model,"x[$l]_1"))
            unfix(variable_by_name(model,"x[$l]_1"))
        end
        set_lower_bound(variable_by_name(model,"x[$l]_1"),0)
        set_upper_bound(variable_by_name(model,"x[$l]_1"),1)
    end
end

function fix_integer_vals(model_info)
    relax_integrality(model_info)
    for l in model_info["branch_ids"]
        if is_fixed(variable_by_name(model_info["model"],"x[$l]_1"))
            unfix(variable_by_name(model_info["model"],"x[$l]_1"))
        end
        fix(variable_by_name(model_info["model"],"x[$l]_1"), model_info["x_soln"][l]; force=true)
    end
end

function relax_integrality(model_dict)
    for l in model_dict["branch_ids"]
	    if JuMP.is_fixed(variable_by_name(model_dict["model"],"x[$l]_1"))
            unfix(variable_by_name(model_dict["model"],"x[$l]_1"))
        end
	    if JuMP.is_integer(variable_by_name(model_dict["model"],"x[$l]_1"))
            JuMP.unset_integer(variable_by_name(model_dict["model"],"x[$l]_1"))
            JuMP.set_lower_bound(variable_by_name(model_dict["model"],"x[$l]_1"),0)
            JuMP.set_upper_bound(variable_by_name(model_dict["model"],"x[$l]_1"),1)
	    else
            JuMP.set_lower_bound(variable_by_name(model_dict["model"],"x[$l]_1"),0)
            JuMP.set_upper_bound(variable_by_name(model_dict["model"],"x[$l]_1"),1)
	    end
    end
end

### "Copied from current master version of JuMP, as of 7 Feb 2020. At some point, we can use the JuMP stable version."
function JuMP_write_to_file(
    model::Model,
    filename::String;
    format::MOI.FileFormats.FileFormat = MOI.FileFormats.FORMAT_AUTOMATIC
)
    dest = MOI.FileFormats.Model(format = format, filename = filename)
    # We add a `full_bridge_optimizer` here because MOI.FileFormats models may not
    # support all constraint types in a JuMP model.
    bridged_dest = MOI.Bridges.full_bridge_optimizer(dest, Float64)
    MOI.copy_to(bridged_dest, backend(model))
    # `dest` will contain the underlying model, with constraints bridged if
    # necessary.
    MOI.write_to_file(dest, filename)
    return
end

function write_to_cbf(model,fn_base::String)
    JuMP.write_to_file( model, string(fn_base,".cbf"), format = MOI.FileFormats.FORMAT_CBF)
end

function add_artificial_var_bds(model::JuMP.Model; bd_mag=1e3)
    all_vars = JuMP.all_variables(model)
    for vv in all_vars
        if !has_lower_bound(vv)
            JuMP.@constraint(model, vv >= -bd_mag )
        end
        if !has_upper_bound(vv)
            JuMP.@constraint(model, vv <= bd_mag )
        end
    end
end

function write_to_cbf_scip(model,fn_base::String)
    model_psd = convertSOCtoPSD(model)
    #add_artificial_var_bds(model_psd)
    fname=string(fn_base,"_scip",".cbf")
    JuMP_write_to_file( model_psd, fname, format = MOI.FileFormats.FORMAT_CBF)
    ### CHANGING VER 3 to VER 2
    fix_version=`sed -i -z 's/VER\n3/VER\n2/g' $fname`
    run(fix_version)
end


text:
######################################################################
# termination condition stuff
#
# TODO: I probably don't need to pass in step_count, the step
#  threshold could simply keep track of this as a field
#
# Although maybe it is useful to have it passed in?
######################################################################

abstract type TerminationCondition end

# determines if a simulation has reached its
function is_complete(f::AbstractFilter, ::TerminationCondition, ::Int)
    error("Termination condition not implemented for this filter or termination condition")
end

struct StepThreshold <: TerminationCondition
    value::Int
end
function is_complete(f::AbstractFilter, st::StepThreshold, step_count::Int)
    ret_val = false
    if step_count >= st.value
        ret_val = true
    end
    return ret_val
end

struct MaxNormThreshold <: TerminationCondition
    value::Float64
end
function is_complete(f::DF, mnt::MaxNormThreshold, ::Int)
    ret_val = false
    if maximum(f.b) > mnt.value
        ret_val = true
    end
    return ret_val
end

struct EntropyThreshold <: TerminationCondition
    value::Float64
end

export OrTC
struct OrTC <: TerminationCondition
    conditions::Vector{TerminationCondition}
end
function is_complete(f::AbstractFilter, st::OrTC, step_count::Int)
    true_sum = 0
    for tc in st.conditions
        if is_complete(f, tc, step_count)
            return true
        end
    end
    return false
end


text:
using Flux, Distributions, Random
using Flux: onehotbatch, Zygote
using AdvancedHMC
using LinearAlgebra
using Plots
include("turing_data.jl")

struct BayesNNet{N, P, RE}
    nnet::N
    p::P
    re::RE

    function BayesNNet(nnet;p = nothing)
        _p, re = Flux.destructure(nnet)
        if p === nothing
            p = _p
        end
        new{typeof(nnet), typeof(p), typeof(re)}(
            nnet, p, re)
    end
end

(b::BayesNNet)(x, p=b.p) = b.re(p)(x)

net_topology = Chain(Dense(2,3,tanh), Dense(3,2,tanh), Dense(2,2,σ))
bnn = BayesNNet(net_topology)

x_train = Flux.normalise(X,dims=2)
y_train = onehotbatch(Y, 0:1)

# log posterior density
alpha = 0.09
prior_sigma = sqrt(1/alpha)
function ℒ(p)
    -sum(Flux.binarycrossentropy.(bnn(x_train, p), y_train)) -
        0.5*norm(p)^2 / prior_sigma^2
end

# Set the number of samples to draw and warmup iterations
n_samples, n_adapts = 6_000, 2_000

# Define a Hamiltonian system
D = length(bnn.p)
metric = DiagEuclideanMetric(D)
hamiltonian = Hamiltonian(metric, ℒ, Zygote)

Random.seed!(1)
initial_p = randn(D)

# Define a leapfrog solver, with initial step size chosen heuristically
initial_ϵ = find_good_stepsize(hamiltonian, initial_p)
integrator = Leapfrog(initial_ϵ)

# Define an HMC sampler, with the following components
#   - multinomial sampling scheme,
#   - generalised No-U-Turn criteria, and
#   - windowed adaption for step-size and diagonal mass matrix
proposal = NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator)
adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(initial_ϵ, integrator))

# Run the sampler to draw samples from the specified Gaussian, where
#   - `samples` will store the samples
#   - `stats` will store diagnostic statistics for each sample
Random.seed!(123)
samples, stats = sample(hamiltonian, proposal, initial_p, n_samples, adaptor, n_adapts;
                        drop_warmup=true, progress=false)

# define test set
x_range = collect(range(-6,stop=6,length=25))
y_range = collect(range(-6,stop=6,length=25))
n_end = size(samples, 1)
anim = @gif for i=1:20:n_end
    plot_data();
    Z = [bnn([x, y], samples[i])[1] for x=x_range, y=y_range]
    contour!(x_range, y_range, Z, title="Iteration $i", clim = (0,1));
end every 30;


text:
@testset "1576.replace-all-s-to-avoid-consecutive-repeating-characters.jl" begin
    @test modify_string("j?qg??b") == "jaqgacb"
    @test modify_string("ubv?w") == "ubvaw"
    @test modify_string("?zs") == "azs"
    @test modify_string("??yw?ipkj?") == "abywaipkja"
    @test modify_string("b?a") == "bca"
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Rotations from an Earth-Fixed Inertial (ECI) reference frame to an
#   Earth-Fixed, Earth-Centered (ECEF) reference frame.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export r_eci_to_ecef

"""
    r_eci_to_ecef([T,] ECI, ECEF, JD_UTC::Number [, eop_data])

Compute the rotation from an Earth-Centered Inertial (`ECI`) reference frame to
an Earth-Centered, Earth-Fixed (`ECEF`) reference frame at the Julian Day [UTC]
`JD_UTC`. The rotation description that will be used is given by `T`, which can
be `DCM` or `Quaternion`. The ECI frame is selected by the input `ECI` and the
`ECEF` frame is selected by the input `ECEF`. The possible values are listed
below. The model used to compute the rotation is specified by the selection of
the origin and destination frames. Currently, there are two models supported:
IAU-76/FK5 and IAU-2006 with 2010 conventions (CIO and equinox approaches).

# Rotation description

The rotations that aligns the ECI with ECEF can be described by Direction Cosine
Matrices or Quaternions. This is selected by the parameter `T`. The possible
values are:

- `DCM`: The rotation will be described by a Direction Cosine Matrix.
- `Quaternion`: The rotation will be described by a Quaternion.

If no value is specified, then it falls back to `DCM`.

# Conversion model

The model that will be used to compute the rotation is automatically inferred
given the selection of the origin and destination frames. **Notice that mixing
IAU-76/FK5 and IAU-2006/2010 frames is not supported.**

# ECI Frame

The ECI frame is selected by the parameter `ECI`. The possible values are:

- `TEME()`: ECI will be selected as the True Equator Mean Equinox (TEME)
    reference frame.
- `TOD()`: ECI will be selected as the True of Date (TOD).
- `MOD()`: ECI will be selected as the Mean of Date (MOD).
- `J2000()`: ECI will be selected as the J2000 reference frame.
- `GCRF()`: ECI will be selected as the Geocentric Celestial Reference Frame
    (GCRF).
- `CIRS()`: ECEF will be selected as the Celestial Intermediate Reference System
    (CIRS).
- `ERS()`: ECI will be selected as the Earth Reference System (ERS).
- `MOD06()`: ECI will be selected as the Mean of Date (MOD) according to the
    definition in IAU-2006/2010 theory.
- `MJ2000()`: ECI will be selected as the J2000 mean equatorial frame (MJ2000).

!!! note
    The frames `MOD()` and `MOD06()` are virtually the same. However, we
    selected different names to make clear which theory are being used since
    mixing transformation between frames from IAU-76/FK5 and IAU-2006/2010 must
    be performed with caution.

# ECEF Frame

The ECEF frame is selected by the parameter `ECEF`. The possible values are:

- `ITRF()`: ECEF will be selected as the International Terrestrial Reference
    Frame (ITRF).
- `PEF()`: ECEF will be selected as the Pseudo-Earth Fixed (PEF) reference
    frame.
- `TIRS()`: ECEF will be selected as the Terrestrial Intermediate Reference
    System (TIRS).

# EOP Data

The conversion between the frames depends on EOP Data (see
[`get_iers_eop`](@ref) and [`read_iers_eop`](@ref)). If IAU-76/FK5 model is
used, then the type of `eop_data` must be [`EOPData_IAU1980`](@ref). Otherwise,
if IAU-2006/2010 model is used, then the type of `eop_data` must be
[`EOPData_IAU2000A`](@ref). The following table shows the requirements for EOP
data given the selected frames.

|   Model                     |   ECI    |  ECEF  |    EOP Data     |
|:----------------------------|:---------|:-------|:----------------|
| IAU-76/FK5                  | `GCRF`   | `ITRF` | EOP IAU1980     |
| IAU-76/FK5                  | `J2000`  | `ITRF` | EOP IAU1980     |
| IAU-76/FK5                  | `MOD`    | `ITRF` | EOP IAU1980     |
| IAU-76/FK5                  | `TOD`    | `ITRF` | EOP IAU1980     |
| IAU-76/FK5                  | `TEME`   | `ITRF` | EOP IAU1980     |
| IAU-76/FK5                  | `GCRF`   | `PEF`  | EOP IAU1980     |
| IAU-76/FK5                  | `J2000`  | `PEF`  | Not required¹   |
| IAU-76/FK5                  | `MOD`    | `PEF`  | Not required¹   |
| IAU-76/FK5                  | `TOD`    | `PEF`  | Not required¹   |
| IAU-76/FK5                  | `TEME`   | `PEF`  | Not required¹   |
| IAU-2006/2010 CIO-based     | `CIRS`   | `ITRF` | EOP IAU2000A    |
| IAU-2006/2010 CIO-based     | `GCRF`   | `ITRF` | EOP IAU2000A    |
| IAU-2006/2010 CIO-based     | `CIRS`   | `TIRS` | Not required¹   |
| IAU-2006/2010 CIO-based     | `GCRF`   | `TIRS` | Not required¹ ² |
| IAU-2006/2010 Equinox-based | `ERS`    | `TIRS` | EOP IAU2000A    |
| IAU-2006/2010 Equinox-based | `MOD06`  | `ITRF` | EOP IAU2000A    |
| IAU-2006/2010 Equinox-based | `MJ2000` | `ITRF` | EOP IAU2000A    |
| IAU-2006/2010 Equinox-based | `ERS`    | `TIRS` | Not required¹ ³ |
| IAU-2006/2010 Equinox-based | `MOD06`  | `TIRS` | Not required¹ ³ |
| IAU-2006/2010 Equinox-based | `MJ2000` | `TIRS` | Not required¹ ³ |

`¹`: In this case, UTC will be assumed equal to UT1 to compute the Greenwich
Mean Sidereal Time. This is an approximation, but should be sufficiently
accurate for some applications. Notice that, if EOP Data is provided, UT1 will
be accurately computed.

`²`: In this case, the terms that account for the free-core nutation and time
dependent effects of the Celestial Intermediate Pole (CIP) position with respect
to the GCRF will not be available, reducing the precision.

## MOD and TOD

In this function, if EOP corrections are not provided, then MOD and TOD frames
will be computed considering the original IAU-76/FK5 theory. Otherwise, the
corrected frame will be used.

# Returns

The rotation description represented by `T` that rotates the ECI reference frame
into alignment with the ECEF reference frame.

# Examples

```julia-repl
julia> eop_IAU1980 = get_iers_eop(Val(:IAU1980));

julia> r_eci_to_ecef(DCM, GCRF(), ITRF(), date_to_jd(1986, 06, 19, 21, 35, 0), eop_IAU1980)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.619267    -0.78518     -0.000797312
  0.78518     -0.619267     0.00106478
 -0.00132979   3.33509e-5   0.999999

julia> r_eci_to_ecef(GCRF(), ITRF(), date_to_jd(1986, 06, 19, 21, 35, 0), eop_IAU1980)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.619267    -0.78518     -0.000797312
  0.78518     -0.619267     0.00106478
 -0.00132979   3.33509e-5   0.999999

julia> r_eci_to_ecef(J2000(), PEF(), date_to_jd(1986, 06, 19, 21, 35, 0))
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.619271    -0.785177    -0.000796885
  0.785176    -0.619272     0.00106622
 -0.00133066   3.45854e-5   0.999999

julia> r_eci_to_ecef(J2000(), PEF(), date_to_jd(1986, 06, 19, 21, 35, 0), eop_IAU1980)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.619267    -0.78518     -0.000796879
  0.78518     -0.619267     0.00106623
 -0.00133066   3.45854e-5   0.999999

julia> r_eci_to_ecef(Quaternion, GCRF(), ITRF(), date_to_jd(1986, 06, 19, 21, 35, 0), eop_IAU1980)
Quaternion{Float64}:
  + 0.43631 + 0.000590997⋅i - 0.000305106⋅j - 0.000305106⋅k

julia> eop_IAU2000A = get_iers_eop(Val(:IAU2000A));

julia> r_eci_to_ecef(GCRF(), ITRF(), date_to_jd(1986, 06, 19, 21, 35, 0), eop_IAU2000A)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.619267    -0.78518     -0.000797311
  0.78518     -0.619267     0.00106478
 -0.00132979   3.33516e-5   0.999999

julia> r_eci_to_ecef(GCRF(), TIRS(), date_to_jd(1986, 06, 19, 21, 35, 0))
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.619271    -0.785177    -0.000796885
  0.785176    -0.619272     0.00106623
 -0.00133066   3.45884e-5   0.999999

julia> r_eci_to_ecef(Quaternion, GCRF(), ITRF(), date_to_jd(1986, 06, 19, 21, 35, 0), eop_IAU2000A)
Quaternion{Float64}:
  + 0.43631 + 0.000590997⋅i - 0.000305106⋅j - 0.000305106⋅k
```
"""
@inline function r_eci_to_ecef(T_ECI::T_ECIs,
    T_ECEF::T_ECEFs,
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    return r_eci_to_ecef(DCM, T_ECI, T_ECEF, JD_UTC, eop_data)
end

@inline function r_eci_to_ecef(
    T::T_ROT,
    T_ECI::T_ECIs,
    T_ECEF::T_ECEFs,
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    return inv_rotation(r_ecef_to_eci(T, T_ECEF, T_ECI, JD_UTC, eop_data))
end

@inline function r_eci_to_ecef(
    T_ECI::T_ECIs_IAU_2006,
    T_ECEF::T_ECEFs_IAU_2006,
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    return r_eci_to_ecef(DCM, T_ECI, T_ECEF, JD_UTC, eop_data)
end

@inline function r_eci_to_ecef(
    T::T_ROT,
    T_ECI::T_ECIs_IAU_2006,
    T_ECEF::T_ECEFs_IAU_2006,
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    return inv_rotation(r_ecef_to_eci(T, T_ECEF, T_ECI, JD_UTC, eop_data))
end

# Specializations for those cases that EOP Data is not needed.
@inline function r_eci_to_ecef(
    T_ECI::Union{Val{:J2000}, Val{:TOD}, Val{:MOD}, Val{:TEME}},
    T_ECEF::Val{:PEF},
    JD_UTC::Number
)
    return r_eci_to_ecef(DCM, T_ECI, T_ECEF, JD_UTC)
end

@inline function r_eci_to_ecef(
    T::T_ROT,
    T_ECI::Union{Val{:J2000}, Val{:TOD}, Val{:MOD}, Val{:TEME}},
    T_ECEF::Val{:PEF},
    JD_UTC::Number
)
    return inv_rotation(r_ecef_to_eci(T, T_ECEF, T_ECI, JD_UTC))
end

@inline function r_eci_to_ecef(
    T_ECI::T_ECIs_IAU_2006,
    T_ECEF::Val{:TIRS},
    JD_UTC::Number
)
    return r_eci_to_ecef(DCM, T_ECI, T_ECEF, JD_UTC)
end

@inline function r_eci_to_ecef(
    T::T_ROT,
    T_ECI::T_ECIs_IAU_2006,
    T_ECEF::Val{:TIRS},
    JD_UTC::Number
)
    return inv_rotation(r_ecef_to_eci(T, T_ECEF, T_ECI, JD_UTC))
end


text:
# Load a common workspace and solve a lasso problem
export jl_complex

function jl_complex()
    exampledir = Pkg.dir("GenSPGL")*"/src/Examples/compare/"
    file = MAT.matopen(exampledir*"vars_complex.mat")
    b = read(file,"b")
    tau = read(file,"tau")
    test = read(file,"test")
    xinit = read(file,"xinit")
    sigma = read(file,"sigma")
    close(file)

    # Choose options for GenSPGL
    opts = spgOptions(  optTol = 1e-5,
                        bpTol = 1e-5,
                        decTol = 1e-4,
                        project = TraceNorm_project,
                        primal_norm = TraceNorm_primal,
                        dual_norm = TraceNorm_dual,
                        proxy = true,
                        ignorePErr = true,
                        iterations = 150,
                        verbosity = 1)

    # Create Params Dict
    afunT(x) = reshape(x,355,709)
    params = Dict("nr"=> 40,
                    "Ind"=> vec(b) .== 0,
                    "numr"=> 355,
                    "numc"=> 709,
                    "funForward"=> NLfunForward,
                    "afunT"=> afunT,
                    "afun"=> afun,
                    "mode"=> 1,
                    "ls"=> 1,
                    "logical"=> 0,
                    "funPenalty"=> funLS)

    # Set threads and solve problems
    #BL.set_num_threads(8)
    @time xLS_jl, r, g, info = spgl1(NLfunForward, vec(b),
                                                   x = vec(xinit),
                                                   tau = tau,
                                                   sigma = sigma,
                                                   options = opts,
                                                   params = params)

    # Write output for MATLAB comparison
    fw = MAT.matopen(exampledir*"xLS_jl.mat", "w")
    write(fw, "xLS_jl", xLS_jl)
    close(fw)

    return xLS_jl, r, g, info
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct OnEventCall <: ProtoType
    _type::Int32
    arguments::EventArguments
    OnEventCall(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct OnEventCall

mutable struct GetTeamInfoCall <: ProtoType
    server_version::Int32
    GetTeamInfoCall(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetTeamInfoCall

mutable struct GetTeamInfoResult <: ProtoType
    team_info::TeamInfo
    GetTeamInfoResult(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetTeamInfoResult

mutable struct GetInstructionCall <: ProtoType
    field::Field
    GetInstructionCall(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetInstructionCall

mutable struct GetInstructionResult <: ProtoType
    wheels::Base.Vector{Wheel}
    command::ControlInfo
    GetInstructionResult(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetInstructionResult

mutable struct GetPlacementCall <: ProtoType
    field::Field
    GetPlacementCall(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetPlacementCall

mutable struct RPCCall <: ProtoType
    on_event::OnEventCall
    get_team_info::GetTeamInfoCall
    get_instruction::GetInstructionCall
    get_placement::GetPlacementCall
    RPCCall(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct RPCCall
const __oneofs_RPCCall = Int[1,1,1,1]
const __oneof_names_RPCCall = [Symbol("method")]
meta(t::Type{RPCCall}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, __oneofs_RPCCall, __oneof_names_RPCCall, ProtoBuf.DEF_FIELD_TYPES)

mutable struct GetPlacementResult <: ProtoType
    placement::Placement
    GetPlacementResult(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetPlacementResult

export RPCCall, OnEventCall, GetTeamInfoCall, GetTeamInfoResult, GetInstructionCall, GetInstructionResult, GetPlacementCall, GetPlacementResult


text:
# ---
# title: 1616. Split Two Strings to Make Palindrome
# id: problem1616
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Two Pointers, String, Greedy
# link: <https://leetcode.com/problems/split-two-strings-to-make-palindrome/description/>
# hidden: true
# ---
#
# You are given two strings `a` and `b` of the same length. Choose an index and
# split both strings **at the same index** , splitting `a` into two strings:
# `aprefix` and `asuffix` where `a = aprefix + asuffix`, and splitting `b` into
# two strings: `bprefix` and `bsuffix` where `b = bprefix + bsuffix`. Check if
# `aprefix + bsuffix` or `bprefix + asuffix` forms a palindrome.
#
# When you split a string `s` into `sprefix` and `ssuffix`, either `ssuffix` or
# `sprefix` is allowed to be empty. For example, if `s = "abc"`, then `"" +
# "abc"`, `"a" + "bc"`, `"ab" + "c"` , and `"abc" + ""` are valid splits.
#
# Return `true` _if it is possible to form_ _a palindrome string, otherwise
# return_`false`.
#
# **Notice** that `x + y` denotes the concatenation of strings `x` and `y`.
#
#
#
# **Example 1:**
#
#
#
#     Input: a = "x", b = "y"
#     Output: true
#     **Explaination:** If either a or b are palindromes the answer is true since you can split in the following way:
#     aprefix = "", asuffix = "x"
#     bprefix = "", bsuffix = "y"
#     Then, aprefix + bsuffix = "" + "y" = "y", which is a palindrome.
#
#
# **Example 2:**
#
#
#
#     Input: a = "abdef", b = "fecab"
#     Output: true
#
#
# **Example 3:**
#
#
#
#     Input: a = "ulacfd", b = "jizalu"
#     Output: true
#     **Explaination:** Split them at index 3:
#     aprefix = "ula", asuffix = "cfd"
#     bprefix = "jiz", bsuffix = "alu"
#     Then, aprefix + bsuffix = "ula" + "alu" = "ulaalu", which is a palindrome.
#
#
# **Example 4:**
#
#
#
#     Input: a = "xbdef", b = "xecab"
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `1 <= a.length, b.length <= 105`
#   * `a.length == b.length`
#   * `a` and `b` consist of lowercase English letters
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using DifferentiableStateSpaceModels, Symbolics, LinearAlgebra, Test
using DifferentiableStateSpaceModels.Examples

@testset "FVGQ20 First Order" begin
    isdefined(Main, :FVGQ20) || include(joinpath(pkgdir(DifferentiableStateSpaceModels), "test/generated_models/FVGQ20.jl"))
    m = PerturbationModel(Main.FVGQ20)
    p_d = (β=0.998, h=0.97, ϑ=1.17, κ=9.51, α=0.21, θp=0.82, χ=0.63, γR=0.77, γy=0.19, γΠ=1.29, Πbar=1.01, ρd=0.12, ρφ=0.93, ρg=0.95, g_bar=0.3, σ_A=exp(-3.97), σ_d=exp(-1.51), σ_φ=exp(-2.36), σ_μ=exp(-5.43), σ_m=exp(-5.85), σ_g=exp(-3.0), Λμ=3.4e-3, ΛA=2.8e-3)
    p_f = (δ=0.025, ε=10, ϕ=0, γ2=0.001, Ω_ii=sqrt(1e-5))

    c = SolverCache(m, Val(1), p_d)
    sol = generate_perturbation(m, p_d, p_f; cache=c)
    generate_perturbation_derivatives!(m, p_d, p_f, c)

    @test sol.retcode == :Success
    @test sol.A ≈ [0.9597066320040756 -0.0016717402006548444 -0.0028476359898478745 2.150192831740316e-18 -0.006077464817194368 0.00022438953543368282 -0.002315915210730716 0.0002480013604186878 0.027777206268578707 -0.00383188379414908 -0.209970722216548 -0.9904138263557944 -0.008023324298765148 -0.004801556142159064; -0.06116209317503432 0.9081857564695927 -0.07117720022389305 1.4035175039973e-17 -0.1589622523374888 0.005869135737739553 -0.04337176006047146 0.001967204841407287 -0.0020170583816015585 -0.05634485080808087 -0.07734742427383892 -0.2888481452521955 -0.20985817937725523 -0.06009042851246476; 0.047149679377472 0.05989700727193322 0.5175441967636624 -8.785633176355972e-18 -0.19497474630453365 0.007198773511681691 0.06605149935298901 -0.001731691775730205 0.0014797381625939485 0.05674820197564877 0.0010549545833333184 -0.06512643814483171 -0.2574010161679107 0.05252057306976145; 0.525972081627061 0.7067647675166118 -0.8250514428090694 -2.0410009966095642e-16 -3.139152046485734 0.1159024182868347 1.0968606775309409 -0.017539660702373366 0.013562717425228943 0.8004565581415306 -0.024947971816113274 -0.8294252329988305 -4.144233763530546 0.6080461766818812; 0.049282152560429886 0.05468765564912119 0.14582797660460228 5.256507747973834e-17 0.6831385774943626 -0.02522254799519321 0.03467394774591567 -0.00174994781502999 0.0014404882645437057 0.020696434081023268 0.014595652596892201 -0.001347799955404436 0.9018632790315961 0.05276005651146111; 1.2382273640942494 1.294467918405045 -0.30675891690633167 -2.1373930578377843e-17 -1.0073328208733665 0.037192308059684084 0.5258480374286546 -0.04336430305608645 0.03512566747437948 0.13191493606079913 0.0883642693379878 -1.3357085443757348 -1.3298567974907014 1.3041462486643989; 0.009599564765915364 0.012194890997879147 -0.0228957392867415 -1.3387600064770917e-17 -0.03969640365736839 0.0014656534991444518 0.8641998463024083 -0.0003525684071545177 0.00030127123908204754 0.011553801582813085 0.0002147862930288743 -0.013259590928356465 -0.05240624662057069 0.01069306619655468; -0.0611620931749605 0.9081857564697143 -0.0711772002238785 -8.141584778182808e-17 -0.15896225233760544 0.005869135737619167 -0.043371760060328696 0.9693452541489845 -0.002017058381581333 -0.05634485080799174 -10.48861951082978 -10.706368869464443 -0.2098581793773417 -0.060090428512097956; -1.2274387057378921e-18 1.129716330863846e-18 -2.6452513161575937e-18 4.303366036318858e-20 1.932605834302807e-18 -4.9112240214915196e-20 -5.5502118613030306e-18 3.2874122628061265e-19 0.12 -4.230985837363517e-18 3.0274528284264972e-18 3.2292323574232476e-18 1.7423640484797458e-18 -1.6502917887199003e-17; 6.5743501365246386e-18 1.69567809349733e-17 1.8453482603434553e-18 1.4104454171393302e-17 1.0721975244983874e-16 -5.853396590979275e-18 2.3016001541456748e-17 -4.972737819827135e-19 -3.546634649666224e-18 0.9300000000000004 3.521544976745966e-16 4.400525620152486e-16 2.113466586628283e-16 2.707831559404897e-16; 2.2568972501382555e-17 1.1941207210767652e-16 -7.886281870178793e-17 -8.262052872093175e-31 -1.7082021990020085e-16 6.306950502338531e-18 1.0385973714711463e-16 -9.933762539483766e-18 7.997775589584522e-19 6.113195184305982e-17 8.986093821921679e-17 9.333325263923877e-17 -6.301399549422952e-17 -6.637664205478931e-17; -2.3295865035336392e-17 -3.648847602877501e-18 2.5686152263757588e-17 -9.1576245100504e-31 5.249861165160181e-17 -1.9383311031746743e-18 -2.9498949743095986e-17 2.4517184677316216e-19 -7.10836567443952e-19 -2.014313976204385e-17 3.4655363692211897e-18 1.5444545238035702e-17 1.0724413225403257e-17 8.242097119011588e-17; -1.8469211863662103e-16 -1.6043964703540862e-16 6.123339946925812e-18 1.6457798003133015e-30 1.012214486991828e-16 -3.737254681960731e-18 -4.189986148448415e-17 6.012143234755371e-18 -5.662286392385015e-18 -3.254482025735823e-17 -1.1428499435437241e-17 1.891055157616275e-16 4.069714234131922e-17 -6.806721471364512e-17; -7.065490002090046e-18 -4.1123352461888506e-18 -3.786545338319949e-17 8.380886195309565e-18 -2.7207396233759008e-17 -2.88508375354103e-19 -9.860474762111539e-18 3.1037019916073786e-19 -6.578099119381701e-19 -1.0956403414738988e-17 -2.7651598835201586e-17 3.323414050181414e-17 -1.8308120469745517e-16 0.9499999999999996]
    @test sol.C ≈ [0.046682850885188414 0.05930396759586929 0.5124199967950123 0.0 -0.19304430327244906 0.007127498526411165 0.06539752411292642 -0.0017145463126460162 0.0014650872902314323 0.056186338595056756 0.0010445094854684119 -0.0644816219419623 -0.2548524912562035 0.052000567405129226; 0.04848049555975208 0.05379806906178239 0.1434558395931211 0.0 0.6720261809473577 -0.024812260881409483 0.03410991773797231 -0.0017214819741272426 0.0014170563030341846 0.020359771812883194 0.014358229788893229 -0.0013258757233242705 0.8871929577850087 0.05190182556133141; 0.4733058027331458 0.6359954782844232 -0.7424379526200766 -0.8998686798130806 -2.8248246078053247 0.10429695613070526 0.9870305698304098 -0.015783391320695707 0.012204664624805324 0.7203057862292772 -0.02244989846562421 -0.7463737894434745 -3.729266165627129 0.5471617102883314; 0.7924415036295586 0.8284343678742737 -0.19631975872429247 0.0 -0.6446734730418766 -0.616178253872422 0.33653254769120866 -0.027752313116301518 0.0224797460923251 0.08442300123162487 0.056551418981011355 -0.854827568836762 -0.851082564294199 0.8346283115774529; 0.7095449970937021 0.7149458978760908 -0.06799364123668006 0.0 -0.10822645419616717 0.003995890475364564 1.0031833138000825 -0.025718053394227534 0.020408148007625304 -0.051400061887996056 0.06913294700293195 -0.7124691120650971 -0.14287798709487276 0.7379111193434185; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.9966057734548978 0.0 0.0 0.0]
    @test sol.Q ≈ [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]
    @test sol.y ≈ [0.8154780241332105, 0.2783051143056207, 1.01, 1.1112732584577985, 1.0165356601604938, 1.562547340626902, 1.0022208173087115, 1.0186852828361062, 1.0, 1.0, 1.2953143777555123, 1.0044580087526465, 0.03489877588693624, 1.0, 0.8982498823496541, 12.044079264159267, 13.382310293510294, 0.004448101265822784, 0.009950330853168092, 0.01640043478966387, 0.0, 0.0, 0.0, 0.0034]
    @test sol.x ≈ [0.8154780241332105, 0.2783051143056207, 1.01, 1.1112732584577985, 1.0165356601604938, 1.562547340626902, 1.0022208173087115, 8.53122233338444, 1.0, 1.0, 1.0034057865562385, 1.0028039236612292, 1.0, 0.4687642021880705]
    @test sol.B ≈ [0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.2209099779593782 0.0 0.0 0.0 0.0 0.0; 0.0 0.09442022319630236 0.0 0.0 0.0 0.0; 0.0 0.0 0.004383095802668776 0.0 0.0 0.0; 0.0 0.0 0.0 0.018873433135151486 0.0 0.0; 0.0 0.0 0.0 0.0 0.002879899158088243 0.0; 0.0 0.0 0.0 0.0 0.0 0.049787068367863944]
    @test sol.η ≈ [0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 1.0 0.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0]
    @test sol.Γ ≈ [0.2209099779593782 0.0 0.0 0.0 0.0 0.0; 0.0 0.09442022319630236 0.0 0.0 0.0 0.0; 0.0 0.0 0.004383095802668776 0.0 0.0 0.0; 0.0 0.0 0.0 0.018873433135151486 0.0 0.0; 0.0 0.0 0.0 0.0 0.002879899158088243 0.0; 0.0 0.0 0.0 0.0 0.0 0.049787068367863944]
end


text:
using MOCNeutronTransport
using Plots
pyplot()
T = Float64
Nsoln = 79
Nset = [12, 27, 48, 79]
f = area

p₁ = Point_2D(T, 0)
p₂ = Point_2D(T, 2)
p₃ = Point_2D(T, 2, 2)
p₄ = Point_2D(T, 3//2, 1//4)
p₅ = Point_2D(T, 3, 1)
p₆ = Point_2D(T, 1, 1)
tri6 = Triangle6_2D((p₁, p₂, p₃, p₄, p₅, p₆))

soln = f(tri6; N=Nsoln)
err = [abs(soln - f(tri6; N = n)) for n = Nset]
plot(Nset[1:length(Nset)-1],
     err[1:length(Nset)-1],
     title = "Quadrature Points vs Triangle6_2D area error (Nsoln = $Nsoln)")
xlabel!("Quadrature Points")
ylabel!("area error")

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
"""
# AsyncLog.jl

`@asynclog` is `@async` with error logging.

Based on: https://github.com/JuliaLang/julia/issues/7626#issuecomment-756395909

e.g.

```
julia> @asynclog "my background loop" while true
           sleep(1)
           @assert false
       end

┌ Error: Error in: "my background loop"
│   exception =
│    AssertionError: false
│    Stacktrace:
│     ...
└ @ Main ~/git/AsyncLog/src/AsyncLog.jl:20
```

"""
module AsyncLog

export @asynclog, @asynclog_loop, @errorlog

"""
    @errorlog "label" expression...

Run `expression` in a `try/catch` block.
Log errors though `@error`.
Use `"label"` to identify error logs.
"""
macro errorlog(label, expr)
    quote
        try
            $(esc(expr))
        catch err
            exception=(err, catch_backtrace())
            label = $(esc(label))
            @error "Error in: \"$label\"" exception
            rethrow(err)
        end
    end
end

"""
    @asynclog "label" expression...

Run `@async expression` in a `try/catch` block.
Log errors though `@error`.
Use `"label"` to identify error logs.
"""
macro asynclog(label, expr)
    esc(quote
        @async @errorlog $label $expr
    end)
end

"""
    @asynclog_loop "label" expression...

Run `@async expression` in a `try/catch` block.
Log errors though `@error`.
Use `"label"` to identify error logs.
Restart `expression` after errors.
"""
macro asynclog_loop(label, expr)
    esc(quote
        @async while true
            try
                @errorlog $label $expr
            catch
            end
        end
    end)
end

# Documentation.

readme() = join([
    Docs.doc(@__MODULE__),
    Docs.@doc @asynclog
   ], "\n\n")

end # module


text:
"""
    module Hamiltonians

This module defines Hamiltonian types and functions for working with
Hamiltonians.

## Exported concrete Hamiltonian types include
- [`MatrixHamiltonian`](@ref)
- [`HubbardReal1D`](@ref)
- [`ExtendedHubbardReal1D`](@ref)
- [`HubbardMom1D`](@ref)
- [`BoseHubbardMom1D2C`](@ref)
- [`BoseHubbardReal1D2C`](@ref)

## Wrappers
- [`GutzwillerSampling`](@ref)
- [`GuidingVectorSampling`](@ref)

## Other
- [`G2Correlator`](@ref)

## Interface for working with Hamiltonians
- [`AbstractHamiltonian`](@ref): defined in the module [`Interfaces`](@ref)
"""
module Hamiltonians

using Parameters, StaticArrays, LinearAlgebra, SparseArrays
using FFTW
using Setfield

using ..StochasticStyles
using ..DictVectors
using ..BitStringAddresses
using ..ConsistentRNG

using ..Interfaces
import ..Interfaces: diagonal_element, num_offdiagonals, get_offdiagonal, starting_address,
    offdiagonals, random_offdiagonal, LOStructure

export AbstractHamiltonian, TwoComponentBosonicHamiltonian
export dimension, rayleigh_quotient, momentum
export BasisSetRep

export MatrixHamiltonian
export HubbardReal1D, HubbardMom1D, ExtendedHubbardReal1D, HubbardRealSpace
export HubbardReal1DEP, shift_lattice, shift_lattice_inv
export HubbardMom1DEP
export BoseHubbardMom1D2C, BoseHubbardReal1D2C
export GutzwillerSampling, GuidingVectorSampling
export Transcorrelated1D
export hubbard_dispersion, continuum_dispersion

export G2Correlator

export LatticeGeometry, PeriodicBoundaries, HardwallBoundaries, LadderBoundaries
export num_neighbours, neighbour_site

export sparse # from SparseArrays

include("abstract.jl")
include("offdiagonals.jl")
include("operations.jl")
include("geometry.jl")
include("excitations.jl")

include("MatrixHamiltonian.jl")

include("HubbardReal1D.jl")
include("HubbardReal1DEP.jl")
include("HubbardMom1D.jl")
include("HubbardMom1DEP.jl")
include("HubbardRealSpace.jl")
include("ExtendedHubbardReal1D.jl")

include("BoseHubbardReal1D2C.jl")
include("BoseHubbardMom1D2C.jl")

include("GutzwillerSampling.jl")
include("GuidingVectorSampling.jl")

include("Transcorrelated1D.jl")

include("TwoBodyCorrelation.jl")

end


text:
#=
Usage: V_i = interp2d_cubic(GCOORD,EL2NOD,els,lc,V,OPTS)

Purpose: Interpolates variables in a triangular finite element mesh using
         quasi-cubic interpolation functions. Based on code developed by
         Chao Shi & Jason Phipps Morgan, 2010.

Input
  GCOORD : [matrix]    : coordinates of all nodes in mesh
  EL2NOD : [matrix]    : finite element connectivity matrix (nnodel x nel)
  els    : [rowvector] : element in which each point is located
  lc     : [matrix]    : local coordinates in each element
  V      : [matrix]    : variables to be interpolated (nnod x nvar, i.e.
                         one variable field in each column)
  OPTS   : [structure] : options for this functions

Output
  V_i    : [matrix]    : interpolated values (npt x nvar, i.e. one
                         interpolated field in each column)
=#

import LinearAlgebra: det

struct QCubic
    x_el::SubArray{Float64, 2, Vector{Float64}, Tuple{Matrix{Int64}}, false}
    z_el::SubArray{Float64, 2, Vector{Float64}, Tuple{Matrix{Int64}}, false}
end

@inline function local_coordinates(vc, pc::Point2D{T}) where {T}
    # Coordinates of triangle's vertices
    x1,x2,x3 = vc[1].x, vc[2].x, vc[3].x
    z1,z2,z3 = vc[1].z, vc[2].z, vc[3].z
    # Coordinates of query point
    x, z = pc.x, pc.z

    denom = -x1*z3 + x1*z2 - x2*z1 + x2*z3 + x3*z1 - x3*z2

    λ2 = -( -x1*z + x1*z3 - x3*z1 + x*z1 +x3*z -x*z3 )/
          denom

    λ3 = (x1*z2 - x1*z - x2*z1 + x2*z + x*z1 -x*z2) /
          denom

    λ1 = 1 - (λ3+λ2)

    return @SVector [λ1, λ2, λ3]

end

@inline quasicubic_SF(λ) = @SVector [λ[1],
                                     λ[2],
                                     λ[3],
                                     λ[1]*λ[2],
                                     √2*λ[3]*λ[2],
                                     λ[3]*λ[1],
                                     λ[1]*λ[2]*(1-λ[3]-2*λ[2]),
                                     √2*λ[3]*λ[2]*(λ[3]-λ[2]),
                                     λ[3]*λ[1]*(1-λ[2]-2*λ[3])
                                     ]

@inline det(v::Vector{Point2D{T}}) where {T} =
    v[2].x*v[3].z+
    v[3].x*v[1].z+
    v[1].x*v[2].z-
    v[2].x*v[1].z-
    v[3].x*v[2].z-
    v[1].x*v[3].z

@inline element_center(v::Vector{Point2D{T}}) where {T} = ( (v[1].x+v[2].x+v[3].x)/3, (v[1].z+v[2].z+v[3].z)/3)

@inline function local_weights(vc)
    center = element_center(vc)
    @SVector [√( (vc[i].x - center[1])^2 + (vc[i].z - center[2])^2) for i in 1:3]
end

function local_derivatives(gr::Grid, vc, field, iel)
    # unpacks
    e2n = gr.e2n_p1
    θ, r = gr.x, gr.z
    # cache indices
    idx = @views e2n[:, iel]
    x_el, z_el = @views (θ[idx]),  @views (r[idx])

    vertices!(vc, x_el, z_el)

    ω = local_weights(vc) # weights

    invdetM = 1/det(vc)

    invMx = @SVector [(vc[2].z - vc[3].z)*invdetM,
                      (vc[3].z - vc[1].z)*invdetM,
                      (vc[1].z - vc[2].z)*invdetM]
    invMz = @SVector [(vc[3].x - vc[2].x)*invdetM,
                      (vc[1].x - vc[3].x)*invdetM,
                      (vc[2].x - vc[1].x)*invdetM]

    B = @views field[idx]
    dBdx = mydot(B, invMx)
    dBdz = mydot(B, invMz)

    return dBdx, dBdz, ω
end

function element_derivatives(gr, vc, fields)
    e2n_p1 = gr.e2n_p1
    nnod = maximum(e2n_p1)
    nels = size(e2n_p1,2)
    dBdx_nod = zeros(nnod)
    dBdz_nod = zeros(nnod)
    ω_nod = zeros(nnod)
    dBdx_el = Vector{Float64}(undef, nels)
    dBdz_el = similar(dBdx_el)

    @inbounds for iel in 1:nels
        dBdx_ipart, dBdz_ipart, ω = local_derivatives(gr, vc, fields, iel)

        dBdx_el[iel] = dBdx_ipart
        dBdz_el[iel] = dBdz_ipart

        # reduction over nodes
        @inbounds for i in 1:3
            el_nod = e2n_p1[i, iel]
            dBdx_nod[el_nod] += dBdx_ipart*ω[i]
            dBdz_nod[el_nod] += dBdz_ipart*ω[i]
            ω_nod[el_nod] += ω[i]
        end

    end

    dBdx, dBdz = similar(dBdx_nod), similar(dBdx_nod)
    @turbo for i in eachindex(dBdx)
        dBdx[i] = dBdx_nod[i]/ω_nod[i]
        dBdz[i] = dBdz_nod[i]/ω_nod[i]
    end

    dBdx_el2nod = view(dBdx, e2n_p1) .- dBdx_el' # remove background slope x
    dBdz_el2nod = view(dBdz, e2n_p1) .- dBdz_el' # remove background slope z

    return dBdx_el2nod, dBdz_el2nod
end

function slope(gr, dBdx_el2nod, dBdz_el2nod, nels)
    slope_info = Array{Float64,2}(undef, 6, nels)
    sqrthalf = √0.5
    @inbounds for i in 1:nels
        idx = @views gr.e2n_p1[:,i]
        x_el, z_el = @views (gr.x[idx]),  @views (gr.z[idx])
        c1x = -x_el[1] + x_el[2]
        c1z = -x_el[1] + x_el[3]
        c2x = -z_el[1] + z_el[2]
        c2z = -z_el[1] + z_el[3]

        slope_xi1 = c1x * dBdx_el2nod[1, i] + c2x.* dBdz_el2nod[1, i]
        slope_xi2 = c1x * dBdx_el2nod[2, i] + c2x.* dBdz_el2nod[2, i]
        slope_xi3 = c1x * dBdx_el2nod[3, i] + c2x.* dBdz_el2nod[3, i]
        slope_zi1 = c1z * dBdx_el2nod[1, i] + c2z.* dBdz_el2nod[1, i]
        slope_zi2 = c1z * dBdx_el2nod[2, i] + c2z.* dBdz_el2nod[2, i]
        slope_zi3 = c1z * dBdx_el2nod[3, i] + c2z.* dBdz_el2nod[3, i]

        dBdSC = sqrthalf*slope_xi3 - sqrthalf*slope_zi3 # dT/dS at point C (slope along side BC at C(0,1)):  dT/dS=dT/dxi*dxi/dS + dT/deta*deta/dS
        dBdSB = sqrthalf*slope_xi2 - sqrthalf*slope_zi2 # dT/dS at point B (slope along side BC at B(1,0))   the 'S+' is from C(0,1) to B(1,0)

        slope_info[1,i] = (slope_xi1-slope_xi2)*0.5; # this is to feed to N4
        slope_info[2,i] = (dBdSC-dBdSB)*0.5;         # N5
        slope_info[3,i] = (slope_zi1-slope_zi3)*0.5; # N6
        slope_info[4,i] = (slope_xi1+slope_xi2)*0.5; # N7
        slope_info[5,i] = (dBdSC+dBdSB)*0.5;         # N8
        slope_info[6,i] = (slope_zi1+slope_zi3)*0.5; # N9
    end
    slope_info
end

function cubic_particle_interpolation!(fields_particles, particle_info, slope_info, gr, fields, vc)
    for (i,info) in  enumerate(particle_info)
        iel = info.t
        pc = info.CCart
        x_el, z_el = @views (gr.x[gr.e2n_p1[:,iel]]),  @views (gr.z[gr.e2n_p1[:,iel]])
        els = view(gr.e2n_p1, :, iel)

        # compute shape functions
        vertices!(vc, x_el, z_el)
        λ = local_coordinates(vc, pc)
        N = quasicubic_SF(λ)

        # # interpolate field to particle
        Fel = view(fields, els)
        slope_info_part = view(slope_info, :, iel)

        fields_particles[i] = Fel[1]*N[1] + Fel[2]*N[2] + Fel[3]*N[3] +
                              slope_info_part[1] * N[4] +
                              slope_info_part[2] * N[5] +
                              slope_info_part[3] * N[6] +
                              slope_info_part[4] * N[7] +
                              slope_info_part[5] * N[8] +
                              slope_info_part[6] * N[9]
    end
end

function vertices!(v::Array{Point2D{T},1}, x, z) where {T}
    @inbounds for i = 1:3
        v[i].x = x[i]
        v[i].z = z[i]
    end
end

function quasicubic_interpolation(gr, particle_info, fields)

    np = length(particle_info)
    nels = size(gr.e2n_p1,2)

    vc = [Point2D{Cartesian}(0.0, 0.0) for _  in 1:3] # element vertices

    dBdx_el2nod, dBdz_el2nod = element_derivatives(gr, vc, fields)

    # also need to map the slope information to the xi, yi (ξ, η) axes
    slope_info = slope(gr, dBdx_el2nod, dBdz_el2nod, nels)

    fields_particles = Vector{Float64}(undef, np)
    cubic_particle_interpolation!(fields_particles, particle_info, slope_info, gr, fields, vc)
    fields_particles
end

function velocity_error(Upolar, particle_info, gr, output_path)
    Uth, Ur = analytic_flow(gr.θ, gr.r)
    # Allocate a couple of buffer arrays
    Uth_node = [u.x for u in Upolar]
    Ur_node = [u.z for u in Upolar]
    θp, rp = particle_coordinates(particle_info, coordinates_system = "polar") # initial particle position

    Uth_an, Ur_an = analytic_flow(θp, rp)
    speed_an = @. √(Uth_an^2 + Ur_an^2)

    Uth_particle = quasicubic_interpolation(gr, particle_info, Uth_node)
    Ur_particle = quasicubic_interpolation(gr, particle_info, Ur_node)
    speed_particle = @. √(Uth_particle^2 + Ur_particle^2)

    error_th = norm(@. Uth_an-Uth_particle)
    error_r = norm(@. Ur_an-Ur_particle)
    error_speed = norm(@. speed_an-speed_particle)

    fname = joinpath(output_path, "velocity_error.txt")
    open(fname, "a") do io
        writedlm(io, [error_th error_r error_speed])
    end

end

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
@testset "Halley" begin
    for n in [1,3,10]

A = rand(n, n)

r = polarfact(A, alg = :halley);

# Test unitary matrix U

U = r.U
H = r.H

@test U'*U ≈ Matrix(I,n,n) atol=1e-7

# Test Hermitian positive semifefinite matrix H

@test issymmetric(H)

for i in eigvals(H)
    @test i >= 0.
end

@test A ≈ U*H

    end
end

##########################################################
@testset "QDWH" begin
    for n in [1,3,10]

m = rand(1:10)

B = rand(m, m)

r2 = polarfact(B, alg = :qdwh);

# Test unitary matrix U

U2 = r2.U
H2 = r2.H

@test U2'*U2 ≈ Matrix(I,m,m) atol=1e-7

# Test Hermitian positive semifefinite matrix H

@test issymmetric(H2)

for i in eigvals(H2)
    @test i >= 0.
end

@test B ≈ U2*H2

    end
end


text:
# this user supplied file creates a dictionary of arguments
# if a key is repeated, the last use of the key is used
# it is a little bit dangerous letting the user run arbitrary code
# as part of the solver
# now that this file is read inside a function, it is better encapsulated

arg_dict = Dict{String, Any}(
"physics" => "Advection",
"run_type" => 1,
"jac_type" => 2,
"order" => 4,
"real_time" => false,
"t_max" => 1000.0,
"IC_name" => "ICp4",
"use_DG" => true,
"Flux_name" => "LFFlux",
"LFalpha" => 0,
"numBC" => 1,
"BC1" => [0, 1, 2, 3],
"BC1_name" => "p4BC",
"SRCname" => "SRCp4",
#"BC2" => [4, 10],
#"BC2_name" => "noPenetrationBC",
"delta_t" => 2*pi/5000,
#"smb_name" => "SRCMESHES/channel.smb",
#"smb_name" => "/users/creanj/test2/p1/m1/parallel2.smb",
#"dmg_name" => "parallel2.dmg",
"smb_name" => "SRCMESHES/serial2.smb",
"dmg_name" => "SRCMESHES/serial2.dmg",
#"smb_name" => "SRCMESHES/channel_verylarge.smb",
#"dmg_name" => "../../mesh_files/vortex.dmg",
"dmg_name" => ".null",
"res_abstol" => 1e-10,
# "step_tol" => 1e-6,
"itermax" => 100,
"res_abstol" => 1e-9,
#"res_reltol" => 1e-9,
#"step_tol" => 1e-10,
"writeq" => false,
"write_edge_vertnums" => false,
"write_face_vertnums" => false,
"write_qic" => false,
"writeboundary" => false,
"write_res" => false,
"print_cond" => false,
"write_counts" => false,
"write_vis" => true,
"output_freq" => 100,
"write_vis" => false,
"solve" => false,
"do_postproc" => true,
"exact_soln_func" => "ICp4",
)


text:
using jInv.Mesh
using jInv.LinearSolvers
using Test
using KrylovMethods
using Multigrid
using LinearAlgebra
using SparseArrays

println("===  Example 2D DivSigGrad ====");

domain = [0.0, 1.0, 0.0, 1.0];
n      = [50,50];
Mr     = getRegularMesh(domain,n)
G      = getNodalGradientMatrix(Mr);
Ar     = G'*G;
Ar     = Ar + 1e-2*norm(Ar,1)*sparse(1.0I,size(Ar,2),size(Ar,2));
N      = size(Ar,2);
B      = Ar*rand(N,4);

levels      = 5;
numCores 	= 8;
maxIter     = 15;
relativeTol = 1e-2;
relaxType   = "SPAI";
relaxParam  = 1.0;
relaxPre 	= 2;
relaxPost   = 2;
cycleType   ='V';
coarseSolveType = "Julia";

MG = getMGparam(Float64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType);

sMGPCG   = getMGsolver(MG,Mr, 1,"GMRES",out=1);
X,  = solveLinearSystem(Ar,B,sMGPCG);
@test norm(Ar*X-B)/norm(B) < sMGPCG.tol

sMGPCG   = getMGsolver(MG,Mr, 1,"PCG",out=1);
X,  = solveLinearSystem(Ar,B,sMGPCG);
@test norm(Ar*X-B)/norm(B) < sMGPCG.tol


text:
struct DynamicNUTS{T} <: Hamiltonian
    n_iters   ::  Integer   # number of samples
    space     ::  Set{T}    # sampling space, emtpy means all
    gid       ::  Integer   # group ID
end

"""
    DynamicNUTS(n_iters::Integer)

Dynamic No U-Turn Sampling algorithm provided by the DynamicHMC package.
To use it, make sure you have the DynamicHMC package installed.

```julia
# Import Turing and DynamicHMC.
using DynamicHMC, Turing

# Model definition.
@model gdemo(x, y) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0,sqrt(s))
  x ~ Normal(m, sqrt(s))
  y ~ Normal(m, sqrt(s))
  return s, m
end

# Pull 2,000 samples using DynamicNUTS.
chn = sample(gdemo(1.5, 2.0), DynamicNUTS(2000))
```
"""
function DynamicNUTS(n_iters::Integer, space...)
    _space = isa(space, Symbol) ? Set([space]) : Set(space)
    DynamicNUTS(n_iters, _space, 0)
end

function Sampler(alg::DynamicNUTS{T}) where T <: Hamiltonian
  return Sampler(alg, Dict{Symbol,Any}())
end

function sample(model::Model,
                alg::DynamicNUTS,
                chunk_size=CHUNKSIZE[]
                ) where T <: Hamiltonian

    if ADBACKEND[] == :forward_diff
        default_chunk_size = CHUNKSIZE[]  # record global chunk size
        setchunksize(chunk_size)        # set temp chunk size
    end

    spl = Sampler(alg)

    n = alg.n_iters
    samples = Array{Sample}(undef, n)
    weight = 1 / n
    for i = 1:n
        samples[i] = Sample(weight, Dict{Symbol, Any}())
    end

    vi = VarInfo()
    model(vi, HamiltonianRobustInit())

    if spl.alg.gid == 0
        link!(vi, spl)
        runmodel!(model, vi, spl)
    end

    function _lp(x)
        value, deriv = gradient(x, vi, model, spl)
        return ValueGradient(-value, -deriv)
    end

    chn_dynamic, _ = NUTS_init_tune_mcmc(FunctionLogDensity(length(vi[spl]), _lp), alg.n_iters)

    for i = 1:alg.n_iters
        vi[spl] = chn_dynamic[i].q
        samples[i].value = Sample(vi, spl).value
    end

    if ADBACKEND[] == :forward_diff
        setchunksize(default_chunk_size)      # revert global chunk size
    end

    return Chain(0, samples)
end


text:
using OmegaCore, Distributions

# A model where x --causes-> y
function m1()
  x = 1 ~ Normal(0, 1)
  y(ω) = (2 ~ Normal(x(ω), 1))(ω)
  (x = x, y = y)
end

# A model where y --causes-> x
function m2()
  x = 3 ~ Normal(0, 1)
  y(ω) = (4 ~ Normal(x(ω), 1))(ω)
  (x = y, y = x)
end

is_model_1 = 3 ~ Bernoulli(0.5)
model = ifelseₚ(is_model_1, m1(), m2())
mx(ω) = model(ω).x
mxω(ω) = mx(ω)(ω)
my(ω) = model(ω).y

x_ = 1.3
y_ = 2.0

# This won't do what we want!
condition1(ω) = (mxω |ᵈ (my => 100.0))(ω) >= 101.0

condition(ω) = (mxω |ᵈ (my =>ˡ y_))(ω) >= x_

# condition = (mxω |ᵈ (my =>ˡ y_)) >=ₚ x_

prob = mean(randsample(is_model_1 |ᶜ condition, 1000))
println("""The model1 is correct (x causes y), given after forcing y to 2 x,
we observed that x was 1.3 is""", prob)


text:
export compute_fft, compute_cc

@deprecate compute_fft() rfft()
@deprecate compute_cc() correlate()

"""

  compute_fft(R)

Computes windowed rfft of ambient noise data. Returns FFTData structure.

# Arguments
- `R::RawData`: RawData structure
"""
function compute_fft(R::RawData)
    FFT = rfft(R.x,1)
    return FFTData(R.name, R.id,R.loc, R.fs, R.gain, R.freqmin, R.freqmax,
                   R.cc_len, R.cc_step, R.whitened, R.time_norm, R.resp,
                   R.misc, R.notes, R.t, FFT)
end

"""
    compute_cc(FFT1::FFTData, FFT2::FFTData, maxlag::Float64;
               corr_type::String="cross-correlation")

Cross-correlate ambient noise data in the frequency domain.

Cross-correlation can be done using one of three options:

- Cross-correlation: ``C_{AB}(ω) = u_A(ω) u^∗_B(ω)``
- Coherence: ``C_{AB}(ω) = \frac{u_A(ω) u^∗_B(ω)}{∣ u_A(ω) ∣  ∣ u_B(ω) ∣}``
- Deconvolution: ``C_{AB}(ω) = \frac{u_A(ω) u^∗_B(ω)}{∣ u_B(ω) ∣^2}``

# Arguments
- `FFT1::FFTData`: FFTData object of fft'd ambient noise data.
- `FFT2::FFTData`: FFTData object of fft'd ambient noise data.
- `maxlag::Float64`: Maximum lag time (in seconds) in cross-correlation to save,
                     e.g. `maxlag = 20.` will save lag times = -20.:20. s.
- `corr_type::String`: Type of correlation: `cross-correlation`, `coherence` or
                       `deconv`.
"""
function compute_cc(FFT1::FFTData, FFT2::FFTData, maxlag::Float64;
                    corr_type::String="CC")

    comp = FFT1.name[end] * FFT2.name[end]
    # get intersect of dates; return nothing if no intersect
    inter = intersect(FFT1.t,FFT2.t)
    if length(inter) == 0
        throw(ArgumentError("No common windows for $(FFT1.name)-$(FFT2.name) $(FFT1.id)"))
    end

    N = convert(Int,round(FFT1.cc_len * FFT1.fs)) # number of data points
    ind1 = findall(x -> x ∈ inter, FFT1.t)
    ind2 = findall(x -> x ∈ inter, FFT2.t)
    corr = correlate(@views(FFT1.fft[:,ind1]), @views(FFT2.fft[:,ind2]), N,
                     convert(Int,round(maxlag * FFT1.fs)))
    rotated = false

    return CorrData(FFT1, FFT2, comp, rotated, corr_type,
                    maxlag, inter, corr)

end