text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -23.82106505746807 <= q <= 20.87497636142748)

                     add_NL_constraint(m, :((-0.7432570628691102 + 0.9028765283094948*(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))*(1 + erf((-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))/sqrt(2)))/2 + -0.9177729976456397*(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))*(1 + erf((-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))/sqrt(2)))/2 + -0.6091404444372248*(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))*(1 + erf((-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))/sqrt(2)))/2 + 0.7075296204914552*(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))*(1 + erf((-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))/sqrt(2)))/2 + -0.9115635329002854*(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))*(1 + erf((-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))/sqrt(2)))/2)*(1 + erf((-0.7432570628691102 + 0.9028765283094948*(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))*(1 + erf((-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))/sqrt(2)))/2 + -0.9177729976456397*(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))*(1 + erf((-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))/sqrt(2)))/2 + -0.6091404444372248*(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))*(1 + erf((-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))/sqrt(2)))/2 + 0.7075296204914552*(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))*(1 + erf((-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))/sqrt(2)))/2 + -0.9115635329002854*(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))*(1 + erf((-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))/sqrt(2)))/2)/sqrt(2)))/2 + (-0.2671327260424561 + -0.018499351907995898*(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))*(1 + erf((-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))/sqrt(2)))/2 + 0.12529725503678657*(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))*(1 + erf((-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))/sqrt(2)))/2 + 0.9839396030143157*(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))*(1 + erf((-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))/sqrt(2)))/2 + 0.28500560463058955*(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))*(1 + erf((-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))/sqrt(2)))/2 + 0.970334521625404*(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))*(1 + erf((-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))/sqrt(2)))/2)*(1 + erf((-0.2671327260424561 + -0.018499351907995898*(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))*(1 + erf((-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))/sqrt(2)))/2 + 0.12529725503678657*(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))*(1 + erf((-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))/sqrt(2)))/2 + 0.9839396030143157*(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))*(1 + erf((-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))/sqrt(2)))/2 + 0.28500560463058955*(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))*(1 + erf((-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))/sqrt(2)))/2 + 0.970334521625404*(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))*(1 + erf((-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))/sqrt(2)))/2)/sqrt(2)))/2 + (0.40040910594674495 + -0.7454849743358074*(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))*(1 + erf((-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))/sqrt(2)))/2 + 0.7108986625228169*(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))*(1 + erf((-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))/sqrt(2)))/2 + -0.30349624151173904*(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))*(1 + erf((-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))/sqrt(2)))/2 + -0.26057109207128715*(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))*(1 + erf((-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))/sqrt(2)))/2 + 0.7858005505148937*(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))*(1 + erf((-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))/sqrt(2)))/2)*(1 + erf((0.40040910594674495 + -0.7454849743358074*(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))*(1 + erf((-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))/sqrt(2)))/2 + 0.7108986625228169*(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))*(1 + erf((-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))/sqrt(2)))/2 + -0.30349624151173904*(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))*(1 + erf((-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))/sqrt(2)))/2 + -0.26057109207128715*(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))*(1 + erf((-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))/sqrt(2)))/2 + 0.7858005505148937*(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))*(1 + erf((-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))/sqrt(2)))/2)/sqrt(2)))/2 + (-0.7120753920031948 + -0.005433419997351141*(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))*(1 + erf((-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))/sqrt(2)))/2 + -0.8839514190404207*(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))*(1 + erf((-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))/sqrt(2)))/2 + -0.2933551592941672*(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))*(1 + erf((-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))/sqrt(2)))/2 + -0.5418193939730438*(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))*(1 + erf((-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))/sqrt(2)))/2 + 0.8978172795430157*(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))*(1 + erf((-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))/sqrt(2)))/2)*(1 + erf((-0.7120753920031948 + -0.005433419997351141*(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))*(1 + erf((-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))/sqrt(2)))/2 + -0.8839514190404207*(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))*(1 + erf((-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))/sqrt(2)))/2 + -0.2933551592941672*(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))*(1 + erf((-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))/sqrt(2)))/2 + -0.5418193939730438*(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))*(1 + erf((-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))/sqrt(2)))/2 + 0.8978172795430157*(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))*(1 + erf((-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))/sqrt(2)))/2)/sqrt(2)))/2 + (0.7950909223530527 + -0.6380481687930124*(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))*(1 + erf((-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))/sqrt(2)))/2 + -0.8467501205760102*(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))*(1 + erf((-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))/sqrt(2)))/2 + 0.5001198683847332*(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))*(1 + erf((-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))/sqrt(2)))/2 + -0.03855272462979853*(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))*(1 + erf((-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))/sqrt(2)))/2 + 0.3301524112226937*(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))*(1 + erf((-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))/sqrt(2)))/2)*(1 + erf((0.7950909223530527 + -0.6380481687930124*(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))*(1 + erf((-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))/sqrt(2)))/2 + -0.8467501205760102*(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))*(1 + erf((-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))/sqrt(2)))/2 + 0.5001198683847332*(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))*(1 + erf((-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))/sqrt(2)))/2 + -0.03855272462979853*(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))*(1 + erf((-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))/sqrt(2)))/2 + 0.3301524112226937*(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))*(1 + erf((-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))/sqrt(2)))/2)/sqrt(2)))/2 - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module HMCExamples
using Base: String
using ArgParse, CSV, DataFrames, DifferenceEquations, DelimitedFiles, DifferentiableStateSpaceModels, JLSO, JSON, LinearAlgebra, Logging, Random, Serialization, Turing, TuringCallbacks, TensorBoardLogger, Zygote
using Turing: @addlogprob!
import ArgParse.parse_item

# auto-generated code from /deps/generate_models.jl
include("generated_models/rbc.jl")
include("generated_models/FVGQ20.jl")

# # Utilities
include("utilities.jl")
include("experiment_results.jl")

# # Specific Models and CLI drivers
include("turing_models.jl")
include("estimate_rbc_1_kalman.jl")
include("estimate_rbc_1_joint.jl")
include("estimate_rbc_2_joint.jl")
include("estimate_FVGQ_1_kalman.jl")
include("estimate_FVGQ_1_joint.jl")
include("estimate_FVGQ_2_joint.jl")

end # module

ext: jl

max_stars_repo_name: HighDimensionalEconLab/HMCExamples

max_issues_repo_name: HighDimensionalEconLab/HMCExamples

max_forks_repo_name: HighDimensionalEconLab/HMCExamples


text:
using LinearAlgebra
using HOODESolver

function testcoeffloat()
    tab=[ (big"0.1", big"0.0001"), (big"0.00001", big"0.01")]

    for (bigeps, bigdt) in tab
        epsilon = convert(Float64,bigeps)
        dt = convert(Float64, bigeps)
        bigceab = CoefExpAB(10,bigeps,32,bigdt)
        ceab = CoefExpAB(10,epsilon,32,bigdt)
        print("epsilon=$epsilon norm=$(norm(bigceab.tab_coef-ceab.tab_coef, Inf))")
    end
 end
 testcoeffloat()


text:
module GreatValueMath

using Statistics
using DataFrames

# min-max feature scaling
function normalize(data::DataFrame, field::String)::DataFrame
    ma = maximum(data[!, field])
    mi = minimum(data[!, field])
    for (index, i) in enumerate(eachrow(data))
        data[index, field] = (i[field] - mi)/(ma - mi)
    end
    return data
end

function theil_sen_linear_regression(data::DataFrame, x::String, y::String)
    𝐧 = nrow(data)
    if 𝐧 > 1
      # independent variable vector
      xᵢ = data[!, x]
      # dependent variable vector
      yᵢ = data[!, y]
      # mean of dependent vector
      yₘ = (1/𝐧)*sum(yᵢ)
      # https://en.wikipedia.org/wiki/Theil%E2%80%93Sen_estimator
      m = median([(yᵢ[i+1]-yᵢ[i])/(xᵢ[i+1]-xᵢ[i]) for i ∈ 1:𝐧-1])
      𝑏 = median([(yᵢ[i]-m*xᵢ[i]) for i ∈ 1:𝐧])
      # https://en.wikipedia.org/wiki/Coefficient_of_determination
      # predicted vector
      fᵢ = [(m*x+𝑏) for x ∈ xᵢ]
      # residual vector
      𝖾ᵢ = [(yᵢ[i]-fᵢ[i]) for i ∈ 1:𝐧]
      SStot = sum([abs2(yᵢ[i]-yₘ) for i ∈ 1:𝐧])
      SSres = sum([abs2(𝖾ᵢ[i]) for i ∈ 1:𝐧])
      𝐑² = 1 - (SSres/SStot)
      return Dict("𝐑²" => 𝐑², "m" => m, "𝑏" => 𝑏)
    else
      return false
    end
end

export normalize
export theil_sen_linear_regression

end # module


text:
import StatsBase.sample, StatsBase.wsample
export init_kmeanspp!, init_svd!, init_nndsvd!
import NMF.nndsvd

# kmeans++ initialization, but with missing data
# we make sure never to look at "unobserved" entries in A
# so that models can be honestly cross validated, for example
function init_kmeanspp!(glrm::GLRM)
	m,n = size(glrm.A)
	k = glrm.k
	possible_centers = Set(1:m)
	glrm.Y = randn(k,n)
	# assign first center randomly
	i = sample(1:m)
	setdiff!(possible_centers, i)
	glrm.Y[1,glrm.observed_features[i]] = glrm.A[i,glrm.observed_features[i]]
	# assign next centers one by one
	for l=1:k-1
		min_dists_per_obs = zeros(m)
		for i in possible_centers
			d = zeros(l)
			for j in glrm.observed_features[i]
				for ll=1:l
					d[ll] += evaluate(glrm.losses[j], glrm.Y[ll,j], glrm.A[i,j])
				end
			end
			min_dists_per_obs[i] = minimum(d)/length(glrm.observed_features[i])
		end
        furthest_index = wsample(1:m,min_dists_per_obs)
		glrm.Y[l+1,glrm.observed_features[furthest_index]] = glrm.A[furthest_index,glrm.observed_features[furthest_index]]
	end
	return glrm
end

function init_svd!(glrm::GLRM; offset=true, scale=true, TOL = 1e-10)
    # only offset if the glrm model is offset
    offset = offset && typeof(glrm.rx) == lastentry1
    # only scale if we also offset
    scale = scale && offset
    m,n = size(glrm.A)
    k = glrm.k

    # find spans of loss functions (for multidimensional losses)
    yidxs = get_yidxs(glrm.losses)
    d = maximum(yidxs[end])

    # create a matrix representation of A with the same dimensions as X*Y
    # by expanding out all data types with embedding dimension greater than 1
    if all(map(length, yidxs) .== 1)
        Areal = glrm.A # save time, but in this case we'll still have a DataFrame
    else
        Areal = zeros(m, d)
        for f=1:n
            if length(yidxs[f]) == 1
                Areal[glrm.observed_examples[f], yidxs[f]] =
                    glrm.A[glrm.observed_examples[f], f]
            else
                if isa(glrm.losses[f].domain, CategoricalDomain)
                    levels = datalevels(glrm.losses[f])
										for e in glrm.observed_examples[f]
											for ilevel in 1:length(levels)
	                        Areal[e, yidxs[f][ilevel]] =
	                            (glrm.A[e, f] == levels[ilevel] ? 1 : -1)
	                    end
									  end
                elseif isa(glrm.losses[f].domain, OrdinalDomain)
                    embed_dim = embedding_dim(glrm.losses[f])
                    mymean = mean(glrm.A[glrm.observed_examples[f], f])
                    levels = datalevels(glrm.losses[f])
										for e in glrm.observed_examples[f]
											for ilevel in 1:(length(levels)-1)
	                        Areal[e, yidxs[f][ilevel]] =
	                            (glrm.A[e, f] > levels[ilevel] ? 1 : -1)
	                    end
									  end
                else
                    error("No default mapping to real valued matrix for domains of type $type(glrm.losses[f].domain)")
                end
            end
        end
    end

    # standardize A, respecting missing values
    means = zeros(d)
    stds  = zeros(d)
    Astd = zeros(m, d)
    for f in 1:n
        for j in yidxs[f]
            nomissing = Areal[glrm.observed_examples[f],j]
            means[j] = mean(nomissing)
            if isnan(means[j])
                means[j] = 1
            end
            stds[j] = std(nomissing)
            if stds[j] < TOL || isnan(stds[j])
                stds[j] = 1
            end
            Astd[glrm.observed_examples[f],j] = Areal[glrm.observed_examples[f],j] - means[j]
        end
    end
    if offset
        k -= 1
        glrm.X[end,:] = 1
        glrm.Y[end,:] = means
        if scale
            Astd = Astd*diagm(1./stds)
        end
        if k <= 0
            warn("Using an offset on a rank 1 model fits *only* the offset. To fit an offset + 1 low rank component, use k=2.")
            return glrm
        end
    end
    # options for rescaling:
    # 1) scale Astd so its mean is the same as the mean of the observations
    Astd *= m*n/sum(map(length, glrm.observed_features))
    # 2) scale columns inversely proportional to number of entries in them & so that column mean is same as mean of observations in it
    # intuition: noise in a dense column is low rank, so downweight dense columns
    # Astd *= diagm(m./map(length, glrm.observed_examples))
    # 3) scale columns proportional to scale of regularizer & so that column mean is same as mean of observations in it
    # Astd *= diagm(m./map(scale, glrm.ry))
		# ASVD = rsvd(Astd, k) - slower than built-in svds, and fails for sparse matrices
		ASVD = svds(Astd, nsv = k)[1]
    # initialize with the top k components of the SVD,
    # rescaling by the variances
    @assert(size(glrm.X, 1) >= k)
    @assert(size(glrm.X, 2) >= m)
    @assert(size(glrm.Y, 1) >= k)
    @assert(size(glrm.Y, 2) >= d)
    glrm.X[1:k,1:m] = Diagonal(sqrt.(ASVD[:S]))*ASVD[:U]' # recall X is transposed as per column major order.
    glrm.Y[1:k,1:d] = Diagonal(sqrt.(ASVD[:S]))*ASVD[:Vt]*Diagonal(stds)
    return glrm
