text:
module Synapses
using RealNeuralNetworks.Utils.SynapseTables
using RealNeuralNetworks.Utils.BoundingBoxes
using DataFrames

export Synapse

mutable struct Synapse{T}
    psdSegmentationId               ::Int
    psdCoordinate                   ::NTuple{3,T} # nm
    psdBoundingBox                  ::BoundingBox
    psdSize                         ::Int
    preSynapticSegmentationId       ::Int
    preSynapticCoordinate           ::NTuple{3,T} # nm
    preSynapticWeight               ::Float32
    postSynapticSegmentationId      ::Int
    postSynapticCoordinate          ::NTuple{3,T} # nm
    postSynapticWeight              ::Float32
end

function Synapse( row::DataFrameRow )
    psdBoundingBox = BoundingBox(   SynapseTables.get_coordinate(row, "BBOX_b"),
                                    SynapseTables.get_coordinate(row, "BBOX_e") )
    psdCoordinate = SynapseTables.get_coordinate(row, "psd_")
    preSynapticCoordinate = SynapseTables.get_coordinate(row, "presyn_")
    postSynapticCoordinate = SynapseTables.get_coordinate(row, "postsyn_")

    Synapse{Float32}(row[:psd_segid], psdCoordinate, psdBoundingBox, row[:size],
            row[:presyn_segid],  preSynapticCoordinate,  row[:presyn_wt],
            row[:postsyn_segid], postSynapticCoordinate, row[:postsyn_wt])
end

function Synapse( df::DataFrame )
    @assert DataFrames.nrow( df ) == 1
    Synapse( df[1,:] )
end

############### Base functions ################
function Base.show(self::Synapse)
    println("psd segmentation id: ", psdSegmentationId)
    println("psd coordinate: ", psdCoordinate)
    println("psd bounding box: ", psdBoundingBox)
    println("psd size: ", psdSize)
    println("\npre synaptic segmentation id: ", preSynapticSegmentationId)
    println("pre synaptic coordinate: ", preSynapticCoordinate)
    println("pre synaptic weight: ", preSynapticWeight)
    println("\npost synaptic segmentation id: ", postSynapticSegmentationId)
    println("post synaptic coordinate: ", postSynapticCoordinate)
    println("post synaptic weight: ", postSynapticWeight)
    nothing
end

@inline function Base.zero(self::Synapse)
    nothing
end

############### properties ####################
"""
Note that the psd size was measured using voxel numbers.
"""
@inline function get_psd_size(self::Synapse) self.psdSize end
@inline function get_psd_segmentation_id(self::Synapse) self.psdSegmentationId end
@inline function get_psd_coordinate( self::Synapse ) self.psdCoordinate end
@inline function get_psd_bounding_box( self::Synapse ) self.psdBoundingBox end

@inline function get_pre_synaptic_segmentation_id( self::Synapse ) self.preSynapticSegmentationId end
@inline function get_pre_synaptic_coordinate( self::Synapse ) self.preSynapticCoordinate end
@inline function get_pre_synaptic_weight( self::Synapse ) self.preSynapticWeight end

@inline function get_post_synaptic_segmentation_id( self::Synapse )
    self.postSynapticSegmentationId
end
@inline function get_post_synaptic_coordinate( self::Synapse )
    self.postSynapticCoordinate
end
@inline function get_post_synaptic_weight( self::Synapse ) self.postSynapticWeight end

#@inline function isbutton(self::Synapse, neuronId::Int)
#    neuronId == get_presynaptic_segmentation_id(self)
#end

end # module


text:
using Test

@testset "Tweaks.jl" begin
    # Run all test in src/*_test.jl
    cd(joinpath(@__DIR__, ".."))
    for (root, dirs, files) in walkdir("src")
        for file in files
            endswith(file, "_test.jl") || continue
            fpath = joinpath(root, file)
            @testset "$fpath" begin
                include(joinpath("..", fpath))
            end
        end
    end
end

ext: jl


text:
using CSV
using DataFrames

cd("/Users/kiristern/Documents/GitHub/SDM/data/prey_sp")

# Import and read all saved GBIF csv files from directory
file_list = readdir()
#remove .DS_Store
file_list = file_list[2:end]
#create an empty array of dimension Any
files = Any[]
#fill array with data from downloaded csv files
for i in 1:length(file_list)
    #see which data file is causing the error problems
    global j = i
    #push! into files, info from file_list
    push!(files, CSV.read(file_list[i], delim="\t"))
end
#create one dataframe
df_prey = vcat(files)

# Select specific columns only
df = df_prey[:, [:species, :taxonRank, :decimalLatitude, :decimalLongitude, :year]]

#save csv file
CSV.write("prey_sp_onedf.csv", df)

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM


text:
using Base.Test
using MathProgComplex

## GLOBAL OPTIMIZATION WITH POLYNOMIALS AND THE PROBLEM OF MOMENTS
# JEAN B. LASSERRE, 2001
# http://www.ii.uib.no/~lennart/drgrad/Lasserre2001.pdf

@testset "Lasserre2001 small real problems" begin

    @testset "Example 5 (p. 811)" begin
        problem = Problem()
        x1 = Variable("x1", Real); add_variable!(problem, x1)
        x2 = Variable("x2", Real); add_variable!(problem, x2)

        set_objective!(problem, -(x1-1)^2 -(x1-x2)^2 -(x2-3)^2)
        add_constraint!(problem, "crt1", (1-(x1-1)^2) >> 0)
        add_constraint!(problem, "crt2", (1-(x1-x2)^2) >> 0)
        add_constraint!(problem, "crt3", (1-(x2-3)^2) >> 0)

        logpath = joinpath(Pkg.dir("MathProgComplex"), "Mosek_runs", "testfolder")

        ## Order 1
        relax_ctx = set_relaxation(problem; hierarchykind=:Real,
                                            d = 1,
                                            params = Dict(:opt_outlev=>0,
                                                          :opt_logpath=>logpath,
                                                          :opt_solver=>testsolver))

        primobj, dualobj = run_hierarchy(problem, relax_ctx)

        @show primobj, dualobj, -3
        @test primobj ≈ -3 atol=1e-4
        @test dualobj ≈ primobj atol=mosek_optgap*min(abs(primobj), abs(dualobj))

        ## Order 2
        relax_ctx = set_relaxation(problem; hierarchykind=:Real,
                                            d = 2,
                                            params = Dict(:opt_outlev=>0,
                                                          :opt_logpath=>logpath,
                                                          :opt_solver=>testsolver))

        primobj, dualobj = run_hierarchy(problem, relax_ctx)

        @show primobj, dualobj, -2
        @test primobj ≈ -2 atol=1e-4
        @test dualobj ≈ primobj atol=mosek_optgap*min(abs(primobj), abs(dualobj))
    end
end


text:
using IntroML
using Test

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


text:
# This is used in Jupyter notebooks
Base.show(io::IO, mime::MIME"image/png", img::AstroImage; kwargs...) =
    show(io, mime, render(img), kwargs...)


text:
"""
    display_matrix(M, rownames, colnames)

Display matrix with rownames and colnames.
"""
function display_matrix(M, rownames, colnames; nrow=6, ncol=6)
    M_part = M[1:nrow, 1:ncol]
    rownames_part = rownames[1:nrow]
    colnames_part = colnames[1:ncol]

    df = DataFrame()
    df[!, :Data] = rownames_part
    for (i,coln) in enumerate(colnames_part)
        df[!, coln] = M_part[:,i]
    end

    display(df)
end

"""
    display_matrix(data, target_col, cue_obj, M, M_type)

Display matrix with rownames and colnames.

# Obligatory Arguments
- `data::DataFrame`: the dataset
- `target_col::Union{String, Symbol}`: the target column name
- `cue_obj::Cue_Matrix_Struct`: the cue matrix structure
- `M::Union{SparseMatrixCSC, Matrix}`: the matrix
- `M_type::Union{String, Symbol}`: the type of the matrix, currently support :C, :S, :F, :G, :Chat, :Shat, :A and :R

# Optional Arguments
- `nrow::Int64 = 6`: the number of rows to display
- `ncol::Int64 = 6`: the number of columns to display

# Examples
```julia
JudiLing.display_matrix(latin, :Word, cue_obj, cue_obj.C, :C)
JudiLing.display_matrix(latin, :Word, cue_obj, S, :S)
JudiLing.display_matrix(latin, :Word, cue_obj, G, :G)
JudiLing.display_matrix(latin, :Word, cue_obj, Chat, :Chat)
JudiLing.display_matrix(latin, :Word, cue_obj, F, :F)
JudiLing.display_matrix(latin, :Word, cue_obj, Shat, :Shat)
JudiLing.display_matrix(latin, :Word, cue_obj, A, :A)
JudiLing.display_matrix(latin, :Word, cue_obj, R, :R)
```
"""
function display_matrix(
    data,
    target_col,
    cue_obj,
    M,
    M_type;
    nrow = 6,
    ncol = 6
    )

    if M_type == :C || M_type == "C"
        rownames = data[:,target_col]
        colnames = [cue_obj.i2f[i] for i in 1:size(M,2)]
    elseif M_type == :Chat || M_type == "Chat"
        rownames = data[:,target_col]
        colnames = [cue_obj.i2f[i] for i in 1:size(M,2)]
    elseif M_type == :S || M_type == "S"
        rownames = data[:,target_col]
        colnames = ["S$i" for i in 1:size(M,2)]
    elseif M_type == :Shat || M_type == "Shat"
        rownames = data[:,target_col]
        colnames = ["S$i" for i in 1:size(M,2)]
    elseif M_type == :F || M_type == "F"
        rownames = colnames = [cue_obj.i2f[i] for i in 1:size(M,2)]
        colnames = ["S$i" for i in 1:size(M,2)]
    elseif M_type == :G || M_type == "G"
        rownames = ["S$i" for i in 1:size(M,2)]
        colnames = [cue_obj.i2f[i] for i in 1:size(M,2)]
    elseif M_type == :A || M_type == "A"
        rownames = [cue_obj.i2f[i] for i in 1:size(M,2)]
        colnames = [cue_obj.i2f[i] for i in 1:size(M,2)]
    elseif M_type == :R || M_type == "R"
        rownames = data[:,target_col]
        colnames = data[:,target_col]
    else
        throw(ArgumentError("type is incorrect, using :C," *
            " :S, :Chat, :Shat, :F, :G, :A and :R"))
    end

    display_matrix(M, rownames, colnames, nrow=nrow, ncol=ncol)
end


text:
mutable struct DDEStats
  nf::Int
  nf2::Int
  nw::Int
  nsolve::Int
  njacs::Int
  nnonliniter::Int
  nnonlinconvfail::Int
  nfpiter::Int
  nfpconvfail::Int
  ncondition::Int
  naccept::Int
  nreject::Int
  maxeig::Float64
end

DDEStats(x::Int = -1) = DDEStats(x, x, x, x, x, x, x, x, x, x, x, x, 0.0)

function Base.show(io::IO, s::DDEStats)
  println(io, summary(s))
  @printf io "%-60s %-d\n" "Number of function 1 evaluations:" s.nf
  @printf io "%-60s %-d\n" "Number of function 2 evaluations:" s.nf2
  @printf io "%-60s %-d\n" "Number of W matrix evaluations:" s.nw
  @printf io "%-60s %-d\n" "Number of linear solves:" s.nsolve
  @printf io "%-60s %-d\n" "Number of Jacobians created:" s.njacs
  @printf io "%-60s %-d\n" "Number of nonlinear solver iterations:" s.nnonliniter
  @printf io "%-60s %-d\n" "Number of nonlinear solver convergence failures:" s.nnonlinconvfail
  @printf io "%-60s %-d\n" "Number of fixed-point solver iterations:" s.nfpiter
  @printf io "%-60s %-d\n" "Number of fixed-point solver convergence failures:" s.nfpconvfail
  @printf io "%-60s %-d\n" "Number of rootfind condition calls:" s.ncondition
  @printf io "%-60s %-d\n" "Number of accepted steps:"  s.naccept
  @printf io "%-60s %-d" "Number of rejected steps:" s.nreject
  iszero(s.maxeig) || @printf io "\n%-60s %-d" "Maximum eigenvalue recorded:" s.maxeig
end

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
rm("../output/"; recursive = true, force = true)
mkdir("../output/")
include("calculated/model.jl")
include("simulated/sampleMaker.jl")
include("scatt.jl")

cd("calculated/")

using ..Model
make_model_data(0.005, 1.5)

cd("..")

using ..Scatt.Finitizer
finitize_data("calculated/calculated.data", 10 .^ (-1:0.01:1), 10.0)
targetIntensity = intensity_at_zero("../output/finitized.data")

cd("simulated/")

using ..SampleMaker.Tools
make_hard_sphere_samples(100, 0.005, 1.5)
set_sld_offsets(targetIntensity)

cd("..")

using ..Scatt.Scatterer
measure("simulated/sample/")

using ..Scatt.Plotter
make_plot()

max_stars_repo_name: JonPallbo/Scatt

max_issues_repo_name: JonPallbo/Scatt

max_forks_repo_name: JonPallbo/Scatt


text:
"""
	struct AuxiliaryExtractor <: AbstractExtractor
		extractor::AbstractExtractor
		extract_fun::Function
	end

Universal extractor for applying any function,
which lets you ambed any transformation into the AbstractExtractor machinery.
Useful e.g. for extractors accompanying trained models, where you need to apply yet another transformation.

```jldocstest
julia> e1 = ExtractDict(Dict(:a=>ExtractString(), :b=>ExtractString()));

julia> e2 = AuxiliaryExtractor(e1, (e,x)->e[:a](x["a"]))
Auxiliary extractor with
  └── Dict
        ├── a: String
        └── b: String

julia> e2(Dict("a"=>"Hello", "b"=>"World"))
ArrayNode{NGramMatrix{String,Array{String,1},Int64},Nothing}:
 "Hello"
```

"""
struct AuxiliaryExtractor <: AbstractExtractor
	extractor::AbstractExtractor
	extract_fun::Function
