text:
Base.run(`clear`)
println("Starting")
println("----------------------------")
using Leapfrog
println("JIT compilation")
mp, cost = init_all(1);
stags, ess, number_of_equilibria = solve(cost, mp);
println("----------------------------")
println("Checking # of equilibria")
println(number_of_equilibria)
println("----------------------------")
println("Timing code")
times = [@elapsed solve(cost,mp) for i = 1:10];
times_mean = mean(times)
times_min = minimum(times)
times_max = maximum(times)
times_var = std(times)
println("----------------------------")
@printf("Mean    Min    Max    Std\n")
@printf("%4.4f  %4.4f %4.4f %4.4f\n", times_mean, times_min, times_max, times_var)
println("----------------------------")


text:
using ReinforcementLearningCore

using AbstractTrees
using TensorBoardLogger: TBLogger

RLCore.is_expand(::TBLogger) = false

ext: jl


text:
function is_ci()
    get(ENV, "TRAVIS", "") == "true" ||
    get(ENV, "APPVEYOR", "") == "true" ||
    get(ENV, "CI", "") == "true"
end

using GLAbstraction, GeometryTypes, ModernGL, FileIO, GLWindow
using ColorTypes
using Base.Test
import GLAbstraction: N0f8

include("macro_test.jl")

if !is_ci() # only do test if not CI... this is for automated testing environments which fail for OpenGL stuff, but I'd like to test if at least including works

window = create_glcontext("test", resolution=(500,500))

include("accessors.jl")
include("uniforms.jl")
include("texture.jl")
include("macro_test.jl")

# Test for creating a Buffer with a 1D Julia Array
# You need to supply the cardinality, as it can't be inferred
# indexbuffer is a shortcut for Buffer(GLUInt[0,1,2,2,3,0], 1, buffertype = GL_ELEMENT_ARRAY_BUFFER)
indexes = indexbuffer(GLuint[0,1,2])
# Test for creating a Buffer with a 1D Julia Array of Vectors
#v = Vec2f[Vec2f(0.0, 0.5), Vec2f(0.5, -0.5), Vec2f(-0.5,-0.5)]

v = Vec2f0[Vec2f0(0.0, 0.5), Vec2f0(0.5, -0.5), Vec2f0(-0.5,-0.5)]

verts = Buffer(v)
@test size(verts, 1) == 3
@test size(verts, 2) == 1

# lets define some uniforms
# uniforms are shader variables, which are supposed to stay the same for an entire draw call
cd(dirname(@__FILE__))
const triangle = std_renderobject(
    Dict(
        :vertex => verts,
        :name_doesnt_matter_for_indexes => indexes
    ),
    LazyShader("test.vert", "test.frag")
)

glClearColor(0,0,0,1)
i = 1
while isopen(window) && i < 20
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    GLAbstraction.render(triangle)
    swapbuffers(window)
    poll_glfw()
    sleep(0.01)
    i += 1
end
GLFW.DestroyWindow(window)

end


text:
include("equilibrium.jl")
using ImpvolEquilibrium
include("calibrate_params.jl")
using .CalibrateParameters
using Base.Test
using Logging
Logging.configure(level=INFO)

function test_data!(parameters, data)
    N, J, T, S = parameters[:N], parameters[:J], parameters[:T], parameters[:S]
    parameters[:beta_j] = zeros(1,1,J,1)
    parameters[:beta_j][1,1,:,1] = 0.5

    parameters[:gamma_jk] = [0.25 0.125 0.125; 0.125 0.25 0.125; 0.125 0.125 0.25]
    @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

    parameters[:S_nt] = zeros(1,N,1,T)
    parameters[:S_nt_data] = zeros(parameters[:S_nt])
    @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

    parameters[:d] = zeros(N,N,J,T)
    parameters[:d][:,:,1,1] = [0.9 0.1; 0.3 0.7]
    parameters[:d][:,:,2,1] = [0.7 0.3; 0.1 0.9]
    parameters[:d][:,:,3,1] = [1.0 0.0; 0.0 1.0]
    @test size(parameters[:d]) == (N,N,J,T)
    @test sum(parameters[:d], 2) ≈ ones(N,1,J,T) atol=1e-9

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

    final_expenditure_shares = zeros(1,1,J,1)
    final_expenditure_shares[1,1,:,1] = [0.3 0.1 0.6]
    @test sum(final_expenditure_shares, 3) ≈ ones(1,1,1,1) atol=1e-9

    # broad country weights for final expenditure
    country_weights = ones(1,N,1,1)
    country_weights = country_weights ./ sum(country_weights, 2)

    CalibrateParameters.calculate_p_and_nu!(parameters, data, final_expenditure_shares, country_weights)
    @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

    display(parameters[:p_sectoral][1,:,:,1])
    parameters[:w_njt] = CalibrateParameters.calculate_nominal_wages(parameters, data)

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

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

    parameters[:A] = CalibrateParameters.calculate_A(parameters, data)
    @test size(parameters[:A]) == (1, N, J, T)
    display(parameters[:A][1,:,:,1])

    info("US wage rate: ", sum(data["va"], 3)[1,end,1,1])

    # total world expenditure in the data - needed to get reasonable starting values
    parameters[:nominal_world_expenditure] = sum(data["va"] ./ parameters[:beta_j], (1,2,3))
    @test size(parameters[:nominal_world_expenditure]) == (1,1,1,T)

    # global, all-time average of sector final expenditure shares
    importance_weight = mean(parameters[:nu_njt], (1, 2, 4))
    parameters[:importance_weight] = importance_weight
    parameters[:A_njs] = Array{Array{Float64, 4}}(T)
    parameters[:A_njs][1] = parameters[:A] .* exp.(0.2*randn(1,N,J,S))
    @test size(parameters[:A_njs]) == (T,)
    @test size(parameters[:A_njs][1]) == (1,N,J,S)

end

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 init_data(parameters)
    N, J, T = parameters[:N], parameters[:J], parameters[:T]
    data = Dict{String,Any}()

    data["p_sectoral_data"] = ones(1,N,J,T)
    data["pwt"] = ones(1,N,1,T)
    data["va"] = ones(1,N,J,T)

    return data
end

srand(7094)

parameters = init_parameters()
data = init_data(parameters)
test_data!(parameters, data)

results = ImpvolEquilibrium.period_wrapper(parameters, 1)


text:
include_dir = joinpath(homedir(), "mbedtls-2.1.0/include/mbedtls")
headers = filter(x->splitext(x)[2]==".h", readdir(include_dir))
r=r"#define (\w+)\s+(-?(0x)?[\dA-F]+)"

c = Dict{ASCIIString, Cint}()

for header in headers
    s = readstring(joinpath(include_dir, header))
    for m in eachmatch(r, s)
        c[m.captures[1]] = parse(Cint, m.captures[2])
    end
end

outfile = joinpath(dirname(@__FILE__), "constants.jl")

open(outfile, "w") do f
    today = Dates.today()
    date = join([Dates.day(today), Dates.month(today), Dates.year(today)], "/")
    println(f, "# Auto-generated by gen_constants.jl on $date")
    for (key,val) in c
        println(f, "const $key=$val")
    end
end


text:
using Test
using ModelingToolkit, OrdinaryDiffEq

# Basic electric components
const t = Sym{ModelingToolkit.Parameter{Real}}(:t)
function Pin(;name)
    @variables v(t) i(t)
    ODESystem(Equation[], t, [v, i], [], name=name, defaults=[v=>1.0, i=>1.0])
end

function Ground(name)
    @named g = Pin()
    eqs = [g.v ~ 0]
    ODESystem(eqs, t, [], [], systems=[g], name=name)
end

function ConstantVoltage(name; V = 1.0)
    val = V
    @named p = Pin()
    @named n = Pin()
    @parameters V
    eqs = [
           V ~ p.v - n.v
           0 ~ p.i + n.i
          ]
    ODESystem(eqs, t, [], [V], systems=[p, n], defaults=Dict(V => val), name=name)
end

function Resistor(name; R = 1.0)
    val = R
    @named p = Pin()
    @named n = Pin()
    @variables v(t)
    @parameters R
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           v ~ p.i * R
          ]
    ODESystem(eqs, t, [v], [R], systems=[p, n], defaults=Dict(R => val), name=name)
end

function Capacitor(name; C = 1.0)
    val = C
    @named p = Pin()
    @named n = Pin()
    @variables v(t)
    @parameters C
    D = Differential(t)
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           D(v) ~ p.i / C
          ]
    ODESystem(eqs, t, [v], [C], systems=[p, n], defaults=Dict(C => val), name=name)
end

R = 1.0
C = 1.0
V = 1.0
resistor = Resistor(:resistor, R=R)
capacitor = Capacitor(:capacitor, C=C)
source = ConstantVoltage(:source, V=V)
ground = Ground(:ground)

function connect(ps...)
    eqs = [
           0 ~ sum(p->p.i, ps) # KCL
          ]
    # KVL
    for i in 1:length(ps)-1
        push!(eqs, ps[i].v ~ ps[i+1].v)
    end

    return eqs
end
rc_eqs = [
          connect(source.p, resistor.p)
          connect(resistor.n, capacitor.p)
          connect(capacitor.n, source.n, ground.g)
         ]

rc_model = ODESystem(rc_eqs, t, systems=[resistor, capacitor, source, ground], name=:rc)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
Return the LOCUS entry of the header.
"""
function parseheader(header::String)
    lines = split(header, "\n")
	firstline = lines[findfirst(line -> occursin("LOCUS", line), lines)]
    locus = split(firstline, r" +")[2]
    return locus
end

"""
Parse lines encoding genomic position, returning the feature as a `Symbol`, and an instance of `Locus`.
"""
function parseposition(line::String)
    feature, posstring = split(strip(line), r" +")
    if occursin(r"(\.\.|\^)", posstring)
        position = UnitRange(parse.(Int, filter.(c->isnumeric(c), split(posstring, r"(\.\.|\^)(.*(\.\.|\^))?")))...)
        strand = occursin("complement", posstring) ? '-' : '+'
    else
        position = parse(Int, posstring)
        position = position:position
        strand = '.'
    end
    complete_left = !occursin('<', posstring)
    complete_right = !occursin('>', posstring)
    order = Vector{UnitRange{Int}}()
	join = occursin("join", posstring)
    if join || occursin("order", posstring)
        for m in eachmatch(r"\d+(\.\.|\^)\d+", posstring)
            r = Meta.parse.(split(m.match, r"(\.\.|\^)"))
            push!(order, r[1]:r[2])
        end
    end
    return Symbol(feature), Locus(position, strand, complete_left, complete_right, order, join)
end

"""
Parse footer (sequence) portion of a GenBank file, returning a `String`.
"""
function filterseq(io::IOBuffer)
    line = String(take!(io))
    String(replace(filter(x -> !(isspace(x) || isnumeric(x)), line), r"[^ATGCNatgcn]" => "n"))
end

"""
    parsechromosome(lines)