end

function init_nndsvd!(glrm::GLRM; scale::Bool=true, zeroh::Bool=false,
                      variant::Symbol=:std, max_iters::Int=0)
    # NNDSVD initialization:
    #    Boutsidis C, Gallopoulos E (2007). SVD based initialization: A head
    #    start for nonnegative matrix factorization. Pattern Recognition
    m,n = size(glrm.A)

    # only initialize based on observed entries
    A_init = zeros(m,n)
    for i = 1:n
        A_init[glrm.observed_examples[i],i] = glrm.A[glrm.observed_examples[i],i]
    end

    # scale all columns by the Loss.scale parameter
    if scale
        for i = 1:n
            A_init[:,i] .*= glrm.losses[i].scale
        end
    end

    # run the first nndsvd initialization
    W,H = nndsvd(A_init, glrm.k, zeroh=zeroh, variant=variant)
    glrm.X = W'
    glrm.Y = H

    # If max_iters>0 do a soft impute for the missing entries of A.
    #   Iterate: Estimate missing entries of A with W*H
    #            Update (W,H) nndsvd estimate based on new A
    for iter = 1:max_iters
        # Update missing entries of A_init
        for j = 1:n
            for i = setdiff(1:m,glrm.observed_examples[j])
                A_init[i,j] = dot(glrm.X[:,i],glrm.Y[:,j])
            end
        end
        # Re-estimate W and H
        W,H = nndsvd(A_init, glrm.k, zeroh=zeroh, variant=variant)
        glrm.X = W'
        glrm.Y = H
    end
end


text:
# global cache
# ============

# XXX `JET_REPORT_CACHE` isn't synced with `JET_CODE_CACHE`, and this may lead to a problem ?

"""
    JET_REPORT_CACHE::$(typeof(JET_REPORT_CACHE))

Keeps JET report cache for a `MethodInstance`.
Reports are cached when `AbstractAnalyzer` exits from `_typeinf`.
"""
const JET_REPORT_CACHE = IdDict{UInt, IdDict{MethodInstance,Vector{InferenceErrorReportCache}}}()

"""
    JET_CODE_CACHE::$(typeof(JET_CODE_CACHE))

Keeps `CodeInstance` cache associated with `mi::MethodInstace` that represent the result of
  an inference on `mi` performed by `AbstractAnalyzer`.
This cache is completely separated from the `NativeInterpreter`'s global cache, so that
  JET analysis never interacts with actual code execution.
"""
const JET_CODE_CACHE = IdDict{UInt, IdDict{MethodInstance,CodeInstance}}()

# just used for interactive developments
__clear_caches!() = (empty!(JET_REPORT_CACHE); empty!(JET_CODE_CACHE))

function CC.code_cache(analyzer::AbstractAnalyzer)
    cache  = JETGlobalCache(analyzer)
    worlds = WorldRange(get_world_counter(analyzer))
    return WorldView(cache, worlds)
end

struct JETGlobalCache
    analyzer::AbstractAnalyzer
end

# cache existence for this `analyzer` is ensured on its construction
jet_report_cache(analyzer::AbstractAnalyzer)     = JET_REPORT_CACHE[get_cache_key(analyzer)]
jet_report_cache(wvc::WorldView{JETGlobalCache}) = jet_report_cache(wvc.cache.analyzer)
jet_code_cache(analyzer::AbstractAnalyzer)       = JET_CODE_CACHE[get_cache_key(analyzer)]
jet_code_cache(wvc::WorldView{JETGlobalCache})   = jet_code_cache(wvc.cache.analyzer)

CC.haskey(wvc::WorldView{JETGlobalCache}, mi::MethodInstance) = haskey(jet_code_cache(wvc), mi)

function CC.typeinf_edge(analyzer::AbstractAnalyzer, method::Method, @nospecialize(atypes), sparams::SimpleVector, caller::InferenceState)
    # NOTE enable the report cache restoration at `code = get(code_cache(interp), mi, nothing)`
    set_cache_enabled!(analyzer, true)
    return @invoke typeinf_edge(analyzer::AbstractInterpreter, method::Method, @nospecialize(atypes), sparams::SimpleVector, caller::InferenceState)
end

function CC.get(wvc::WorldView{JETGlobalCache}, mi::MethodInstance, default)
    ret = get(jet_code_cache(wvc), mi, default) # will ignore native code cache for a `MethodInstance` that is not analyzed by JET yet

    analyzer = wvc.cache.analyzer

    # XXX this relies on a very dirty analyzer state manipulation, the reason for this is
    # that this method (and `code_cache(::AbstractAnalyzer)`) can be called from multiple
    # contexts including edge inference, constant prop' heuristics and inlining, where we
    # want to use report cache only in edge inference, but we can't tell which context is
    # the caller of this specific method call here and thus can't tell whether we should
    # enable report cache reconstruction without the information
    if get_cache_enabled(analyzer)
        if isa(ret, CodeInstance)
            # cache hit, now we need to append cached reports associated with this `MethodInstance`
            global_cache = get(jet_report_cache(wvc), mi, nothing)
            if isa(global_cache, Vector{InferenceErrorReportCache})
                for cached in global_cache
                    restored = restore_cached_report!(cached, analyzer)
                    push!(get_to_be_updated(analyzer), restored) # should be updated in `abstract_call` (after exiting `typeinf_edge`)
                    # TODO make this holds when the `analyzer` hooks into `finish` or `optimize`
                    # more generally, handle cycles correctly
                    @static JET_DEV_MODE && if isa(analyzer, JETAnalyzer)
                        actual, expected = first(restored.vst).linfo, mi
                        @assert actual === expected "invalid global cache restoration, expected $expected but got $actual"
                    end
                end
            end
        end
    end

    set_cache_enabled!(analyzer, false)
    return ret
end

function CC.getindex(wvc::WorldView{JETGlobalCache}, mi::MethodInstance)
    r = CC.get(wvc, mi, nothing)
    r === nothing && throw(KeyError(mi))
    return r::CodeInstance
end

function CC.setindex!(wvc::WorldView{JETGlobalCache}, ci::CodeInstance, mi::MethodInstance)
    setindex!(jet_code_cache(wvc), ci, mi)
    add_jet_callback!(mi) # register the callback on invalidation
    return nothing
end

function add_jet_callback!(linfo)
    @static BACKEDGE_CALLBACK_ENABLED || return nothing

    if !isdefined(linfo, :callbacks)
        linfo.callbacks = Any[invalidate_jet_cache!]
    else
        if !any(@nospecialize(cb)->cb===invalidate_jet_cache!, linfo.callbacks)
            push!(linfo.callbacks, invalidate_jet_cache!)
        end
    end
    return nothing
end

function invalidate_jet_cache!(replaced, max_world, depth = 0)
    for cache in values(JET_REPORT_CACHE); delete!(cache, replaced); end
    for cache in values(JET_CODE_CACHE); delete!(cache, replaced); end

    if isdefined(replaced, :backedges)
        for mi in replaced.backedges
            mi = mi::MethodInstance
            if !(any(cache->haskey(cache, mi), values(JET_REPORT_CACHE)) ||
                 any(cache->haskey(cache, mi), values(JET_CODE_CACHE)))
                continue # otherwise fall into infinite loop
            end
            invalidate_jet_cache!(mi, max_world, depth+1)
        end
    end
    return nothing
end

# local cache
# ===========

struct JETLocalCache
    analyzer::AbstractAnalyzer
    cache::Vector{InferenceResult}
end

CC.get_inference_cache(analyzer::AbstractAnalyzer) = JETLocalCache(analyzer, get_inference_cache(get_native(analyzer)))

function CC.cache_lookup(linfo::MethodInstance, given_argtypes::Vector{Any}, cache::JETLocalCache)
    inf_result = cache_lookup(linfo, given_argtypes, cache.cache)

    isa(inf_result, InferenceResult) || return inf_result

    # constant prop' hits a cycle (recur into same non-constant analysis), we just bail out
    isa(inf_result.result, InferenceState) && return inf_result

    # cache hit, try to restore local report caches
    analyzer = cache.analyzer
    sv = get_current_frame(analyzer)::InferenceState

    analysis_result = jet_cache_lookup(linfo, given_argtypes, get_cache(analyzer))

    # corresponds to report throw away logic in `_typeinf(analyzer::AbstractAnalyzer, frame::InferenceState)`
    filter!(!is_from_same_frame(sv.linfo, linfo), get_reports(analyzer))

    if isa(analysis_result, AnalysisResult)
        for cached in analysis_result.cache
            restored = restore_cached_report!(cached, analyzer)
            push!(get_to_be_updated(analyzer), restored) # should be updated in `abstract_call_method_with_const_args`
            # TODO make this holds when the `analyzer` hooks into `finish` or `optimize`
            # more generally, handle cycles correctly
            @static JET_DEV_MODE && if isa(analyzer, JETAnalyzer)
                actual, expected = first(restored.vst).linfo, linfo
                @assert actual === expected "invalid local cache restoration, expected $expected but got $actual"
            end
        end
    end

    return inf_result
end

# should sync with the implementation of `cache_lookup(linfo::MethodInstance, given_argtypes::Vector{Any}, cache::Vector{InferenceResult})`
function jet_cache_lookup(linfo::MethodInstance, given_argtypes::Vector{Any}, cache::Vector{AnalysisResult})
    method = linfo.def::Method
    nargs::Int = method.nargs
    method.isva && (nargs -= 1)
    length(given_argtypes) >= nargs || return nothing
    for cached_result in cache
        cached_result.linfo === linfo || continue
        cache_match = true
        cache_argtypes = cached_result.argtypes
        cache_overridden_by_const = cached_result.overridden_by_const
        for i in 1:nargs
            if !is_argtype_match(given_argtypes[i],
                                 cache_argtypes[i],
                                 CC.getindex(cache_overridden_by_const, i))
                cache_match = false
                break
            end
        end
        if method.isva && cache_match
            cache_match = is_argtype_match(tuple_tfunc(given_argtypes[(nargs + 1):end]),
                                           cache_argtypes[end],
                                           CC.getindex(cache_overridden_by_const, CC.length(cache_overridden_by_const)))
        end
        cache_match || continue
        return cached_result
    end
    return nothing