end

(e::AuxiliaryExtractor)(x::HierarchicType; store_input=false) = e.extract_fun(e.extractor, x; store_input)


text:
module TestCancellation

using Julio
using Test
using ..Utils: @test_error, ⊏

function test_onclose()
    ncalls = Ref(0)
    err = @test_error Julio.withtaskgroup() do tg
        Julio.oncancel() do
            ncalls[] += 1
        end
        error("_CANCELLING_")
    end
    @test "_CANCELLING_" ⊏ sprint(showerror, err)
    @test ncalls[] == 1
end

function test_shield()
    ie1, oe1 = Julio.channel()
    ie2, oe2 = Julio.channel()
    ie3, oe3 = Julio.channel()
    local task
    Julio.withtaskgroup() do tg
        task = Julio.spawn!(tg) do
            @test_error Julio.withtaskgroup() do tg
                Julio.spawn!(tg) do
                    Julio.shield() do
                        put!(ie2, :enter_shield)
                        put!(ie2, :leave_shield)
                    end
                end
                Julio.spawn!(tg) do
                    try
                        put!(ie3, :waiting)
                    finally
                        Julio.shield() do
                            put!(ie1, :thrown)
                        end
                    end
                end
                error("_CANCELLING_")
            end
        end
        @test take!(oe2) === :enter_shield
        @test take!(oe1) === :thrown
        @test take!(oe2) === :leave_shield
    end
    @test "_CANCELLING_" ⊏ sprint(showerror, fetch(task))
end

function spin()
    while true
        Julio.yield()
    end
end

function test_yield()
    Julio.withtaskgroup() do tg
        Julio.spawn!(spin, tg)
        Julio.spawn!(spin, tg)
        sleep(0.01)
        Julio.cancel!(tg)
    end
    @test true
end

function test_nested_cancel()
    timeout = Ref(false)
    Julio.withtaskgroup() do tg0
        Julio.withtaskgroup() do tg1
            Julio.spawn!(tg1) do
                Julio.sleep(10)
                timeout[] = true
            end
            Julio.cancel!(tg0)
        end
    end
    @test !timeout[]
end

function test_iscancelled()
    local a, b, c
    ie, oe = Julio.channel()
    Julio.withtaskgroup() do tg
        Julio.spawn!(tg) do
            a = Julio.iscancelled()
            put!(ie, 111)
            try
                put!(ie, 222)
            finally
                b = Julio.iscancelled()
                Julio.shield() do
                    c = Julio.iscancelled()
                end
            end
        end
        @test take!(oe) == 111
        Julio.cancel!(tg)
    end
    @test !a
    @test b
    @test !c
end

function test_close_cancelled()
    sept, rept = Julio.channel()
    Julio.withtaskgroup() do tg
        Julio.spawn!(tg) do
            put!(sept, 111)
            try
                put!(sept, 222)   # does not succeed
            finally
                close(sept)
            end
        end
        @test take!(rept) == 111
        Julio.cancel!(tg)
    end
    ans = :notset
    Julio.withtimeout(1) do
        ans = iterate(rept)
    end
    @test ans === nothing
end

end  # module


text:
using Dolo

# model = yaml_import("examples/models/rbc_mc.yaml")
model = yaml_import("examples/models/rbc.yaml")
# model = yaml_import("examples/models/rbc_iid.yaml")

sol = Dolo.improved_time_iteration(model)

G = Dolo.distG(model, sol)

Π, dΠ = Dolo.transition_matrix(G; diff=true)

DΠ = Π .!= 0.0
DdΠ = [ maximum(abs,DΠ[c]) for c in CartesianIndices(dΠ)]

Π = Dolo.transition_matrix(G)
sum(Π; dims=2) # should be only ones

x0 = G.x0

# G contains all informations to compute the ergodic_distribution
m = Dolo.ergodic_distribution(G)

Dolo.ergodic_distribution(model, sol)

# we can compute distribution transitions

@time μ1 = G(G.μ0, G.x0);
@time μ1 = G(G.μ0, G.x0);

x0 = G.x0
μ0 = G.μ0

# we compute differential operators w.r.t. arguments
μ1, ∂G_∂μ, ∂G_∂x = G(μ0, x0, diff=true);

# we can pass vectors as arguments
x0_flat = cat(G.x0.data...; dims=1)
μ1, ∂G_∂μ, ∂G_∂x = G(μ0, x0_flat; diff=true)

# differential operate *like* matrices and can be composed like usual matrices
∂G_∂μ(μ0*0.01)
∂G_∂μ*(μ0*0.01) # equivalent

∂G_∂x*(x0_flat*0.01)

[∂G_∂μ ∂G_∂μ]
[∂G_∂μ ∂G_∂x]([μ0; x0_flat]) # yes, it is kind of beautiful!

using FiniteDiff

Jμ_num = FiniteDiff.finite_difference_jacobian(mu->G(mu, x0_flat), μ0)
Jμ_exact = convert(Matrix, ∂G_∂μ)
maximum(abs, Jμ_num-Jμ_exact)
@assert maximum(abs, Jμ_num-Jμ_exact)<1e-9  # This works

@time Jx_num = FiniteDiff.finite_difference_jacobian(x->G(μ0, x), x0_flat)
@time Jx_exact = convert(Matrix,∂G_∂x)
maximum(abs, Jx_num-Jx_exact)

@assert maximum(abs, Jx_num-Jx_exact)<1e-9

using Plots

pl1 = spy(abs.(Jx_num).>1e-8, title="Numerical")
pl2 = spy(abs.(Jx_exact).>1e-8, title="Exact")
pl3 = spy(abs.(Jx_exact - Jx_num).>1e-8, title="Diff")
plot(pl1,pl2,pl3)

# M = reshape(Jx_exact, 5, 10, 5, 20)
# MM = reshape(permutedims(M, [1, 2, 4, 3]), 50, 100)

pl1 = spy(abs.(Jx_num).>1e-8, title="Numerical")
pl2 = spy(abs.(Jx_exact).>1e-8, title="Exact")
pl3 = spy(abs.(Jx_exact - Jx_num).>1e-8, title="Diff")
plot(pl1,pl2,pl3)


text:
##-------------------------------------------##
##先定义形状喝流行阵列参数,然后在调用函数生成对应A
##-------------------------------------------##
##name = "线路输入 (LIS_MIIS_8X8_192K)"
export linear1,
        linear2,
        circle,
        A_matrix

abstract type ArrayShape end

struct linear1{T1<:Float64, T2<:Int} <: ArrayShape #一维直线阵列
    d::T1  #阵列间距
    n::T2  #阵列个数
end
struct linear2{T1<:Float64, T2<:Int} <: ArrayShape   #二维阵列
    d1::T1
    d2::T1 #阵列间距
    n1::T2 #长度阵列个数
    n2::T2 #宽度阵列个数
end
struct circle{T1<:Float64, T2<:Int} <:ArrayShape #圆形阵列
    n::T2  #阵列个数
    r::T1 #阵列半径
end
struct A_matrix{N,T<:Float64,F<:ArrayShape,T1<:NTuple{N, Float64},T2<:NTuple{N, Float64},I<:Int }
    angle::T1  #角度范围,单位为°
    length::T2  #声源波长,初始化是用Vector声明，不然无法修改数值
    step1::T    #角度步进值1,单位为弧度制
    shape::F    #排列形状
    dims::I     #空间角度划分维度
end

function (A::A_matrix)()
    if typeof(A.shape) == linear1
        result = linear1_methon(A)
    elseif typeof(A.shape) == linear2
        result = linear2_methon(A)
    elseif typeof(A.shape) == circle
        result = circle_methon(A)
    else
        error("shape shouldn't be $(typeof(A.shape))")
    end
    return result
end