Parse and return one chromosome entry, and the line number that it ends at.
"""
function parsechromosome(lines, G::Type = Gene)
    iobuffer = IOBuffer()
	isheader = true
    isfooter = false
    spanning = false
    position_spanning = false
    qualifier = String("")
    content = String("")
    header = ""

    feature = :source
    locus = Locus()

    chromosome = Chromosome{G}()

    linecount = 0
    for line in lines
        linecount += 1

		if length(line) == 0
			continue
		end

        # Catch cases where there's no header
        if linecount == 1 && occursin(r"gene", line)
            isheader = false
        end

        ### HEADER
		if isheader && occursin(r"FEATURES", line)
            chromosome.header = String(take!(iobuffer))
			isheader = false

        elseif isheader
            linecount == 1 ? print(iobuffer, line) : print(iobuffer, '\n', line)

        # Check if the footer has been reached
		elseif !isheader && !isfooter && (occursin(r"^BASE COUNT", line) || occursin(r"ORIGIN", line))
			# Stop parsing the file when the list of genes is over
            isfooter = true
			iobuffer = IOBuffer()

        ### BODY
		elseif !isheader && !isfooter
            if position_spanning && occursin(r"  /", line)
                position_spanning = false
                spanningline = filter(x -> x != ' ', String(take!(iobuffer)))
                try
                    feature, locus = parseposition(spanningline)
                catch
                    println(spanningline)
                    println(line)
                    @error "parseposition(spanningline) failed at line $linecount"
                end
            elseif position_spanning
                print(iobuffer, line)
            end
            if occursin(r"^ {5}\S", line)
                spanning = false
                try
                    feature, locus = parseposition(line)
                catch
                    println(line)
                    @error "parseposition(line) failed at line $linecount"
                end
                addgene!(chromosome, feature, locus)
            elseif !spanning && occursin(r"^ +/", line)
                if occursin(r"=", line)
                    if occursin("=\"", line)
                        (qualifier, content) = match(r"^ +/([^=]+)=\"?([^\"]*)\"?$", line).captures
                        content = String(content)
                    else
                        (qualifier, content) = match(r"^ +/(\S+)=(\S+)$", line).captures
                        try
                            tmpcontent = Meta.parse(content)
                            tmpcontent isa Expr && throw(Meta.ParseError)
							content = tmpcontent
                        catch
                            content = Symbol(content)
                        end
                    end

                    if occursin(r"\".*[^\"]$", line)
                        spanning = true
                    end

					isempty(names(chromosome.genedata)) ?
						(chromosome.genedata[!, Symbol(qualifier)] = Union{Missing, typeof(content)}[content]) :
	                    pushproperty!(chromosome.genes[end], Symbol(qualifier), content)

                else
                    # Qualifiers without a value assigned to them end up here
                    qualifier = split(line, '/')[2]
                    pushproperty!(chromosome.genes[end], Symbol(qualifier), true)
                end
            elseif spanning
                try
                    content = match(r" {21}([^\"]*)\"?$", line)[1]
                catch
                    @warn "Couldn't read content (line $linecount)"
                end
                if line[end] == '"'
                    spanning = false
                end
                if eltype(chromosome.genedata[!, Symbol(qualifier)]).b <: AbstractArray
                    i = index(chromosome.genes[end])
                    chromosome.genedata[!, Symbol(qualifier)][end][end] = Base.getproperty(chromosome.genes[end], Symbol(qualifier))[end] * "\n" * content
                else
                    Base.setproperty!(chromosome.genes[end], Symbol(qualifier), Base.getproperty(chromosome.genes[end], Symbol(qualifier)) * "\n" * content)
                end
            end

        ### FOOTER
        elseif isfooter
            if line == "//"
                break
            end
            if occursin(r"^ ", line)
                print(iobuffer, line)
            end
        end
    end
    chromosome.name = parseheader(chromosome.header)
    chromosome.sequence = LongDNASeq(filterseq(iobuffer))
    return linecount, chromosome
end

"""
    readgbk(input, G::Type = Gene; gunzip = false)

Parse GenBank-formatted file, returning a `Vector{Chromosome}`. `input` can be a file path or an `IOStream`. File names ending in ".gz" are assumed to be gzipped and are decompressed. Setting `gunzip` to `true` forces this behaviour.
The type of `AbstractGene` to be used can be specified with `G`, though currently the only option is `Gene`.
"""
function readgbk(filename::AbstractString, G::Type = Gene; gunzip = false)
    gz = filename[end-2:end] == ".gz"
    if gz || gunzip
        GZip.open(f -> readgbk(f, G), filename)
    else
        open(f -> readgbk(f, G), filename)
    end
end

function readgbk(input::IO, G::Type = Gene; gunzip = false)
	finished = false
	chrs = Chromosome{G}[]
	if gunzip
		lines = readlines(gzdopen(input))
	else
		lines = readlines(input)
	end
	currentline = 1
	while !finished
		if currentline >= length(lines)
			break
		end
		i, chr = parsechromosome(lines[currentline:end], G)
		currentline += i
		push!(chrs, chr)
	end
	chrs
end

## Printing

"""
    printgbk(io::IO, chr)
    printgbk(path::AbstractString, chr)

Print `chr` in GenBank format.
"""
function printgbk(path::AbstractString, chrs)
	open(path, "w") do f
		printgbk(f, chrs)
	end
end
function printgbk(chrs::AbstractVector{C}) where {C <: Chromosome}
    io = IOBuffer()
    printgbk(io, chrs)
end
function printgbk(io::IO, chrs::AbstractVector{C}) where {C <: Chromosome}
    for chr in chrs
        printgbk(io, chr)
    end
    return io
end
function printgbk(chr::C) where {C <: Chromosome}
    io = IOBuffer()
    printgbk(io, chr)
end
function printgbk(io::IO, chr::C) where {C <: Chromosome}
	if chr.header[1] != '#'
	    println(io, chr.header)
	else
		### GFF3 header
		println(io, "LOCUS       unknown", lpad(length(chr.sequence), 10, ' '), " bp     DNA")
	end
    println(io, rpad("FEATURES", 21, ' '), "Location/Qualifiers")
    println(io, chr.genes)
    println(io, "ORIGIN")
    formatsequence(chr.sequence, io)
    println(io)
    println(io, "//")
    return io
end

function formatsequence(sequence, io = IOBuffer)
    p = length(string(length(sequence))) + 2
    if length(sequence) > 60
        intervals = [i:i+60 for i in range(1; step = 60, stop = length(sequence)-60)]
        for interval in intervals
            println(io, lpad(string(first(interval)), p, ' '), " ", sequence[interval[1:10]],
                " ", sequence[interval[11:20]], " ", sequence[interval[21:30]],
                " ", sequence[interval[31:40]], " ", sequence[interval[41:50]],
                " ", sequence[interval[51:60]])
        end
    else
        intervals = [1:1]
    end
    i = intervals[end].stop
    if i <= length(sequence)
        print(io, lpad(i, p, ' '), " ")
        j = 0
        while i+j <= length(sequence)
            print(io, sequence[i+j])
            (j+1) % 10 == 0 && print(io, " ")
            j += 1
        end
    end
    return io
end


text:
function initialize_closest_weight_vectors!(parameters)
    distances = zeros(parameters.N, parameters.N)
    λ = parameters.λ
    for i in 1:parameters.N
        for j in (i+1):parameters.N
            distances[i, j] = norm(λ[i] - λ[j])
            distances[j, i] = distances[i, j]
        end
        I = sortperm(distances[i, :])
        push!(parameters.B, I[2:parameters.T+1])
    end
end

function update_reference_point!(z::Vector{Float64}, F::Vector{Float64})
    for i in 1:length(z)
        if z[i] > F[i]
            z[i] = F[i]
        end
    end
end

function update_reference_point!(z::Vector{Float64}, sol::xFgh_indiv)
    update_reference_point!(z, sol.f)
end

function update_reference_point!(z::Vector{Float64}, population)
    for sol in population
        update_reference_point!(z, sol)
    end
end

function g(fx, λ, z)
    m = -Inf
    for i in eachindex(z)
        @inbounds v = λ[i] * abs(fx[i] - z[i])
        if m < v
            m = v
        end
    end

    return m
end


text: has442A() = true


text:
using Test
using TypeClasses
using DataTypesBasic
using Suppressor
using Documenter

# Test Utils
splitln(str) = split(strip(str), "\n")

if v"1.6" <= VERSION < v"1.7"
  # somehow only Julia 1.6 does this correctly
  @test isempty(detect_ambiguities(TypeClasses))
  # doctests are super instable, hence we only do it for a specific Julia Version
  doctest(TypeClasses)
end

@testset "TypeClasses" begin
  @testset "MonoidAlternative" begin
    include("TypeClasses/MonoidAlternative.jl")
  end
  @testset "FunctorApplicativeMonad" begin
    include("TypeClasses/FunctorApplicativeMonad.jl")
  end
  @testset "FlipTypes" begin
    include("TypeClasses/FlipTypes.jl")
  end
end

@testset "TypeInstances" begin
  @testset "AbstractDictionary" begin
    include("TypeInstances/AbstractDictionary.jl")
  end
  @testset "AbstractVector" begin
    include("TypeInstances/AbstractVector.jl")
  end
  @testset "Callable" begin
    include("TypeInstances/Callable.jl")
  end
  @testset "Dict" begin
    include("TypeInstances/Dict.jl")
  end
  @testset "Future" begin
    include("TypeInstances/Future.jl")
  end
  @testset "Iterable" begin
    include("TypeInstances/Iterable.jl")
  end
  @testset "Pair" begin
    include("TypeInstances/Pair.jl")
  end
  @testset "State" begin
    include("TypeInstances/State.jl")
  end
  @testset "Task" begin
    include("TypeInstances/Task.jl")
  end
  @testset "Tuple" begin
    include("TypeInstances/Tuple.jl")
  end
  @testset "Writer" begin
    include("TypeInstances/Writer.jl")
  end
end

@testset "TypeInstances_DataTypesBasic" begin
  @testset "Const" begin
    include("TypeInstances_DataTypesBasic/Const.jl")
  end
  @testset "Identity" begin
    include("TypeInstances_DataTypesBasic/Identity.jl")
  end
  @testset "Either" begin
    include("TypeInstances_DataTypesBasic/Either.jl")
  end
  @testset "Option" begin
    include("TypeInstances_DataTypesBasic/Option.jl")
  end
  @testset "Try" begin
    include("TypeInstances_DataTypesBasic/Try.jl")
  end
  @testset "ContextManager" begin
    include("TypeInstances_DataTypesBasic/ContextManager.jl")
  end
end


text:
# setup callback functions for Ipopt
function eval_f_par(pcof::Vector{Float64}, params:: Juqbox.objparams, wa::Working_Arrays)
    #@show(pcof)
    f =Juqbox.traceobjgrad(pcof,params,wa,false,false)
    # @show(f)
    return f[1]
  end

function eval_g(pcof)
    return 0.0
end

function eval_grad_f_par(pcof::Vector{Float64}, grad_f::Vector{Float64}, params:: Juqbox.objparams, wa::Working_Arrays)
    objf, Gtemp, primaryobjf, secondaryobjf, traceinfid = Juqbox.traceobjgrad(pcof,params,wa,false, true)

    Gtemp = vcat(Gtemp...)
    for i in 1:length(Gtemp)
        grad_f[i] = Gtemp[i]
    end
    # remember the value of the primary obj func (for termination in intermediate_par)
    params.lastTraceInfidelity = traceinfid
    params.lastLeakIntegral = secondaryobjf

    # Save intermediate parameter vectors
    if params.save_pcof_hist
        push!(params.pcof_hist, copy(pcof)) #pcof_hist is and Array of Vector{Float64}
    end
end

function eval_jac_g(
    x::Vector{Float64},         # Current solution
    mode,                       # Either :Structure or :Values
    rows::Vector{Int32},        # Sparsity structure - row indices
    cols::Vector{Int32},        # Sparsity structure - column indices
    values::Vector{Float64})    # The values of the Hessian

    if mode == :Structure
        # rows[...] = ...
        # ...
        # cols[...] = ...
    else
        # values[...] = ...
    end
end

function intermediate_par(
    alg_mod::Int,
    iter_count::Int,
    obj_value::Float64,
    inf_pr::Float64, inf_du::Float64,
    mu::Float64, d_norm::Float64,
    regularization_size::Float64,
    alpha_du::Float64, alpha_pr::Float64,
    ls_trials::Int,
    params:: Juqbox.objparams)
  # ...
    if params.saveConvHist
        push!(params.objHist, obj_value)
        push!(params.dualInfidelityHist, inf_du)
        push!(params.primaryHist, params.lastTraceInfidelity)
        push!(params.secondaryHist,  params.lastLeakIntegral)
    end
    if params.lastTraceInfidelity < params.traceInfidelityThreshold
        println("Stopping because trace infidelity = ", params.lastTraceInfidelity,
                " < threshold = ", params.traceInfidelityThreshold)
        return false
    else
        return true  # Keep going
    end
end

"""
    prob = setup_ipopt_problem(params, wa, nCoeff, minCoeff, maxCoeff [, maxIter=50,
                            lbfgsMax=10, startFromScratch=true, ipTol=1e-5, acceptTol=1e-5, acceptIter=15])