end

CC.push!(cache::JETLocalCache, inf_result::InferenceResult) = CC.push!(cache.cache, inf_result)


text:
function cp_file(file, srcDir, dstDir)
  srcFile = joinpath(srcDir, file)
  dstFile = joinpath(dstDir, file)
  !isfile(dstFile) && cp(srcFile, dstFile; remove_destination = true)
end

function create_cfg()
  srcDir = dirname(@__FILE__)
  cfgDir = joinpath(srcDir, "cfg")
  !isdir(cfgDir) && mkdir(cfgDir)
  cp_file("appinfo.json", srcDir, cfgDir)
  cp_file("config.json", srcDir, cfgDir)
end

create_cfg()


text:
# calculate probability by pdf
# covariance is plused 20
# plot as 2d-contour

module ContourPdfPlus20
    using DataFrames, CSV, Statistics, LinearAlgebra
    using Base.Iterators, Plots
    pyplot()

    function pdf(x, mu, cov_mat)
        # inverse matrix
        inv_mat = inv(cov_mat)

        # determinant
        det_mat = det(cov_mat)

        diff_x = x - mu
        return exp(-((diff_x'/cov_mat)*diff_x)/2)/(2*pi*sqrt(det_mat))
    end

    function main(is_test=false)
        # input
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_700.txt")
        df_org = CSV.read(data_path, DataFrame,
                      header=["date", "time", "ir", "lidar"],
                      delim=' ')

        # extract between 12 and 16 o'clock
        df_ext = df_org[(df_org.time.>=120000).&(df_org.time.<160000),:]

        # variance
        var_ir = var(df_ext.ir)
        var_lidar = var(df_ext.lidar)

        #covariance
        diff_ir = df_ext.ir .- mean(df_ext.ir)
        diff_lidar = df_ext.lidar .- mean(df_ext.lidar)
        a = diff_ir .* diff_lidar
        covar = cov(df_ext.ir, df_ext.lidar, corrected=false)

        # covariance matrix
        cov_mat = [var_ir covar; covar var_lidar]
        cov_mat[1, 2] += 20
        cov_mat[2, 1] += 20

        # mean(mu)
        mu = [mean(df_ext.ir); mean(df_ext.lidar)]

        # plot contour
        vx = 0:40
        vy = 710:750
        z = [pdf([x; y], mu, cov_mat) for x in vx, y in vy]
        contour(vx, vy, z', label="contour", c=:haline, aspect_ratio=:equal)

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/multi_dim_gauss_dist/contour_pdf/contour_pdf_plus20.png")
            savefig(save_path)
        end
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
include("equilibrium.jl")
include("calibrate_params.jl")
using CalibrateParameters
using Base.Test
using Logging

function init_parameters()
    parameters = Dict{Symbol, Any}()
    N, J, T = 2, 3, 1
    parameters[:N], parameters[:J], parameters[:T] = N, J, T
    # CES parameters
	parameters[:sigma] = 0.999
	parameters[:theta] = 4.0
	parameters[:eta] = 4.0

	########## parameters common across scenarios
	## these are function of data
	# inverse of adjustment cost, 0 if cannot readjust
	parameters[:one_over_rho] = 0.1
    parameters[:S] = 100

    parameters[:numerical_zero] = 1e-12

    parameters[:bp_weights] = [0.774074394803123; -0.201004684236153; -0.135080548288772; -0.0509519648766360]

    # adaptive step size. large lambda means large steps
    parameters[:inner_step_size] = exp(-0.10*(parameters[:J]-1)^0.75)
    # large substitution needs more dampening
    parameters[:middle_step_size] = exp(-0.275*max(1,parameters[:sigma]))
    parameters[:adjustment_step_size] = 0.25
    # any deviation from sigma=1 needs more dampening
    parameters[:outer_step_size] = exp(-0.5*abs(log(parameters[:sigma])))
    # this is log points of average input price differences
    parameters[:inner_tolerance]  = 0.001
    parameters[:middle_tolerance] = 0.001
    parameters[:adjustment_tolerance] = 0.0005
    parameters[:outer_tolerance] = 0.001

    # maximum number of iterations in each loop
    parameters[:max_iter_inner] = 1000
    parameters[:max_iter_middle] = 50
    parameters[:max_iter_adjustment] = 50
    parameters[:max_iter_outer] = 50
    return parameters
end

function test_parameters(parameters)
    function positive(args...)
        for item in args
            @test minimum(item) > 0.0
        end
    end
    @testset "CalibrateParameters tests" begin
    N, J, T = parameters[:N], parameters[:J], parameters[:T]
    @test size(parameters[:beta_j]) == (1,1,J,1)
    @test size(parameters[:gamma_jk]) == (J,J)

    C = diagm(parameters[:beta_j][1,1,:,1])+parameters[:gamma_jk]'-eye(J)
    @test sum(C, 2) ≈ zeros(J,1) atol=1e-9

    @test size(parameters[:S_nt_data]) == (1,N,1,T)
    @test sum(parameters[:S_nt_data], 2) ≈ zeros(1,1,1,T) atol=1e-9

    @test size(parameters[:d]) == (N,N,J,T)
    @test sum(parameters[:d], 2)[:,:,1:end-1,:] ≈ ones(N,1,J-1,T) atol=1e-9

    @test size(parameters[:kappa_mnjt]) == (N,N,J,T)

    @test size(parameters[:final_expenditure_shares]) == (1,N,J,T)
    @test sum(parameters[:final_expenditure_shares], 3) ≈ ones(1,N,1,T) atol=1e-9

    @test size(parameters[:nu_njt]) == (1,1,J,T)
    @test sum(parameters[:nu_njt], 3) ≈ ones(1,1,1,T) atol=1e-9
    @test size(parameters[:p_sectoral]) == (1,N,J,T)
    @test any(isnan.(parameters[:p_sectoral])) == false
    @test parameters[:p_sectoral][1,end,:,1] ≈ ones(J) atol=1e-9

    @test size(parameters[:w_njt]) == (1,N,J,T)

    @test size(parameters[:B_j]) == (1,1,J,1)

    @test typeof(parameters[:xi]) == Float64

    @test size(parameters[:A]) == (1, N, J, T)
    # productivity does not change more than 50% per year
    logA = log.(parameters[:A])
    @test maximum(abs.(logA[:,:,:,T] .- logA[:,:,:,1])) < log(1.50)*T

    @test size(parameters[:nominal_world_expenditure]) == (1,1,1,T)

    @test size(parameters[:A_njs]) == (T,)
    @test size(parameters[:A_njs][1]) == (1,N,J,1)
    @test size(parameters[:A_njs][2]) == (1,N,J,parameters[:S])

    for t=1:T
        @test parameters[:A_njs][t][:,:,:,1] ≈ parameters[:A][:,:,:,t]
    end

    positive(parameters[:final_expenditure_shares], parameters[:nu_njt], parameters[:A_njs][2], parameters[:p_sectoral], parameters[:A])
    end
end

srand(7094)

parameters = init_parameters()
CalibrateParameters.calibrate_parameters!(parameters, "data/impvol_data.jld2")
test_parameters(parameters)


text:
include("LinearModel.jl")
include("BoltzmannModel.jl")
include("PowerLawModel.jl")
include("VacuumModel.jl")

show(io::IO, tm::AbstractTemperatureModel) = print(io, tm)
show(io::IO,::MIME"text/plain", tm::AbstractTemperatureModel) = show(io, tm)


text:
function generateMassMatrix(FEMmesh)
    M = Mass2DMatrix(FEMmesh)
    #Mlu = lu(M)
    return M
end

function generateMassMatrixWScalar(FEMmesh,f)
    M = Mass2DMatrix(FEMmesh,f)
    #Mlu = lu(M)
    return M
end


text:
using RRTMGP
using DocStringExtensions
using RRTMGP.OpticalProps
using RRTMGP.Utilities
using RRTMGP.Gases
using RRTMGP.GasOptics
using RRTMGP.GasConcentrations
using RRTMGP.RTESolver
using RRTMGP.Fluxes
using RRTMGP.SourceFunctions
using RRTMGP.CloudOptics

import RRTMGP.OpticalProps: validate!, subset_range!, delta_scale!, get_nmom

include(joinpath("..", "ReadInputData", "ReadInputs.jl"))
include(joinpath("..", "ReadInputData", "LoadCoefficients.jl"))
include(joinpath("..", "ReadInputData", "LoadCloudCoefficients.jl"))

include("OpticalPropsNStream.jl")
include("OpticalPropsNStream_kernels.jl")
include("ExtractSubset.jl")

max_stars_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
struct Numeric end

# scalar ops

@pure Numeric !, !=, &, *, +, -, /, \, <, <=, >, >=, ^, |, abs, abs2, acos, acosd,
              acosh, acot, acotd, acoth, acsc, acscd, acsch, asec, asecd, asech,
              asin, asind, asinh, atan, atand, atanh, cbrt, conj, cos, cosd, cosh,
              cospi, cot, cotd, coth, csc, cscd, csch, deg2rad, exp, exp10, exp2,
              expm1, float, inv, log, log10, log1p, log2, rad2deg, sec, secd, sech,
              sin, sind, sinh, sinpi, sqrt, tan, tand, tanh, transpose, trailing_zeros,
              >>, <<, unsigned, rem, //

@abstract Numeric rand() = Float64
@abstract Numeric randn() = Float64
@abstract Numeric rand(::Type{Bool}) where T = Bool

@abstract Numeric sum(xs::Array{T,N}; dims = :) where {T,N} =
  dims == Const(:) ? T : Array{T,N}


text:
using Dierckx

# Structure Refinement
function Refine(data::Matrix{T}, tmatch::Matrix{T}, ensize::Int64, qstart::T; population::Int64=10, maxiter::Int64=100000, reports::Int64=100) where T<:Real

    # Initialize
    local NGENE = size(tmatch, 2) - 1;
    local GENES = collect(1:NGENE);
    chromosomes = reshape(sample(GENES, ensize*population, replace=true), ensize, population);

    # Pre-process the data according to qstart
    id = findall(data[:,1] .> qstart)[1];
    expdata = data[id:end, :];

    # filter out Problematic data
    expdata = expdata[(expdata[:, 3] .!= 0), :]; # error = 0.0
    expdata[:, 2] = abs.(expdata[:, 2]);

    t = zeros(size(expdata,1), size(tmatch,2));
    [t[:,i] = Spline1D(tmatch[:, 1], tmatch[:, i]; k=1)(expdata[:,1]) for i = 1:size(tmatch, 2)];
    t[:,1] = expdata[:,1];

    # Fitness Value
    function fitness(data::Matrix{T}, t::Matrix{T}, chromosomes::Matrix{Int64}, population::Int64) where T<:Real
        ft = zeros(population);

        function logfit(iexp::Matrix{T}, ifit::Vector{T}) where T<:Real

            #σh = iexp[:,2] + iexp[:,3];
            #σl = iexp[:,2] - iexp[:,3];

            σ = iexp[:,3]./(iexp[:,2] * log(10));
            #σ[index_pos] = 0.5 * (log10.(σh[index_pos]) - log10.(σl[index_pos]));
            return mean(((log10.(iexp[:,2]) - log10.(ifit))./σ).^2);
        end

        [ft[i] = logfit(expdata, mean(t[:, chromosomes[:,i] .+ 1], dims=2)[:,1]) for i = 1:population];
        return ft;
    end

    # Selection
    function selection(chromosomes::Matrix{Int64}, ft::Vector{T}, population::Int64) where T<:Real
        m = Int64(floor(population/3));
        sele = sample(collect(1:population), ProbabilityWeights(exp.(-3.0*ft)), m, replace=false);
        parents = chromosomes[:, sele];
        return parents;
    end

    # Crossover
    function crossover(parents::Matrix{Int64}, ensize::Int64, population::Int64)
        child = Matrix{Int64}(undef, ensize, population);
        nparents = size(parents, 2);
        for i = 1:population
            p = sample(collect(1:nparents), 2);
            father, mother = parents[:, p[1]], parents[:, p[2]];
            sameindex = findall(indexin(father, mother) .!== nothing);
            m = length(sameindex);
            if m == 0
                child[:, i] = sample([father; mother], ensize, replace=false);
            else
                child[1:m, i] = father[sameindex];
                child[m+1:end, i] = sample(unique([father; mother]), ensize-m);
            end
        end
        return child;
    end

    # Mutation
    function mutation(child::Matrix{Int64}, ensize::Int64, population::Int64; rate::Float64=0.3)
        mutated_child = copy(child);
        for i = 1:population
            if rand() <= rate
                nmutation = sample(collect(1:ensize), ProbabilityWeights(1.0 ./collect(1:ensize)));
                mutated_child[sample(collect(1:ensize), nmutation, replace=false), i] = sample(GENES, nmutation, replace=false);
            end
        end
        return mutated_child;
    end

    # Main process
    best_fitness = Inf;
    best_chromosome = Vector{Int64}(undef, size(chromosomes, 1));

    tick = Int64(round(maxiter/reports));

    for i = 1:maxiter
        ft = fitness(expdata, t, chromosomes, population);
        parents = selection(chromosomes, ft, population);
        child = crossover(parents, ensize, population);
        mutated_child = mutation(child, ensize, population);

        m = minimum(ft);

        if m <= best_fitness
            best_fitness = m;
            best_chromosome = chromosomes[:, findall(ft.==m)];
        end

        chromosomes = mutated_child;
        i % tick == 0 ? println("Generation #$i: fitness = $(round(m, digits=5)), best fit so far = $(round(best_fitness, digits=5)).") : nothing;
    end
    println("Finished Genetic Algorithm, with best fitness = $(round(best_fitness, digits=5)).");
    return best_fitness, best_chromosome;
end


text:
"""
    struct SolidStateDetector{T,SC,CT,PT,VDM} <: AbstractConfig{T}

Struct to describe all parts of a solid state detector, i.e.
the [`Semiconductor`](@ref), a set of [`Contact`](@ref) and (optionally) [`Passive`](@ref) and virtual drift volumes.

The properties of the parts (charge densities, fixed potentials, relative permittivity
of the materials) will be used as input to the calculation of [`ElectricPotential`](@ref) and
[`WeightingPotential`](@ref) in the [`Simulation`](@ref).

## Parametric types
* `T`: Precision type.
* `SC`: Type of the `semiconductor`.
* `CT`: Type of the `contacts`.
* `PT`: Type of the `passives`.
* `VDM`: Type of the `virtual_drift_volumes`.

## Fields
* `name::String`: Name of the detector.
* `semiconductor::SC`: [`Semiconductor`](@ref) of the detector.
* `contacts::CT`: Vector of [`Contact`](@ref) of the detector.
* `passives::PT`: Vector of [`Passive`](@ref) objects, e.g. holding structures around the detector.
* `virtual_drift_volumes::VDM`: Vector of virtual drift volumes in which the drift can be modulated
    by user-defined methods for `modulate_driftvector`, e.g. [`DeadVolume`](@ref).

See also [`Semiconductor`](@ref), [`Contact`](@ref), [`Passive`](@ref) and [`DeadVolume`](@ref).
"""
struct SolidStateDetector{T,SC,CT,PT,VDM} <: AbstractConfig{T}
    name::String  # optional
    semiconductor::SC
    contacts::CT
    passives::PT
    virtual_drift_volumes::VDM

    SolidStateDetector{T}(n::AbstractString,s::SC,c::C,p::P,v::VDM) where {T,SC,C,P,VDM}= new{T,SC,C,P,VDM}(n,s,c,p,v)
end

function SolidStateDetector(det::SolidStateDetector{T,SC,CT,PT,VDM}, impurity_density::AbstractImpurityDensity{T}) where {T,SC,CT,PT,VDM}
    sc = Semiconductor(det.semiconductor, impurity_density)
    SolidStateDetector{T}(
        det.name, sc, det.contacts, det.passives, det.virtual_drift_volumes
    )
end
function SolidStateDetector(det::SolidStateDetector{T,SC,CT,PT,VDM}, chargedriftmodel::AbstractChargeDriftModel{T}) where {T,SC,CT,PT,VDM}
    sc = Semiconductor(det.semiconductor, chargedriftmodel)
    SolidStateDetector{T}(
        det.name, sc, det.contacts, det.passives, det.virtual_drift_volumes
    )
end
function SolidStateDetector(det::SolidStateDetector{T,SC,CT,PT,VDM}; contact_id::Int, contact_potential::Real) where {T,SC,CT,PT,VDM}
    oc = det.contacts[contact_id]
    nc = Contact(T(contact_potential), oc.material, oc.id, oc.name, oc.geometry )
    contacts = [c.id == contact_id ? nc : c for c in det.contacts]
    SolidStateDetector{T}( det.name, det.semiconductor, contacts, det.passives, det.virtual_drift_volumes )
end

get_precision_type(::SolidStateDetector{T}) where {T} = T

function get_world_limits_from_objects(::Type{Cylindrical}, det::SolidStateDetector{T}) where {T <: SSDFloat}
    ax1l::T, ax1r::T, ax2l::T, ax2r::T, ax3l::T, ax3r::T = 0, 1, 0, 1, 0, 1
    t::Array{T, 2} = hcat(vcat(Vector{CylindricalPoint{T}}.(ConstructiveSolidGeometry.extreme_points.(ConstructiveSolidGeometry.surfaces(det.semiconductor.geometry)))...)...)
    (ax1l, ax1r), (ax3l, ax3r) = broadcast(i -> extrema(t[i,:]), 1:2:3)
    for c in det.contacts
        t = hcat([ax1l ax1r; ax2l ax2r; ax3l ax3r], hcat.(Vector{CylindricalPoint{T}}.(ConstructiveSolidGeometry.extreme_points.(ConstructiveSolidGeometry.surfaces(c.geometry)))...)...)
        (ax1l, ax1r), (ax3l, ax3r) = broadcast(i -> extrema(t[i,:]), 1:2:3)
    end
    if !ismissing(det.passives)
        for p in det.passives
            t = hcat([ax1l ax1r; ax2l ax2r; ax3l ax3r], hcat.(Vector{CylindricalPoint{T}}.(ConstructiveSolidGeometry.extreme_points.(ConstructiveSolidGeometry.surfaces(p.geometry)))...)...)
            (ax1l, ax1r), (ax3l, ax3r) = broadcast(i -> extrema(t[i,:]), 1:2:3)
        end
    end
    return ax1l, ax1r, ax2l, ax2r, ax3l, ax3r
end

function get_world_limits_from_objects(::Type{Cartesian}, det::SolidStateDetector{T}) where {T <: SSDFloat}
    ax1l::T, ax1r::T, ax2l::T, ax2r::T, ax3l::T, ax3r::T = 0, 1, 0, 1, 0, 1
    t::Array{T, 2} = hcat(vcat(ConstructiveSolidGeometry.extreme_points.(ConstructiveSolidGeometry.surfaces(det.semiconductor.geometry))...)...)
    (ax1l, ax1r), (ax2l, ax2r), (ax3l, ax3r) = broadcast(i -> extrema(t[i,:]), 1:3)
    for c in det.contacts
        t = hcat([ax1l ax1r; ax2l ax2r; ax3l ax3r], hcat.(ConstructiveSolidGeometry.extreme_points.(ConstructiveSolidGeometry.surfaces(c.geometry))...)...)
        (ax1l, ax1r), (ax2l, ax2r), (ax3l, ax3r) = broadcast(i -> extrema(t[i,:]), 1:3)
    end
    if !ismissing(det.passives)
        for p in det.passives
            t = hcat([ax1l ax1r; ax2l ax2r; ax3l ax3r], hcat.(ConstructiveSolidGeometry.extreme_points.(ConstructiveSolidGeometry.surfaces(p.geometry))...)...)
            (ax1l, ax1r), (ax2l, ax2r), (ax3l, ax3r) = broadcast(i -> extrema(t[i,:]), 1:3)
        end
    end
    return ax1l, ax1r, ax2l, ax2r, ax3l, ax3r
end

function SolidStateDetector{T}(config_file::Dict, input_units::NamedTuple) where {T <: SSDFloat}

    @assert !haskey(config_file, "objects") "Configuration file deprecation.\n
        The configuration file format was updated in v0.6.0.
        However, this configuration file still seems to be in the old format.\n
        To update your configuration file to the new format (v0.6.0 and newer),
        open a new Julia session and load the following file:\n
        \tinclude(\"<path_to_SolidStateDetectors.jl>/test/update_config_files.jl\")\n
        Afterwards, run\n
        \tupdate_config_file(\"<path_to_configuration_file>\")\n
        This method returns the file name of the updated configuration file.
        Please close the Julia session after updating the configuration files, as some
        parsing methods are overridden with old methods.\n
        Note that if your old geometries were defined using `difference`, you might
        need to increase the dimensions of the subtracted geometry as it is now treated as
        open primitive. Please have a look at the documentation:\n
        \thttps://juliaphysics.github.io/SolidStateDetectors.jl/stable/\n"

    @assert haskey(config_file, "detectors") "Config file needs an entry `detectors` that defines the detector(s)."
    config_detector = config_file["detectors"][1] # still only one detector

    transformations = parse_CSG_transformation(T, config_detector, input_units)

    @assert haskey(config_detector, "semiconductor") "Each detector needs an entry `semiconductor`. Please define the semiconductor."
    semiconductor = Semiconductor{T}(config_detector["semiconductor"], input_units, transformations)

    @assert haskey(config_detector, "contacts") "Each detector needs at least two contacts. Please define the them in the configuration file."
    contacts = broadcast(c -> Contact{T}(c, input_units, transformations), config_detector["contacts"])

    passives = []
    if haskey(config_detector, "passives") # "passives" as entry of "detectors"
        append!(passives, broadcast(p -> Passive{T}(p, input_units, transformations), config_detector["passives"]))
    end
    if haskey(config_file, "surroundings") # "surroundings" as entry in the configuration file
        append!(passives, broadcast(p -> Passive{T}(p, input_units, parse_CSG_transformation(T, p, input_units)), config_file["surroundings"]))
    end
    if isempty(passives)
        passives = missing
    end

    virtual_drift_volumes = if haskey(config_detector, "virtual_drift_volumes")
        broadcast(v -> construct_virtual_volume(T, v, input_units, transformations), config_detector["virtual_drift_volumes"])
    else
        missing
    end

    name = haskey(config_file, "name") ? config_file["name"] : "NoNameDetector"
    SolidStateDetector{T}( name, semiconductor, contacts, passives, virtual_drift_volumes )
end

function SolidStateDetector(dict::Dict)
    SolidStateDetector{Float32}(dict)
end

function in(pt::AbstractCoordinatePoint{T}, det::SolidStateDetector{T})::Bool where T
    in(pt, det.semiconductor) || reduce((x, contact) -> x || in(pt, contact), det.contacts, init = false)
end

function println(io::IO, det::SolidStateDetector{T}) where {T <: SSDFloat}
    println("________"*det.name*"________\n")
    println("---General Properties---")
    println("- Precision type: $(T)")
    println()
    println("\t_____Semiconductor_____\n")
    println(det.semiconductor)
    println()
    println("# Contacts: $(length(det.contacts))")
    if length(det.contacts)<=5
        for contact in det.contacts
            println(contact)
        end
    end
    if !ismissing(det.passives)
        println()
        println("# Passives: $(length(det.passives))")
        if length(det.passives) <= 5
            for passive in det.passives
                println(passive)
            end
        end
    end
end

function show(io::IO, det::SolidStateDetector{T}) where {T <: SSDFloat} println(io, det) end
function print(io::IO, det::SolidStateDetector{T}) where {T <: SSDFloat} println(io, det) end
function show(io::IO,::MIME"text/plain", det::SolidStateDetector) where {T <: SSDFloat} show(io, det) end


text: struct SaturnType <: abstractPlanet end


text:
using FractionalCalculus
using SymbolicUtils
using Test

macro eqtest(expr)
    @assert expr.head == :call && expr.args[1] in [:(==), :(!=)]
    if expr.args[1] == :(==)
        :(@test isequal($(expr.args[2]), $(expr.args[3])))
    else
        :(@test !isequal($(expr.args[2]), $(expr.args[3])))
    end |> esc
end

@testset "Test INT CONSTANT AND POWERS" begin
    @syms x
    @eqtest semiint(x^3) == 0.51583047638652*(x^3.5)
    @eqtest semiint(log(x)) == 2*log(4*x - 2)*sqrt(0.3183098861837907*x)
    #@eqtest semiint(sqrt(1+x)) == sqrt(x/pi) + ((1+x)*atan(sqrt(x)))/sqrt(pi)
    #@eqtest semiint(sqrt(1-x)) == sqrt(x/pi) + ((1-x)*atanh(sqrt(x)))/sqrt(pi)
    @eqtest semiint(1/sqrt(1+x)) == 2/sqrt(pi)*atan(sqrt(x))
    #@eqtest semiint(1/sqrt(1-x)) == 2/sqrt(pi)*atanh(sqrt(x))
end


text:
module PFFRGSolver

    # load external modules
    using LinearAlgebra
    using StaticArrays
    using HDF5
    using QuadGK
    using Cubature
    using LoopVectorization
    using Dates
    using Test
    using TimerOutputs

    # load source code
    include("Lattice/Lattice.jl")
    include("Frequency/Frequency.jl")
    include("Action/Action.jl")
    include("Flow/Flow.jl")
    include("Observable/Observable.jl")
    include("Launcher/Launcher.jl")
    include("Timers/Timers.jl")

    # export types, structs and functions
    export
        # from Lattice/unitcell.jl
        Unitcell,
        lattice_avail,
        get_unitcell,

        # from Lattice/site.jl
        Site,
        get_metric,
        get_nbs,

        # from Lattice/bond.jl
        Bond,

        # from Lattice/build.jl
        Lattice,
        get_lattice,
        model_avail,
        init_model!,
        get_site,
        get_bond,

        # from Lattice/reduced.jl
        Reduced_lattice,
        get_trafos_orig,
        get_trafos_uc,
        get_reduced_lattice,

        # from Lattice/disk.jl
        read_lattice,

        # from Lattice/test.jl
        test_lattice,

        # from Lattice/timers.jl
        get_lattice_timers,

        # from Frequency/param.jl
        Param,
        get_param,

        # from Frequency/mesh.jl
        Mesh,
        get_mesh,

        # from Frequency/buffer.jl
        Buffer,

        # from Frequency/test.jl
        test_frequencies,

        # from Frequency/timers.jl
        get_frequency_timers,

        # from Action/channel.jl
        Channel,
        get_abs_max,

        # from Action/vertex.jl
        Vertex,

        # from Action/Action.jl
        Action,
        Action_su2,
        Action_u1_dm,
        read_checkpoint,

        # from Action/disk.jl
        read_self,

        # from Action/test.jl
        test_action,

        # from Action/timers.jl
        get_action_timers,

        # from Flow/test.jl
        test_flow,

        # from Flow/timers.jl
        get_flow_timers,

        # from Observable/momentum.jl
        get_momenta,
        get_path,
        compute_structure_factor,

        # from Observable/disk.jl
        read_χ_labels,
        read_χ,
        read_χ_all,
        read_χ_flow_at_site,
        compute_structure_factor_flow!,
        compute_structure_factor_flow_all!,
        read_structure_factor_labels,
        read_structure_factor,
        read_structure_factor_all,
        read_structure_factor_flow_at_momentum,
        read_reference_momentum,

        # from Observable/fluctuations.jl
        compute_fluctuations,
        compute_fluctuations_flow,

        # from Observable/test.jl
        test_observable,

        # from Observable/timers.jl
        get_observable_timers,

        # from Launcher/Launcher.jl
        save_launcher!,
        make_job!,
        make_repository!,
        collect_repository!,
        launch!,

        # from Timers/Timers.jl
        get_PFFRG_timers
end


text:
# This file contains all the functions related to the obstacles

# this function computes the position of the obstacle
function obstaclePosition(obs_now::Array{Float64},modelParams::ModelParams,obstacle::Obstacle,posInfo::PosInfo)

    n_obs = obstacle.n_obs
    obs_next=zeros(1,3,n_obs)::Array{Float64}    # initialize the array that will contain the next states of the obstacle

    dt = modelParams.dt

    obs_next[:,1,:] = (obs_now[:,1,:] + obs_now[:,3,:]*dt)%posInfo.s_target # s
    obs_next[:,2,:] = obs_now[:,2,:]                     # ey
    obs_next[:,3,:] = obs_now[:,3,:]                     # v

    # if any(x->x>posInfo.s_target, obs_next[1,1,:]) == true

    #     index = find(obs_next[1,1,:].>posInfo.s_target)

    #     obs_next[1,1,index] = obs_next[1,1,index] - posInfo.s_target
    # end

    return obs_next

end

max_stars_repo_name: RoyaFiroozi/barc

max_issues_repo_name: RoyaFiroozi/barc

max_forks_repo_name: RoyaFiroozi/barc


text:
################################################################################
#
#  Basic field access
#
################################################################################

@doc Markdown.doc"""
    nf(x::NumFieldOrdIdl) -> AnticNumberField

Returns the number field, of which $x$ is an integral ideal.
"""
nf(x::NumFieldOrdIdl) = nf(order(x))

@doc Markdown.doc"""
    order(I::NumFieldOrdIdl) -> NfOrd

Returns the order of $I$.
"""
@inline order(a::NumFieldOrdIdl) = a.order

_algebra(a::NumFieldOrdIdl) = nf(a)

################################################################################
#
#   Degree and ramification index
#
################################################################################

function degree(P::NfRelOrdIdl)
  @assert isprime(P)
  return P.splitting_type[2]*degree(minimum(P))
end

function ramification_index(P::NfRelOrdIdl)
  @assert isprime(P)
  return P.splitting_type[1]
end

function absolute_ramification_index(P::NfRelOrdIdl)
  @assert isprime(P)
  return P.splitting_type[1]*absolute_ramification_index(minimum(P))
end

function absolute_ramification_index(P::NfAbsOrdIdl)
  @assert isprime(P)
  return ramification_index(P)
end

################################################################################
#
#   Absolute basis
#
################################################################################

function absolute_basis(I::T) where T <: Union{NfAbsOrdIdl, NfAbsOrdFracIdl}
  return basis(I)
end

function absolute_basis(I::T) where T <: Union{NfRelOrdIdl, NfRelOrdFracIdl}
  res = elem_type(nf(order(I)))[]
  pb = pseudo_basis(I)
  pbb = pseudo_basis(order(I))
  for i in 1:length(pb)
    (e, I) = pb[i]
    for b in absolute_basis(I)
      push!(res, e * b)
    end
  end
  return res
end

################################################################################
#
#   Absolute Minimum
#
################################################################################

@doc Markdown.doc"""
    absolute_minimum(I::NumFieldOrdIdl) -> fmpz

Given an ideal $I$, returns a generator of the ideal $I \cap \mathbb Z$.
"""
absolute_minimum(::NumFieldOrdIdl)

function absolute_minimum(I::NfAbsOrdIdl)
  return minimum(I)
end

function absolute_minimum(I::NfRelOrdIdl)
  return absolute_minimum(minimum(I))::fmpz
end

################################################################################
#
#   Absolute norm
#
################################################################################

@doc Markdown.doc"""
    absolute_norm(I::NumFieldOrdIdl) -> fmpz

Given an ideal $I$, returns its norm over $\mathbb Q$.
"""
absolute_norm(::NumFieldOrdIdl)

function absolute_norm(x::NfAbsOrdIdl)
  return norm(x)
end

function absolute_norm(a::NfRelOrdIdl)
  return norm(a, FlintQQ)
end

function norm(I::NumFieldOrdIdl, K::NumField)
  O = order(I)
  if K == base_field(O)
    return norm(I)
  else
    return norm(norm(I), K)
  end
end

function norm(I::NumFieldOrdIdl, ::FlintRationalField)
  return absolute_norm(I)
end

################################################################################
#
#  Uniformizer
#
################################################################################

@doc Markdown.doc"""
    uniformizer(P::NumFieldOrdIdl) -> NumFieldOrdElem

Given a prime ideal $P$, returns an element $u \in P$ with valuation(u, P) == 1.
"""
function uniformizer(P::NfRelOrdIdl)
  @assert isprime(P)

  if isone(ramification_index(P))
    return order(P)(uniformizer(minimum(P, copy = false)))
  end

  r = 500 # hopefully enough
  z = rand(P, r)
  while true
    if !iszero(z) && valuation(z, P) == 1
      break
    end
    z = rand(P, r)
  end
  return z
end

function uniformizer(P::NfAbsOrdIdl)
  @assert isprime(P)
  p = minimum(P)
  if isdefined(P, :gens_normal) && P.gens_normal == p
    return P.gen_two
  elseif isone(ramification_index(P))
    return order(P)(p)
  else
    if p > 250
      r = 500  # should still have enough elements...
    else
      r = Int(div(p, 2))
    end
    z = rand(P, r)
    while true
      if !iszero(z) && valuation(z, P) == 1
        break
      end
      z = rand(P, r)
    end
    return z
  end
end

################################################################################
#
#  p-uniformizer
#
################################################################################

@doc Markdown.doc"""
    p_uniformizer(P::NumFieldOrdIdl) -> NumFieldOrdElem

Given a prime ideal P, returns an element $u \in P$ with valuation(u, P) == 1 and valuation 0 at all
other prime ideals lying over minimum(P).
"""
p_uniformizer(::NumFieldOrdIdl)

function p_uniformizer(P::NfAbsOrdIdl)
  assure_2_normal(P)
  return P.gen_two
end

function p_uniformizer(P::NfRelOrdIdl{S, T, U}) where {S, T, U}
  @assert isprime(P)

  if isdefined(P, :p_uniformizer)
    return P.p_uniformizer::elem_type(order(P))
  end

  p = minimum(P, copy = false)
  prime_dec = prime_decomposition(order(P), p)
  primes = Vector{typeof(P)}()
  for (PP, e) in prime_dec
    if PP != P
      push!(primes, PP)
    end
  end
  P2 = P^2
  r = 500
  z = rand(P, r)
  while !_is_p_uniformizer(z, P, P2, primes)
    z = rand(P, r)
  end
  P.p_uniformizer = z
  return z
end

################################################################################
#
#  Anti uniformizer
#
################################################################################

@doc Markdown.doc"""
    anti_uniformizer(P::NumFieldOrdIdl) -> NumFieldElem

Given a prime ideal $P$, returns an element $a$ in the number field containing $P$
with valuation(a, P) == -1, valuation(a, Q) = 0 at the prime conjugate to $P$
and non-negative valuation at all other prime ideals.
"""
anti_uniformizer(::NumFieldOrdIdl)

################################################################################
#
#   Absolute anti uniformizer
#
################################################################################

@doc Markdown.doc"""
    absolute_anti_uniformizer(P::NumFieldOrdIdl) -> NumFieldElem

Given a prime ideal $P$, this function returns an element $x$ with valuation(x, P) = -1$,
valuation(x, Q) = 0$ for every ideal Q conjugate to $P$ and non-negative valuation
at any other prime ideal.
"""
absolute_anti_uniformizer(::NumFieldOrdIdl)

function absolute_anti_uniformizer(P::NfAbsOrdIdl)
  return anti_uniformizer(P)
end

function absolute_anti_uniformizer(P::NfRelOrdIdl)
  OL = order(P)
  L = nf(OL)
  A = absolute_basis(OL)
  d = absolute_degree(nf(OL))
  OLmat = zero_matrix(FlintQQ, d, d)
  Lbas = absolute_basis(L)
  for i in 1:d
    c = elem_in_nf(A[i], copy = false)
    coords = absolute_coordinates(c)
    for j in 1:d
      OLmat[i, j] = coords[j]
    end
  end

  OLmatinv = inv(OLmat)

  u = elem_in_nf(p_uniformizer(P))

  @show isintegral(u)

  umat = zero_matrix(FlintQQ, d, d)

  for i in 1:d
    c = u * Lbas[i]
    coordsc = absolute_coordinates(c)
    for j in 1:d
      umat[i, j] = coordsc[j]
    end
  end

  N = OLmat * umat * OLmatinv

  p = absolute_minimum(P)

  z = zero_matrix(GF(p, cached = false), d, d)

  for i in 1:d
    for j in 1:d
      z[i, j] = FlintZZ(N[i, j])
    end
  end

  K = left_kernel_basis(z)

  k = K[1]
  return inv(L(p)) * elem_in_nf(sum(elem_type(OL)[A[i] * lift(k[i]) for i in 1:d]))
end

################################################################################
#
#  Prime number in a prime ideal
#
################################################################################

@doc Markdown.doc"""
    prime_number(P::NumFieldOrdIdl) -> fmpz

Given a prime ideal $P$, returns the unique prime number $p$ contained in $P$.
"""
function prime_number(P::NumFieldOrdIdl; check::Bool = true)
  if check
    @assert isprime(P)
  end
  return prime_number(minimum(P), check = false)
end

function prime_number(P::NfAbsOrdIdl; check::Bool = true)
  if check
    @assert isprime(P)
  end
  return minimum(P)
end

################################################################################
#
#   Support
#
################################################################################

function support(I::T) where T <: Union{NumFieldOrdIdl, NumFieldOrdFracIdl}
  lp = factor(I)
  return collect(keys(lp))
end

################################################################################
#
#   Is integral
#
################################################################################

isintegral(I::NumFieldOrdIdl) = true

function isintegral(I::NfOrdFracIdl)
  @assert ismaximal(order(I))
  simplify(I)
  return denominator(I) == 1
end

function isintegral(a::NfRelOrdFracIdl)
  @assert ismaximal(order(a))
  return defines_ideal(order(a), basis_pmatrix(a, copy = false))
end

################################################################################
#
#   Trace ideal
#
################################################################################

function tr(I::NfAbsOrdIdl)
  E = nf(order(I))
  K = base_field(E)
  return gcd(fmpz[tr(x) for x in basis(I)])
end

function tr(I::NfAbsOrdFracIdl)
  E = nf(order(I))
  K = base_field(E)
  traces = fmpq[trace(b) for b in basis(I)]
  #TODO: This is deeply wrong.
  return reduce(gcd, traces; init = fmpq(0))
end

function tr(I::T) where T <: Union{NfRelOrdIdl, NfRelOrdFracIdl}
  E = nf(order(I))
  K = base_field(E)
  return fractional_ideal(maximal_order(K), elem_type(K)[trace(b) for b in absolute_basis(I)])
end

################################################################################
#
#   Gens
#
################################################################################

gens(I::NumFieldOrdIdl) = small_generating_set(I)

function gens(I::NumFieldOrdFracIdl)
  K = nf(order(I))
  nI = numerator(I)
  dI = denominator(I)
  gnI = gens(nI)
  return elem_type(K)[elem_in_nf(x)//dI for x in gnI]
end

function small_generating_set(I::NfAbsOrdIdl)
  OK = order(I)
  if isone(I)
    return elem_type(OK)[one(OK)]
  end
  if has_2_elem(I)
    return elem_type(OK)[OK(I.gen_one), OK(I.gen_two)]
  end
  if ismaximal_known_and_maximal(OK)
    _assure_weakly_normal_presentation(I)
    return elem_type(OK)[OK(I.gen_one), OK(I.gen_two)]
  end
  id_gen = zero_matrix(FlintZZ, 2*n, n)
  m = minimum(I, copy = false)
  B = basis(I, copy = false)
  gens = NfOrdElem[]
  for i = 1:length(B)
    if i != 1
      c = matrix(FlintZZ, 1, n, coordinates(B[i]))
      reduce_mod_hnf_ll!(c, id_gen)
      if iszero(c)
        continue
      end
    end
    M = representation_matrix_mod(B[i], modu)
    _copy_matrix_into_matrix(id_gen, 1, 1, M)
    hnf_modular_eldiv!(id_gen, m, :lowerleft)
    push!(gens, B[i])
    if view(id_gen, n+1:2*n, 1:n) == basis_matrix(a, copy = false)
      break
    end
  end
  return gens
end

function small_generating_set(I::NfRelOrdIdl)
  OK = order(I)
  K = nf(OK)
  B = pseudo_basis(I, copy = false)
  starting_gens = elem_type(OK)[]
  for i = 1:length(B)
    gensI = small_generating_set(numerator(B[i][2], copy = false))
    for x in gensI
      push!(starting_gens, OK(divexact(elem_in_nf(x, copy = false)*B[i][1], denominator(B[i][2], copy = false))))
    end
  end
  #Now, I have a set of generators as a OK-module.
  #Let's discard the non relevant elements
  indices = Int[]
  Id = ideal(OK, 0)
  id_gen = pseudo_matrix(zero_matrix(base_field(K), 0, degree(OK)))
  for i = 1:length(starting_gens)
    if i != 1
      if starting_gens[i] in Id
        continue
      end
    end
    Id = Id + ideal(OK, starting_gens[i])
    push!(indices, i)
    if Id == I
      break
    end
  end
  return starting_gens[indices]
end

function isramified(O::NumFieldOrd, P::NumFieldOrdIdl)
  OK = order(P)
  d = discriminant(O, nf(OK))
  return !iscoprime(P, d)
end


text:
# Algorithm 3.3
# Hermite Interpolation

function hermite_interpolate(x, y, yp)
	n = length(x)
	Q = zeros(2n + 1, 2n + 1)
	z = zeros[1, 2n + 1]
	for i ∈ 1:n
		z[2i] = x[i]
		z[2i + 1] = x[i]
		Q[2i, 1] = y[i]
		Q[2i + 1, 2] = yp[i]

		if i ≠ 1
			Q[2i, 2] = (Q[2i, 1] - Q[2i - 1, 1])/(z[2i] - z[2i - 1])
		end
	end
	for i ∈ 3:(2n + 1)
		for j ∈ 3:i
			Q[i, j] = (Q[i, j - 1] - Q[i - 1, j - 1])/(z[i] - z[i - j])
		end
	end
	return diag(Q)
end

print(hermite_interpolate([1, 2, 3], [0.84, 0,91, 0,14], [0.54, -0.41, -0.98]))


text:
#---------------------------------------------------------------------#
#This module stores everything needed to build different qudrature points
#in 1D.
#Written by F.X. Giraldo on 4/10/2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey, CA 93943-5216
#---------------------------------------------------------------------#
module QuadraturePoints
    include("chebyshev_gauss.jl")
    include("equispaced_points.jl")
    include("legendre_poly.jl")
    include("legendre_gauss.jl")
    include("lobatto_gauss.jl")
    include("lagrange_basis.jl")
end


text:
module steadystate

using ..states, ..operators, ..operators_dense, ..superoperators
using ..timeevolution

"""
    steadystate.master(H, J; <keyword arguments>)

Calculate steady state using long time master equation evolution.

# Arguments
* `H`: Arbitrary operator specifying the Hamiltonian.
* `J`: Vector containing all jump operators which can be of any arbitrary
        operator type.
* `rho0=dm(basisstate(b))`: Initial density operator. If not given the
        ``|0⟩⟨0|`` state in respect to the choosen basis is used.
* `tol=1e-3`: Tracedistance used as termination criterion.
* `hmin=1e-7`: Minimal time step used in the time evolution.
* `rates=ones(N)`: Vector or matrix specifying the coefficients for the
        jump operators.
* `Jdagger=dagger.(Jdagger)`: Vector containing the hermitian conjugates of the
        jump operators. If they are not given they are calculated automatically.
* `fout=nothing`: If given this function `fout(t, rho)` is called every time an
        output should be displayed. To limit copying to a minimum the given
        density operator `rho` is further used and therefore must not be changed.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master(H::Operator, J::Vector;
                rho0::DenseOperator=tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1))),
                hmin=1e-7, tol=1e-3,
                rates::Union{Vector{Float64}, Matrix{Float64}, Void}=nothing,
                Jdagger::Vector=dagger.(J),
                fout::Union{Function,Void}=nothing,
                kwargs...)
    t,u = timeevolution.master([0., Inf], rho0, H, J; rates=rates, Jdagger=Jdagger,
                        hmin=hmin, hmax=Inf,
                        display_initialvalue=false,
                        display_finalvalue=false,
                        display_intermediatesteps=true,
                        fout=fout,
                        steady_state = true,
                        tol = tol, kwargs...)
end

"""
    steadystate.eigenvector(L)
    steadystate.eigenvector(H, J)

Find steady state by calculating the eigenstate of the Liouvillian matrix `l`.
"""
function eigenvector(L::DenseSuperOperator; kwargs...)
    d, v = Base.eig(L.data; kwargs...)
    index = findmin(abs.(d))[2]
    data = reshape(v[:,index], length(L.basis_r[1]), length(L.basis_r[2]))
    op = DenseOperator(L.basis_r[1], L.basis_r[2], data)
    return op/trace(op)
end

function eigenvector(L::SparseSuperOperator; sigma::Float64=1e-30, kwargs...)
    d, v, nconv, niter, nmult, resid = try
      Base.eigs(L.data; nev=1, sigma=sigma, kwargs...)
    catch err
      if isa(err, LinAlg.SingularException)
        error("Base.LinAlg.eigs() algorithm failed; try using DenseOperators")
      else
        rethrow(err)
      end
    end
    data = reshape(v[:,1], length(L.basis_r[1]), length(L.basis_r[2]))
    op = DenseOperator(L.basis_r[1], L.basis_r[2], data)
    return op/trace(op)
end

eigenvector(H::Operator, J::Vector; rates::Union{Vector{Float64}, Matrix{Float64}}=ones(Float64, length(J)), kwargs...) = eigenvector(liouvillian(H, J; rates=rates); kwargs...)

end # module


text:
module SpiderType
    @enum SType Z X H In Out
end  # module SpiderType

"""
    ZXDiagram{T, P}
This is the type for representing ZX-diagrams.
"""
struct ZXDiagram{T<:Integer, P} <: AbstractZXDiagram{T, P}
    mg::Multigraph{T}

    st::Dict{T, SpiderType.SType}
    ps::Dict{T, P}

    layout::ZXLayout{T}
    phase_ids::Dict{T, Tuple{T, Int}}

    scalar::Scalar{P}
    inputs::Vector{T}
    outputs::Vector{T}

    function ZXDiagram{T, P}(mg::Multigraph{T}, st::Dict{T, SpiderType.SType}, ps::Dict{T, P},
        layout::ZXLayout{T}, phase_ids::Dict{T, Tuple{T, Int}} = Dict{T, Tuple{T, Int}}(),
        s::Scalar{P} = Scalar{P}(),
        inputs::Vector{T} = Vector{T}(), outputs::Vector{T} = Vector{T}()) where {T<:Integer, P}
        if nv(mg) == length(ps) && nv(mg) == length(st)
            if length(phase_ids) == 0
                for v in vertices(mg)
                    if st[v] in (SpiderType.Z, SpiderType.X)
                        phase_ids[v] = (v, 1)
                    end
                end
            end
            if length(inputs) == 0
                for v in vertices(mg)
                    if st[v] == SpiderType.In
                        push!(inputs, v)
                    end
                end
                if layout.nbits > 0
                    sort!(inputs, by = (v -> qubit_loc(layout, v)))
                end
            end
            if length(outputs) == 0
                for v in vertices(mg)
                    if st[v] == SpiderType.Out
                        push!(outputs, v)
                    end
                end
                if layout.nbits > 0
                    sort!(outputs, by = (v -> qubit_loc(layout, v)))
                end
            end
            zxd = new{T, P}(mg, st, ps, layout, phase_ids, s, inputs, outputs)
            round_phases!(zxd)
            return zxd
        else
            error("There should be a phase and a type for each spider!")
        end
    end
end

"""
    ZXDiagram(mg::Multigraph{T}, st::Dict{T, SpiderType.SType}, ps::Dict{T, P},
        layout::ZXLayout{T} = ZXLayout{T}(),
        phase_ids::Dict{T,Tuple{T, Int}} = Dict{T,Tuple{T,Int}}()) where {T, P}
    ZXDiagram(mg::Multigraph{T}, st::Vector{SpiderType.SType}, ps::Vector{P},
        layout::ZXLayout{T} = ZXLayout{T}()) where {T, P}

Construct a ZXDiagram with all information.

```jldoctest
julia> using LightGraphs, Multigraphs, ZXCalculus;

julia> using ZXCalculus.SpiderType: In, Out, H, Z, X;

julia> mg = Multigraph(5);

julia> for i = 1:4
           add_edge!(mg, i, i+1)
       end;

julia> ZXDiagram(mg, [In, Z, H, X, Out], [0//1, 1, 0, 1//2, 0])
ZX-diagram with 5 vertices and 4 multiple edges:
(S_1{input} <-1-> S_2{phase = 1//1⋅π})
(S_2{phase = 1//1⋅π} <-1-> S_3{H})
(S_3{H} <-1-> S_4{phase = 1//2⋅π})
(S_4{phase = 1//2⋅π} <-1-> S_5{output})

```
"""
ZXDiagram(mg::Multigraph{T}, st::Dict{T, SpiderType.SType}, ps::Dict{T, P},
    layout::ZXLayout{T} = ZXLayout{T}(),
    phase_ids::Dict{T,Tuple{T, Int}} = Dict{T,Tuple{T,Int}}()) where {T, P} = ZXDiagram{T, P}(mg, st, ps, layout, phase_ids)
ZXDiagram(mg::Multigraph{T}, st::Vector{SpiderType.SType}, ps::Vector{P},
    layout::ZXLayout{T} = ZXLayout{T}()) where {T, P} =
    ZXDiagram(mg, Dict(zip(sort!(vertices(mg)), st)), Dict(zip(sort!(vertices(mg)), ps)), layout)

"""
    ZXDiagram(nbits)

Construct a ZXDiagram of a empty circuit with qubit number `nbit`

```jldoctest; setup = :(using ZXCalculus)
julia> zxd = ZXDiagram(3)
ZX-diagram with 6 vertices and 3 multiple edges:
(S_1{input} <-1-> S_2{output})
(S_3{input} <-1-> S_4{output})
(S_5{input} <-1-> S_6{output})

```
"""
function ZXDiagram(nbits::T) where {T<:Integer}
    mg = Multigraph(2*nbits)
    st = [SpiderType.In for _ = 1:2*nbits]
    ps = [Phase(0//1) for _ = 1:2*nbits]
    spider_q = Dict{T, Rational{Int}}()
    spider_col = Dict{T, Rational{Int}}()
    for i = 1:nbits
        add_edge!(mg, 2*i-1, 2*i)
        @inbounds st[2*i] = SpiderType.Out
        spider_q[2*i-1] = i
        spider_col[2*i-1] = 2
        spider_q[2*i] = i
        spider_col[2*i] = 1
    end
    layout = ZXLayout(nbits, spider_q, spider_col)
    return ZXDiagram(mg, st, ps, layout)
end

Base.copy(zxd::ZXDiagram{T, P}) where {T, P} = ZXDiagram{T, P}(copy(zxd.mg), copy(zxd.st), copy(zxd.ps), copy(zxd.layout),
    deepcopy(zxd.phase_ids), copy(zxd.scalar), copy(zxd.inputs), copy(zxd.outputs))

"""
    spider_type(zxd, v)

Returns the spider type of a spider.
"""
spider_type(zxd::ZXDiagram{T, P}, v::T) where {T<:Integer, P} = zxd.st[v]

"""
    phase(zxd, v)

Returns the phase of a spider. If the spider is not a Z or X spider, then return 0.
"""
phase(zxd::ZXDiagram{T, P}, v::T) where {T<:Integer, P} = zxd.ps[v]

"""
    set_phase!(zxd, v, p)

Set the phase of `v` in `zxd` to `p`.
"""
function set_phase!(zxd::ZXDiagram{T, P}, v::T, p::P) where {T, P}
    if has_vertex(zxd.mg, v)
        while p < 0
            p += 2
        end
        p = rem(p, 2)
        zxd.ps[v] = p
        return true
    end
    return false
end

"""
    nqubits(zxd)

Returns the qubit number of a ZX-diagram.
"""
nqubits(zxd::ZXDiagram) = zxd.layout.nbits

function print_spider(io::IO, zxd::ZXDiagram{T, P}, v::T) where {T<:Integer, P}
    st_v = spider_type(zxd, v)
    if st_v == SpiderType.Z
        printstyled(io, "S_$(v){phase = $(zxd.ps[v])"*(zxd.ps[v] isa Phase ? "}" : "⋅π}"); color = :green)
    elseif st_v == SpiderType.X
        printstyled(io, "S_$(v){phase = $(zxd.ps[v])"*(zxd.ps[v] isa Phase ? "}" : "⋅π}"); color = :red)
    elseif st_v == SpiderType.H
        printstyled(io, "S_$(v){H}"; color = :yellow)
    elseif st_v == SpiderType.In
        print(io, "S_$(v){input}")
    elseif st_v == SpiderType.Out
        print(io, "S_$(v){output}")
    end
end

function Base.show(io::IO, zxd::ZXDiagram{T, P}) where {T<:Integer, P}
    println(io, "ZX-diagram with $(nv(zxd.mg)) vertices and $(ne(zxd.mg)) multiple edges:")
    for v1 in sort!(vertices(zxd.mg))
        for v2 in neighbors(zxd.mg, v1)
            if v2 >= v1
                print(io, "(")
                print_spider(io, zxd, v1)
                print(io, " <-$(mul(zxd.mg, v1, v2))-> ")
                print_spider(io, zxd, v2)
                print(io, ")\n")
            end
        end
    end
end

"""
    nv(zxd)

Returns the number of vertices (spiders) of a ZX-diagram.
"""
LightGraphs.nv(zxd::ZXDiagram) = nv(zxd.mg)

"""
    ne(zxd; count_mul = false)

Returns the number of edges of a ZX-diagram. If `count_mul`, it will return the
sum of multiplicities of all multiple edges. Otherwise, it will return the
number of multiple edges.
"""
LightGraphs.ne(zxd::ZXDiagram; count_mul::Bool = false) = ne(zxd.mg, count_mul = count_mul)

LightGraphs.outneighbors(zxd::ZXDiagram, v; count_mul::Bool = false) = outneighbors(zxd.mg, v, count_mul = count_mul)
LightGraphs.inneighbors(zxd::ZXDiagram, v; count_mul::Bool = false) = inneighbors(zxd.mg, v, count_mul = count_mul)

LightGraphs.degree(zxd::ZXDiagram, v::Integer) = degree(zxd.mg, v)
LightGraphs.indegree(zxd::ZXDiagram, v::Integer) = degree(zxd, v)
LightGraphs.outdegree(zxd::ZXDiagram, v::Integer) = degree(zxd, v)

"""
    neighbors(zxd, v; count_mul = false)

Returns a vector of vertices connected to `v`. If `count_mul`, there will be
multiple copy for each vertex. Otherwise, each vertex will only appear once.
"""
LightGraphs.neighbors(zxd::ZXDiagram, v; count_mul::Bool = false) = neighbors(zxd.mg, v, count_mul = count_mul)
function LightGraphs.rem_edge!(zxd::ZXDiagram, x...)
    rem_edge!(zxd.mg, x...)
end
function LightGraphs.add_edge!(zxd::ZXDiagram, x...)
    add_edge!(zxd.mg, x...)
end

"""
    rem_spiders!(zxd, vs)

Remove spiders indexed by `vs`.
"""
function rem_spiders!(zxd::ZXDiagram{T, P}, vs::Vector{T}) where {T<:Integer, P}
    if rem_vertices!(zxd.mg, vs)
        for v in vs
            delete!(zxd.ps, v)
            delete!(zxd.st, v)
            delete!(zxd.phase_ids, v)
            rem_vertex!(zxd.layout, v)
        end
        return true
    end
    return false
end

"""
    rem_spider!(zxd, v)

Remove a spider indexed by `v`.
"""
rem_spider!(zxd::ZXDiagram{T, P}, v::T) where {T<:Integer, P} = rem_spiders!(zxd, [v])

"""
    add_spider!(zxd, spider_type, phase = 0, connect = [])

Add a new spider which is of the type `spider_type` with phase `phase` and
connected to the vertices `connect`.
"""
function add_spider!(zxd::ZXDiagram{T, P}, st::SpiderType.SType, phase::P = zero(P), connect::Vector{T}=T[]) where {T<:Integer, P}
    v = add_vertex!(zxd.mg)[1]
    set_phase!(zxd, v, phase)
    zxd.st[v] = st
    if st in (SpiderType.Z, SpiderType.X)
        zxd.phase_ids[v] = (v, 1)
    end
    if all(has_vertex(zxd.mg, c) for c in connect)
        for c in connect
            add_edge!(zxd.mg, v, c)
        end
    end
    return v
end

"""
    insert_spider!(zxd, v1, v2, spider_type, phase = 0)

Insert a spider of the type `spider_type` with phase = `phase`, between two
vertices `v1` and `v2`. It will insert multiple times if the edge between
`v1` and `v2` is a multiple edge. Also it will remove the original edge between
`v1` and `v2`.
"""
function insert_spider!(zxd::ZXDiagram{T, P}, v1::T, v2::T, st::SpiderType.SType, phase::P = zero(P)) where {T<:Integer, P}
    mt = mul(zxd.mg, v1, v2)
    vs = Vector{T}(undef, mt)
    for i = 1:mt
        v = add_spider!(zxd, st, phase, [v1, v2])
        @inbounds vs[i] = v
        rem_edge!(zxd, v1, v2)
    end
    return vs
end

"""
    round_phases!(zxd)

Round phases between [0, 2π).
"""
function round_phases!(zxd::ZXDiagram{T, P}) where {T<:Integer, P}
    ps = zxd.ps
    for v in keys(ps)
        while ps[v] < 0
            ps[v] += 2
        end
        ps[v] = rem(ps[v], 2)
    end
    return
end

spiders(zxd::ZXDiagram) = vertices(zxd.mg)
qubit_loc(zxd::ZXDiagram{T, P}, v::T) where {T, P} = qubit_loc(zxd.layout, v)
function column_loc(zxd::ZXDiagram{T, P}, v::T) where {T, P}
    c_loc = column_loc(zxd.layout, v)
    return c_loc
end

"""
    push_gate!(zxd, ::Val{M}, locs...[, phase]; autoconvert=true)

Push an `M` gate to the end of qubit `loc` where `M` can be `:Z`, `:X`, `:H`, `:SWAP`, `:CNOT` and `:CZ`.
If `M` is `:Z` or `:X`, `phase` will be available and it will push a
rotation `M` gate with angle `phase * π`.
If `autoconvert` is `false`, the input `phase` should be a rational numbers.
"""
function push_gate!(zxd::ZXDiagram{T, P}, ::Val{:Z}, loc::T, phase = zero(P); autoconvert::Bool=true) where {T, P}
    @inbounds out_id = get_outputs(zxd)[loc]
    @inbounds bound_id = neighbors(zxd, out_id)[1]
    rphase = autoconvert ? safe_convert(P, phase) : phase
    insert_spider!(zxd, bound_id, out_id, SpiderType.Z, rphase)
    return zxd
end

function push_gate!(zxd::ZXDiagram{T, P}, ::Val{:X}, loc::T, phase = zero(P); autoconvert::Bool=true) where {T, P}
    @inbounds out_id = get_outputs(zxd)[loc]
    @inbounds bound_id = neighbors(zxd, out_id)[1]
    rphase = autoconvert ? safe_convert(P, phase) : phase
    insert_spider!(zxd, bound_id, out_id, SpiderType.X, rphase)
    return zxd
end

function push_gate!(zxd::ZXDiagram{T, P}, ::Val{:H}, loc::T) where {T, P}
    @inbounds out_id = get_outputs(zxd)[loc]
    @inbounds bound_id = neighbors(zxd, out_id)[1]
    insert_spider!(zxd, bound_id, out_id, SpiderType.H)
    return zxd
end

function push_gate!(zxd::ZXDiagram{T, P}, ::Val{:SWAP}, locs::Vector{T}) where {T, P}
    q1, q2 = locs
    push_gate!(zxd, Val{:Z}(), q1)
    push_gate!(zxd, Val{:Z}(), q2)
    push_gate!(zxd, Val{:Z}(), q1)
    push_gate!(zxd, Val{:Z}(), q2)
    v1, v2, bound_id1, bound_id2 = (sort!(spiders(zxd)))[end-3:end]
    rem_edge!(zxd, v1, bound_id1)
    rem_edge!(zxd, v2, bound_id2)
    add_edge!(zxd, v1, bound_id2)
    add_edge!(zxd, v2, bound_id1)
    return zxd
end

function push_gate!(zxd::ZXDiagram{T, P}, ::Val{:CNOT}, loc::T, ctrl::T) where {T, P}
    push_gate!(zxd, Val{:Z}(), ctrl)
    push_gate!(zxd, Val{:X}(), loc)
    @inbounds v1, v2 = (sort!(spiders(zxd)))[end-1:end]
    add_edge!(zxd, v1, v2)
    add_power!(zxd, 1)
    return zxd
end

function push_gate!(zxd::ZXDiagram{T, P}, ::Val{:CZ}, loc::T, ctrl::T) where {T, P}
    push_gate!(zxd, Val{:Z}(), ctrl)
    push_gate!(zxd, Val{:Z}(), loc)
    @inbounds v1, v2 = (sort!(spiders(zxd)))[end-1:end]
    add_edge!(zxd, v1, v2)
    insert_spider!(zxd, v1, v2, SpiderType.H)
    add_power!(zxd, 1)
    return zxd
end

"""
    pushfirst_gate!(zxd, ::Val{M}, loc[, phase])

Push an `M` gate to the beginning of qubit `loc` where `M` can be `:Z`, `:X`, `:H`, `:SWAP`, `:CNOT` and `:CZ`.
If `M` is `:Z` or `:X`, `phase` will be available and it will push a
rotation `M` gate with angle `phase * π`.
"""
function pushfirst_gate!(zxd::ZXDiagram{T, P}, ::Val{:Z}, loc::T, phase::P = zero(P)) where {T, P}
    @inbounds in_id = get_inputs(zxd)[loc]
    @inbounds bound_id = neighbors(zxd, in_id)[1]
    insert_spider!(zxd, in_id, bound_id, SpiderType.Z, phase)
    return zxd
end

function pushfirst_gate!(zxd::ZXDiagram{T, P}, ::Val{:X}, loc::T, phase::P = zero(P)) where {T, P}
    @inbounds in_id = get_inputs(zxd)[loc]
    @inbounds bound_id = neighbors(zxd, in_id)[1]
    insert_spider!(zxd, in_id, bound_id, SpiderType.X, phase)
    return zxd
end

function pushfirst_gate!(zxd::ZXDiagram{T, P}, ::Val{:H}, loc::T) where {T, P}
    @inbounds in_id = get_inputs(zxd)[loc]
    @inbounds bound_id = neighbors(zxd, in_id)[1]
    insert_spider!(zxd, in_id, bound_id, SpiderType.H)
    return zxd
end

function pushfirst_gate!(zxd::ZXDiagram{T, P}, ::Val{:SWAP}, locs::Vector{T}) where {T, P}
    q1, q2 = locs
    pushfirst_gate!(zxd, Val{:Z}(), q1)
    pushfirst_gate!(zxd, Val{:Z}(), q2)
    pushfirst_gate!(zxd, Val{:Z}(), q1)
    pushfirst_gate!(zxd, Val{:Z}(), q2)
    @inbounds v1, v2, bound_id1, bound_id2 = (sort!(spiders(zxd)))[end-3:end]
    rem_edge!(zxd, v1, bound_id1)
    rem_edge!(zxd, v2, bound_id2)
    add_edge!(zxd, v1, bound_id2)
    add_edge!(zxd, v2, bound_id1)
    return zxd
end

function pushfirst_gate!(zxd::ZXDiagram{T, P}, ::Val{:CNOT}, loc::T, ctrl::T) where {T, P}
    pushfirst_gate!(zxd, Val{:Z}(), ctrl)
    pushfirst_gate!(zxd, Val{:X}(), loc)
    @inbounds v1, v2 = (sort!(spiders(zxd)))[end-1:end]
    add_edge!(zxd, v1, v2)
    add_power!(zxd, 1)
    return zxd
end

function pushfirst_gate!(zxd::ZXDiagram{T, P}, ::Val{:CZ}, loc::T, ctrl::T) where {T, P}
    pushfirst_gate!(zxd, Val{:Z}(), ctrl)
    pushfirst_gate!(zxd, Val{:Z}(), loc)
    @inbounds v1, v2 = (sort!(spiders(zxd)))[end-1:end]
    add_edge!(zxd, v1, v2)
    insert_spider!(zxd, v1, v2, SpiderType.H)
    add_power!(zxd, 1)
    return zxd
end

"""
    tcount(zxd)

Returns the T-count of a ZX-diagram.
"""
tcount(cir::ZXDiagram) = sum([phase(cir, v) % 1//2 != 0 for v in spiders(cir)])

"""
    get_inputs(zxd)

Returns a vector of input ids.
"""
get_inputs(zxd::ZXDiagram) = zxd.inputs

"""
    get_outputs(zxd)

Returns a vector of output ids.
"""
get_outputs(zxd::ZXDiagram) = zxd.outputs

"""
    scalar(zxd)

Returns the scalar of `zxd`.
"""
scalar(zxd::ZXDiagram) = zxd.scalar

function add_global_phase!(zxd::ZXDiagram{T, P}, p::P) where {T, P}
    add_phase!(zxd.scalar, p)
    return zxd
end

function add_power!(zxd::ZXDiagram, n)
    add_power!(zxd.scalar, n)
    return zxd
end

function spider_sequence(zxd::ZXDiagram{T, P}) where {T, P}
    seq = []
    generate_layout!(zxd, seq)
    return seq
end

function generate_layout!(zxd::ZXDiagram{T, P}, seq::Vector{Any} = []) where {T, P}
    layout = zxd.layout
    nbits = length(zxd.inputs)
    vs_frontier = copy(zxd.inputs)
    vs_generated = Set(vs_frontier)
    frontier_col = [1//1 for _ = 1:nbits]
    frontier_active = [true for _ = 1:nbits]
    for i = 1:nbits
        set_qubit!(layout, vs_frontier[i], i)
        set_column!(layout, vs_frontier[i], 1//1)
    end

    while !(zxd.outputs ⊆ vs_frontier)
        while any(frontier_active)
            for q in 1:nbits
                if frontier_active[q]
                    v = vs_frontier[q]
                    nb = neighbors(zxd, v)
                    if length(nb) <= 2
                        set_loc!(layout, v, q, frontier_col[q])
                        push!(seq, v)
                        push!(vs_generated, v)
                        q_active = false
                        for v1 in nb
                            if !(v1 in vs_generated)
                                vs_frontier[q] = v1
                                frontier_col[q] += 1
                                q_active = true
                                break
                            end
                        end
                        frontier_active[q] = q_active
                    else
                        frontier_active[q] = false
                    end
                end
            end
        end
        for q = 1:nbits
            v = vs_frontier[q]
            nb = neighbors(zxd, v)
            isupdated = false
            for v1 in nb
                if !(v1 in vs_generated)
                    q1 = findfirst(isequal(v1), vs_frontier)
                    if q1 !== nothing
                        col = maximum(frontier_col[min(q, q1):max(q, q1)])
                        set_loc!(layout, v, q, col)
                        set_loc!(layout, v1, q1, col)
                        push!(vs_generated, v, v1)
                        push!(seq, (v, v1))
                        nb_v1 = neighbors(zxd, v1)
                        new_v1 = nb_v1[findfirst(v -> !(v in vs_generated), nb_v1)]
                        new_v = nb[findfirst(v -> !(v in vs_generated), nb)]
                        vs_frontier[q] = new_v
                        vs_frontier[q1] = new_v1
                        for i in min(q, q1):max(q, q1)
                            frontier_col[i] = col + 1
                        end
                        frontier_active[q] = true
                        frontier_active[q1] = true
                        isupdated = true
                        break
                    elseif spider_type(zxd, v1) == SpiderType.H && degree(zxd, v1) == 2
                        nb_v1 = neighbors(zxd, v1)
                        v2 = nb_v1[findfirst(!isequal(v), nb_v1)]
                        q2 = findfirst(isequal(v2), vs_frontier)
                        if q2 !== nothing
                            col = maximum(frontier_col[min(q, q2):max(q, q2)])
                            set_loc!(layout, v, q, col)
                            set_loc!(layout, v2, q2, col)
                            q1 = (q + q2)//2
                            denominator(q1) == 1 && (q1 += 1//2)
                            set_loc!(layout, v1, q1, col)
                            push!(vs_generated, v, v1, v2)
                            push!(seq, (v, v1, v2))
                            nb_v2 = neighbors(zxd, v2)
                            new_v = nb[findfirst(v -> !(v in vs_generated), nb)]
                            new_v2 = nb_v2[findfirst(v -> !(v in vs_generated), nb_v2)]
                            vs_frontier[q] = new_v
                            vs_frontier[q2] = new_v2
                            for i in min(q, q2):max(q, q2)
                                frontier_col[i] = col + 1
                            end
                            frontier_active[q] = true
                            frontier_active[q2] = true
                            isupdated = true
                            break
                        end
                    end
                end
                isupdated && break
            end
        end
    end
    for q = 1:length(zxd.outputs)
        set_loc!(layout, zxd.outputs[q], q, maximum(frontier_col))
    end
    return layout
end

"""
    continued_fraction(ϕ, n::Int) -> Rational

Obtain `s` and `r` from `ϕ` that satisfies `|s//r - ϕ| ≦ 1/2r²`
"""
function continued_fraction(fl, n::Int)
    if n == 1 || abs(mod(fl, 1)) < 1e-10
        Rational(floor(Int, fl), 1)
    else
        floor(Int, fl) + 1//continued_fraction(1/mod(fl, 1), n-1)
    end
end

safe_convert(::Type{T}, x) where T = convert(T, x)
safe_convert(::Type{T}, x::T) where T<:Rational = x
function safe_convert(::Type{T}, x::Real) where T<:Rational
    local fr
    for n=1:16 # at most 20 steps, otherwise the number may overflow.
        fr = continued_fraction(x, n)
        abs(fr - x) < 1e-12 && return fr
    end
    @warn "converting phase to rational, but with rounding error $(abs(fr-x))."
    return fr
end


text:
module TestBaseline

using Test
using TSML

function test_baseline()
    Random.seed!(123)
    iris=getiris()
    instances=iris[:,1:4]
    labels=iris[:,5] |> collect
    bl = Baseline()
    fit!(bl,instances,labels)
    @test bl.model == "setosa"
    @test sum(transform!(bl,instances) .== repeat(["setosa"],nrow(iris))) == nrow(iris)
    idy = Identity()
    fit!(idy,instances,labels)
    @test idy.model == nothing
    @test (transform!(idy,instances) .== instances) |> Matrix |> sum == 150*4
    @test idy.args == Dict()
end
@testset "Baseline Tests" begin
  test_baseline()
end

end


text:
using RRModels

using Test

include("test_LinearBucket.jl")


text:
using ComponentArrays
using OrdinaryDiffEq
using LinearAlgebra
using Unitful

function newton(du, u, p, t)
    mu = 398600.4418u"km^3/s^2"
    r = norm(u.r)
    du.r = u.v
    du.v = -mu .* u.r / r^3
end

r0 = [1131.340, -2282.343, 6672.423]u"km"
v0 = [-5.64305, 4.30333, 2.42879]u"km/s"
Δt = 86400.0*365u"s"
rv0 = ComponentArray(r=r0, v=v0)

prob = ODEProblem(newton, rv0, (0.0u"s", Δt))
@time sol = solve(prob, Vern8(), dt=100u"s", adaptive=false)
#sol |> display
println()

@time sol2 = solve(prob, Tsit5(), dt=100u"s",adaptive=false)

sol2.u[1].r |> display
sol2.u[1].v |> display