#一维直线阵列只能用于测量一维空间角度
function linear1_methon(A::A_matrix)
    angle, bolength, step1, shape = (A.angle .*pi ./ 180), A.length[1], A.step1, A.shape
    d, n = shape.d, shape.n
    @assert A.dims == 1

    d = 0: d: (n-1)*d
    theta = range(angle[1], stop=angle[2], step=step1)
    Am = Array{Complex, 2}(undef, n, length(theta))
    @. Am = exp(-im*2*pi/bolength*d*theta')
    return @view Am[:,:] #返回结果
end

#适用于空间网格二维划分，并且对应空间二维ULA阵列形式
function linear2_methon(A::A_matrix)
    angle, bolength, step1, shape = (A.angle .*pi ./ 180), A.length[1], A.step1, A.shape
    d1, d2, n1, n2 = shape.d1, shape.d2, shape.n1, shape.n2

    @assert length(angle) == 4
    @assert length(step1) == 2

    theta = range(angle[1], step=step1[1], stop=angle[2])
    beta = range(angle[3], step=step1[2], stop=angle[4])
    l1, l2 = length(theta), length(beta)
    array_dim1 = 0: d1: (n1-1)*d1
    array_dim2 = 0: d2: (n2-1)*d2
    Am = Array{Complex,2}(undef, n1*n2, l1*l2)
    #生成d的序列,具体X，Y，Z轴的赋值与麦克风安装位置有关
    d_seq_buff = Array{Float64,2}(undef, n1*n2, 3)
    @views Threads.@threads for i in 1:length(array_dim1)
        @. d_seq_buff[(i-1)*n1+1 : (i-1)*n1+n2, 1] = array_dim1[i]
        d_seq_buff[(i-1)*n1+1 : (i-1)*n1+n2, 3] = array_dim2[1:end]
    end
    @. d_seq_buff[:, 2] = 0
    #生成theta序列
    Trigonometric_buff1 = Array{Float64, 2}(undef, 2, l1) #存储三角函数值
    @. Trigonometric_buff1[1, :] = sin(theta)
    @. Trigonometric_buff1[2, :] = cos(theta)
    Trigonometric_buff2 = Array{Float64, 2}(undef, 2, l2) #存储三角函数值
    @. Trigonometric_buff2[1, :] = sin(beta)
    @. Trigonometric_buff2[2, :] = cos(beta)
    #序列相乘得到Am
    @views Threads.@threads for i in 1: l1
        @. Am[:, (i-1)*l1+1:(i-1)*l1+l2] = (d_seq_buff[:, 1] * Trigonometric_buff2[2, :]' * Trigonometric_buff1[1, i])
                        +(d_seq_buff[:, 3] * Trigonometric_buff1[2, i])
    end

    @. Am = exp(-im * 2 * pi / bolength * Am)

    return @view Am[:, :]
end

#圆形阵列，对应空间二维划分
function circle_methon(A::A_matrix)
    angle, bolength, step1, shape = (A.angle .*pi ./ 180), A.length[1], A.step1, A.shape
    n, r = shape.n, shape.r
    @assert length(angle) == 4
    @assert length(step1) == 2

    theta = range(angle[1], step=step1[1], stop=angle[2])
    beta = range(angle[3], step=step1[2], stop=angle[4])
    l1, l2 = length(theta), length(beta)
    Am = Array{Complex, 2}(undef, n, l1*l2)
    #获得阵列位置序列
    d_seq_buff = Array{Float64, 2}(undef, n, 3)
    P_buff = (0:n-1) * 2 * pi / n
    @. d_seq_buff[1:n, 1] = sin(P_buff) * r
    @. d_seq_buff[1:n, 2] = 0
    @. d_seq_buff[1:n, 3] = cos(P_buff) * r
    #获得theta序列
    Trigonometric_buff1 = Array{Float64, 2}(undef, 2, l1) #存储三角函数值
    @. Trigonometric_buff1[1, :] = sin(theta)
    @. Trigonometric_buff1[2, :] = cos(theta)
    Trigonometric_buff2 = Array{Float64, 2}(undef, 2, l2) #存储三角函数值
    @. Trigonometric_buff2[1, :] = sin(beta)
    @. Trigonometric_buff2[2, :] = cos(beta)
    #序列相乘得到Am
    @views Threads.@threads for i in 1: l1
        @. Am[:, (i-1)*l1+1:(i-1)*l1+l2] = (d_seq_buff[:, 1] * Trigonometric_buff2[2, :]' * Trigonometric_buff1[1, i])
                        +(d_seq_buff[:, 3] * Trigonometric_buff1[2, i])
    end
    @. Am = exp(-im * 2 * pi / bolength * Am)

    return @view Am[:, :]
end


text:
# -----------------------------------------------------------------------------
# Implementation of Atomistic MolecularDynamicsSimulator API
# -----------------------------------------------------------------------------

"""
    NBSimulator <: MolecularDynamicsSimulator

A wrapper around NBodySimulator to implement the Atomistic API.

**Field descriptions**
- `Δt::Real` the time between timesteps, assumed to be in atomic units
- `steps::Integer` the number of timesteps for the simulation
- `t₀::Real` the starting time of the simulation, assumed to be in atomic units;
    defaults to 0
- `thermostat::Thermostat` the thermostat for the simulation;
    many options are defined by `NBodySimulator`, but a user could also define a custom thermostat;
    defaults to the `NullThermostat`
- `simulator::OrdinaryDiffEqAlgorithm` the algorithm to be used for the ODE;
    defaults to VelocityVerlet
- `potentials::Dict{Symbol,PotentialParameters}` dictionary of potentials;
    shouldn't be manipulated directly by the user
"""
@kwdef struct NBSimulator <: MolecularDynamicsSimulator
    Δt::Real       # in TIME_UNIT
    steps::Integer
    t₀::Real = 0.0 # in TIME_UNIT
    thermostat::Thermostat = NullThermostat()
    simulator::OrdinaryDiffEqAlgorithm = VelocityVerlet()
    potentials::Dict{Symbol,PotentialParameters} = Dict{Symbol,PotentialParameters}()
end
function NBSimulator(
    Δt::Unitful.Time,
    steps::Integer;
    t₀::Unitful.Time = 0.0 * TIME_UNIT,
    thermostat::Thermostat = NullThermostat(),
    simulator::OrdinaryDiffEqAlgorithm = VelocityVerlet(),
    potentials::Dict{Symbol,PotentialParameters} = Dict{Symbol,PotentialParameters}()
)
    NBSimulator(austrip(Δt), steps, austrip(t₀), thermostat, simulator, potentials)
end

# Extract the tuple of start_time, end_time from the simulator
time_range(simulator::NBSimulator) = (simulator.t₀, simulator.t₀ + simulator.steps * simulator.Δt)

function simulate(system::AbstractSystem, simulator::NBSimulator, potential::ArbitraryPotential)
    simulator.potentials[:custom] = InteratomicPotentialParameters(potential)
    simulate(system, simulator)
end
function simulate(system::AbstractSystem, simulator::NBSimulator, potential::LennardJonesParameters)
    simulator.potentials[:lennard_jones] = potential
    simulate(system, simulator)
end
function simulate(system::AbstractSystem, simulator::NBSimulator)
    nb_system = PotentialNBodySystem{ElementMassBody}(bodies(system), simulator.potentials)
    simulation = NBodySimulation(nb_system, time_range(simulator), nbs_boundary_conditions(system), simulator.thermostat, austrip(u"k"))
    NBSResult(run_simulation(simulation, simulator.simulator, dt = simulator.Δt))
end

# -----------------------------------------------------------------------------
# Integration with InteratomicPotentials
# -----------------------------------------------------------------------------

# Internal struct that represents a potential in the format accepted by NBodySimulator
# Wraps the underlying InteratomicPotential with a cache of the forces for the current timestep
struct InteratomicPotentialParameters <: PotentialParameters
    potential::ArbitraryPotential
    timestep_cache::RefValue{Real}
    force_cache::RefValue{Vector{SVector{3,Real}}}
    InteratomicPotentialParameters(potential::ArbitraryPotential) = new(potential, Ref{Real}(), Ref{Vector{SVector{3,Real}}}())
end

function NBodySimulator.get_accelerating_function(parameters::InteratomicPotentialParameters, simulation::NBodySimulation)
    bounding_box = ab_bounding_box(simulation.boundary_conditions)
    boundary_conditions = ab_boundary_conditions(simulation.boundary_conditions)
    symbols = get_atomic_symbols(simulation.system)
    numbers = get_atomic_numbers(simulation.system)
    masses = get_masses(simulation.system)
    (dv, u, v, t, i) -> begin
        if !isassigned(parameters.timestep_cache) || t != parameters.timestep_cache[]
            positions = [u[:, j] for j ∈ 1:size(u, 2)] * LENGTH_UNIT
            system = FastSystem(bounding_box, boundary_conditions, positions, symbols, numbers, masses * MASS_UNIT)
            wrapped_system = DynamicSystem(system, t * TIME_UNIT)
            parameters.timestep_cache[] = t
            parameters.force_cache[] = force(wrapped_system, parameters.potential)
        end
        dv .+= parameters.force_cache[][i] / masses[i]
    end
end

# -----------------------------------------------------------------------------
# Unitful Support for Built-in LennardJonesParameters
# -----------------------------------------------------------------------------

function NBodySimulator.LennardJonesParameters(ϵ::Unitful.Energy, σ::Unitful.Length, R::Unitful.Length)
    LennardJonesParameters(austrip(ϵ), austrip(σ), austrip(R))
end


text:
import Bcrypt

import Printf: @sprintf
using Test

@testset "test bcrypt is easy" begin
	pass = Array{UInt8,1}("mypassword")
	hashPass = Bcrypt.GenerateFromPassword(pass, 0)
	@test Bcrypt.CompareHashAndPassword(hashPass, pass)

	notPass = Array{UInt8,1}("notthepass")
	@test Bcrypt.CompareHashAndPassword(hashPass, notPass) == false
end

@testset "test bcrypt is correct" begin
	pass = Array{UInt8, 1}("allmine")
	salt = Array{UInt8, 1}("XajjQvNhvvRt5GSeFk1xFe")
	expect = Array{UInt8, 1}("\$2a\$10\$XajjQvNhvvRt5GSeFk1xFeyqRrsxkhBkUiQeg0dt.wU1qD4aFDcga")

	hsh = Bcrypt.bcrypt(pass[1:end], 10, salt[1:end])
	@test endswith(String(expect[1:end]), String(hsh[1:end]))

	h = Bcrypt.newFromHash(expect)
	@test expect == hash(h)
end

@testset "test very short passwords" begin
	key = UInt8[0x6b]
	salt = Array{UInt8, 1}("XajjQvNhvvRt5GSeFk1xFe")
	Bcrypt.bcrypt(key, 10, salt)
	@test true
end

@testset "test too long passwords work" begin
	salt = Array{UInt8, 1}("XajjQvNhvvRt5GSeFk1xFe")
	longPass = Array{UInt8, 1}("012345678901234567890123456789012345678901234567890123456")
	longExpect = Array{UInt8, 1}("\$2a\$10\$XajjQvNhvvRt5GSeFk1xFe5l47dONXg781AmZtd869sO8zfsHuw7C")
	hsh = Bcrypt.bcrypt(longPass[1:end], 10, salt[1:end])
	@test endswith(String(longExpect[1:end]), String(hsh[1:end]))
end

struct invalidHashTest
	err::Type{<:Exception}
	hsh::Array{UInt8, 1}
end

const invalidTests = invalidHashTest[
	invalidHashTest(Bcrypt.HashTooShort, Array{UInt8,1}("\$2a\$10\$fooo")),
	invalidHashTest(Bcrypt.HashTooShort, Array{UInt8,1}("\$2a")),
	invalidHashTest(Bcrypt.HashVersionTooNew, Array{UInt8,1}("\$3a\$10\$sssssssssssssssssssssshhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh")),
	invalidHashTest(Bcrypt.InvalidHashPrefix, Array{UInt8,1}("%2a\$10\$sssssssssssssssssssssshhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh")),
	invalidHashTest(Bcrypt.InvalidCost, Array{UInt8,1}("\$2a\$32\$sssssssssssssssssssssshhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh")),
]

@testset "test invalid hash exceptions" begin
	for case in invalidTests
		try
			Bcrypt.newFromHash(case.hsh)
		catch ex
			@test isa(ex, case.err)
		end
		try
			Bcrypt.CompareHashAndPassword(case.hsh, Array{UInt8,1}("rando"))
		catch ex
			@test isa(ex, case.err)
		end
	end
end

@testset "test unpadded b64 encoding" begin
	original = UInt8[101, 201, 101, 75, 19, 227, 199, 20, 239, 236, 133, 32, 30, 109, 243, 30]
	encodedOriginal = Array{UInt8, 1}("XajjQvNhvvRt5GSeFk1xFe")

	encoded = Bcrypt.bcrypt_b64_encode(original[1:end])
	@test encoded == encodedOriginal

	decoded = Bcrypt.bcrypt_b64_decode(encodedOriginal[1:end])
	@test decoded == original
end

@testset "test cost" begin
	suffix = "XajjQvNhvvRt5GSeFk1xFe5l47dONXg781AmZtd869sO8zfsHuw7C"
	for version in ["2a", "2"]
		for cost in [4, 10]
			s = @sprintf "\$%s\$%02d\$%s" version cost suffix
			h = Array{UInt8, 1}(s)
			actual = Bcrypt.Cost(h)
			@test actual == cost
		end
	end
	try
		Bcrypt.Cost(Array{UInt8, 1}("\$a\$a\$$suffix"))
	catch ex
		@test isa(ex, Bcrypt.HashTooShort)
	end
end


text:
include("writeMat.jl")

nxs = 10
SP, X = PolynomialRing(Nemo.ZZ,["X$i" for i in 1:nxs])

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       X[2]+X[6] X[10];
       X[7] X[1]+X[4];
       X[8] X[2]+X[4];
       X[9] X[1]+X[6];
       X[10] X[2]+X[6]];

toMaple(net, X, "M2.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       0 0;
       0 0;
       0 0];

toMaple(net, X, "M2nonrev0000.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       0 0;
       0 0;
       X[2]+X[6] X[10]];

toMaple(net, X, "M2nonrev0001.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       0 0;
       X[1]+X[6] X[9];
       0 0];

toMaple(net, X, "M2nonrev0010.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       0 0;
       X[1]+X[6] X[9];
       X[2]+X[6] X[10]];

toMaple(net, X, "M2nonrev0011.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       0 0;
       0 0];

toMaple(net, X, "M2nonrev0100.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       0 0;
       X[2]+X[6] X[10]];

toMaple(net, X, "M2nonrev0101.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       0 0];

toMaple(net, X, "M2nonrev0110.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       X[2]+X[6] X[10]];

toMaple(net, X, "M2nonrev0111.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       0 0;
       0 0];

toMaple(net, X, "M2nonrev1000.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       0 0;
       X[2]+X[6] X[10]];

toMaple(net, X, "M2nonrev1001.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       X[1]+X[6] X[9];
       0 0];

toMaple(net, X, "M2nonrev1010.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       X[1]+X[6] X[9];
       X[2]+X[6] X[10]];

toMaple(net, X, "M2nonrev1011.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       0 0;
       0 0];

toMaple(net, X, "M2nonrev1100.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       0 0;
       X[2]+X[6] X[10]];

toMaple(net, X, "M2nonrev1101.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       0 0];

toMaple(net, X, "M2nonrev1110.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       X[2]+X[6] X[10]];

toMaple(net, X, "M2nonrev1111.mpl")

max_stars_repo_name: LauraBMo/CRNTscripts

max_issues_repo_name: LauraBMo/CRNTscripts

max_forks_repo_name: LauraBMo/CRNTscripts


text:
using Distributions, Plots
using FluPredictibility.BioTools

export plot_all_trajectories
export pfix, pfix_v_freq, pfix_v_freq_positivederivative, fitness_plot, meanfreq

function plot_all_trajectories(ph; yearticks=false, label=true, lw=2.5)
	p = plot(size=(900,600))
	for z in ph
		X,Y,tmp,α = frequency_series(z)
		if yearticks
			X = [year(x) + month(x) /12. for x in X]
		end
		for a in 1:size(Y,2)
			plot!(p, X, Y[:,a], label=label ? "$(α[a])" : "", linewidth=lw)
		end
	end
	plot!(p, ylim=(-0.05,1.05))
	return p
end

function plot_all_trajectories(trajectories::Array{FrequencyTraj{A},1}; label=false, lw = 1, yearticks=false) where A
	p = plot(size=(900,600))
	dmin=missing; dmax=missing
	for traj in trajectories
		X = traj.date .+ traj.t
		if yearticks
			X = [year(x) + month(x) /12. for x in X]
			dmin = ismissing(dmin) ? minimum(X) : min(dmin, minimum(X))
			dmax = ismissing(dmax) ? maximum(X) : max(dmax, maximum(X))
		end
		Y = traj.freq
		lab = "$(traj.i) - $(traj.val)"
		if traj.fixation == :fixed
			st = (lw, :blue, 0.5)
		elseif traj.fixation == :lost
			st = (lw, :red, 0.5)
		else
			st = (lw, :black, :dashdot, 0.5)
		end
		plot!(p, X, Y, label = label ? lab : "", line=st)
	end
	plot!(p, frame=:box, xticks=floor(dmin):ceil(dmax))
	return p
end

plot_trajectory(traj::FrequencyTraj) = plot(traj.date .+ traj.t, traj.freq, label="$(traj.i) - $(traj.val)")
plot_trajectory!(p, traj::FrequencyTraj) = plot!(p, traj.date .+ traj.t, traj.freq, label="$(traj.i) - $(traj.val)")

function trajectory_freqbin(traj, alphabins)
	freqtraj_cf = Dict()
	for (α,dα) in alphabins
	    freqtraj_cf[α] = frequency_condition(traj, α, dα=dα)
	end
	return freqtraj_cf
end

function bernoulli_estimator(x,n)
	naivemean = x/n
	estmean = (x+1/2) / (n+1)
	α = x + 1/2
	β = n - x + 1/2
	P = Distributions.Beta(α,β)
	#
	lowerbound = 0.
	while cdf(P, lowerbound) < 0.05
		lowerbound += 0.01
	end
	#
	higherbound = 1.
	while cdf(P, higherbound) > 0.95
		higherbound -= 0.01
	end
	#
	return estmean, estmean-lowerbound, higherbound-estmean
end

# pfix(traj) = sum([x.fixation==:fixed for x in traj]) / (sum([x.fixation==:fixed for x in traj]) + sum([x.fixation==:lost for x in traj]));
meanfreq(traj) = mean(t.freq[t.index[:active]] for t in traj)
"""
"""
function pfix(trajectories)
	n = length(trajectories)
	y = count(t->t.fixation==:fixed, trajectories)
	yf, errdown, errup = bernoulli_estimator.(y,n)
	#
	xf = meanfreq(trajectories)
	err = (errdown, errup)
	return xf, yf, err
end

"""
	pfix_v_freq(ph, alphabins; v=false, effective=false, effective_kwargs=(clustering=:naive, distance=:strains))

Used for clustering trajectories. Needs `fp` as an argument since trajectory distance is by default based on common strains.
"""
function pfix_v_freq(ph, fp::FluPop, alphabins;
	v=false, effective=false, effective_kwargs=(clustering=:naive, distance=:strains), outn=false)
	trajectories = all_trajectories(ph, keep_unfinished=false)
	trajectories = previous_state_condition(trajectories, :lost)
	get_strains!(trajectories, fp)
	# Binning by frequency
	traj_fb = sort(OrderedDict(trajectory_freqbin(trajectories, alphabins)));
	# Keeping only trajectories that have a frequency backed by 20 strains at the time where it is binned.
	for (k,trajs) in traj_fb
	    traj_fb[k] = population_size_condition(trajs, 20, mode=:active)
	    if effective
	    	n0 = length(traj_fb[k])
			traj_fb[k] = effective_trajectories(traj_fb[k]; effective_kwargs...)
			v && println("Freq. bin $k: $n0 --> $(length(traj_fb[k])) effective trajectories")
		end
	end
	#
	# v && println("Based on a total of $(sum(length(v) for v in values(traj_fb))) trajectories")
	#
	n = [length(traj_fb[x]) for x in keys(traj_fb)] # For error bars
	x = [count(t->t.fixation==:fixed, traj_fb[x]) for x in keys(traj_fb)]
	out = bernoulli_estimator.(x,n)
	#
	xf = [meanfreq(traj_fb[x]) for x in keys(traj_fb)]
	yf = [x[1] for x in out];
	# yf = [count(t->t.fixation==:fixed, traj_fb[x])/length(traj_fb[x]) for x in keys(traj_fb)]
	errup = [x[3] for x in out]; errdown = [x[2] for x in out]
	#
	err = (errdown, errup)
	# return sort(collect(keys(traj_fb))), yf, err
	# println(xf)
	if outn
		return xf, yf, err, n
	else
		return xf, yf, err
	end
end
"""
	pfix_v_freq(ph, alphabins; v=false)

Usual function used for main text plots.
"""
function pfix_v_freq(ph, alphabins; v=false, outn=false)
	trajectories = all_trajectories(ph, keep_unfinished=false)
	trajectories = previous_state_condition(trajectories, :lost)
	# Binning by frequency
	traj_fb = sort(OrderedDict(trajectory_freqbin(trajectories, alphabins)));
	# Keeping only trajectories that have a frequency backed by 20 strains at the time where it is binned.
	for (k,trajs) in traj_fb
	    traj_fb[k] = population_size_condition(trajs, 20, mode=:active)
	end
	#
	# v && println("Based on a total of $(sum(length(v) for v in values(traj_fb))) trajectories")
	#
	n = [length(traj_fb[x]) for x in keys(traj_fb)] # For error bars
	x = [count(t->t.fixation==:fixed, traj_fb[x]) for x in keys(traj_fb)]
	out = bernoulli_estimator.(x,n)
	#
	xf = [meanfreq(traj_fb[x]) for x in keys(traj_fb)]
	yf = [x[1] for x in out];
	# yf = [count(t->t.fixation==:fixed, traj_fb[x])/length(traj_fb[x]) for x in keys(traj_fb)]
	errup = [x[3] for x in out]; errdown = [x[2] for x in out]
	#
	err = (errdown, errup)
	# return sort(collect(keys(traj_fb))), yf, err
	# println(xf)
	if outn
		return xf, yf, err, n
	else
		return xf, yf, err
	end
end

"""
	fitness_plot(trajectories, field, alphabins; verbose=false, dq = 0., alldat=false)

Trajectories higher or lower than median fitness.
"""
function fitness_plot(trajectories, field, alphabins; verbose=false, dq = 0., alldat=false)
    traj_fb = trajectory_freqbin(trajectories, alphabins);
    #
    for (k,v) in traj_fb
        traj_fb[k] = population_size_condition(v, 20, mode=:active)
    end
    # High and low fitnesses
    high_fit = Dict()
    low_fit = Dict()
    medfit = Dict()
    all_fit = Dict()
    for (k,v) in traj_fb
        fvalues = collect(skipmissing([x.data[field][x.index[:active]] for x in v]))
        medfit[k] = median(fvalues)
        # println(fvalues)
        verbose && println("Frequency $k -- median fitness $(medfit[k])")
        fitness(x) = x.data[field][x.index[:active]]
        all_fit[k] = v[findall(x->!ismissing(fitness(x)), v)]
        high_fit[k] = v[findall(x->!ismissing(fitness(x)) && fitness(x) > quantile(fvalues, 0.5 + dq), v)]
        low_fit[k] = v[findall(x->!ismissing(fitness(x)) && fitness(x) <= quantile(fvalues, 0.5 - dq), v)]
    end
    # Arrays
    # dat = vcat([reshape(collect(pfix(traj_fb[x[1]])), 1, 3) for x in alphabins]...)
    dat = vcat([reshape(collect(pfix(all_fit[x[1]])), 1, 3) for x in alphabins]...)
    dat_low = vcat([reshape(collect(pfix(low_fit[x[1]])), 1, 3) for x in alphabins]...)
    dat_high = vcat([reshape(collect(pfix(high_fit[x[1]])), 1, 3) for x in alphabins]...)
	#
    alldat && return dat, dat_low, dat_high, traj_fb, high_fit, low_fit, medfit
    return dat, dat_low, dat_high
end
"""
	fitness_plot(trajectories, field, alphabins, fitval; verbose=false alldat=false)

Trajectories higher or lower than `fitval` for field `field`.
"""
function fitness_plot(trajectories, field, alphabins, fitval; verbose=false, alldat=false)
    traj_fb = trajectory_freqbin(trajectories, alphabins);
    #
    for (k,v) in traj_fb
        traj_fb[k] = population_size_condition(v, 20, mode=:active)
    end
    # High and low fitnesses
    high_fit = Dict()
    low_fit = Dict()
    medfit = Dict()
    for (k,v) in traj_fb
        high_fit[k] = v[findall(x->x.data[field][x.index[:active]] > fitval, v)]
        low_fit[k] = v[findall(x->x.data[field][x.index[:active]] <= fitval, v)]
    end
    # Arrays
    dat = vcat([reshape(collect(pfix(traj_fb[x[1]])), 1, 3) for x in alphabins]...)
    dat_low = vcat([reshape(collect(pfix(low_fit[x[1]])), 1, 3) for x in alphabins]...)
    dat_high = vcat([reshape(collect(pfix(high_fit[x[1]])), 1, 3) for x in alphabins]...)
	#
    alldat && return dat, dat_low, dat_high, traj_fb, high_fit, low_fit, medfit
    return dat, dat_low, dat_high
end

"""
"""
function pfix_v_freq_positivederivative(ph, alphabins)
	trajectories = all_trajectories(ph, keep_unfinished=false)
	trajectories = previous_state_condition(trajectories, :lost)
	# Binning by frequency
	traj_fb = sort(OrderedDict(trajectory_freqbin(trajectories, alphabins)));
	# Keeping only trajectories that have a frequency backed by 50 strains at the time where it is binned.
	for (k,v) in traj_fb
	    traj_fb[k] = population_size_condition(v, 20, mode=:active)
	    traj_fb[k] = derivative_condition(traj_fb[k])
	end
	#
	n = [length(traj_fb[x]) for x in keys(traj_fb)] # For error bars
	x = [count(t->t.fixation==:fixed, traj_fb[x]) for x in keys(traj_fb)]
	out = bernoulli_estimator.(x,n)
	#
	xf = [meanfreq(traj_fb[x]) for x in keys(traj_fb)]
	yf = [x[1] for x in out];
	errup = [x[3] for x in out]; errdown = [x[2] for x in out]
	#
	err = (errdown, errup)
	return xf, yf, err
end

function mean_trajectory()

end

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


text:
"""
	dwdclimateraw(datein,siteid;url)
Read DWD CDC global monthly raw observations

**Input**
* datein: input date in DateTime format (can be a vector for multiple sites)
* siteid: WMO site ID, e.g. 87593 for La Plata, Argentina (one ID only)
* param: parameter to be extracted (one only), e.g. "Rx" (ftp://ftp-cdc.dwd.de/pub/CDC/observations_global/CLIMAT/monthly/raw/readme_RAW_CLIMATs_eng.txt)
* url: URL od local folder link
* downto: download temporary data to this folder (file will be deleted afterwards)

**Ouptut**
See ftp://ftp-cdc.dwd.de/pub/CDC/observations_global/CLIMAT/monthly/raw/readme_RAW_CLIMATs_eng.txt

**Example**
dwd = dwdclimateraw(DateTime(2017,10),87593,"Rx");
# Vector input
timein = collect(DateTime(2017,10):Dates.Month(1):DateTime(2017,12));
dwdvec = dwdclimateraw(timein,87593,"Rx");
"""
function dwdclimateraw(datein,siteid::Int,param::String;
	url::String="ftp://ftp-cdc.dwd.de/pub/CDC/observations_global/CLIMAT/monthly/raw/",
	downto::String="")::DataFrame
	dfo = DataFrame(datetime=datein,paramout=Dates.value.(datein).*0);
	for (i,v) in enumerate(typeof(datein)==DateTime ? [datein] : datein)
		filein = joinpath(url,@sprintf("CLIMAT_RAW_%s.txt",Dates.format(v,"yyyymm")));
		if isdir(url)
			fileout = filein;
		else
			fileout = isdir(downto) ? joinpath(downto,"dwd_month_raw_TEMP.txt") : "dwd_month_raw_TEMP.txt"
			download(filein,fileout);
		end
		d,h = readdlm(fileout,';',header=true)
		r,c = dwdclimateraw_indices(d,h,param,siteid);
		dfo[:paramout][i] = eltype(d[r,c])==Int ? d[r,c] : NA;
		!isdir(url) ? rm(fileout) : nothing
	end
	return rename!(dfo,:paramout=>Symbol(param));
end

"""
Auxiliary function to find indices of the site and parameter
"""
function dwdclimateraw_indices(datain::Array{Any,2},headerin::Array{AbstractString,2},
								paramin::String,sitein::Int)::Tuple{Int,Int}
	r,c = 1,1;
	for i in headerin
		i == paramin ? break : c += 1
	end
	for i in datain[:,3]
		i == sitein ? break : r += 1
	end
	return r,c
end

"""
    readgssd(file_in)
Read global surface summary of day data
See https://www1.ncdc.noaa.gov/pub/data/gsod/readme.txt

**Input**
* `file_in`: input file name (full)

**Output**
* DataFrame containin all columns of the input file

**Example**
```
file_in = joinpath(dirname(@__DIR__),"test","input","nndc_climate_cdo.txt");
df = readgssd(file_in);
```
"""
function readgssd(file_in)::DataFrame
    # see https://www1.ncdc.noaa.gov/pub/data/gsod/readme.txt
    par = [:STN,:WBAN,:YEARMODA,:TEMP,:Count_TEMP,:DEWP,:Count_DEWP,:SLP,:Count_SLP,:STP,
        :Count_STP,:VISIB,:Count_VISIB,:WDSP,:Count_WDSP,:MXSPD,:GUST,:MAX,:Flag_MAX,:MIN,
        :Flag_MIN,:PRCP,:Flag_PRCP,:SNDP,:FRSHTT]
    i1 = [1,8,15,25,32,36,43,47,54,58,65,69,75,79,85,89,96,103,109,111,
            117,119,124,126,133];
    i2 = [6,12,22,30,33,41,44,52,55,63,66,73,76,83,86,93,100,108,109,116,
            117,123,124,130,138];
    typei = [Int,Int,Int,Float64,Int,Float64,Int,Float64,Int,Float64,
                Int,Float64,Int,Float64,Int,Float64,Float64,Float64,String,
                Float64,String,Float64,String,Float64,Int];
    function parse_col(s::String,i1,i2,typei)::Vector{Any}
        out_vec = Vector{Any}(undef,0);
        for i in 1:length(i1)
            push!(out_vec,
                typei[i]!=String ? Base.parse(typei[i],s[i1[i]:i2[i]]) : s[i1[i]:i2[i]])
        end
        return out_vec
    end

    data_out = DataFrame([Vector{x}(undef,0) for x in typei],par);
    open(file_in,"r") do fid
        headtxt = readline(fid);
        while !eof(fid)
            push!(data_out,parse_col(readline(fid),i1,i2,typei));
        end
    end
    return data_out;
end

"""
    convertgssd(df)
Convert global surface summary of day data to SI units and set flags to NaN

**Input**
* `df`: readgssd output

**Output**
* dataframe with converted units and only essential columns (for me :-)
* temperature in degrees C, windspeed in m/s, pressure in hPa, precipitation in mm

**Example**
```
file_in = joinpath(dirname(@__DIR__),"test","input","nndc_climate_cdo.txt");
df = readgssd(file_in) |> convertgssd
```
"""
function convertgssd(df::DataFrame)::DataFrame
    data_out = DataFrame(datetime=FileTools.pattern2time.(df[:YEARMODA],"day"));
    d = copy(df);
    # flag => NaN;
    flags = [9999.9,9999.9,9999.9,9999.9,999.9,999.9,99.99,9999.9];
    for (i,v) in enumerate([:TEMP,:DEWP,:MAX,:MIN,:WDSP,:MXSPD,:PRCP,:STP])
        d[v][d[v].==flags[i]] .= NaN;
    end
    # temperature
    for i in [:TEMP,:DEWP,:MAX,:MIN]
        data_out[i] = (d[i].-32).*5/9;
    end
    # wind speed
    for i in [:WDSP,:MXSPD]
        data_out[i] = d[i]./10.0.*0.514444444;
    end
    # precipitation (mm)
    data_out[:PRCP] = d[:PRCP].*25.4
    data_out[:STP] = d[:STP];
    return data_out;
end


text:
module Parse

using ExprTools
using Markdown
using Pkg
using CRC32c
using MatchCore
using ..Comonicon
using ..Comonicon.Types
using ..Comonicon.CodeGen

export @cast,
    @main,
    read_doc,
    command,
    rm_lineinfo,
    default_name,
    get_version,
    iscached,
    cachefile,
    create_cache,
    enable_cache,
    disable_cache

include("cast.jl")
include("markdown.jl")
include("runtime.jl")
include("utils.jl")

end


text:
struct ODEBacksolveSensitivityFunction{rateType,uType,uType2,UF,PF,G,JC,GC,DG,TJ,PJT,PJC,CP,SType,CV} <: SensitivityFunction
  uf::UF
  pf::PF
  g::G
  J::TJ
  pJ::PJT
  dg_val::uType
  jac_config::JC
  g_grad_config::GC
  paramjac_config::PJC
  sensealg::BacksolveAdjoint
  f_cache::rateType
  discrete::Bool
  y::uType2
  sol::SType
  dg::DG
  checkpoints::CP
  colorvec::CV
end

@noinline function ODEBacksolveSensitivityFunction(g,u0,p,sensealg,discrete,sol,dg,checkpoints,tspan,colorvec)
  numparams = length(p)
  numindvar = length(u0)
  # if there is an analytical Jacobian provided, we are not going to do automatic `jac*vec`
  f = sol.prob.f
  isautojacvec = DiffEqBase.has_jac(f) ? false : get_jacvec(sensealg)
  J = isautojacvec ? nothing : similar(u0, numindvar, numindvar)

  if !discrete
    if dg != nothing || isautojacvec
      pg = nothing
      pg_config = nothing
    else
      pg = UGradientWrapper(g,tspan[2],p)
      pg_config = build_grad_config(sensealg,pg,u0,p)
    end
  else
    pg = nothing
    pg_config = nothing
  end

  if DiffEqBase.has_jac(f) || isautojacvec
    jac_config = nothing
    uf = nothing
  else
    uf = DiffEqDiffTools.UJacobianWrapper(f,tspan[2],p)
    jac_config = build_jac_config(sensealg,uf,u0)
  end

  y = copy(sol(tspan[1])) # TODO: Has to start at interpolation value!

  if DiffEqBase.has_paramjac(f) || isautojacvec
    paramjac_config = nothing
    pf = nothing
  else
    pf = DiffEqDiffTools.ParamJacobianWrapper(f,tspan[1],y)
    paramjac_config = build_param_jac_config(sensealg,pf,y,p)
  end

  pJ = isautojacvec ? nothing : similar(sol.prob.u0, numindvar, numparams)

  dg_val = similar(u0, numindvar) # number of funcs size
  f_cache = deepcopy(u0)

  return ODEBacksolveSensitivityFunction(uf,pf,pg,J,pJ,dg_val,
                               jac_config,pg_config,paramjac_config,
                               sensealg,f_cache,
                               discrete,y,sol,dg,checkpoints,colorvec)
end

# u = λ'
function (S::ODEBacksolveSensitivityFunction)(du,u,p,t)
  @unpack y, sol, J, uf, sensealg, f_cache, jac_config, discrete, dg, dg_val, g, g_grad_config = S
  idx = length(y)
  f = sol.prob.f
  isautojacvec = DiffEqBase.has_jac(f) ? false : get_jacvec(sensealg)

  λ     = @view u[1:idx]
  dλ    = @view du[1:idx]
  grad  = @view u[idx+1:end-idx]
  dgrad = @view du[idx+1:end-idx]
  _y    = @view u[end-idx+1:end]
  dy    = @view du[end-idx+1:end]
  copyto!(vec(y), _y)
  isautojacvec || f(dy, _y, p, t)

  if !isautojacvec
    if DiffEqBase.has_jac(f)
      f.jac(J,y,p,t) # Calculate the Jacobian into J
    else
      uf.t = t
      jacobian!(J, uf, y, f_cache, sensealg, jac_config)
    end
    mul!(dλ',λ',J)
  else
    _dy, back = Tracker.forward(y, sol.prob.p) do u, p
      if DiffEqBase.isinplace(sol.prob)
        out_ = map(zero, u)
        f(out_, u, p, t)
        Tracker.collect(out_)
      else
        vec(f(u, p, t))
      end
    end
    dλ[:], dgrad[:] = Tracker.data.(back(λ))
    dy[:] = vec(Tracker.data(_dy))
  end

  dλ .*= -one(eltype(λ))

  if !discrete
    if dg != nothing
      dg(dg_val,y,p,t)
    else
      g.t = t
      gradient!(dg_val, g, y, sensealg, g_grad_config)
    end
    dλ .+= dg_val
  end

  if !isautojacvec
    @unpack pJ, pf, paramjac_config = S
    if DiffEqBase.has_paramjac(f)
      f.paramjac(pJ,y,sol.prob.p,t) # Calculate the parameter Jacobian into pJ
    else
      jacobian!(pJ, pf, sol.prob.p, f_cache, sensealg, paramjac_config)
    end
    mul!(dgrad',λ',pJ)
  end
  nothing
end

# g is either g(t,u,p) or discrete g(t,u,i)
@noinline function ODEAdjointProblem(sol,sensesensealg::BacksolveAdjoint,
                                     g,t=nothing,dg=nothing;
                                     checkpoints=sol.t,
                                     callback=CallbackSet())
  f = sol.prob.f
  tspan = (sol.prob.tspan[2],sol.prob.tspan[1])
  discrete = t != nothing

  p = sol.prob.p
  p === DiffEqBase.NullParameters() && error("Your model does not have parameters, and thus it is impossible to calculate the derivative of the solution with respect to the parameters. Your model must have parameters to use parameter sensitivity calculations!")

  u0 = zero(sol.prob.u0)

  len = length(u0)+length(p)
  λ = similar(u0, len)
  sense = ODEBacksolveSensitivityFunction(g,u0,
                                        p,sensesensealg,discrete,
                                        sol,dg,checkpoints,tspan,f.colorvec)

  init_cb = t !== nothing && sol.prob.tspan[2] == t[end]
  cb = generate_callbacks(sense, g, λ, t, callback, init_cb)
  z0 = [vec(zero(λ)); vec(sense.y)]
  ODEProblem(sense,z0,tspan,p,callback=cb)
end


text:
#-------------------------------------------------------------------------------
# File: constants.jl
# Description:
# Date: October 31, 2019
# Author: Jonathan Fontaine, Killian Fretaud, Rémi Garcia,
#         Boualem Lamraoui, Benoît Le Badezet, Benoit Loger
#-------------------------------------------------------------------------------

const WEIGHTS_OBJECTIVE_FUNCTION = [1000, 1, 1000000]

# TIME
const TIME_LIMIT = 60.0
const TIME_PART_ILS = 0.9
const TIME_PART_VNS_LPRC = 0.96
const TIME_PART_VNS_PCC = 1 # Unused for now

# ILS-HPRC
const ALPHA_ILS = 25
const BETA_ILS = 50
const STOPPING_CRITERIA_ILS_HPRC = 3
const NBCAR_PERTURBATION = 5
const NBCAR_DIVERSIFICATION = 50

# VNS-LPRC
const VNS_LPRC_MAX_NON_IMPROVEMENT = 150
const VNS_LPRC_MIN_INSERT = 3
const VNS_LPRC_MAX_INSERT = 8
const VNS_LPRC_MIN_EXCHANGE = 5
const VNS_LPRC_MAX_EXCHANGE = 12

# VNS-PCC
const VNS_PCC_MIN_INSERT = 2
const VNS_PCC_MAX_INSERT = 8
const VNS_PCC_MIN_EXCHANGE = 2
const VNS_PCC_MAX_EXCHANGE = 10
const VNS_PCC_MINMAX = [(VNS_PCC_MIN_INSERT, VNS_PCC_MAX_INSERT), (VNS_PCC_MIN_EXCHANGE, VNS_PCC_MAX_EXCHANGE)]

# Files names
const INSTANCES = Dict{String, Array{String, 1}}(
    "A" => [
        "022_3_4_EP_RAF_ENP",
        "022_3_4_RAF_EP_ENP",
        "024_38_3_EP_ENP_RAF",
        "024_38_3_EP_RAF_ENP",
        "024_38_5_EP_ENP_RAF",
        "024_38_5_EP_RAF_ENP",
        "025_38_1_EP_ENP_RAF",
        "025_38_1_EP_RAF_ENP",
        "039_38_4_EP_RAF_ch1",
        "039_38_4_RAF_EP_ch1",
        "048_39_1_EP_ENP_RAF",
        "048_39_1_EP_RAF_ENP",
        "064_38_2_EP_RAF_ENP_ch1",
        "064_38_2_EP_RAF_ENP_ch2",
        "064_38_2_RAF_EP_ENP_ch1",
        "064_38_2_RAF_EP_ENP_ch2"
    ],
    "B" => [
        "022_EP_ENP_RAF_S22_J1",
        "022_EP_RAF_ENP_S22_J1",
        "022_RAF_EP_ENP_S22_J1",
        "023_EP_ENP_RAF_S23_J3",
        "023_EP_RAF_ENP_S23_J3",
        "023_RAF_EP_ENP_S23_J3",
        "024_V2_EP_ENP_RAF_S22_J1",
        "024_V2_EP_RAF_ENP_S22_J1",
        "024_V2_RAF_EP_ENP_S22_J1",
        "025_EP_ENP_RAF_S22_J3",
        "025_EP_RAF_ENP_S22_J3",
        "025_RAF_EP_ENP_S22_J3",
        "028_ch1_EP_ENP_RAF_S22_J2",
        "028_ch1_EP_RAF_ENP_S22_J2",
        "028_ch1_RAF_EP_ENP_S22_J2",
        "028_ch2_EP_ENP_RAF_S23_J3",
        "028_ch2_EP_RAF_ENP_S23_J3",
        "028_ch2_RAF_EP_ENP_S23_J3",
        "029_EP_ENP_RAF_S21_J6",
        "029_EP_RAF_ENP_S21_J6",
        "029_RAF_EP_ENP_S21_J6",
        "035_ch1_EP_ENP_RAF_S22_J3",
        "035_ch1_EP_RAF_ENP_S22_J3",
        "035_ch1_RAF_EP_ENP_S22_J3",
        "035_ch2_EP_ENP_RAF_S22_J3",
        "035_ch2_EP_RAF_ENP_S22_J3",
        "035_ch2_RAF_EP_ENP_S22_J3",
        "039_ch1_EP_ENP_RAF_S22_J4",
        "039_ch1_EP_RAF_ENP_S22_J4",
        "039_ch1_RAF_EP_ENP_S22_J4",
        "039_ch3_EP_ENP_RAF_S22_J4",
        "039_ch3_EP_RAF_ENP_S22_J4",
        "039_ch3_RAF_EP_ENP_S22_J4",
        "048_ch1_EP_ENP_RAF_S22_J3",
        "048_ch1_EP_RAF_ENP_S22_J3",
        "048_ch1_RAF_EP_ENP_S22_J3",
        "048_ch2_EP_ENP_RAF_S22_J3",
        "048_ch2_EP_RAF_ENP_S22_J3",
        "048_ch2_RAF_EP_ENP_S22_J3",
        "064_ch1_EP_ENP_RAF_S22_J3",
        "064_ch1_EP_RAF_ENP_S22_J3",
        "064_ch1_RAF_EP_ENP_S22_J3",
        "064_ch2_EP_ENP_RAF_S22_J4",
        "064_ch2_EP_RAF_ENP_S22_J4",
        "064_ch2_RAF_EP_ENP_S22_J4"
    ],
    "X" => [
        "022_RAF_EP_ENP_S49_J2",
        "023_EP_RAF_ENP_S49_J2",
        "024_EP_RAF_ENP_S49_J2",
        "025_EP_ENP_RAF_S49_J1",
        "028_CH1_EP_ENP_RAF_S50_J4",
        "028_CH2_EP_ENP_RAF_S51_J1",
        "029_EP_RAF_ENP_S49_J5",
        "034_VP_EP_RAF_ENP_S51_J1_J2_J3",
        "034_VU_EP_RAF_ENP_S51_J1_J2_J3",
        "035_CH1_RAF_EP_S50_J4",
        "035_CH2_RAF_EP_S50_J4",
        "039_CH1_EP_RAF_ENP_S49_J1",
        "039_CH3_EP_RAF_ENP_S49_J1",
        "048_CH1_EP_RAF_ENP_S50_J4",
        "048_CH2_EP_RAF_ENP_S49_J5",
        "064_CH1_EP_RAF_ENP_S49_J1",
        "064_CH2_EP_RAF_ENP_S49_J4",
        "655_CH1_EP_RAF_ENP_S51_J2_J3_J4",
        "655_CH2_EP_RAF_ENP_S52_J1_J2_S01_J1"
    ]
)


text:
struct NNRODE{C,O} <: NeuralNetDiffEqAlgorithm
    chain::C
    opt::O
end
NNRODE(chain;opt=Flux.ADAM(0.1)) = NNRODE(chain,opt)

function DiffEqBase.solve(
    prob::DiffEqBase.AbstractRODEProblem,
    alg::NeuralNetDiffEqAlgorithm,
    args...;
    dt,
    timeseries_errors = true,
    save_everystep=true,
    adaptive=false,
    abstol = 1f-6,
    verbose = false,
    maxiters = 100)

    DiffEqBase.isinplace(prob) && error("Only out-of-place methods are allowed!")

    u0 = prob.u0
    tspan = prob.tspan
    f = prob.f
    p = prob.p
    t0 = tspan[1]

    #hidden layer
    chain  = alg.chain
    opt    = alg.opt
    ps     = Flux.params(chain)
    data   = Iterators.repeated((), maxiters)

    #train points generation
    ts = tspan[1]:dt:tspan[2]

    #The phi trial solution
    ϕ(t,W) = u0 .- (t.-tspan[1]).*v(Tracker.collect([t,W]))

    dfdx = (t,W) -> Tracker.gradient((t,W) -> sum(ϕ(t,W)), t,W; nest = true)[1]
    loss = () -> sum(abs2,sum(abs2,dfdx(t,W).-f(ϕ(t,W)[1],p,t,W)[1]) for (t,W) in zip(ts,randn(length(ts))))

    cb = function ()
        l = loss()
        verbose && println("Current loss is: $l")
        l < abstol && Flux.stop()
    end
    Flux.train!(loss, ps, data, opt; cb = cb)

    u = [ϕ(t,W)[1].data for (t,W) in zip(ts,randn(length(ts)))]
    u
end #solve


text:
include("./produce_data.jl")
using NPZ

# generate all data sets for paper
function main()
    h=1e-7
    ms=[1.09, 1.10, 1.20]
    for m in ms
        println("m=",m)
        data=produce_data_box_scaling_branching_network(m, h)
        npzwrite(@sprintf("%s/data/data_box_scaling_m%06.4f_h%.2e.npz", path_project,m,h),data)
    end
end

# call function
main()


text:
@testset "distributions" begin
    for mtype in (:dgmesh,:frapmesh,:fvmesh)
        for closing_vec in ( :normalised_closing_operator,
                :naive_normalised_closing_operator)#, :unnormalised_closing_operator)

            # pdfs
            f(x,i) = (i-1)/12.0./sum(1:3)
            msh = @eval $mtype
            dq = DiscretisedFluidQueue(am,msh)
            d = SFMDistribution(f,dq)#point_mass(11.990,2,am,$i)#
            if mtype!=:fvmesh
                @test sum(d)≈1.0
            end
            f_rec = (mtype!==:frapmesh) ? pdf(d) : pdf(d,@eval($(Symbol(closing_vec,"_pdf"))))
            x = (0.01:0.4:11.99)'
            if mtype!=:frapmesh
                @test all(isapprox.( f_rec.(x,1:4)-f.(x,1:4), 0, atol=sqrt(eps()) ))
            end

            # cdfs
            cdf_rec = (mtype!==:frapmesh) ? cdf(d) : cdf(d,@eval($(Symbol(closing_vec,"_cdf"))))
            f_cdf(x,i) = f(x,i)*x
            @test sum(cdf_rec.(msh.nodes[end],1:4))≈1.0
            @test sum(cdf_rec.(msh.nodes[end]+1.0,1:4))≈sum(f_cdf.(msh.nodes[end],1:4))
            x = (-1:0.4:12.99)'
            @test all(isapprox.( f_cdf.(x,1:4), cdf_rec.(x,1:4), atol=5e-2 ))

            # point masses at boundaries
            d = d * 0.5
            d[1] = 0.25
            d[end] = 0.25
            f_cdf_2(x,i) = 0.25*(i∈(2))*(x>=msh.nodes[1])+0.5*f_cdf(x,i)*(x>=msh.nodes[1])+0.25*(x>=msh.nodes[end])*(i∈(3))
            cdf_rec_2 = (mtype!==:frapmesh) ? cdf(d) : cdf(d,@eval($(Symbol(closing_vec,"_cdf"))))
            @test sum(cdf_rec_2.(msh.nodes[end],1:4))≈1.0
            @test sum(cdf_rec_2.(msh.nodes[end]+1.0,1:4))≈sum(f_cdf.(msh.nodes[end],1:4))
            @test all(isapprox.( f_cdf_2.(x,1:4), cdf_rec_2.(x,1:4), atol=5e-2 ))

            @test_throws DomainError left_point_mass(1,dq)
            tst = zeros(n_phases(am)*n_bases_per_phase(msh)+N₋(am.S)+N₊(am.S))
            tst[1] = 1.0
            @test left_point_mass(2,dq)==tst

            @test_throws DomainError right_point_mass(2,dq)
            tst2 = zeros(n_phases(am)*n_bases_per_phase(msh)+N₋(am.S)+N₊(am.S))
            tst2[end] = 1.0
            @test right_point_mass(3,dq)==tst2
        end
        # test point masses not at boundaries...
    end
end


text:
# Load modules
using JuMP
using SpineModel
using Cbc

# Custon constraints required for case study A4
function add_constraints(m::Model)

    # Cyclic storage bounds on the first solve
    @fetch stor_state = m.ext[:variables]
    cons = m.ext[:constraints][:stor_cyclic] = Dict()
    filtered_storage__commodity = filter(stor_c -> stor_c.storage in SpineModel.indices(cyclic), storage__commodity())
    stor_start = [first(stor_state_indices(storage=stor, commodity=c)) for (stor, c) in filtered_storage__commodity]
    stor_end = [last(stor_state_indices(storage=stor, commodity=c)) for (stor, c) in filtered_storage__commodity]
    for (stor, c, t_first) in stor_start
        for (stor, c, t_last) in stor_end
            (stor, c, t_first, t_last) in keys(cons) && continue
            cons[stor, c, t_first, t_last] = @constraint(
                m,
                stor_state[stor, c, t_first]
                ==
                stor_state[stor, c, t_last]
            )
        end
    end

    # THIS ISN'T NECESSARY FOR THE BUILDING MODEL!
    # Input-output conversion with online consumption
    @fetch flow, units_on = m.ext[:variables]
    cons = m.ext[:constraints][:conversion] = Dict()
    for (u, c1, c2) in SpineModel.indices(idle_heat_rate)
        for t in SpineModel.t_lowest_resolution(map(x -> x.t, flow_indices(unit=u, commodity=[c1, c2])))
            cons[u, c1, c2, t] = @constraint(
                m,
                + reduce(
                    +,
                    flow[u_, n, c1_, d, t_] * SpineModel.duration(t_)
                    for (u_, n, c1_, d, t_) in flow_indices(
                        unit=u, commodity=c1, direction=direction(:from_node), t=t_in_t(t_long=t)
                    );
                    init=0
                )
                ==
                + variable_heat_rate(unit=u, commodity1=c1, commodity2=c2, t=t)
                * reduce(
                    +,
                    flow[u_, n, c2_, d, t_] * SpineModel.duration(t_)
                    for (u_, n, c2_, d, t_) in flow_indices(
                        unit=u, commodity=c2, direction=direction(:to_node), t=t_in_t(t_long=t)
                    );
                    init=0
                )
                + idle_heat_rate(unit=u, commodity1=c1, commodity2=c2, t=t)
                #* unit_capacity(unit=u, commodity=c2, direction=direction(:to_node))
                * reduce(
                    +,
                    units_on[u_, t_] * SpineModel.duration(t_)
                    for (u_, t_) in units_on_indices(unit=u, t=t_in_t(t_long=t));
                    init=0
                )
            )
        end
    end

end

function update_constraints(m::Model)

    # Remove cyclic storage bounds from subsequent solves
    cons = pop!(m.ext[:constraints], :stor_cyclic, nothing)
    cons === nothing && return
    delete.(m, values(cons))

    # Update input-output convetrsion equation
    @fetch flow, units_on = m.ext[:variables]
    cons = m.ext[:constraints][:conversion]
    for (u, c1, c2) in SpineModel.indices(idle_heat_rate)
        for t in SpineModel.t_lowest_resolution(map(x -> x.t, flow_indices(unit=u, commodity=[c1, c2])))
            # Update idle heat rates
            set_normalized_coefficient(
                cons[u, c1, c2, t],
                units_on[u,t],
                (
                    - idle_heat_rate(unit=u, commodity1=c1, commodity2=c2, t=t)
                    #* unit_capacity(unit=u, commodity=c2, direction=direction(:to_node))
                    * SpineModel.duration(t)
                )
            )
            # Update variable heat rates
            for (u_, n, c2_, d, t_) in flow_indices(unit=u, commodity=c2, direction=direction(:to_node), t=t_in_t(t_long=t))
                set_normalized_coefficient(
                    cons[u_, c1, c2_, t_],
                    flow[u_, n, c2_, d, t_],
                    - variable_heat_rate(unit=u_, commodity1=c1, commodity2=c2_, t=t) * SpineModel.duration(t_)
                )
            end
        end
    end

end

# Catch an error if no db_url is provided
try
    ARGS[1]
catch e
    if isa(e, BoundsError)
        println("!!! Database url required as a command line argument !!!")
    end
end

# Run the model from the chosen database
m = run_spinemodel(
    ARGS[1];
    with_optimizer=optimizer_with_attributes(Cbc.Optimizer, "logLevel"=>1, "allowableGap"=>0, "ratioGap"=>0),
    add_constraints=m->add_constraints(m),
    update_constraints=m->update_constraints(m),
)


text:
function quoteid(conn, str)
    if conn.quoteidentifierchar == '\0'
        conn.quoteidentifierchar = API.getinfostring(conn.dbc, API.SQL_IDENTIFIER_QUOTE_CHAR)[1]
    end
    # avoid double quoting
    if str[1] == conn.quoteidentifierchar && str[end] == conn.quoteidentifierchar
        return str
    else
        return string(conn.quoteidentifierchar, str, conn.quoteidentifierchar)
    end
end

sqltype(conn, ::Type{Union{T, Missing}}) where {T} = sqltype(conn, T)

function sqltype(conn, T)
    if isempty(conn.types)
        types = Tables.columntable(Cursor(API.gettypes(conn.dbc)))
        conn.alltypes = types
        i = findfirst(==(API.SQL_VARCHAR), types.DATA_TYPE)
        if i === nothing
            defaultT = "VARCHAR(255)"
        else
            defaultT = types.TYPE_NAME[i] * "($(types.COLUMN_SIZE[i]))"
        end
        for jlT in BINDTYPES
            _, sqlT = bindtypes(jlT)
            i = findfirst(==(sqlT), types.DATA_TYPE)
            nm = i !== nothing ? types.TYPE_NAME[i] : defaultT
            if i !== nothing && types.CREATE_PARAMS[i] !== missing && nm != "DOUBLE" && nm != "FLOAT"
                nm *= occursin(',', types.CREATE_PARAMS[i]) ? "($(typeprecision(jlT)),$(typescale(jlT)))" : "($(types.COLUMN_SIZE[i]))"
            end
            conn.types[jlT] = nm
        end
    end
    return conn.types[T]
end

checkdupnames(names) = length(unique(map(x->lowercase(String(x)), names))) == length(names) || error("duplicate case-insensitive column names detected; sqlite doesn't allow duplicate column names and treats them case insensitive")

function createtable(conn::Connection, nm::AbstractString, sch::Tables.Schema; debug::Bool=false, quoteidentifiers::Bool=true, createtableclause::AbstractString="CREATE TABLE", columnsuffix=Dict())
    names = sch.names
    checkdupnames(names)
    types = [sqltype(conn, T) for T in sch.types]
    columns = (string(quoteidentifiers ? quoteid(conn, String(names[i])) : names[i], ' ', types[i], ' ', get(columnsuffix, names[i], "")) for i = 1:length(names))
    debug && @info "executing create table statement: `$createtableclause $nm ($(join(columns, ", ")))`"
    return DBInterface.execute(conn, "$createtableclause $nm ($(join(columns, ", ")))")
end

"""
    ODBC.load(table, conn, name; append=true, quoteidentifiers=true, limit=typemax(Int64), createtableclause=nothing, columnsuffix=Dict(), debug=false)
    table |> ODBC.load(conn, name; append=true, quoteidentifiers=true, limit=typemax(Int64), createtableclause=nothing, columnsuffix=Dict(), debug=false)

Attempts to take a Tables.jl source `table` and load into the database represented by `conn` with table name `name`.

It first detects the `Tables.Schema` of the table source and generates a `CREATE TABLE` statement
with the appropriate column names and types. If no table name is provided, one will be autogenerated, like `odbcjl_xxxxx`.
The `CREATE TABLE` clause can be provided manually by passing the `createtableclause` keyword argument, which
would allow specifying a temporary table or `if not exists`.
Column definitions can also be enhanced by providing arguments to `columnsuffix` as a `Dict` of
column name (given as a `Symbol`) to a string of the enhancement that will come after name and type like
`[column name] [column type] enhancements`. This allows, for example, specifying the charset of a string column
by doing something like `columnsuffix=Dict(:Name => "CHARACTER SET utf8mb4")`.

Do note that databases vary wildly in requirements for `CREATE TABLE` and column definitions
so it can be extremely difficult to load data generically. You may just need to tweak some of the provided
keyword arguments, but you may also need to execute the `CREATE TABLE` and `INSERT` statements
yourself. If you run into issues, you can [open an issue](https://github.com/JuliaDatabases/ODBC.jl/issues) and
we can see if there's something we can do to make it easier to use this function.
"""
function load end

load(conn::Connection, table::AbstractString="odbcjl_"*Random.randstring(5); kw...) = x->load(x, conn, table; kw...)

function load(itr, conn::Connection, name::AbstractString="odbcjl_"*Random.randstring(5); append::Bool=true, quoteidentifiers::Bool=true, debug::Bool=false, limit::Integer=typemax(Int64), kw...)
    # get data
    rows = Tables.rows(itr)
    sch = Tables.schema(rows)
    if sch === nothing
        # we want to ensure we always have a schema, so materialize if needed
        rows = Tables.rows(columntable(rows))
        sch = Tables.schema(rows)
    end
    # ensure table exists
    if quoteidentifiers
        name = quoteid(conn, name)
    end
    try
        createtable(conn, name, sch; quoteidentifiers=quoteidentifiers, debug=debug, kw...)
    catch e
        @warn "error creating table" (e, catch_backtrace())
    end
    if !append
        DBInterface.execute(conn, "DELETE FROM $name")
    end
    # start a transaction for inserting rows
    transaction(conn) do
        params = chop(repeat("?,", length(sch.names)))
        stmt = DBInterface.prepare(conn, "INSERT INTO $name VALUES ($params)")
        for (i, row) in enumerate(rows)
	    i > limit && break
            debug && @info "inserting row $i; $(Tables.Row(row))"
            DBInterface.execute(stmt, Tables.Row(row); debug=debug)
        end
    end

    return name
end

function transaction(f::Function, conn)
    API.setcommitmode(conn.dbc, false)
    try
        f()
        API.endtran(conn.dbc, true)
    catch
        API.endtran(conn.dbc, false)
        rethrow()
    finally
        API.setcommitmode(conn.dbc, true)
    end
end


text:
abstract type AbstractCMPS{T,N} end

scalartype(ψ::AbstractCMPS{T}) where T = scalartype(T)

function leftreducedoperator(op::FieldOperator, Ψ::AbstractCMPS, ρL = nothing; kwargs...)
    if isnothing(ρL)
        ρL = leftenv(Ψ; kwargs...)
    end
    _leftreducedoperator(op, Ψ.Q, Ψ.Rs, ρL)
end
function rightreducedoperator(op::FieldOperator, Ψ::AbstractCMPS, ρR = nothing; kwargs...)
    if isnothing(ρR)
        ρR = rightenv(Ψ; kwargs...)
    end
    _rightreducedoperator(op, Ψ.Q, Ψ.Rs, ρR)
end

function _leftreducedoperator(op::FieldOperator, Q, Rs, ρL)
    A, B = _ketbrafactors(op, Q, Rs)
    return B'*(ρL*A)
end

function _rightreducedoperator(op::FieldOperator, Q, Rs, ρR)
    A, B = _ketbrafactors(op, Q, Rs)
    return (A*ρR)*B'
end

function _expval(op::FieldOperator, Q, Rs, ρL, ρR)
    A, B = _ketbrafactors(op, Q, Rs)
    return localdot(ρL*B, A*ρR)
end

function leftreducedoperator(ops::LocalOperator, Ψ::AbstractCMPS, ρL = nothing; kwargs...)
    if isnothing(ρL)
        ρL, = leftenv(Ψ; kwargs...)
    end
    hL = zero(ρL)
    for (coeff, op) in zip(coefficients(ops), operators(ops))
        if coeff isa Number
            axpy!(coeff, leftreducedoperator(op, Ψ, ρL), hL)
        else
            mul!(hL, coeff, leftreducedoperator(op, Ψ, ρL), true, true)
        end
    end
    return hL
end

function rightreducedoperator(ops::LocalOperator, Ψ::AbstractCMPS, ρR = nothing; kwargs...)
    if isnothing(ρR)
        ρR, = rightenv(Ψ; kwargs...)
    end
    hR = zero(ρR)
    for (coeff, op) in zip(coefficients(ops), operators(ops))
        if coeff isa Number
            axpy!(coeff, rightreducedoperator(op, Ψ, ρR), hR)
        else
            mul!(hR, coeff, rightreducedoperator(op, Ψ, ρR), true, true)
        end
    end
    return hR
end

function expval(ops::LocalOperator, Ψ::AbstractCMPS, ρL = nothing, ρR = nothing; kwargs...)
    if isnothing(ρL)
        ρL, = leftenv(Ψ; kwargs...)
    end
    if isnothing(ρR)
        ρR, = rightenv(Ψ; kwargs...)
    end
    (a, b) = domain(Ψ)
    if !isfinite(a)
        a = zero(a)
    end
    if !isfinite(b)
        b = oneunit(b)
    end
    ZL = localdot(ρL, ρR)
    Z = real(ZL(a))
    Zb = ZL(b)
    Z ≈ Zb || @warn "error in computing normalisation: Za = $Z, Zb = $Zb"
    evs = _expval.(operators(ops), Ref(Ψ.Q), Ref(Ψ.Rs), Ref(ρL), Ref(ρR))
    ev = sum(coefficients(ops) .* evs)
    return ev/Z
end

function expval(H::LocalHamiltonian, Ψ::AbstractCMPS, ρL = nothing, ρR = nothing; kwargs...)
    @assert domain(H) == domain(Ψ)
    return integrate(expval(H.h, Ψ, ρL, ρR; kwargs...), domain(H))
end


text:
"""model definition of line2d"""

struct StateLine2D <: AbsState
    x::Float64
    y::Float64
    theta::Float64
end

function get_mid_status(p::StateLine2D, q::StateLine2D)::StateLine2D
    return StateLine2D(
        (p.x + q.x) / 2,
        (p.y + q.y) / 2,
        diff_angles(q.theta, p.theta) / 2 + p.theta,
    )
end

function dist(q1::StateLine2D, q2::StateLine2D)::Float64
    # normalize
    return norm([q1.x - q2.x, q1.y - q2.y, diff_angles(q1.theta, q2.theta) / π])
end

function gen_connect(
    q::StateLine2D,  # to identify type
    obstacles::Vector{CircleObstacle2D},
    rads::Vector{Float64};
    step_dist::Float64 = STEP_DIST,
    max_dist::Union{Nothing,Float64} = nothing,
)::Function

    # check: q \in C_free
    f(q::StateLine2D, i::Int64)::Bool = begin
        a = [q.x, q.y]
        b = [cos(q.theta), sin(q.theta)] * rads[i] + a

        if any(x -> (x < rads[i] || 1 - rads[i] < x), vcat(a, b))
            return false
        end

        if any(o -> dist(a, b, [o.x, o.y]) < o.r, obstacles)
            return false
        end

        return true
    end

    f(q_from::StateLine2D, q_to::StateLine2D, i::Int64)::Bool = begin
        D = dist(q_from, q_to)
        !isnothing(max_dist) && D > max_dist && return false

        dt = diff_angles(q_to.theta, q_from.theta)
        for e in vcat(collect(0:step_dist:D) / D, 1.0)
            # intermediate: root
            a = (1 - e) * [q_from.x, q_from.y] + e * [q_to.x, q_to.y]
            # intermediate: theta
            t = q_from.theta + e * dt
            # intermediate: tip
            b = a + rads[i] * [cos(t), sin(t)]

            # outside
            any(x -> (x < 0 || 1 < x), vcat(a, b)) && return false

            # obstacles
            any(o -> dist(a, b, [o.x, o.y]) < o.r, obstacles) && return false
        end

        return true
    end

    return f
end

function gen_collide(
    q::StateLine2D,
    rads::Vector{Float64};
    step_dist::Float64 = STEP_DIST,
    safety_dist::Float64 = SAFETY_DIST_LINE,
)::Function

    N = length(rads)

    f(
        q_i_from::StateLine2D,
        q_i_to::StateLine2D,
        q_j_from::StateLine2D,
        q_j_to::StateLine2D,
        i::Int64,
        j::Int64,
    ) = begin
        # check each pair of step
        D_i = dist(q_i_from, q_i_to)
        D_j = dist(q_j_from, q_j_to)

        dt_i = diff_angles(q_i_to.theta, q_i_from.theta)
        dt_j = diff_angles(q_j_to.theta, q_j_from.theta)

        for e_i in vcat(collect(0:step_dist:D_i) / D_i, 1.0)
            # root, angle, tip
            a_i = (1 - e_i) * [q_i_from.x, q_i_from.y] + e_i * [q_i_to.x, q_i_to.y]
            t_i = q_i_from.theta + e_i * dt_i
            b_i = rads[i] * [cos(t_i), sin(t_i)] + a_i

            for e_j in vcat(collect(0:step_dist:D_j) / D_j, 1.0)
                # root, angle, tip
                a_j = (1 - e_j) * [q_j_from.x, q_j_from.y] + e_j * [q_j_to.x, q_j_to.y]
                t_j = q_j_from.theta + e_j * dt_j
                b_j = rads[j] * [cos(t_j), sin(t_j)] + a_j

                # check collision
                if dist(a_i, b_i, a_j, b_j) < safety_dist
                    return true
                end
            end
        end
        return false
    end

    f(q_i::StateLine2D, q_j::StateLine2D, i::Int64, j::Int64) = begin
        a_i = [q_i.x, q_i.y]
        b_i = [cos(q_i.theta), sin(q_i.theta)] * rads[i] + a_i
        a_j = [q_j.x, q_j.y]
        b_j = [cos(q_j.theta), sin(q_j.theta)] * rads[j] + a_j
        return dist(a_i, b_i, a_j, b_j) < safety_dist
    end

    f(Q_from::Vector{Node{StateLine2D}}, Q_to::Vector{Node{StateLine2D}}) = begin
        for i = 1:N, j = i+1:N
            if f(Q_from[i].q, Q_to[i].q, Q_from[j].q, Q_to[j].q, i, j)
                return true
            end
        end
        return false
    end

    f(Q::Vector{Node{StateLine2D}}, q_i_to::StateLine2D, i::Int64) = begin
        q_i_from = Q[i].q
        D_i = dist(q_i_from, q_i_to)

        dt_i = diff_angles(q_i_to.theta, q_i_from.theta)

        for e_i in vcat(collect(0:step_dist:D_i) / D_i, 1.0)
            # root, angle, tip
            a_i = (1 - e_i) * [q_i_from.x, q_i_from.y] + e_i * [q_i_to.x, q_i_to.y]
            t_i = q_i_from.theta + e_i * dt_i
            b_i = rads[i] * [cos(t_i), sin(t_i)] + a_i

            for j in filter(k -> k != i, 1:N)
                q_j = Q[j].q
                a_j = [q_j.x, q_j.y]
                b_j = rads[j] * [cos(q_j.theta), sin(q_j.theta)] + a_j

                # check collision
                if dist(a_i, b_i, a_j, b_j) < safety_dist
                    return true
                end
            end
        end

        return false
    end

    return f
end

function gen_uniform_sampling(q::StateLine2D)::Function
    () -> begin
        return StateLine2D(rand(2)..., rand() * 2π)
    end
end

function plot_motion!(q_from::StateLine2D, q_to::StateLine2D, rad::Float64, params)
    plot!([q_from.x, q_to.x], [q_from.y, q_to.y]; params...)
end

function plot_start_goal!(q_init::StateLine2D, q_goal::StateLine2D, rad::Float64, params)
    plot!([q_init.x], [q_init.y]; markershape = :hex, params...)
    plot!([q_goal.x], [q_goal.y]; markershape = :star, params...)
end

function plot_agent!(q::StateLine2D, rad::Float64, color::String)
    plot!(
        [q.x, rad * cos(q.theta) + q.x],
        [q.y, rad * sin(q.theta) + q.y],
        color = color,
        lw = 5,
        label = nothing,
    )
end

function gen_random_instance_StateLine2D(; params...)
    return gen_random_instance(StateLine2D(0, 0, 0); params...)
end


text:
# Compound Layers with combination of primitive, recurrent and convolutional layers
"""
MLP is a bunch of stacked Dense Layers, i.e. MLP(;input::Int, hiddens::Array{Int}, output::Int)
"""
struct MLP
    layers::Tuple{Vararg{Layer}}
end

function MLP(;hiddens::Array, input::Int, output::Int, o...)
    ls = Layer[]
    for h in hiddens
        l = Dense(;input=input, output=h, activation=relu)
        input = h
        push!(ls, l)
    end
    push!(ls, Dense(;input=input, output=output)) # no need non-linear activity for final layer
    return MLP(Tuple(ls))
end

function (m::MLP)(x)
    for l in m.layers
        x = l(x)
    end
    return x
end


text:
include("vanderpol.jl")

using DynamicPolynomials, SemialgebraicSets, MathOptInterfaceMosek

# TODO: do we return the triple (𝑃, 𝑂, 𝑂_XFZ18) or rather a problem instance 𝑃
# generic options 𝑂 and a continuous post solver instance with the given specific
# options, XFZ18(𝑂_XFZ18)

function vanderpol_XFZ18(; k=6)
    @polyvar x₁ x₂

    # define the set of initial states X₀ = {x: V₀(x) <= 0}
    V₀ = (x₁ - 1.4)^2 + (x₂ - 2.4)^2  - 0.15
    X0 = @set V₀ <= 0

    (𝑃, 𝑂) = vanderpol(X0=X0, variables = (x₁, x₂))

    # algorithm-specific options
    𝑂_XFZ18 = Options()

    # constraints Y = {x: g(x) >= 0} compact search space Y x [0, T]
    g = 25 - x₁^2 - x₂^2
    𝑂_XFZ18[:search_space] = @set g >= 0

    # degree of the relaxation
    𝑂_XFZ18[:relaxation_degree] = k

    # define the optimization solver
    𝑂_XFZ18[:solver] = MosekOptimizer

    return (𝑃, 𝑂, 𝑂_XFZ18)
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
export Vground, head_wind, cross_wind

"""
    Vground(Vtas::Float64, Vwind::Float64, ∠wind::Float64, course::Float64)

Return the ground speed `Vground` [m/s] given the `course` [deg], airspeed
`Vtas` [m/s], wind speed `Vwind`[m/s] and wind direction (from) `∠wind` [deg].

source: edwilliams.org/avform.htm
"""
function Vground(Vtas::Float64, Vwind::Float64, ∠wind::Float64, course::Float64)
    swc = (Vwind / Vtas) * sind(∠wind - course)
    if abs(swc) > 1.0
        return NaN
    end
    return Vtas * √(1-swc^2) - Vwind * cosd(∠wind - course)
end

"""
    head_wind(Vwind::Float64, ∠wind::Float64, course::Float64)

Return the `head_wind` [m/s] component for a given wind speed `Vwind` [m/s],
`course` [deg], and wind direction (from) `∠wind` [deg]. A tail-wind is a
negative head-wind.

source: edwilliams.org/avform.htm
"""
function head_wind(Vwind::Float64, ∠wind::Float64, course::Float64)
    return Vwind*cosd(∠wind - course)
end

"""
    cross_wind(Vwind::Float64, ∠wind::Float64, course::Float64)

Return the `cross_wind` [m/s] component for a given wind speed `Vwind` [m/s],
`course` [deg], and wind direction (from) `∠wind` [deg]. A positive cross-wind
component indicates a wind from the right.

source: edwilliams.org/avform.htm
"""
function cross_wind(Vwind::Float64, ∠wind::Float64, course::Float64)
    return Vwind*sind(∠wind - course)
end


text:
using EarthDataLab
using Test
using Dates
using IterTools
using WeightedOnlineStats
using Tables

c=Cube()
d = subsetcube(c,variable=["air_temperature_2m", "gross_primary_productivity", "land_surface_temperature"],lon=(30,31),lat=(50,51),
              time=(Date("2002-01-01"),Date("2008-12-31")))
mytable = CubeTable(value=d, expandaxes=("var",))

mytable2 = CubeTable(data=d)

@testset "cubefittable and WeightedCovMatrix fittable" begin
    covmCube = cubefittable(mytable, WeightedCovMatrix, :value, weight=(x->cosd(x.lat)),showprog=false)
    covmVal = WeightedCovMatrix()
    for t in mytable
        for row in Tables.rows(t)
            obs = row.value
            obslat = row.lat
            if !any(ismissing.(obs))
                fit!(covmVal, obs, cosd(obslat))
            end
        end
    end
    @test all(isapprox.(covmCube.data, value(covmVal)))
end

@testset "remaining functions in tablestats.jl" begin
    meanCube = cubefittable(mytable2, WeightedMean, :data, weight=(x->cosd(x.lat)),
                            by=(:Variable,))
    weightedTableMeans = fittable(mytable2, WeightedMean, :data, by=(:Variable,),
                                    weight=(x -> cosd(x.lat)))
    foreach(enumerate(["air_temperature_2m", "gross_primary_productivity", "land_surface_temperature"])) do (i,v)
      @test isapprox(value(weightedTableMeans)[(v,)], meanCube.data[i])
    end
end


text:
# # 215: 2D Nonlinear Poisson with boundary reaction
# ([source code](SOURCE_URL))

module Example215_NonlinearPoisson2D_BoundaryReaction

using Printf
using VoronoiFVM
using ExtendableGrids
using GridVisualize
using ExtendableSparse

function main(;n=10,Plotter=nothing,verbose=false, unknown_storage=:sparse,tend=100,max_lureuse=0,
              factorization=LUFactorization())
    h=1.0/convert(Float64,n)
    X=collect(0.0:h:1.0)
    Y=collect(0.0:h:1.0)

    grid=VoronoiFVM.Grid(X,Y)

    eps=1.0e-2
    physics=VoronoiFVM.Physics(
        num_species=2,
        breaction=function(f,u,node)
        if  node.region==2
            f[1]=1*(u[1]-u[2])
            f[2]=1*(u[2]-u[1])
        else
            f[1]=0
            f[2]=0
        end
        end,

    flux=function(f,u0,edge)
        u=unknowns(edge,u0)
        f[1]=eps*(u[1,1]-u[1,2])
        f[2]=eps*(u[2,1]-u[2,2])
    end,

    storage=function(f,u,node)
        f[1]=u[1]
        f[2]=u[2]
    end)

    sys=VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)
    enable_species!(sys,1,[1])
    enable_species!(sys,2,[1])

    inival=unknowns(sys)
    U=unknowns(sys)
    inival[1,:].=map((x,y)->exp(-5.0*((x-0.5)^2+(y-0.5)^2)),grid)
    inival[2,:].=0

    control=VoronoiFVM.NewtonControl()
    control.verbose=verbose
    control.tol_linear=1.0e-5
    control.max_lureuse=max_lureuse
    control.factorization=factorization

    tstep=0.01
    time=0.0
    istep=0
    u25=0

    p=GridVisualizer(Plotter=Plotter,layout=(2,1))
    while time<tend
        time=time+tstep
        solve!(U,inival,sys,control=control,tstep=tstep)
        inival.=U
        if verbose
            @printf("time=%g\n",time)
        end
        I=integrate(sys,physics.storage,U)
        Uall=sum(I)
        tstep*=1.2
        istep=istep+1
        u25=U[25]
        scalarplot!(p[1,1],grid,U[1,:],title=@sprintf("U1: %.3g U1+U2:%8.3g",I[1,1],Uall),flimits=(0,1))
        scalarplot!(p[2,1],grid,U[2,:],title=@sprintf("U2: %.3g",I[2,1]),flimits=(0,1))
        reveal(p)
    end
    return u25
end

function test()
    testval=0.2760603343272377
    main(unknown_storage=:sparse) ≈ testval &&
        main(unknown_storage=:dense) ≈ testval
end
end


text:
import Dates

function calcage(input_day, input_month, input_value)
    """
        Method to to either calculate a person's age or year of birth.
    """
    if typeof(input_day) != Int64
        input_day = parse(Int, input_day)
    end
    if typeof(input_month) != Int64
        input_month = parse(Int, input_month)
    end
    if typeof(input_value) != Int64
        input_value = parse(Int, input_value)
    end

    today = split(string(Dates.today()), "-")
    current_year = parse(Int, today[1])
    current_month = parse(Int, today[2])
    current_day = parse(Int, today[3])

    value = current_year - input_value
    if current_month < input_month
        value -= 1
    end
    if current_month == input_month
        if current_day < input_day
            value -= 1
        end
    end

    return value
end

max_stars_repo_name: urbanware-org/snippets

max_issues_repo_name: urbanware-org/snippets

max_forks_repo_name: urbanware-org/snippets


text:
abstract type AbstractCycle end

""
@inline function _descend!(problem, level, cycle_t)
    coarsen!(problem, level)
    _cycle!(problem, level, cycle_t)
    uncoarsen!(problem, level)
    return nothing
end

@inline function _cycle!(problem, level, cycle_t)
    if doinitial(problem, level)
        initial!(problem, level)
        process_initial!(problem, level)
    else
        _core_cycle!(problem, level, cycle_t)
    end
    return nothing
end

function _cycle!(problem, level, cycle_t, numcycle)
    for i in 1:numcycle
        precycle!(problem, level)
        _cycle!(problem, level, cycle_t)
        postcycle!(problem, level)
    end
    return nothing
end

"""
"""
function cycle!(problem, level, cycle_t = VCycle(), numcycle = 1)
    preprocessing!(problem, level)
    _cycle!(problem, level, cycle_t, numcycle)
    postprocessing!(problem, level)
    return nothing
end


text:
using YAML

type SettingsYAMLFile <: SettingsSource
    data::Dict
    filename::AbstractString
    watched::Bool
    file_watcher::Task

    function SettingsYAMLFile(filename::AbstractString, watched::Bool=true)
        data = YAML.load(open(filename))

        # Now get our settings source object so we can update the data
        # asynchroniously if desired
        settings_source = new(data, filename, watched, @async begin end)

        # Now set up the real watcher
        if watched
            settings_source.file_watcher = @async begin
                while settings_source.watched
                    name, status = watch_file(settings_source.filename)
                    if status.changed
                        settings_source.data = YAML.load(open(settings_source.filename))
                    end
                end
            end
        end

        return settings_source
    end
end

Base.getindex(item::SettingsYAMLFile, key) = item.data[key]

Base.keys(item::SettingsYAMLFile) = keys(item.data)


text:
using Documenter, ModiaLang

makedocs(
  #modules  = [Modia],
  sitename = "ModiaLang",
  authors  = "Hilding Elmqvist (Mogram) and Martin Otter (DLR-SR)",
  format = Documenter.HTML(prettyurls = false),
  pages    = [
     "Home"      => "index.md",
  ]
)

max_stars_repo_name: ModiaSim/ModiaLang

max_issues_repo_name: ModiaSim/ModiaLang

max_forks_repo_name: ModiaSim/ModiaLang


text:
# code from https://github.com/Gnimuc/Nuklear.jl/blob/master/gen/generator.jl

using Clang

const h3_include_dir = normpath(@__DIR__, "..", "deps", "usr", "include", "h3")
((root, dirs, files),) = walkdir(h3_include_dir)
headers = normpath.(root, files)

# create a work context
ctx = DefaultContext()

# parse headers
parse_headers!(ctx, headers, includes=[LLVM_INCLUDE])

# settings
ctx.libname = "libh3"
ctx.options["is_function_strictly_typed"] = false
ctx.options["is_struct_mutable"] = false  # for nested struct

# write output
api_file = normpath(@__DIR__, "libh3_functions.jl")
api_stream = open(api_file, "w")
for trans_unit in ctx.trans_units
    root_cursor = getcursor(trans_unit)
    push!(ctx.cursor_stack, root_cursor)
    header = spelling(root_cursor)
    @info "wrapping header: $header ..."
    # loop over all of the child cursors and wrap them, if appropriate.
    ctx.children = children(root_cursor)
    for (i, child) in enumerate(ctx.children)
        child_name = name(child)
        child_header = filename(child)
        ctx.children_index = i
        # choose which cursor to wrap
        startswith(child_name, "__") && continue  # skip compiler definitions
        child_name in keys(ctx.common_buffer) && continue  # already wrapped
        child_header != header && continue  # skip if cursor filename is not in the headers to be wrapped

        wrap!(ctx, child)
    end
    @info "writing $(api_file)"
    println(api_stream, "# Julia wrapper for header: $header")
    println(api_stream, "# Automatically generated using Clang.jl")
    print_buffer(api_stream, ctx.api_buffer)
    println(api_stream)
    empty!(ctx.api_buffer)  # clean up api_buffer for the next header
end
close(api_stream)

common_file = normpath(@__DIR__, "libh3_types.jl")
expr_buffer = dump_to_buffer(ctx.common_buffer)
io = IOBuffer()
print_buffer(io, expr_buffer)
body = String(take!(io))
open(common_file, "w") do f
    println(f, """
# Automatically generated using Clang.jl

using CEnum
""")
    for line in split(body, '\n')
        !startswith(line, "#") && occursin("MASK_NEGATIVE", line) && print(f, "# ")
        println(f, line)
    end
end


text:
module TestExamples2DMHD

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2d_dgsem")

@testset "MHD" begin
  @trixi_testset "elixir_mhd_alfven_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [0.00011149543672225127, 5.888242524520296e-6, 5.888242524510072e-6, 8.476931432519067e-6, 1.3160738644036652e-6, 1.2542675002588144e-6, 1.2542675002747718e-6, 1.8705223407238346e-6, 4.651717010670585e-7],
      linf = [0.00026806333988971254, 1.6278838272418272e-5, 1.627883827305665e-5, 2.7551183488072617e-5, 5.457878055614707e-6, 8.130129322880819e-6, 8.130129322769797e-6, 1.2406302192291552e-5, 2.373765544951732e-6])
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl with flux_derigs_etal" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [1.7201098719531215e-6, 8.692057393373005e-7, 8.69205739320643e-7, 1.2726508184718958e-6, 1.040607127595208e-6, 1.07029565814218e-6, 1.0702956581404748e-6, 1.3291748105236525e-6, 4.6172239295786824e-7],
      linf = [9.865325754310206e-6, 7.352074675170961e-6, 7.352074674185638e-6, 1.0675656902672803e-5, 5.112498347226158e-6, 7.789533065905019e-6, 7.789533065905019e-6, 1.0933531593274037e-5, 2.340244047768378e-6],
      volume_flux = (flux_derigs_etal, flux_nonconservative_powell))
  end

  @trixi_testset "elixir_mhd_alfven_wave_mortar.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave_mortar.jl"),
      l2   = [3.7762324533854616e-6, 1.5534623833573546e-6, 1.4577234868196855e-6, 1.7647724628707057e-6, 1.4831911814574333e-6, 1.456369119716533e-6, 1.4115666913995062e-6, 1.804758237422838e-6, 8.320469738087189e-7],
      linf = [3.670661330201774e-5, 1.530289442645827e-5, 1.3592183785327006e-5, 1.5173897443654383e-5, 9.43771379136038e-6, 1.0906323046233624e-5, 1.0603954940346938e-5, 1.5900499596113726e-5, 5.978772247650426e-6],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_mhd_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
      l2   = [0.03637302248881514, 0.043002991956758996, 0.042987505670836056, 0.02574718055258975, 0.1621856170457943, 0.01745369341302589, 0.017454552320664566, 0.026873190440613117, 5.336243933079389e-16],
      linf = [0.23623816236321427, 0.3137152204179957, 0.30378397831730597, 0.21500228807094865, 0.9042495730546518, 0.09398098096581875, 0.09470282020962917, 0.15277253978297378, 4.307694418935709e-15])
  end

  @trixi_testset "elixir_mhd_orszag_tang.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_orszag_tang.jl"),
      l2   = [0.21967600768935716, 0.2643126515795721, 0.31488287201980875, 0.0, 0.5160141621186931, 0.23028914748088603, 0.34413527376463915, 0.0, 0.003178793090381426],
      linf = [1.2749969218080568, 0.6737013368774057, 0.8604154399895696, 0.0, 2.799342099887639, 0.6473347557712643, 0.9691773375490476, 0.0, 0.05729832038724348],
      tspan = (0.0, 0.09))
  end

  @trixi_testset "elixir_mhd_orszag_tang.jl with flux_hll" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_orszag_tang.jl"),
      l2   = [0.10806619664693064, 0.20199136742199922, 0.22984589847526207, 0.0, 0.29950152196422647, 0.15688413207147794, 0.24293641543490646, 0.0, 0.003246181006326598],
      linf = [0.560316034595759, 0.5095520363866776, 0.6536748458764621, 0.0, 0.9627447086204038, 0.3981375420906146, 0.673472146198816, 0.0, 0.04879208429337193],
      tspan = (0.0, 0.06), surface_flux = (flux_hll, flux_nonconservative_powell))
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl one step with initial_condition_constant" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [7.144325530681224e-17, 2.123397983547417e-16, 5.061138912500049e-16, 3.6588423152083e-17, 8.449816179702522e-15, 3.9171737639099993e-16, 2.445565690318772e-16, 3.6588423152083e-17, 9.971153407737885e-17],
      linf = [2.220446049250313e-16, 8.465450562766819e-16, 1.8318679906315083e-15, 1.1102230246251565e-16, 1.4210854715202004e-14, 8.881784197001252e-16, 4.440892098500626e-16, 1.1102230246251565e-16, 4.779017148551244e-16],
      maxiters = 1,
      initial_condition = initial_condition_constant,
      atol = 2.0e-13)
  end

  @trixi_testset "elixir_mhd_rotor.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"),
      l2   = [1.2623319195262743, 1.8273050553090515, 1.7004151198284634, 0.0, 2.2978570581460818, 0.2147235065899803, 0.23558337696054493, 0.0, 0.0032515115395693483],
      linf = [11.003677581472843, 14.70614192714736, 15.687648666952708, 0.0, 17.098104835553823, 1.3283750501377847, 1.4365828094434892, 0.0, 0.07886241196068537],
      tspan = (0.0, 0.05))
  end

  @trixi_testset "elixir_mhd_blast_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_blast_wave.jl"),
      l2   = [0.17646728395490927, 3.866230215339417, 2.4867304651291255, 0.0, 355.4562971958441, 2.359493623565687, 1.4030741420730297, 0.0, 0.029613599942667133],
      linf = [1.581630420824181, 44.15725488910748, 13.056964982196554, 0.0, 2244.875490238186, 13.07679044647926, 9.14612176426092, 0.0, 0.5154756722488522],
      tspan = (0.0, 0.003),
      # Calling the AnalysisCallback before iteration 9 causes the interpolation
      # of this IC to have negative density/pressure values, crashing the simulation.
      coverage_override = (maxiters=9,))
  end
end

end # module