Setup structure containing callback functions and convergence criteria for
optimization via IPOPT.

# Arguments
- `params:: objparams`: Struct with problem definition
- `wa::Working_Arrays`: Struct containing preallocated working arrays
- `nCoeff:: Int64`: Number of parameters in optimization
- `minCoeff:: Array{Float64, 1}`: Minimum allowable value for each parameter
- `maxCoeff:: Array{Float64, 1}`: Maximum allowable value for each parameter
- `maxIter:: Int64`: Maximum number of iterations to be taken by optimizer
- `lbfgsMax:: Int64`: Maximum number of past iterates for Hessian approximation by L-BFGS
- `startFromScratch:: Bool`: Specify whether the optimization is starting from file or not
- `ipTol:: Float64`: Desired convergence tolerance (relative)
- `acceptTol:: Float64`: Acceptable convergence tolerance (relative)
- `acceptIter:: Int64`: Number of acceptable iterates before triggering termination
"""
function setup_ipopt_problem(params:: Juqbox.objparams, wa::Working_Arrays, nCoeff:: Int64, minCoeff:: Array{Float64, 1}, maxCoeff:: Array{Float64, 1}, maxIter:: Int64=50, lbfgsMax:: Int64=10, startFromScratch:: Bool=true, ipTol:: Float64=1e-5, acceptTol:: Float64=1e-5, acceptIter:: Int64=15)
    # callback functions need access to the params object
    eval_f(pcof) = eval_f_par(pcof, params, wa)
    eval_grad_f(pcof, grad_f) = eval_grad_f_par(pcof, grad_f, params, wa)
    intermediate(alg_mod, iter_count, obj_value, inf_pr, inf_du, mu,
                 d_norm, regularization_size, alpha_du, alpha_pr, ls_trials) =
                     intermediate_par(alg_mod, iter_count, obj_value, inf_pr, inf_du, mu,
                                      d_norm, regularization_size, alpha_du, alpha_pr, ls_trials, params)
    # setup the Ipopt data structure
    mNLconstraints = 0;
    nEleJac = 0;
    nEleHess = 0;
    dum0 = zeros(0);
    dum1 = zeros(0);
    # tmp
    #    println("setup_ipopt_problem: nCoeff = ", nCoeff, " length(minCoeff) = ", length(minCoeff))
    prob = createProblem( nCoeff, minCoeff, maxCoeff, mNLconstraints, dum0, dum1, nEleJac, nEleHess, eval_f, eval_g, eval_grad_f, eval_jac_g);
    addOption( prob, "hessian_approximation", "limited-memory");
    addOption( prob, "limited_memory_max_history", lbfgsMax);
    addOption( prob, "max_iter", maxIter);
    addOption( prob, "tol", ipTol);
    addOption( prob, "acceptable_tol", acceptTol);
    addOption( prob, "acceptable_iter", acceptIter);
    if !startFromScratch # enable warm start of Ipopt
        addOption( prob, "warm_start_init_point", "yes")
        #        addOption( prob, "mu_init", 1e-6) # not sure how to set this parameter
        #        addOption( prob, "nlp_scaling_method", "none") # what about scaling?
        #
        # the following settings prevent the initial parameters to be pushed away from their bounds
        addOption( prob, "warm_start_bound_push", 1e-16)
        addOption( prob, "warm_start_bound_frac", 1e-16)
        addOption( prob, "warm_start_slack_bound_frac", 1e-16)
        addOption( prob, "warm_start_slack_bound_push", 1e-16)

        if !params.quiet
            println("Ipopt: Enabling warm start option")
        end
    end

    # intermediate callback function
    setIntermediateCallback(prob, intermediate)

# output some of the settings
    if !params.quiet
        println("Ipopt parameters: max # iterations = ", maxIter)
        println("Ipopt parameters: max history L-BFGS = ", lbfgsMax)
        println("Ipopt parameters: tol = ", ipTol)
        println("Ipopt parameters: atol = ", acceptTol)
        println("Ipopt parameters: accept # iter = ", acceptIter)
    end

    return prob
end

"""
    pcof = run_optimizer(prob, pcof0 [, baseName:: String=""])

Call IPOPT to  optimizize the control functions.

# Arguments
- `prob:: IpoptProblem`: Struct containing Ipopt parameters callback functions
- `pcof0:: Array{Float64, 1}`: Initial guess for the parameter values
- `baseName:: String`: Name of file for saving the optimized parameters; extension ".jld2" is appended
"""
function run_optimizer(prob:: IpoptProblem, pcof0:: Array{Float64, 1}, baseName:: String="")
    # takes at most max_iter >= 0 iterations. Set with addOption(prob, "max_iter", nIter)

    # initial guess for IPOPT
    prob.x = pcof0;

    # Ipopt solver
    println("*** Starting the optimization ***")
    @time solveProblem(prob);
    pcof = prob.x;

    #save the b-spline coeffs on a JLD2 file
    if length(baseName)>0
        fileName = baseName * ".jld2"
        save_pcof(fileName, pcof)
        println("Saved B-spline parameters on binary jld2-file '", fileName, "'");
    end

    return pcof

end # run_optimizer


text:
using Pkg.TOML, Libdl

config_toml = joinpath(first(DEPOT_PATH), "prefs", "PAPI.toml")
mkpath(dirname(config_toml))

if !isfile(config_toml)
    touch(config_toml)
end

config = TOML.parsefile(config_toml)
update_config = false

if haskey(ENV, "JULIA_PAPI_BINARY")
    config["binary"] = ENV["JULIA_PAPI_BINARY"]
    update_config = true
end

if haskey(ENV, "JULIA_PAPI_LIBRARY")
    config["library"] = ENV["JULIA_PAPI_LIBRARY"]
    update_config = true
end

if haskey(ENV, "JULIA_PAPI_PATH")
    config["path"] = ENV["JULIA_PAPI_PATH"]
    update_config = true
end

if update_config
    open(config_toml, write=true) do f
        TOML.print(f, config)
    end
end

binary = get(config, "binary", "")

# 2. generate deps.jl
deps = if binary == "system"
    @info "using system PAPI"
    library = get(config, "library", "libpapi")
    path    = get(config, "path", "")

    const libpapi = find_library(library, path == "" ? [] : [joinpath(path, "lib"), joinpath(path, "lib64")])
    if libpapi == ""
        error("libpapi could not be found")
    end

    papi_version = try
        papi_version_cmd = Cmd([path == "" ? "papi_version" : joinpath(path, "bin", "papi_version")])
        output = chomp(read(papi_version_cmd, String))
        startswith(output, "PAPI Version: ") || error("unexpected output from $papi_version")
        v = output[15:end]
        m = match(r"^(\d+)(?:\.(\d+))?(?:\.(\d+))?(?:\.(\d+))?$", v)
        m === nothing && throw(ArgumentError("invalid version string: $v"))
        major, minor, patch, _ = map(m.captures) do x
            x !== nothing ? parse(Int, x) : 0
        end
        VersionNumber(major, minor, patch)
    catch e
        error("failed to find papi version: $e")
    end

    quote
        const libpapi = $libpapi
        const papi_version = $papi_version
    end
else
    @warn "using PAPI_jll: the use of a system library is probably recommended"
    papi_version = v"6.0.0"

    quote
        using PAPI_jll
        const papi_version = $papi_version
    end
end

remove_line_numbers(x) = x
function remove_line_numbers(ex::Expr)
    if ex.head == :macrocall
        ex.args[2] = nothing
    else
        ex.args = [remove_line_numbers(arg) for arg in ex.args if !(arg isa LineNumberNode)]
    end
    return ex
end

# only update deps.jl if it has changed.
# allows users to call Pkg.build("PAPI") without triggering another round of precompilation
deps_str =
    """
    # This file has been generated automatically.
    # It will be overwritten the next time `Pkg.build("PAPI")` is called.
    """ *
    string(remove_line_numbers(deps))

if !isfile("deps.jl") || deps_str != read("deps.jl", String)
    write("deps.jl", deps_str)
end


text:
# This is a test file

#+MyTag1,MyTag2
#
# My documentedItem
#
function myFunction (;x::Float64=sin(1))
     return 0
end

# final comment (no return)


text:
pynulltuple(len) = pynew(errcheck(C.PyTuple_New(len)))

function pytuple_setitem(xs::Py, i, x)
    x_ = Py(x)
    err = C.PyTuple_SetItem(getptr(xs), i, getptr(x_))
    pystolen!(x_)
    errcheck(err)
    return xs
end

function pytuple_getitem(xs::Py, i)
    pynew(incref(errcheck(C.PyTuple_GetItem(getptr(xs), i))))
end

function pytuple_fromiter(xs)
    sz = Base.IteratorSize(typeof(xs))
    if sz isa Base.HasLength || sz isa Base.HasShape
        # length known, e.g. Tuple, Pair, Vector
        ans = pynulltuple(length(xs))
        for (i, x) in enumerate(xs)
            pytuple_setitem(ans, i-1, x)
        end
        return ans
    else
        # length unknown
        xs_ = pylist_fromiter(xs)
        ans = pylist_astuple(xs_)
        pydel!(xs_)
        return ans
    end
end

@generated function pytuple_fromiter(xs::Tuple)
    n = length(xs.parameters)
    code = []
    push!(code, :(ans = pynulltuple($n)))
    for i in 1:n
        push!(code, :(pytuple_setitem(ans, $(i-1), xs[$i])))
    end
    push!(code, :(return ans))
    return Expr(:block, code...)
end

"""
    pytuple(x=())

Convert `x` to a Python `tuple`.

If `x` is a Python object, this is equivalent to `tuple(x)` in Python.
Otherwise `x` must be iterable.
"""
pytuple() = pynulltuple(0)
pytuple(x) = ispy(x) ? pybuiltins.tuple(x) : pytuple_fromiter(x)
export pytuple

pyistuple(x) = pytypecheckfast(x, C.Py_TPFLAGS_TUPLE_SUBCLASS)


text:
using SossGen
using Test

@testset "SossGen.jl" begin
    # Write your own tests here.
end

max_stars_repo_name: cscherrer/SossGen

max_issues_repo_name: cscherrer/SossGen

max_forks_repo_name: cscherrer/SossGen


text:
## ------------------------------------------------------
# utils
red_resolution(v, n::Int) = v[begin:max(div(end, n), 1):end]
normalizeP!(P::AbstractVector; tol = 1e-5) = (Z = sum(P); (abs(Z - 1.0) > tol) && (P .= P ./ Z); P)

function chuncks(v, c::Int)
    vlen = length(v)
    clen = max(div(vlen, c) + 1, 1)
    chnks =
    idx0 = firstindex(v)
    while true
        idx1 = min(vlen, idx0 + clen)
        push!(chnks, view(v, idx0:idx1))
        (idx1 >= vlen) && return chnks
        idx0 = idx1 + 1
    end
end

function hist(Vi, P)
    hist = Dict{Float64, Float64}()
    for (v, p) in zip(Vi, P)
        get!(hist, v, 0.0)
        hist[v] += p
    end
    hist
end

function n_ave_conv!(v, n; aux = copy(v))
    idx0 = firstindex(v)
    idx1 = lastindex(v)
    for i in eachindex(v)
        ran = max(idx0, i - n):min(idx1, i + n)
        v[i] = sum(view(aux, ran)) / length(ran)
    end
    return v
end

function tail_ave(v, n)
    idx0 = firstindex(v)
    idx1 = lastindex(v)
    tail = view(v, max(idx0, idx1 - n):idx0)
    mean(tail)
end

function check_stst(w, th, vs...)
    for v in vs
        idx1 = lastindex(v)
        (idx1 <= w) && return false
        ran = (idx1 - w):idx1
        vi = view(v, ran)
        rel_var = abs(std(vi) / mean(vi))
        (rel_var > th) && return false
    end
    return true
end

ismul(v, m) = iszero(rem(v, m))

function findapproxi(v, r::AbstractRange)
    v0, v1 = first(r), last(r)
    (v0 > v1) && error("Only increasing ranges are allowed")
    idx = Int(div(v - v0, step(r), RoundNearest)) + 1
    clamp(idx, firstindex(r), lastindex(r))
end

findapproxv(v, r::AbstractRange) = r[findapproxi(v, r)]


text:
# this illustrates effect and detection of influential observations
using Plots, LinearAlgebra
n = 20
x = (1:n-1)/(n-1)
x = [x; 3] # the last observation is an outlying value of x

x = [ones(n,1) x]
P = x*inv(x'*x)*x'

beta = [10, -1]
e = 2*randn(n,1)
y = x*beta + e

# The fit
yhat = P*y

# calculate leverage and influence
leverage = diag(P)
e = y - yhat
influence = (leverage ./ (1.0 .-leverage)) .* e

xlabel = "X"
x = x[:,2]
scatter(x, y, label = "Data points")
plot!(x, yhat, label = "fitted")
plot!(x, leverage, label = "Leverage")
plot!(x, influence, label = "Influence")
gui()
#savefig("InfluentialObservation.svg")

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
# ---
# title: 1658. Minimum Operations to Reduce X to Zero
# id: problem1658
# author: AquaIndigo
# date: 2020-11-23
# difficulty: Medium
# categories: Two Pointers, Binary Search, Greedy
# link: <https://leetcode.com/problems/minimum-operations-to-reduce-x-to-zero/description/>
# hidden: true
# ---
#
# You are given an integer array `nums` and an integer `x`. In one operation,
# you can either remove the leftmost or the rightmost element from the array
# `nums` and subtract its value from `x`. Note that this **modifies** the array
# for future operations.
#
# Return _the **minimum number** of operations to reduce _`x` _to **exactly**_
# `0` _if it 's possible_ _, otherwise, return_`-1`.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,1,4,2,3], x = 5
#     Output: 2
#     Explanation: The optimal solution is to remove the last two elements to reduce x to zero.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [5,6,7,8,9], x = 4
#     Output: -1
#
#
# **Example 3:**
#
#
#
#     Input: nums = [3,2,20,1,1,3], x = 10
#     Output: 5
#     Explanation: The optimal solution is to remove the last three elements and the first two elements (5 operations in total) to reduce x to zero.
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 105`
#   * `1 <= nums[i] <= 104`
#   * `1 <= x <= 109`
#
#
## @lc code=start
using LeetCode

function min_operations(nums::Vector{Int}, x::Int)
    rnums = reverse(nums)
    if nums[1] > x && rnums[1] > x
        return -1
    end
    len = length(nums)
    rlm, llm, res = len, len, len
    for i in 1:(len - 1)
        if ((nums[i + 1] += nums[i]) > x)
            rlm = i
            break
        end
    end
    if rm == len
        if nums[end] < x
            return -1
        else
            return len
        end
    end
    for i in 1:(len - 1)
        if ((rnums[i + 1] += rnums[i]) > x)
            llm = i
            break
        end
    end
    i, j = 1, llm
    while i <= rlm && j > 0
        if nums[i] == x
            res = min(i, res)
            i += 1
        elseif rnums[j] == x
            res = min(j, res)
            j -= 1
        else
            sm = nums[i] + rnums[j]
            ## println(i," ", j," ", sm)
            if sm == x
                res = min(i + j, res)
                i += 1
                j -= 1
            elseif sm > x
                j -= 1
            else
                i += 1
            end
        end
    end
    return res
end
## @lc code=end


text:
using RCall, JLD

R"
library(gridExtra)
library(cowplot)
library(tidyverse)
library(latex2exp)
library(viridis)
"

DF =load("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_7_1.jld")
DF_PY = DF["DF_PY"]

R"
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
df_r= $DF_PY
df_r$beta_scaled= range01(df_r$beta_log)
df_r$color = rgb(df_r$beta_scaled, df_r$sigma,0.5,1)
#df_r$color = rgb(0, df_r$sigma,0.5,1)
#df_r$color_sigma = rgb(0, df_r$sigma,0.5,1)
df_r$color_sigma = rgb(df_r$beta_scaled, df_r$sigma,0.5,1)
df_r$color_beta = rgb(0, 0,0.5,1)
p <- df_r %>% ggplot(aes(x=Exp, y = std, group=sigma))  + geom_line(alpha = 0.8,linetype = 'longdash') + geom_point(aes(colour=color), size=4)
p_py <- p  + geom_line(data =df_r, aes(x=Exp, y = std, group=beta), alpha = 0.8, linetype='dotted') +
  xlim(0, 100)+ ylim(0,16)+ labs(y='Std', x='Expectation')+scale_colour_identity()+theme_classic()+
   theme(plot.title = element_text(hjust = 0.5,size = 15), axis.text.x = element_text(size=15), axis.text.y = element_text(size=15),legend.position='none')
"

R"p_py
pdf('/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/figures/figure_7_1.pdf', width= 4, height = 4)
plot(p_py)
dev.off()"

max_stars_repo_name: konkam/approximatingBNPpriors

max_issues_repo_name: konkam/approximatingBNPpriors

max_forks_repo_name: konkam/approximatingBNPpriors


text:
# The infite loop was patched with a killswitch after certain number of iterations.
# Now the obtained reconstructions are too good to be real and I need to check what
push!(LOAD_PATH, "../models")
push!(LOAD_PATH, "../")
push!(LOAD_PATH, "../SparseInverseProblems/src")
using SparseInverseProblemsMod
using SuperResModels
using TestCases
using Utils
using PyCall
@pyimport numpy as np
@pyimport matplotlib.pyplot as plt

## Parameters

# Static Parameters
f_c = 20
x_max = 1.0
num_x = 10*f_c
filter_x = ones(2*f_c+1)

# Dynamic Parameters
K  = 2
tau = 1.0/(K*1.0)
v_max = 0.5
num_v = 10

# Variability of weights
minWeights = 1
maxWeights = 10

# Variability of number of particles
min_number_part = 15
max_number_part = 15

# Number of generated examples
num_trials = 10

# test case
test_case = () -> TestCases.cloud_1d_full(x_max, v_max, minWeights, maxWeights, K, tau, rand(min_number_part:max_number_part))
# WARNING: changing anything in the test_case will generate troubles if the rejection_sampling algorithm
# in Utils.js is not either updated for it, or disabled.

# Noise parameters
noises_data = linspace(0,0.1,5)
noises_data = noises_data[2:end]
noises_position = linspace(0,0.01,5)
noises_position = noises_position[2:end]

# To discard faulty reconstructions
threshold_weight = 0.1

# Randomness set for replication
srand(1)

model_static = SuperResModels.Fourier1d(f_c, x_max, filter_x, num_x)
model_dynamic = SuperResModels.DynamicFourier1d(model_static, v_max, tau, K, num_v*K)

### Location of data folder
dataFolder = "data/1Dsimulations"

bins = np.load("../"*dataFolder*"/separationDistribBins.npy")
density = np.load("../"*dataFolder*"/separationDistribVal.npy")

### Number of points
@everywhere begin
    model_static = $model_static
    model_dynamic = $model_dynamic
    x_max = $x_max
    v_max = $v_max
    minWeights = $minWeights
    maxWeights = $maxWeights
    K = $K
    tau = $tau
    dataFolder = $dataFolder
    test_case = $test_case
    noises_data = $noises_data
    noises_position = $noises_position
    bins = $bins
    density = $density
end

println("### Beggining round ###")
    K = div(length(model_dynamic.times)-1,2)
    println("K = ", K)
    tau = maximum(model_dynamic.times)/K
    println("tau = ", tau)
    println(" Visualization of generated example" )
    (thetas, weights) =  Utils.Rejection_sampling(test_case, bins, density, K, tau, model_static.x_max)
    # Obtain the separation of the configuration of particles
    separation = Utils.separation_val(thetas,K,tau,model_static.x_max)
    # And the separationg given by the dynamic norm.
    separation_dyn = Utils.separation_norm(thetas,K,tau,model_static.x_max)
    sizeAmp = 50
#    plt.scatter(thetas[1,:], thetas[2,:], s = weights*sizeAmp, alpha = 0.8)
#    plt.xlim((0, x_max))
#    plt.ylim((-v_max, v_max))
#    plt.title(" Sep: "*string(separation)*", SepDyn: "*string(separation_dyn))
#    plt.show()
    d = dim(model_dynamic)
    println("dimension : ",d)
    results = zeros(1, 7)
    ## 	Noiseless case
    # Dynamic case
    # We generate the target measure
    target = Utils.generate_target(model_dynamic, thetas, weights)
    (results[1, 1], results[1, 2], results[1,3]) = Utils.generate_and_reconstruct_dynamic(model_dynamic, thetas, weights, target)
    println(" Results: dx = "*string(results[1,1])*", dv = "*string(results[1,2])*", dw = "*string(results[1,3]))
    # Inspecting the explicit reconstruction
    results2 = zeros(1,7)
    (thetas_est, weights_est) =  Utils.run_simulation_target(model_dynamic, thetas, weights, target)
    # Discard obviously defective reconstructions
    outputSize = length(thetas_est)
    thetas_est = thetas_est[:, weights_est .> threshold_weight]
    weights_est = weights_est[weights_est .> threshold_weight]
    println(" Discarded reconstructions = ", outputSize- length(thetas_est))
    if (length(thetas) == length(thetas_est))
        # The output of the minimization algorithm returns an unsorted vector with the
        # estimated positions ans weights of particles. match_points is a function that
        # returns a vector that tries to match the best possible the corresponding particles
        # to the estimated particles.
        corres = Utils.match_points(thetas, thetas_est)
        dist_x = norm(thetas[1:d, :] - thetas_est[1:d, corres], Inf)
        dist_v = norm(thetas[d+1:end, :] - thetas_est[d+1:end, corres], Inf)
        dist_w = norm(weights[:] - weights_est[corres], Inf)
        # returns the infinite distance between the matched particles, in space and velocity.
        (results2[1,1], results2[1,2], results2[1,3]) = (dist_x, dist_v, dist_w)
    else
         (results2[1,1], results2[1,2], results2[1,3]) = model_dynamic.static.x_max, model_dynamic.v_max, maximum(weights)
    end
    println(" Results2: dx = "*string(results[1,1])*", dv = "*string(results[1,2])*", dw = "*string(results[1,3]))
        # Plotting reconstructions constrasted to the real solutions.
    plt.scatter(thetas[1,:], thetas[2,:], s = weights*sizeAmp, alpha = 0.8)
    plt.scatter(thetas_est[1,:], thetas_est[2,:], s=weights_est*sizeAmp/2, alpha = 0.8)
    plt.xlim((0, x_max))
    plt.ylim((-v_max, v_max))
    plt.title(" Sep: "*string(separation)*", SepDyn: "*string(separation_dyn))
    plt.show()
	# Trying the static reconstructions
    (results[1,4], results[1,5], results[1,6], results[1,7]) = Utils.generate_and_reconstuct_static_best(model_dynamic, model_static, thetas, weights, target)
	# Checking deeply the static reconstructions.
    static_target = Utils.target_to_static(target_dynamic, model_dynamic.times)
        # Projecting the measurements to static
    k = 3
    thetas_t = Utils.to_static(thetas,model_dynamic.times[k],model_static.xmax)
        # Show the static measurements
    plt.scatter(thetas_t, zeros(length(thetas_t)), s = weights*sizeAmp, alpha = 0.8)
    plt.xlim((0,x_max))
    plt.ylim((-v_max,v_max))
    (thetas_t_est, weights_t_est) = Utils.run_simulation_target(model_static, thetas_t, weights, static_target[k])
    plt.scatter(thetas_t_est, zeros(length(thetas_t)), s = weights_t*sizeAmp/2, alpha=0.8)
    plt.show()


text:
using Simplices: subsample_coeffs
using Simplices: simplexintersection
using Simplices: intersectingvertices

function interiorpoints(s::SArray{D, T}, n::Int) where {D, T}
    R = rand(size(s, 2), n)
    convex_coeffs = ((1 ./ sum(R, dims = 1)) .* R)
    s * convex_coeffs
end

function interiorpoints(s, n::Int)
    R = rand(size(s, 2), n)
    convex_coeffs = ((1 ./ sum(R, dims = 1)) .* R)
    s * convex_coeffs
end

connectvertices(s) = hcat(s, hcat([s[:, i]
                    for i in [1, 3, 1, 4, 2]]...,))
splitaxes(x) = ([x[k, :] for k = 1:size(x, 1)]...,)
preparesimplex(s) = splitaxes(connectvertices(s))

"""
    getsimplex(r::Embeddings.AbstractEmbedding,
                t::DelaunayTriangulation,
                i::Int)

Get the vertices of the i-th simplex of the triangulation.
"""
getsimplex(E::Embeddings.AbstractEmbedding, DT::DelaunayTriangulation, i::Int) =
    E.points[:, DT.indices[:, i]]

getsimplex(d::Dataset, DT::DelaunayTriangulation, i::Int) =
    d.points[:, DT.indices[:, i]]

getsimplex(pts, DT::DelaunayTriangulation, i::Int) =
    pts[:, DT.indices[:, i]]

"""
    forwardmap(r::Embeddings.AbstractEmbedding,
                t::DelaunayTriangulation,
                i::Int)

Get the vertices of the i-th simplex of the triangulation projected
one step forward in time.
"""
forwardmap(r::Embeddings.AbstractEmbedding, t::DelaunayTriangulation, i::Int) =
    E.points[:, DT.indices[:, i] .+ 1]

forwardmap(r::Dataset, DT::DelaunayTriangulation, i::Int) =
    E.points[:, DT.indices[:, i] .+ 1]

forwardmap(pts, DT::DelaunayTriangulation, i::Int) =
    pts[:, DT.indices[:, i] .+ 1]

"""
Generate a total of `n_fillpoints` points lying inside the intersection
`sᵢ ∩ sⱼ` between simplices `sᵢ` and `sⱼ`.

This is done by first performing a Delaunay triangulation of the set of vertices
forming the intersecting volume. Then, points are inserted into the subsimplices
of the convex hull of the intersection.

If `sample_randomly = true`, then points are generated randomly inside each
subsimplex. If `sample_randomly = false` (default), then points are generated
as the centroids of the subsimplices resulting from a shape-preserving
refinement of the simplices in the convex hull.

## Arguments
- **``sᵢ``**: An array where each column is a simplex vertex.
- **``sⱼ``**: An array where each column is a simplex vertex.
- **``n_fillpoints``**: The total number of points to fill the convex
    hull of the intersection with. If the intersection is more complex than
    a single simplex, then the convex hull of the intersection is triangulated
    and `n_fillpoints` points are distributed evenly among the resulting
    subsimplices.
- **``sample_randomly``**: Should points be sampled randomly within the
    intersection? Default is `sample_randomly = false`, which inserts points
    inside the intersecting volume in a regular manner.
"""
function generate_fillpoints(sᵢ, sⱼ, n_fillpoints = 100, sample_randomly = true)

    intvol = simplexintersection(sᵢ, sⱼ)

    if intvol > 0 + 1e-9
        intersecting_vertices = intersectingvertices(sᵢ, sⱼ)
        if size(intersecting_vertices, 1) > 4
            triang = hcat(delaunay(intersecting_vertices)...,)
        else
            triang = zeros(Int, 4, 1)
            triang[:, 1] = 1:4
        end

        n_simplices = size(triang, 2)
        if n_simplices == 1
            npts_per_simplex = n_fillpoints
        elseif n_simplices > 1

            npts_per_simplex = ceil(Int, n_fillpoints/n_simplices)
        end
        C = transpose(subsample_coeffs(3, npts_per_simplex, sample_randomly))

        fillpts = Vector{Array{Float64, 2}}(undef, n_simplices)
        for k = 1:n_simplices
            fillpts[k] = transpose(C * intersecting_vertices[triang[:, k], :])
        end
        fillpts = hcat(fillpts...,)
        return intersecting_vertices, fillpts
    else
        return nothing
    end
end

export
interiorpoints,
connectvertices,
splitaxes,
preparesimplex,
getsimplex,
forwardmap,
generate_fillpoints


text:
# Temperature
# ===========

const pt = [
     9.86573e-01, 1.62228e-02, 1.55270e-02,-1.04323e-01,-3.75801e-03,
    -1.18538e-03,-1.24043e-01, 4.56820e-03, 8.76018e-03,-1.36235e-01,
    -3.52427e-02, 8.84181e-03,-5.92127e-03,-8.61650e+00, 0.00000e+00,
     1.28492e-02, 0.00000e+00, 1.30096e+02, 1.04567e-02, 1.65686e-03,
    -5.53887e-06, 2.97810e-03, 0.00000e+00, 5.13122e-03, 8.66784e-02,
     1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00,-7.27026e-06,
     0.00000e+00, 6.74494e+00, 4.93933e-03, 2.21656e-03, 2.50802e-03,
     0.00000e+00, 0.00000e+00,-2.08841e-02,-1.79873e+00, 1.45103e-03,
     2.81769e-04,-1.44703e-03,-5.16394e-05, 8.47001e-02, 1.70147e-01,
     5.72562e-03, 5.07493e-05, 4.36148e-03, 1.17863e-04, 4.74364e-03,
     6.61278e-03, 4.34292e-05, 1.44373e-03, 2.41470e-05, 2.84426e-03,
     8.56560e-04, 2.04028e-03, 0.00000e+00,-3.15994e+03,-2.46423e-03,
     1.13843e-03, 4.20512e-04, 0.00000e+00,-9.77214e+01, 6.77794e-03,
     5.27499e-03, 1.14936e-03, 0.00000e+00,-6.61311e-03,-1.84255e-02,
    -1.96259e-02, 2.98618e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
     6.44574e+02, 8.84668e-04, 5.05066e-04, 0.00000e+00, 4.02881e+03,
    -1.89503e-03, 0.00000e+00, 0.00000e+00, 8.21407e-04, 2.06780e-03,
     0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
    -1.20410e-02,-3.63963e-03, 9.92070e-05,-1.15284e-04,-6.33059e-05,
    -6.05545e-01, 8.34218e-03,-9.13036e+01, 3.71042e-04, 0.00000e+00,
     4.19000e-04, 2.70928e-03, 3.31507e-03,-4.44508e-03,-4.96334e-03,
    -1.60449e-03, 3.95119e-03, 2.48924e-03, 5.09815e-04, 4.05302e-03,
     2.24076e-03, 0.00000e+00, 6.84256e-03, 4.66354e-04, 0.00000e+00,
    -3.68328e-04, 0.00000e+00, 0.00000e+00,-1.46870e+02, 0.00000e+00,
     0.00000e+00, 1.09501e-03, 4.65156e-04, 5.62583e-04, 3.21596e+00,
     6.43168e-04, 3.14860e-03, 3.40738e-03, 1.78481e-03, 9.62532e-04,
     5.58171e-04, 3.43731e+00,-2.33195e-01, 5.10289e-04, 0.00000e+00,
     0.00000e+00,-9.25347e+04, 0.00000e+00,-1.99639e-03, 0.00000e+00,
     0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
     0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
]

# Density
# =======

# He
# --

const pd_He = [
     1.09979E+00,-4.88060E-02,-1.97501E-01,-9.10280E-02,-6.96558E-03,
     2.42136E-02, 3.91333E-01,-7.20068E-03,-3.22718E-02, 1.41508E+00,
     1.68194E-01, 1.85282E-02, 1.09384E-01,-7.24282E+00, 0.00000E+00,
     2.96377E-01,-4.97210E-02, 1.04114E+02,-8.61108E-02,-7.29177E-04,
     1.48998E-06, 1.08629E-03, 0.00000E+00, 0.00000E+00, 8.31090E-02,
     1.12818E-01,-5.75005E-02,-1.29919E-02,-1.78849E-02,-2.86343E-06,
     0.00000E+00,-1.51187E+02,-6.65902E-03, 0.00000E+00,-2.02069E-03,
     0.00000E+00, 0.00000E+00, 4.32264E-02,-2.80444E+01,-3.26789E-03,
     2.47461E-03, 0.00000E+00, 0.00000E+00, 9.82100E-02, 1.22714E-01,
    -3.96450E-02, 0.00000E+00,-2.76489E-03, 0.00000E+00, 1.87723E-03,
    -8.09813E-03, 4.34428E-05,-7.70932E-03, 0.00000E+00,-2.28894E-03,
    -5.69070E-03,-5.22193E-03, 6.00692E-03,-7.80434E+03,-3.48336E-03,
    -6.38362E-03,-1.82190E-03, 0.00000E+00,-7.58976E+01,-2.17875E-02,
    -1.72524E-02,-9.06287E-03, 0.00000E+00, 2.44725E-02, 8.66040E-02,
     1.05712E-01, 3.02543E+04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
    -6.01364E+03,-5.64668E-03,-2.54157E-03, 0.00000E+00, 3.15611E+02,
    -5.69158E-03, 0.00000E+00, 0.00000E+00,-4.47216E-03,-4.49523E-03,
     4.64428E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     4.51236E-02, 2.46520E-02, 6.17794E-03, 0.00000E+00, 0.00000E+00,
    -3.62944E-01,-4.80022E-02,-7.57230E+01,-1.99656E-03, 0.00000E+00,
    -5.18780E-03,-1.73990E-02,-9.03485E-03, 7.48465E-03, 1.53267E-02,
     1.06296E-02, 1.18655E-02, 2.55569E-03, 1.69020E-03, 3.51936E-02,
    -1.81242E-02, 0.00000E+00,-1.00529E-01,-5.10574E-03, 0.00000E+00,
     2.10228E-03, 0.00000E+00, 0.00000E+00,-1.73255E+02, 5.07833E-01,
    -2.41408E-01, 8.75414E-03, 2.77527E-03,-8.90353E-05,-5.25148E+00,
    -5.83899E-03,-2.09122E-02,-9.63530E-03, 9.77164E-03, 4.07051E-03,
     2.53555E-04,-5.52875E+00,-3.55993E-01,-2.49231E-03, 0.00000E+00,
     0.00000E+00, 2.86026E+01, 0.00000E+00, 3.42722E-04, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
]

# O
# -

const pd_O = [
     1.02315E+00,-1.59710E-01,-1.06630E-01,-1.77074E-02,-4.42726E-03,
     3.44803E-02, 4.45613E-02,-3.33751E-02,-5.73598E-02, 3.50360E-01,
     6.33053E-02, 2.16221E-02, 5.42577E-02,-5.74193E+00, 0.00000E+00,
     1.90891E-01,-1.39194E-02, 1.01102E+02, 8.16363E-02, 1.33717E-04,
     6.54403E-06, 3.10295E-03, 0.00000E+00, 0.00000E+00, 5.38205E-02,
     1.23910E-01,-1.39831E-02, 0.00000E+00, 0.00000E+00,-3.95915E-06,
     0.00000E+00,-7.14651E-01,-5.01027E-03, 0.00000E+00,-3.24756E-03,
     0.00000E+00, 0.00000E+00, 4.42173E-02,-1.31598E+01,-3.15626E-03,
     1.24574E-03,-1.47626E-03,-1.55461E-03, 6.40682E-02, 1.34898E-01,
    -2.42415E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 6.13666E-04,
    -5.40373E-03, 2.61635E-05,-3.33012E-03, 0.00000E+00,-3.08101E-03,
    -2.42679E-03,-3.36086E-03, 0.00000E+00,-1.18979E+03,-5.04738E-02,
    -2.61547E-03,-1.03132E-03, 1.91583E-04,-8.38132E+01,-1.40517E-02,
    -1.14167E-02,-4.08012E-03, 1.73522E-04,-1.39644E-02,-6.64128E-02,
    -6.85152E-02,-1.34414E+04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     6.07916E+02,-4.12220E-03,-2.20996E-03, 0.00000E+00, 1.70277E+03,
    -4.63015E-03, 0.00000E+00, 0.00000E+00,-2.25360E-03,-2.96204E-03,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     3.92786E-02, 1.31186E-02,-1.78086E-03, 0.00000E+00, 0.00000E+00,
    -3.90083E-01,-2.84741E-02,-7.78400E+01,-1.02601E-03, 0.00000E+00,
    -7.26485E-04,-5.42181E-03,-5.59305E-03, 1.22825E-02, 1.23868E-02,
     6.68835E-03,-1.03303E-02,-9.51903E-03, 2.70021E-04,-2.57084E-02,
    -1.32430E-02, 0.00000E+00,-3.81000E-02,-3.16810E-03, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00,-9.05762E-04,-2.14590E-03,-1.17824E-03, 3.66732E+00,
    -3.79729E-04,-6.13966E-03,-5.09082E-03,-1.96332E-03,-3.08280E-03,
    -9.75222E-04, 4.03315E+00,-2.52710E-01, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
 ]

# N2
# --

const pd_N2 = [
     1.16112E+00, 0.00000E+00, 0.00000E+00, 3.33725E-02, 0.00000E+00,
     3.48637E-02,-5.44368E-03, 0.00000E+00,-6.73940E-02, 1.74754E-01,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 1.74712E+02, 0.00000E+00,
     1.26733E-01, 0.00000E+00, 1.03154E+02, 5.52075E-02, 0.00000E+00,
     0.00000E+00, 8.13525E-04, 0.00000E+00, 0.00000E+00, 8.66784E-02,
     1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00,-2.50482E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.48894E-03,
     6.16053E-04,-5.79716E-04, 2.95482E-03, 8.47001E-02, 1.70147E-01,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 2.47425E-05, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
]

# TLB
# ---

const pd_TLB = [
     9.44846E-01, 0.00000E+00, 0.00000E+00,-3.08617E-02, 0.00000E+00,
    -2.44019E-02, 6.48607E-03, 0.00000E+00, 3.08181E-02, 4.59392E-02,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 1.74712E+02, 0.00000E+00,
     2.13260E-02, 0.00000E+00,-3.56958E+02, 0.00000E+00, 1.82278E-04,
     0.00000E+00, 3.07472E-04, 0.00000E+00, 0.00000E+00, 8.66784E-02,
     1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 3.83054E-03, 0.00000E+00, 0.00000E+00,
    -1.93065E-03,-1.45090E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00,-1.23493E-03, 1.36736E-03, 8.47001E-02, 1.70147E-01,
     3.71469E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     5.10250E-03, 2.47425E-05, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 3.68756E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
]

# O2
# --

const pd_O2 = [
     1.35580E+00, 1.44816E-01, 0.00000E+00, 6.07767E-02, 0.00000E+00,
     2.94777E-02, 7.46900E-02, 0.00000E+00,-9.23822E-02, 8.57342E-02,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 2.38636E+01, 0.00000E+00,
     7.71653E-02, 0.00000E+00, 8.18751E+01, 1.87736E-02, 0.00000E+00,
     0.00000E+00, 1.49667E-02, 0.00000E+00, 0.00000E+00, 8.66784E-02,
     1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00,-3.67874E+02, 5.48158E-03, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 8.47001E-02, 1.70147E-01,
     1.22631E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     8.17187E-03, 3.71617E-05, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.10826E-03,
    -3.13640E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
    -7.35742E-02,-5.00266E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 1.94965E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
]

# Ar
# --

const pd_Ar = [
     1.04761E+00, 2.00165E-01, 2.37697E-01, 3.68552E-02, 0.00000E+00,
     3.57202E-02,-2.14075E-01, 0.00000E+00,-1.08018E-01,-3.73981E-01,
     0.00000E+00, 3.10022E-02,-1.16305E-03,-2.07596E+01, 0.00000E+00,
     8.64502E-02, 0.00000E+00, 9.74908E+01, 5.16707E-02, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 8.66784E-02,
     1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 3.46193E+02, 1.34297E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-3.48509E-03,
    -1.54689E-04, 0.00000E+00, 0.00000E+00, 8.47001E-02, 1.70147E-01,
     1.47753E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     1.89320E-02, 3.68181E-05, 1.32570E-02, 0.00000E+00, 0.00000E+00,
     3.59719E-03, 7.44328E-03,-1.00023E-03,-6.50528E+03, 0.00000E+00,
     1.03485E-02,-1.00983E-03,-4.06916E-03,-6.60864E+01,-1.71533E-02,
     1.10605E-02, 1.20300E-02,-5.20034E-03, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
    -2.62769E+03, 7.13755E-03, 4.17999E-03, 0.00000E+00, 1.25910E+04,
     0.00000E+00, 0.00000E+00, 0.00000E+00,-2.23595E-03, 4.60217E-03,
     5.71794E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
    -3.18353E-02,-2.35526E-02,-1.36189E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 2.03522E-02,-6.67837E+01,-1.09724E-03, 0.00000E+00,
    -1.38821E-02, 1.60468E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.51574E-02,
    -5.44470E-04, 0.00000E+00, 7.28224E-02, 6.59413E-02, 0.00000E+00,
    -5.15692E-03, 0.00000E+00, 0.00000E+00,-3.70367E+03, 0.00000E+00,
     0.00000E+00, 1.36131E-02, 5.38153E-03, 0.00000E+00, 4.76285E+00,
    -1.75677E-02, 2.26301E-02, 0.00000E+00, 1.76631E-02, 4.77162E-03,
     0.00000E+00, 5.39354E+00, 0.00000E+00,-7.51710E-03, 0.00000E+00,
     0.00000E+00,-8.82736E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
]

# H
# -

const pd_H = [
     1.26376E+00,-2.14304E-01,-1.49984E-01, 2.30404E-01, 2.98237E-02,
     2.68673E-02, 2.96228E-01, 2.21900E-02,-2.07655E-02, 4.52506E-01,
     1.20105E-01, 3.24420E-02, 4.24816E-02,-9.14313E+00, 0.00000E+00,
     2.47178E-02,-2.88229E-02, 8.12805E+01, 5.10380E-02,-5.80611E-03,
     2.51236E-05,-1.24083E-02, 0.00000E+00, 0.00000E+00, 8.66784E-02,
     1.58727E-01,-3.48190E-02, 0.00000E+00, 0.00000E+00, 2.89885E-05,
     0.00000E+00, 1.53595E+02,-1.68604E-02, 0.00000E+00, 1.01015E-02,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.84552E-04,
    -1.22181E-03, 0.00000E+00, 0.00000E+00, 8.47001E-02, 1.70147E-01,
    -1.04927E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00,-5.91313E-03,
    -2.30501E-02, 3.14758E-05, 0.00000E+00, 0.00000E+00, 1.26956E-02,
     8.35489E-03, 3.10513E-04, 0.00000E+00, 3.42119E+03,-2.45017E-03,
    -4.27154E-04, 5.45152E-04, 1.89896E-03, 2.89121E+01,-6.49973E-03,
    -1.93855E-02,-1.48492E-02, 0.00000E+00,-5.10576E-02, 7.87306E-02,
     9.51981E-02,-1.49422E+04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     2.65503E+02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 6.37110E-03, 3.24789E-04,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     6.14274E-02, 1.00376E-02,-8.41083E-04, 0.00000E+00, 0.00000E+00,
     0.00000E+00,-1.27099E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
    -3.94077E-03,-1.28601E-02,-7.97616E-03, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00,-6.71465E-03,-1.69799E-03, 1.93772E-03, 3.81140E+00,
    -7.79290E-03,-1.82589E-02,-1.25860E-02,-1.04311E-02,-3.02465E-03,
     2.43063E-03, 3.63237E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
]

# N
# -

const pd_N = [
     7.09557E+01,-3.26740E-01, 0.00000E+00,-5.16829E-01,-1.71664E-03,
     9.09310E-02,-6.71500E-01,-1.47771E-01,-9.27471E-02,-2.30862E-01,
    -1.56410E-01, 1.34455E-02,-1.19717E-01, 2.52151E+00, 0.00000E+00,
    -2.41582E-01, 5.92939E-02, 4.39756E+00, 9.15280E-02, 4.41292E-03,
     0.00000E+00, 8.66807E-03, 0.00000E+00, 0.00000E+00, 8.66784E-02,
     1.58727E-01, 9.74701E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 6.70217E+01,-1.31660E-03, 0.00000E+00,-1.65317E-02,
     0.00000E+00, 0.00000E+00, 8.50247E-02, 2.77428E+01, 4.98658E-03,
     6.15115E-03, 9.50156E-03,-2.12723E-02, 8.47001E-02, 1.70147E-01,
    -2.38645E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.37380E-03,
    -8.41918E-03, 2.80145E-05, 7.12383E-03, 0.00000E+00,-1.66209E-02,
     1.03533E-04,-1.68898E-02, 0.00000E+00, 3.64526E+03, 0.00000E+00,
     6.54077E-03, 3.69130E-04, 9.94419E-04, 8.42803E+01,-1.16124E-02,
    -7.74414E-03,-1.68844E-03, 1.42809E-03,-1.92955E-03, 1.17225E-01,
    -2.41512E-02, 1.50521E+04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     1.60261E+03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00,-3.54403E-04,-1.87270E-02,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     2.76439E-02, 6.43207E-03,-3.54300E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00,-2.80221E-02, 8.11228E+01,-6.75255E-04, 0.00000E+00,
    -1.05162E-02,-3.48292E-03,-6.97321E-03, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00,-1.45546E-03,-1.31970E-02,-3.57751E-03,-1.09021E+00,
    -1.50181E-02,-7.12841E-03,-6.64590E-03,-3.52610E-03,-1.87773E-02,
    -2.22432E-03,-3.93895E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
]

# Hot O
# -----

const pd_hotO = [
     6.04050E-02, 1.57034E+00, 2.99387E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-1.51018E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00,-8.61650E+00, 1.26454E-02,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 5.50878E-03, 0.00000E+00, 0.00000E+00, 8.66784E-02,
     1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 6.23881E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 8.47001E-02, 1.70147E-01,
    -9.45934E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
]

const pd = [pd_He pd_O pd_N2 pd_TLB pd_O2 pd_Ar pd_H pd_N pd_hotO]'

# S Param
# =======

const ps = [
     9.56827E-01, 6.20637E-02, 3.18433E-02, 0.00000E+00, 0.00000E+00,
     3.94900E-02, 0.00000E+00, 0.00000E+00,-9.24882E-03,-7.94023E-03,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 1.74712E+02, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 2.74677E-03, 0.00000E+00, 1.54951E-02, 8.66784E-02,
     1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00,-6.99007E-04, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 1.24362E-02,-5.28756E-03, 8.47001E-02, 1.70147E-01,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 2.47425E-05, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
]

# Turbo
# =====

const pdl_1 = [
     1.09930E+00, 3.90631E+00, 3.07165E+00, 9.86161E-01, 1.63536E+01,
     4.63830E+00, 1.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 1.28840E+00, 3.10302E-02, 1.18339E-01
]

const pdl_2 = [
     1.00000E+00, 7.00000E-01, 1.15020E+00, 3.44689E+00, 1.28840E+00,
     1.00000E+00, 1.08738E+00, 1.22947E+00, 1.10016E+00, 7.34129E-01,
     1.15241E+00, 2.22784E+00, 7.95046E-01, 4.01612E+00, 4.47749E+00,
     1.23435E+02,-7.60535E-02, 1.68986E-06, 7.44294E-01, 1.03604E+00,
     1.72783E+02, 1.15020E+00, 3.44689E+00,-7.46230E-01, 9.49154E-01
]

const pdl = [pdl_1 pdl_2]'

# Lower Boundary
# ==============

const ptm = [
     1.04130E+03, 3.86000E+02, 1.95000E+02, 1.66728E+01, 2.13000E+02,
     1.20000E+02, 2.40000E+02, 1.87000E+02,-2.00000E+00, 0.00000E+00
]

const pdm_1 = [
     2.45600E+07, 6.71072E-06, 1.00000E+02, 0.00000E+00, 1.10000E+02,
     1.00000E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
]

const pdm_2 = [
     8.59400E+10, 1.00000E+00, 1.05000E+02,-8.00000E+00, 1.10000E+02,
     1.00000E+01, 9.00000E+01, 2.00000E+00, 0.00000E+00, 0.00000E+00
]

const pdm_3 = [
     2.81000E+11, 0.00000E+00, 1.05000E+02, 2.80000E+01, 2.89500E+01,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
]

const pdm_4 = [
     3.30000E+10, 2.68270E-01, 1.05000E+02, 1.00000E+00, 1.10000E+02,
     1.00000E+01, 1.10000E+02,-1.00000E+01, 0.00000E+00, 0.00000E+00
]

const pdm_5 = [
     1.33000E+09, 1.19615E-02, 1.05000E+02, 0.00000E+00, 1.10000E+02,
     1.00000E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
]

const pdm_6 = [
     1.76100E+05, 1.00000E+00, 9.50000E+01,-8.00000E+00, 1.10000E+02,
     1.00000E+01, 9.00000E+01, 2.00000E+00, 0.00000E+00, 0.00000E+00,
]

const pdm_7 = [
     1.00000E+07, 1.00000E+00, 1.05000E+02,-8.00000E+00, 1.10000E+02,
     1.00000E+01, 9.00000E+01, 2.00000E+00, 0.00000E+00, 0.00000E+00
]

const pdm_8 = [
     1.00000E+06, 1.00000E+00, 1.05000E+02,-8.00000E+00, 5.50000E+02,
     7.60000E+01, 9.00000E+01, 2.00000E+00, 0.00000E+00, 4.00000E+03
]

const pdm = [pdm_1 pdm_2 pdm_3 pdm_4 pdm_5 pdm_6 pdm_7 pdm_8]'

# TN1(2)
# ======

const ptl_1 = [
     1.00858E+00, 4.56011E-02,-2.22972E-02,-5.44388E-02, 5.23136E-04,
    -1.88849E-02, 5.23707E-02,-9.43646E-03, 6.31707E-03,-7.80460E-02,
    -4.88430E-02, 0.00000E+00, 0.00000E+00,-7.60250E+00, 0.00000E+00,
    -1.44635E-02,-1.76843E-02,-1.21517E+02, 2.85647E-02, 0.00000E+00,
     0.00000E+00, 6.31792E-04, 0.00000E+00, 5.77197E-03, 8.66784E-02,
     1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00,-8.90272E+03, 3.30611E-03, 3.02172E-03, 0.00000E+00,
    -2.13673E-03,-3.20910E-04, 0.00000E+00, 0.00000E+00, 2.76034E-03,
     2.82487E-03,-2.97592E-04,-4.21534E-03, 8.47001E-02, 1.70147E-01,
     8.96456E-03, 0.00000E+00,-1.08596E-02, 0.00000E+00, 0.00000E+00,
     5.57917E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 9.65405E-03, 0.00000E+00, 0.00000E+00, 2.00000E+00
]

# TN1(3)
# ======

const ptl_2 = [
     9.39664E-01, 8.56514E-02,-6.79989E-03, 2.65929E-02,-4.74283E-03,
     1.21855E-02,-2.14905E-02, 6.49651E-03,-2.05477E-02,-4.24952E-02,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 1.19148E+01, 0.00000E+00,
     1.18777E-02,-7.28230E-02,-8.15965E+01, 1.73887E-02, 0.00000E+00,
     0.00000E+00, 0.00000E+00,-1.44691E-02, 2.80259E-04, 8.66784E-02,
     1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 2.16584E+02, 3.18713E-03, 7.37479E-03, 0.00000E+00,
    -2.55018E-03,-3.92806E-03, 0.00000E+00, 0.00000E+00,-2.89757E-03,
    -1.33549E-03, 1.02661E-03, 3.53775E-04, 8.47001E-02, 1.70147E-01,
    -9.17497E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     3.56082E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00,-1.00902E-02, 0.00000E+00, 0.00000E+00, 2.00000E+00
]

# TN1(4)
# ======

const ptl_3 = [
     9.85982E-01,-4.55435E-02, 1.21106E-02, 2.04127E-02,-2.40836E-03,
     1.11383E-02,-4.51926E-02, 1.35074E-02,-6.54139E-03, 1.15275E-01,
     1.28247E-01, 0.00000E+00, 0.00000E+00,-5.30705E+00, 0.00000E+00,
    -3.79332E-02,-6.24741E-02, 7.71062E-01, 2.96315E-02, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 6.81051E-03,-4.34767E-03, 8.66784E-02,
     1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 1.07003E+01,-2.76907E-03, 4.32474E-04, 0.00000E+00,
     1.31497E-03,-6.47517E-04, 0.00000E+00,-2.20621E+01,-1.10804E-03,
    -8.09338E-04, 4.18184E-04, 4.29650E-03, 8.47001E-02, 1.70147E-01,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
    -4.04337E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-9.52550E-04,
     8.56253E-04, 4.33114E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.21223E-03,
     2.38694E-04, 9.15245E-04, 1.28385E-03, 8.67668E-04,-5.61425E-06,
     1.04445E+00, 3.41112E+01, 0.00000E+00,-8.40704E-01,-2.39639E+02,
     7.06668E-01,-2.05873E+01,-3.63696E-01, 2.39245E+01, 0.00000E+00,
    -1.06657E-03,-7.67292E-04, 1.54534E-04, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
]

# TN1(5) TN2(1)
# =============

const ptl_4 = [
     1.00320E+00, 3.83501E-02,-2.38983E-03, 2.83950E-03, 4.20956E-03,
     5.86619E-04, 2.19054E-02,-1.00946E-02,-3.50259E-03, 4.17392E-02,
    -8.44404E-03, 0.00000E+00, 0.00000E+00, 4.96949E+00, 0.00000E+00,
    -7.06478E-03,-1.46494E-02, 3.13258E+01,-1.86493E-03, 0.00000E+00,
    -1.67499E-02, 0.00000E+00, 0.00000E+00, 5.12686E-04, 8.66784E-02,
     1.58727E-01,-4.64167E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     4.37353E-03,-1.99069E+02, 0.00000E+00,-5.34884E-03, 0.00000E+00,
     1.62458E-03, 2.93016E-03, 2.67926E-03, 5.90449E+02, 0.00000E+00,
     0.00000E+00,-1.17266E-03,-3.58890E-04, 8.47001E-02, 1.70147E-01,
     0.00000E+00, 0.00000E+00, 1.38673E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.60571E-03,
     6.28078E-04, 5.05469E-05, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-1.57829E-03,
    -4.00855E-04, 5.04077E-05,-1.39001E-03,-2.33406E-03,-4.81197E-04,
     1.46758E+00, 6.20332E+00, 0.00000E+00, 3.66476E-01,-6.19760E+01,
     3.09198E-01,-1.98999E+01, 0.00000E+00,-3.29933E+02, 0.00000E+00,
    -1.10080E-03,-9.39310E-05, 1.39638E-04, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
]

const ptl = [ptl_1 ptl_2 ptl_3 ptl_4]'

# TN2(2)
# ======

const pma_1 = [
     9.81637E-01,-1.41317E-03, 3.87323E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-3.58707E-02,
    -8.63658E-03, 0.00000E+00, 0.00000E+00,-2.02226E+00, 0.00000E+00,
    -8.69424E-03,-1.91397E-02, 8.76779E+01, 4.52188E-03, 0.00000E+00,
     2.23760E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00,-7.07572E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00,
    -4.11210E-03, 3.50060E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00,-8.36657E-03, 1.61347E+01, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00,-1.45130E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.24152E-03,
     6.43365E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.33255E-03,
     2.42657E-03, 1.60666E-03,-1.85728E-03,-1.46874E-03,-4.79163E-06,
     1.22464E+00, 3.53510E+01, 0.00000E+00, 4.49223E-01,-4.77466E+01,
     4.70681E-01, 8.41861E+00,-2.88198E-01, 1.67854E+02, 0.00000E+00,
     7.11493E-04, 6.05601E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
]

# TN2(3)
# ======

const pma_2 = [
     1.00422E+00,-7.11212E-03, 5.24480E-03, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-5.28914E-02,
    -2.41301E-02, 0.00000E+00, 0.00000E+00,-2.12219E+01,-1.03830E-02,
    -3.28077E-03, 1.65727E-02, 1.68564E+00,-6.68154E-03, 0.00000E+00,
     1.45155E-02, 0.00000E+00, 8.42365E-03, 0.00000E+00, 0.00000E+00,
     0.00000E+00,-4.34645E-03, 0.00000E+00, 0.00000E+00, 2.16780E-02,
     0.00000E+00,-1.38459E+02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 7.04573E-03,-4.73204E+01, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 1.08767E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-8.08279E-03,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 5.21769E-04,
    -2.27387E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 3.26769E-03,
     3.16901E-03, 4.60316E-04,-1.01431E-04, 1.02131E-03, 9.96601E-04,
     1.25707E+00, 2.50114E+01, 0.00000E+00, 4.24472E-01,-2.77655E+01,
     3.44625E-01, 2.75412E+01, 0.00000E+00, 7.94251E+02, 0.00000E+00,
     2.45835E-03, 1.38871E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
]

# TN2(4) TN3(1)
# =============

const pma_3 = [
     1.01890E+00,-2.46603E-02, 1.00078E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-6.70977E-02,
    -4.02286E-02, 0.00000E+00, 0.00000E+00,-2.29466E+01,-7.47019E-03,
     2.26580E-03, 2.63931E-02, 3.72625E+01,-6.39041E-03, 0.00000E+00,
     9.58383E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00,-1.85291E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 1.39717E+02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 9.19771E-03,-3.69121E+02, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00,-1.57067E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-7.07265E-03,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.92953E-03,
    -2.77739E-03,-4.40092E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.47280E-03,
     2.95035E-04,-1.81246E-03, 2.81945E-03, 4.27296E-03, 9.78863E-04,
     1.40545E+00,-6.19173E+00, 0.00000E+00, 0.00000E+00,-7.93632E+01,
     4.44643E-01,-4.03085E+02, 0.00000E+00, 1.15603E+01, 0.00000E+00,
     2.25068E-03, 8.48557E-04,-2.98493E-04, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
]

# TN3(2)
# ======

const pma_4 = [
     9.75801E-01, 3.80680E-02,-3.05198E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 3.85575E-02,
     5.04057E-02, 0.00000E+00, 0.00000E+00,-1.76046E+02, 1.44594E-02,
    -1.48297E-03,-3.68560E-03, 3.02185E+01,-3.23338E-03, 0.00000E+00,
     1.53569E-02, 0.00000E+00,-1.15558E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 4.89620E-03, 0.00000E+00, 0.00000E+00,-1.00616E-02,
    -8.21324E-03,-1.57757E+02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 6.63564E-03, 4.58410E+01, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00,-2.51280E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 9.91215E-03,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-8.73148E-04,
    -1.29648E-03,-7.32026E-05, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-4.68110E-03,
    -4.66003E-03,-1.31567E-03,-7.39390E-04, 6.32499E-04,-4.65588E-04,
    -1.29785E+00,-1.57139E+02, 0.00000E+00, 2.58350E-01,-3.69453E+01,
     4.10672E-01, 9.78196E+00,-1.52064E-01,-3.85084E+03, 0.00000E+00,
    -8.52706E-04,-1.40945E-03,-7.26786E-04, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
]

# TN3(3)
# ======

const pma_5 = [
     9.60722E-01, 7.03757E-02,-3.00266E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.22671E-02,
     4.10423E-02, 0.00000E+00, 0.00000E+00,-1.63070E+02, 1.06073E-02,
     5.40747E-04, 7.79481E-03, 1.44908E+02, 1.51484E-04, 0.00000E+00,
     1.97547E-02, 0.00000E+00,-1.41844E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 5.77884E-03, 0.00000E+00, 0.00000E+00, 9.74319E-03,
     0.00000E+00,-2.88015E+03, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00,-4.44902E-03,-2.92760E+01, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 2.34419E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 5.36685E-03,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-4.65325E-04,
    -5.50628E-04, 3.31465E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.06179E-03,
    -3.08575E-03,-7.93589E-04,-1.08629E-04, 5.95511E-04,-9.05050E-04,
     1.18997E+00, 4.15924E+01, 0.00000E+00,-4.72064E-01,-9.47150E+02,
     3.98723E-01, 1.98304E+01, 0.00000E+00, 3.73219E+03, 0.00000E+00,
    -1.50040E-03,-1.14933E-03,-1.56769E-04, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
]

# TN3(4)
# ======

const pma_6 = [
     1.03123E+00,-7.05124E-02, 8.71615E-03, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-3.82621E-02,
    -9.80975E-03, 0.00000E+00, 0.00000E+00, 2.89286E+01, 9.57341E-03,
     0.00000E+00, 0.00000E+00, 8.66153E+01, 7.91938E-04, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 4.68917E-03, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 7.86638E-03, 0.00000E+00, 0.00000E+00, 9.90827E-03,
     0.00000E+00, 6.55573E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00,-4.00200E+01, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 7.07457E-03, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 5.72268E-03,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.04970E-04,
     1.21560E-03,-8.05579E-06, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.49941E-03,
    -4.57256E-04,-1.59311E-04, 2.96481E-04,-1.77318E-03,-6.37918E-04,
     1.02395E+00, 1.28172E+01, 0.00000E+00, 1.49903E-01,-2.63818E+01,
     0.00000E+00, 4.70628E+01,-2.22139E-01, 4.82292E-02, 0.00000E+00,
    -8.67075E-04,-5.86479E-04, 5.32462E-04, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
]

# TN3(5) SURFACE TEMP TSL
# =======================

const pma_7 = [
     1.00828E+00,-9.10404E-02,-2.26549E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.32420E-02,
    -9.08925E-03, 0.00000E+00, 0.00000E+00, 3.36105E+01, 0.00000E+00,
     0.00000E+00, 0.00000E+00,-1.24957E+01,-5.87939E-03, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 2.79765E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 2.01237E+03, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00,-1.75553E-02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 3.29699E-03,
     1.26659E-03, 2.68402E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.17894E-03,
     1.48746E-03, 1.06478E-04, 1.34743E-04,-2.20939E-03,-6.23523E-04,
     6.36539E-01, 1.13621E+01, 0.00000E+00,-3.93777E-01, 2.38687E+03,
     0.00000E+00, 6.61865E+02,-1.21434E-01, 9.27608E+00, 0.00000E+00,
     1.68478E-04, 1.24892E-03, 1.71345E-03, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
]

# TGN3(2) SURFACE GRAD TSLG
# =========================

const pma_8 = [
     1.57293E+00,-6.78400E-01, 6.47500E-01, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-7.62974E-02,
    -3.60423E-01, 0.00000E+00, 0.00000E+00, 1.28358E+02, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 4.68038E+01, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00,-1.67898E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 2.90994E+04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 3.15706E+01, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
]

# TGN2(1) TGN1(2)
# ===============

const pma_9 = [
     8.60028E-01, 3.77052E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-1.17570E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 7.77757E-03, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 1.01024E+02, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 6.54251E+02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-1.56959E-02,
     1.91001E-02, 3.15971E-02, 1.00982E-02,-6.71565E-03, 2.57693E-03,
     1.38692E+00, 2.82132E-01, 0.00000E+00, 0.00000E+00, 3.81511E+02,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
]

# TGN3(1) TGN2(2)
# ===============

const pma_10 = [
     1.06029E+00,-5.25231E-02, 3.73034E-01, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 3.31072E-02,
    -3.88409E-01, 0.00000E+00, 0.00000E+00,-1.65295E+02,-2.13801E-01,
    -4.38916E-02,-3.22716E-01,-8.82393E+01, 1.18458E-01, 0.00000E+00,
    -4.35863E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00,-1.19782E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 2.62229E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00,-5.37443E+01, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00,-4.55788E-01, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 3.84009E-02,
     3.96733E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 5.05494E-02,
     7.39617E-02, 1.92200E-02,-8.46151E-03,-1.34244E-02, 1.96338E-02,
     1.50421E+00, 1.88368E+01, 0.00000E+00, 0.00000E+00,-5.13114E+01,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     5.11923E-02, 3.61225E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
     0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
]

const pma = [pma_1 pma_2 pma_3 pma_4 pma_5 pma_6 pma_7 pma_8 pma_9 pma_10]'

# Middle Atmosphere Averages
# ==========================

const pavgm = [
     2.61000E+02, 2.64000E+02, 2.29000E+02, 2.17000E+02, 2.17000E+02,
     2.23000E+02, 2.86760E+02,-2.93940E+00, 2.50000E+00, 0.00000E+00
]


text:
# test with the real data set
# Test Bayesian Detctor
include("./lib/metadata.jl")
include("./2bayes_detector_replica.jl")
include("./lib/SignalGenerator.jl")
include("./lib/roc.jl")

# FIRST OBTAINE SAMPLE DATA USING (TEST_BD_REAL.R -> real-signal-prepare.R) SCRIPT

using SignalGenerator
using BayesChpOrig
using Metadata
using Roc

function test_bd_detector_real_data()
    work_path = "./data/activity-recognition-sample/"
    signal = readdlm(work_path * "sig-activ.dat", ',')
    x = signal[:,1]
    y = signal[:,2]
    # matrix_max_vals, matrix = bayes_detector(x, lambda = 500)
    # matrix_max_vals, matrix = bayes_detector(x, lambda = 1000)
    # matrix_max_vals, matrix = bayes_detector(x, lambda = 2000)
    # matrix_max_vals, matrix = bayes_detector(x, lambda = 2500)
    matrix_max_vals, matrix = bayes_detector(x, lambda = 3000)
    detections = extract_changes(matrix_max_vals)
    #
    # max_delay = 30
    println("... test_bd_real > save detections into ", work_path*"test_bd_real_detections.txt")
    writedlm(work_path*"test_bd_real_detections.txt", detections)
    # writedlm(outdir*"test_bd_signal.txt", sig)
    # writedlm(outdir*"test_bd_changes.txt", chps)
    # writedlm(outdir*"test_bd_matrix.txt", matrix, ",")
    # writedlm(outdir*"test_bd_matrix_max_vals.txt", matrix_max_vals, ",")
end
test_bd_detector_real_data()

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency