text:
struct EOS_LNG <: MultiFluidModel
    components::Vector{String}
    icomponents::UnitRange{Int}
    properties::MultiFluidPropertyParam
    ideal::MultiFluidIdealParam
    single::MultiFluidSingleParam
    pair::MultiFluidPairParam
    references::Vector{String}
end

@registermodel EOS_LNG

function EOS_LNG(components::Vector{String})
    I = getnames_gerg2008(components)

    params = getparams_eos_lng()
    Mw = each_split_model(params[:Mw],I)
    rhoc = each_split_model(params[:rhoc],I)
    vc = SingleParam("critical volume (dm3/mol)",rhoc.components,1 ./ rhoc.values)
    Tc = each_split_model(params[:Tc],I)
    pc = each_split_model(params[:pc],I)
    lb_v = each_split_model(params[:lb_v],I)
    properties = MultiFluidPropertyParam(Mw,rhoc,vc,Tc,pc,lb_v)
    n0 = each_split_model(params[:n0],I)
    theta = each_split_model(params[:theta],I)
    gamma_T = each_split_model(params[:gamma_T],I)
    gamma_v = each_split_model(params[:gamma_v],I)
    beta_T = each_split_model(params[:beta_T],I)
    beta_v = each_split_model(params[:beta_v],I)
    ideal = MultiFluidIdealParam(n0,theta,gamma_T,gamma_v,beta_T,beta_v)
    n = each_split_model(params[:n],I) |> pack_vectors
    t = each_split_model(params[:t],I) |> pack_vectors
    d = each_split_model(params[:d],I) |> pack_vectors
    c = each_split_model(params[:c],I) |> pack_vectors
    single = MultiFluidSingleParam(n,t,d,c)
    nij = each_split_model(params[:nij],I) |> pack_vectors
    tij = each_split_model(params[:tij],I) |> pack_vectors
    dij = each_split_model(params[:dij],I) |> pack_vectors
    Fij = each_split_model(params[:Fij],I)
    eta_ij = each_split_model(params[:eta_ij],I) |> pack_vectors
    epsilon_ij = each_split_model(params[:epsilon_ij],I) |> pack_vectors
    beta_ij = each_split_model(params[:beta_ij],I) |> pack_vectors
    gamma_ij = each_split_model(params[:gamma_ij],I) |> pack_vectors
    pair = MultiFluidPairParam(nij,tij,dij,Fij,beta_ij,gamma_ij,eta_ij,epsilon_ij)
    icomponents = 1:length(I)
    references = ["10.1021/je300655b"]
    return EOS_LNG(components,icomponents,properties,ideal,single,pair,references)
end

export EOS_LNG
#nij tij dij ηij εij βij γij
#C1,C4

function getparams_eos_lng()
    params = getparams_gerg2008()
    γ_T =params[:gamma_T].values
    γ_v = params[:gamma_v].values
    β_T = params[:beta_T].values
    β_v = params[:beta_v].values
    lng_newpairs = [(1,6),(1,7),(1,8),(1,9)]
    _β_T = [0.9421, 0.9405, 0.9082, 0.886]
    _γ_T = [1.0307, 0.9917, 1.03884, 0.993]
    _β_v = [1.035, 1.0434, 1.02874, 1.023]
    _γ_v = [1.118, 1.143, 1.13209, 1.076]

    for (k,tup) in pairs(lng_newpairs)
        i,j = tup
        γ_T[i,j] =  _γ_T[k]
        γ_v[i,j] =  _γ_v[k]
        γ_T[j,i] =  _γ_T[k]
        γ_v[j,i] =  _γ_v[k]

        β_T[i,j] = _β_T[k]
        β_v[i,j] = _β_v[k]
        β_T[j,i] = 1/_β_T[k]
        β_v[j,i] = 1/_β_v[k]
    end

    nij = params[:nij].values
    tij = params[:tij].values
    dij = params[:dij].values
    eta_ij = params[:eta_ij].values
    epsilon_ij =params[:epsilon_ij].values
    beta_ij = params[:beta_ij].values
    gamma_ij = params[:gamma_ij].values
#C1,C4i
    nij[1,6] = [0.7588, -0.4386, -0.02273, 45.05, -2.291, -62.51, 33.32, -12.14]
    tij[1,6] = [1.02, 0.71, 1.57, 3.41, 2.12, 3.28, 3.37, 3.4]
    dij[1,6] = [1,2,3,1,1,1,2,1]
    eta_ij[1,6] = [1.34, 1.45, 0.96, 1.33, 1.9]
    epsilon_ij[1,6] = [0.59, 1.9, 0.87, 1.12, 1.43]
    beta_ij[1,6] = [1.07, 1.06, 1.11, 1.2, 1.23]
    gamma_ij[1,6] = [0.0, 0.0, 0.0, 0.0, 0.0]
#C1,C4i
    nij[1,7] = [0.9396, -0.1439, -0.1413, 35.32, -4.216, 59.17, -76.68, -41.39]
    tij[1,7] = [1.43, 0.3, 1.2, 3.1, 1.78, 3.36, 2.7, 3.7]
    dij[1,7] = [1,2,3,1,1,1,2,1]
    eta_ij[1,7] = [1.87, 1.05, 1.78, 1.19, 2.0]
    epsilon_ij[1,7] = [1.73, 0.78, 1.75, 1.84, 1.71]
    beta_ij[1,7] = [1.67, 1.76, 1.02, 1.76, 1.06]
    gamma_ij[1,7] = [0.0, 0.0, 0.0, 0.0, 0.0]

#C1C5
    nij[1,8] = [0.03711, -0.12154, 27.086, -13.614, -14.45, -0.46867]
    tij[1,8] = [1.54, 0.95, 0.47, 0.9, 0.22, 2.65]
    dij[1,8] = [2,3,1,1,1,2]
    eta_ij[1,8] = [0.6, 0.43, 0.7, 1.4]
    epsilon_ij[1,8] = [0.7, 0.61, 0.7, 0.46]
    beta_ij[1,8] = [0.916, 0.87, 0.86, 2.8]
    gamma_ij[1,8] = [0.5,0.5,0.5,0.5]
#C1C5i
    nij[1,9] =[0.051, -0.158, -67.49, -88.27, 154.9, 3.725]
    tij[1,9] = [0.2, 0.53, 1.79, 2.1, 2.0, 0.2]
    dij[1,9] = [2,3,1,1,1,2]
    eta_ij[1,9] = [0.64, 0.39, 0.48, 1.0]
    epsilon_ij[1,9] = [0.5, 0.5, 0.5, 0.5]
    beta_ij[1,9] = [1.56, 1.33, 1.46, 2.7]
    gamma_ij[1,9] = [0.0, 0.0, 0.0, 0.0]

    return params
end


text:
"""
Computes tendencies to qt and θ_liq_ice due to precipitation formation
"""
function compute_precipitation_formation_tendencies(
    grid,
    state,
    up::UpdraftVariables,
    precip::PrecipVariables,
    dt,
    param_set,
)
    p0_c = center_ref_state(state).p0
    ρ0_c = center_ref_state(state).ρ0
    aux_up = center_aux_updrafts(state)
    aux_bulk = center_aux_bulk(state)
    prog_pr = center_prog_precipitation(state)
    tendencies_pr = center_tendencies_precipitation(state)

    @inbounds for i in 1:(up.n_updrafts)
        @inbounds for k in real_center_indices(grid)
            T_up = aux_up[i].T[k]
            q_tot_up = aux_up[i].q_tot[k]
            ts_up = TD.PhaseEquil_pTq(param_set, p0_c[k], T_up, q_tot_up)

            # autoconversion and accretion
            mph = precipitation_formation(
                param_set,
                precip.precipitation_model,
                prog_pr.q_rai[k],
                prog_pr.q_sno[k],
                aux_up[i].area[k],
                ρ0_c[k],
                dt,
                ts_up,
            )
            aux_up[i].qt_tendency_precip_formation[k] = mph.qt_tendency * aux_up[i].area[k]
            aux_up[i].θ_liq_ice_tendency_precip_formation[k] = mph.θ_liq_ice_tendency * aux_up[i].area[k]

            tendencies_pr.q_rai[k] += mph.qr_tendency * aux_up[i].area[k]
            tendencies_pr.q_sno[k] += mph.qs_tendency * aux_up[i].area[k]
        end
    end
    # TODO - to be deleted once we sum all tendencies elsewhere
    @inbounds for k in real_center_indices(grid)
        aux_bulk.θ_liq_ice_tendency_precip_formation[k] = 0
        aux_bulk.qt_tendency_precip_formation[k] = 0
        @inbounds for i in 1:(up.n_updrafts)
            aux_bulk.θ_liq_ice_tendency_precip_formation[k] += aux_up[i].θ_liq_ice_tendency_precip_formation[k]
            aux_bulk.qt_tendency_precip_formation[k] += aux_up[i].qt_tendency_precip_formation[k]
        end
    end
    return
end

ext: jl


text:
# Exceptions
function exception1()
    try
        fin = open("this_file_does_not_exist.txt")
    catch exc
        println("Oh no! : $exc")
    end
end

function exception2()
    f = open("output.txt")
    try
        line = readline(f)
        println(line)
    finally
        close(f)
    end
end

# exception1()
# exception2()

# Serialization
using Serialization

io = IOBuffer()
a = [1, 2, 3, 4, 5]
serialize(io, a) # Translate object to byte array
s = take!(io)
println(s)
b = deserialize(IOBuffer(s));
println(b)
close(io)
# Note that they are no longer the same object:
println("a === b: $(a === b)")

# Command objects
cmd = `cmd /c echo hello` # In Windows
# cmd = `echo hello` in Unix
# Alternatively: cm = @cmd " echo 1 "
run(cmd)
a = read(cmd, String)
print(a)

filename = ".gitignore"
md5 = `powershell /c Get-FileHash $filename -algorithm md5`
# run(md5)

# Importing Julia files / modules
# include("wc.jl")

s = "1 2\t 3\n 4";
repr(s) # Produces: "\"1 2\\t 3\\n 4\""
dump(s) # Produces: String "1 2\t 3\n 4"


text:
"""
Overview: We here collect all functions and data structures required to manage
solid materials. The package has an internal materials database that is
updateed by the user to insert his own materials.
"""

"""
> struct material
- ρ : density of the material
- Dᵢⱼₖₗ : fourth order elasticity tensor in Voigt notation
"""
struct material
  #
  ρ::Float64
  Dᵢⱼₖₗ::Matrix{Float64}
  #
end

"""
> load_material( name::String )
It loads the data structure from the database
- name : name of the material to load
"""
function MORFE_load_material(name::String)
  pth = joinpath(get_materials_path(),name*".jld2")
  data = load(pth)
  return data["material"]
end

"""
> add_material( name::String, ρ::Float64, E::Float64, ν::Float64 )
It adds an isotropic material to the internal database of the package.
- name = name of the material
- ρ = density of the material
- E = Young's modulus
- ν = Poisson's ratio
**return**  *nothing*
"""
function MORFE_add_material(name::String, ρ::Float64, E::Float64, ν::Float64)
  #
  pth = get_materials_path()
  #
  Dᵢⱼₖₗ = zeros(Float64,(dim*(dim-1),dim*(dim-1)))
  #
  λ = E*ν/((1.0+ν)*(1.0-2.0*ν))
  μ = E/(1+ν)
  #
  Dᵢⱼₖₗ[1,1] = λ+μ
  Dᵢⱼₖₗ[1,2] = λ
  Dᵢⱼₖₗ[2,1] = λ
  Dᵢⱼₖₗ[2,2] = λ+μ
  Dᵢⱼₖₗ[1,3] = λ
  Dᵢⱼₖₗ[3,1] = λ
  Dᵢⱼₖₗ[2,3] = λ
  Dᵢⱼₖₗ[3,2] = λ
  Dᵢⱼₖₗ[3,3] = λ+μ
  Dᵢⱼₖₗ[4,4] = μ/2.0
  Dᵢⱼₖₗ[5,5] = μ/2.0
  Dᵢⱼₖₗ[6,6] = μ/2.0
  #
  mat = material(ρ, Dᵢⱼₖₗ)
  mname = joinpath(pth,name*".jld2")
  save(mname, Dict("material"=>mat))
  #
  return nothing
end

"""
> add_material( name::String, ρ::Float64, Dᵢⱼₖₗ::Matrix{Float64} )
It adds an anisotropic material to the internal database of the package.
- name = name of the material
- ρ = density of the material
- Dᵢⱼₖₗ = fourth order elasticity tensor
**return**  *nothing*
"""
function MORFE_add_material(name::String, ρ::Float64, Dᵢⱼₖₗ::Matrix{Float64})
  #
  pth = get_materials_path()
  #
  mat = material(ρ, Dᵢⱼₖₗ)
  mname = joinpath(pth,name*".jld2")
  save(mname, Dict("material"=>mat))
  #
  return nothing
end

"""
>  get_materials_path()
It returns the path associated to the material database
**return**  pth::String
"""
function get_materials_path()
  tail = length(joinpath("src","MORFEInvariantManifold.jl"))
  pth = chop(pathof(MORFEInvariantManifold), head = 0, tail = tail)*"materials"
  return pth
end

"""
> list_materials()
It lists all materials in the database.
**return**  *nothing*
"""
function MORFE_list_materials()
  pth = get_materials_path()
  tail = length(".jld2")
  files = readdir(pth)
  for f in files
    println(chop(f,head=0,tail=tail))
  end
  return nothing
end

"""
> delete_materials()
It delets a material in the database.
- name = name of the material to delete
**return**  *nothing*
"""
function MORFE_delete_material(name::String)
  pth = joinpath(get_materials_path(),name)*".jld2"
  try
    rm(pth)
  catch
    println("Material not present")
  end
  return nothing
end


text:
##
##!------------------------------------------------------------------------------------------------!
##  graph.jl
##
##  Tests for graph.jl.
##
##  Creation date: 18/08/15
##  Author: Arindam Biswas <arindam dot b at eml dot cc>
##!------------------------------------------------------------------------------------------------!
##

using SpecialGraphs
using Base.Test

for g in (GraphDynamic(UInt32(9)), BareGraphDynamic(UInt32(9)))
    # After initialization, order = 9, size = 0, degree(*) = 0
    # vertices = 1:9, edges = [], neighbours(*) = []
    # exists(*) = true, adjacent(*, *) = false
    @test order(g) == 9
    @test size(g) == 0
    for i in 1:9
        @test degree(g, UInt32(i)) == 0
        @test length(neighbours(g, UInt32(i))) == 0
        @test exists(g, UInt32(i))
    end

    @test length(vertices(g)) == 9
    @test length(edges(g)) == 0

    @test !adjacent(g, UInt32(1), UInt32(2))   # 1 -- 2
    @test !adjacent(g, UInt32(3), UInt32(2))   # 2 -- 3
    @test !adjacent(g, UInt32(5), UInt32(7))   # 5 -- 7
    @test !adjacent(g, UInt32(3), UInt32(6))   # 3 -- 6
    @test !adjacent(g, UInt32(4), UInt32(2))   # 2 -- 4
    @test !adjacent(g, UInt32(5), UInt32(8))   # 5 -- 8
    @test !adjacent(g, UInt32(9), UInt32(3))   # 3 -- 9

    # Add a vertex and edges corresponding to the tests above.
    add_vertex!(g)
    add_edge!(g, UInt32(1), UInt32(2))   # 1 -- 2
    add_edge!(g, UInt32(3), UInt32(2))   # 2 -- 3
    add_edge!(g, UInt32(5), UInt32(7))   # 5 -- 7
    add_edge!(g, UInt32(3), UInt32(6))   # 3 -- 6
    add_edge!(g, UInt32(4), UInt32(2))   # 2 -- 4
    add_edge!(g, UInt32(5), UInt32(8))   # 5 -- 8
    add_edge!(g, UInt32(9), UInt32(3))   # 3 -- 9

    # If g is not a bare graph, add vertex and edge information.
    if (typeof(g) == GraphDynamic)
        vertex_info!(g, UInt32(3), Int64(10))
        vertex_info!(g, UInt32(4), Int64(2))
        vertex_info!(g, UInt32(5), Int64(7))
        edge_info!(g, UInt32(2), UInt32(1), Int64(11))
        edge_info!(g, UInt32(3), UInt32(2), Int64(17))
        edge_info!(g, UInt32(2), UInt32(4), Int64(13))
        edge_info!(g, UInt32(9), UInt32(3), Int64(7))
    end

    # Here, order = 10, size = 7, degree(1) = 1, degree(2) = 3, degree(3) = 3, degree(5) = 2
    # length(vertices) = 10 , length(edges) = 7, sort(neighbours(3)) = [2, 6, 9], neighbours(6) = [3]
    # exists(*) = true, adjacent(1, 2) = true, adjacent(3, 2) = true, adjacent(8, 5) = true
    # adjacent(3, 6) = true, adjacent(3, 9) = true, adjacent(5, 6) = false, adjacent(9, 2) = false
    @test order(g) == 10
    @test size(g) == 7
    @test degree(g, UInt32(1)) == 1
    @test degree(g, UInt32(2)) == 3
    @test degree(g, UInt32(3)) == 3
    @test degree(g, UInt32(5)) == 2

    @test length(vertices(g)) == 10
    @test length(edges(g)) == 7

    @test sort(neighbours(g, UInt32(3))) == [2, 6, 9]
    @test neighbours(g, UInt32(6)) == [3]

    for i in 1:10
        @test exists(g, UInt32(i))
    end

    @test adjacent(g, UInt32(1), UInt32(2))    # 1 -- 2
    @test adjacent(g, UInt32(3), UInt32(2))    # 2 -- 3
    @test adjacent(g, UInt32(8), UInt32(5))    # 5 -- 8
    @test adjacent(g, UInt32(3), UInt32(6))    # 3 -- 6
    @test adjacent(g, UInt32(3), UInt32(9))    # 3 -- 9
    @test !adjacent(g, UInt32(5), UInt32(6))   # 5 -- 6
    @test !adjacent(g, UInt32(9), UInt32(2))   # 2 -- 9

    # If g is not a bare graph, check vertex and edge information.
    if (typeof(g) == GraphDynamic)
        @test vertex_info(g, UInt32(3)) == Int64(10)
        @test vertex_info(g, UInt32(4)) == Int64(2)
        @test vertex_info(g, UInt32(5)) == Int64(7)
        @test edge_info(g, UInt32(2), UInt32(1)) == Int64(11)
        @test edge_info(g, UInt32(3), UInt32(2)) == Int64(17)
        @test edge_info(g, UInt32(2), UInt32(4)) == Int64(13)
        @test edge_info(g, UInt32(9), UInt32(3)) == Int64(7)
    end

    # Remove 2 vertices and 3 edges.
    remove_vertex!(g, UInt32(4))
    remove_vertex!(g, UInt32(7))
    remove_edge!(g, UInt32(1), UInt32(2))      # 1 -- 2
    remove_edge!(g, UInt32(3), UInt32(2))      # 2 -- 3
    remove_edge!(g, UInt32(8), UInt32(5))      # 5 -- 8

    # Here, order = 8, size = 4, degree(1) = 0, degree(2) = 0, degree(3) = 2, degree(5) = 1
    # length(vertices) = 8, length(edges) = 4, sort(neighbours(3)) = [6, 9], neighbours(6) = [3]
    # exists(3) = true, exists(4) = false, exists(6) = true, exists(7) = false
    # adjacent(1, 2) = false, adjacent(3, 2) = false, adjacent(8, 5) = false, adjacent(3, 6) = true,
    # adjacent(3, 9) = true, adjacent(5, 6) = false, adjacent(9, 2) = false
    @test order(g) == 8
    @test size(g) == 2
    @test degree(g, UInt32(1)) == 0
    @test degree(g, UInt32(2)) == 0
    @test degree(g, UInt32(3)) == 2

    @test length(vertices(g)) == 8
    @test length(edges(g)) == 2

    @test sort(neighbours(g, UInt32(3))) == [6, 9]
    @test neighbours(g, UInt32(6)) == [3]

    @test exists(g, UInt32(3))
    @test !exists(g, UInt32(4))
    @test exists(g, UInt32(6))
    @test !exists(g, UInt32(7))

    @test !adjacent(g, UInt32(1), UInt32(2))   # 1 -- 2
    @test !adjacent(g, UInt32(3), UInt32(2))   # 2 -- 3
    @test !adjacent(g, UInt32(8), UInt32(5))   # 5 -- 8
    @test adjacent(g, UInt32(3), UInt32(6))    # 3 -- 6
    @test adjacent(g, UInt32(3), UInt32(9))    # 3 -- 9
    @test !adjacent(g, UInt32(5), UInt32(6))   # 5 -- 6
    @test !adjacent(g, UInt32(9), UInt32(2))   # 2 -- 9

    # If g is not a bare graph, check vertex and edge information.
    if (typeof(g) == GraphDynamic)
        @test vertex_info(g, UInt32(3)) == Int64(10)
        # @test vertex_info(g, UInt32(4)) == Int64(2) deleted
        @test vertex_info(g, UInt32(5)) == Int64(7)
        # @test edge_info(g, UInt32(2), UInt32(1)) == Int64(11) deleted
        # @test edge_info(g, UInt32(3), UInt32(2)) == Int64(17) deleted
        # @test edge_info(g, UInt32(2), UInt32(4)) == Int64(13) deleted
        @test edge_info(g, UInt32(9), UInt32(3)) == Int64(7)
    end

    # Add a deleted vertex, a deleted edge and a new edge.
    add_vertex!(g, UInt32(4))
    add_edge!(g, UInt32(3), UInt32(2))
    add_edge!(g, UInt32(9), UInt32(2))

    # If g is not a bare graph, add vertex and edge information.
    if (typeof(g) == GraphDynamic)
        vertex_info!(g, UInt32(4), Int64(9))
        edge_info!(g, UInt32(3), UInt32(2), Int64(3))
        edge_info!(g, UInt32(2), UInt32(9), Int64(13))
    end

    # Here, order = 9, size = 6, degree(1) = 0, degree(2) = 3, degree(3) = 3, degree(5) = 1
    # length(vertices) = 9, length(edges) = 6, sort(neighbours(3)) = [2, 6, 9], neighbours(6) = [3]
    # exists(3) = true, exists(4) = true, exists(6) = true, exists(7) = false
    # adjacent(1, 2) = false, adjacent(3, 2) = true, adjacent(8, 5) = false, adjacent(3, 6) = true,
    # adjacent(3, 9) = true, adjacent(5, 6) = false, adjacent(9, 2) = true
    @test order(g) == 9
    @test size(g) == 4
    @test degree(g, UInt32(1)) == 0
    @test degree(g, UInt32(2)) == 2
    @test degree(g, UInt32(3)) == 3

    @test length(vertices(g)) == 9
    @test length(edges(g)) == 4

    @test sort(neighbours(g, UInt32(3))) == [2, 6, 9]
    @test neighbours(g, UInt32(6)) == [3]

    @test exists(g, UInt32(3))
    @test exists(g, UInt32(4))
    @test exists(g, UInt32(6))
    @test !exists(g, UInt32(7))

    @test !adjacent(g, UInt32(1), UInt32(2))   # 1 -- 2
    @test adjacent(g, UInt32(3), UInt32(2))    # 2 -- 3
    @test !adjacent(g, UInt32(8), UInt32(5))   # 5 -- 8
    @test adjacent(g, UInt32(3), UInt32(6))    # 3 -- 6
    @test adjacent(g, UInt32(3), UInt32(9))    # 3 -- 9
    @test !adjacent(g, UInt32(5), UInt32(6))   # 5 -- 6
    @test adjacent(g, UInt32(9), UInt32(2))    # 2 -- 9

    # If g is not a bare graph, check vertex and edge information.
    if (typeof(g) == GraphDynamic)
        @test vertex_info(g, UInt32(3)) == Int64(10)
        @test vertex_info(g, UInt32(4)) == Int64(9) # added again
        @test vertex_info(g, UInt32(5)) == Int64(7)
        # @test edge_info(g, UInt32(2), UInt32(1)) == Int64(11) deleted
        @test edge_info(g, UInt32(3), UInt32(2)) == Int64(3)  # added again
        # @test edge_info(g, UInt32(2), UInt32(4)) == Int64(13) deleted
        @test edge_info(g, UInt32(9), UInt32(3)) == Int64(7)
        @test edge_info(g, UInt32(9), UInt32(2)) == Int64(13) # new
    end
end


text:
# Optimization interface

"""
    optimize(f[, F], indiv, algo[, opts])
    optimize(f[, F], constr, algo[, opts])
    optimize(f[, F], constr, indiv, algo[, opts])
    optimize(f[, F], constr, algo, poplt[, opts])

Perform optimization of the function `f` using the alorithm `algo` with the population, composed of
the initial population `poplt`, or individuals similar to the original individual `indiv`,
or generated from the constraints `constr`, with the options `opts`.
- For multi-objective optimization, the objective value `F` *must* be provided.
"""
optimize(f, individual, method::M,
         opts::Options = Options(;default_options(method)...)) where {M<:AbstractOptimizer} =
    optimize(f, NoConstraints(), individual, method, opts)
optimize(f, F::AbstractVector, individual, method::M,
         opts::Options = Options(;default_options(method)...)) where {M<:AbstractOptimizer} =
    optimize(f, F, NoConstraints(), individual, method, opts)
optimize(f, bounds::ConstraintBounds, method::M,
         opts::Options = Options(;default_options(method)...)) where {M<:AbstractOptimizer} =
    optimize(f, BoxConstraints(bounds), method, opts)
optimize(f, F, bounds::ConstraintBounds, method::M,
         opts::Options = Options(;default_options(method)...)) where {M<:AbstractOptimizer} =
    optimize(f, F, BoxConstraints(bounds), method, opts)
function optimize(f, constraints::C, method::M,
                  opts::Options = Options(;default_options(method)...)
                  ) where {M<:AbstractOptimizer, C<:AbstractConstraints}
    population = initial_population(method, bounds(constraints))
    optimize(f, constraints, method, population, opts)
end
function optimize(f, F, constraints::C, method::M,
                  opts::Options = Options(;default_options(method)...)
                  ) where {M<:AbstractOptimizer, C<:AbstractConstraints}
    population = initial_population(method, bounds(constraints))
    optimize(f, F, constraints, method, population, opts)
end
function optimize(f, constraints::C, individual, method::M,
                  opts::Options = Options(;default_options(method)...)
                  ) where {M<:AbstractOptimizer, C<:AbstractConstraints}
    population = initial_population(method, individual)
    optimize(f, constraints, method, population, opts)
end
function optimize(f, F, constraints::C, individual, method::M,
                  opts::Options = Options(;default_options(method)...)
                  ) where {M<:AbstractOptimizer, C<:AbstractConstraints}
    population = initial_population(method, individual)
    optimize(f, F, constraints, method, population, opts)
end
function optimize(f, constraints::C, method::M, population,
                  opts::Options = Options(;default_options(method)...)
                 ) where {M<:AbstractOptimizer, C<:AbstractConstraints}
    @assert length(population) > 0 "Population is empty"
    objfun = EvolutionaryObjective(f, first(population); eval=opts.parallelization)
    optimize(objfun, constraints, method, population, opts)
end
function optimize(f, F, constraints::C, method::M, population,
                  opts::Options = Options(;default_options(method)...)
                 ) where {M<:AbstractOptimizer, C<:AbstractConstraints}
    @assert length(population) > 0 "Population is empty"
    objfun = EvolutionaryObjective(f, first(population), F; eval=opts.parallelization)
    optimize(objfun, constraints, method, population, opts)
end

function optimize(objfun::D, constraints::C, method::M, population::AbstractArray,
                  options::Options = Options(;default_options(method)...),
                  state = initial_state(method, options, objfun, population)
                 ) where {D<:AbstractObjective, C<:AbstractConstraints, M<:AbstractOptimizer}
    # setup trace
    tr = OptimizationTrace{typeof(value(objfun)), typeof(method)}()
    tracing = options.store_trace || options.show_trace || options.callback !== nothing

    # prepare iteration counter (used to make "initial state" trace entry)
    iteration = 0
    t0 = time()
    stopped, stopped_by_callback = false, false
    converged, counter_tol = false, 0 # tolerance convergence
    is_moo = ismultiobjective(objfun)

    options.show_trace && print_header(method)
    trace!(tr, iteration, objfun, state, population, method, options, time()-t0)

    _time = time()
    while !converged && !stopped && iteration < options.iterations
        iteration += 1

        # perform state update
        update_state!(objfun, constraints, state, population, method, options, iteration) && break

        # evaluate convergence
        converged = assess_convergence(objfun, state, method, options)

        # update the function value
        !is_moo && value!(objfun, minimizer(state))

        # check convergence persistence
        counter_tol = converged ? counter_tol+1 : 0
        converged = converged && (counter_tol > options.successive_f_tol)

        # custom state-based termination
        converged = converged || terminate(state)

        if tracing
            # update trace; callbacks can stop routine early by returning true
            stopped_by_callback = trace!(tr, iteration, objfun, state, population, method, options, time()-t0)
        end

        _time = time()
        stopped_by_time_limit = _time-t0 > options.time_limit

        if stopped_by_callback || stopped_by_time_limit
            stopped = true
        end
    end

    after_while!(objfun, state, method, options)

    return EvolutionaryOptimizationResults(
        method,
        minimizer(state),
        is_moo ? NaN : value(objfun),
        iteration,
        iteration == options.iterations,
        converged,
        options.abstol,
        options.reltol,
        abschange(objfun, state),
        relchange(objfun, state),
        tr,
        f_calls(objfun),
        options.time_limit,
        _time-t0,
    )
end

after_while!(objfun, state, method, options) = nothing


text:
using MPI
using Test
using ClimateMachine.ODESolvers: AbstractODESolver
using ClimateMachine.GenericCallbacks

mutable struct PseudoSolver <: AbstractODESolver
    t::Float64
    steps::Int
    PseudoSolver() = new(42.0, 0)
end
gettime(ps::PseudoSolver) = ps.t

mutable struct MyCallback
    initialized::Bool
    calls::Int
    finished::Bool
end
MyCallback() = MyCallback(false, 0, false)

GenericCallbacks.init!(cb::MyCallback, _...) = cb.initialized = true
GenericCallbacks.call!(cb::MyCallback, _...) = (cb.calls += 1; nothing)
GenericCallbacks.fini!(cb::MyCallback, _...) = cb.finished = true

MPI.Init()
mpicomm = MPI.COMM_WORLD
ps = PseudoSolver()

wtcb = GenericCallbacks.EveryXWallTimeSeconds(MyCallback(), 2, mpicomm)
stcb = GenericCallbacks.EveryXSimulationTime(MyCallback(), 0.5)
sscb = GenericCallbacks.EveryXSimulationSteps(MyCallback(), 10)

fn_calls = 0
fncb = () -> (global fn_calls += 1)

wtfn_calls = 0
wtfncb = GenericCallbacks.EveryXWallTimeSeconds(
    AtInit(() -> global wtfn_calls += 1),
    2,
    mpicomm,
)

stfn_calls = 0
stfncb = GenericCallbacks.EveryXSimulationTime(
    AtInitAndFini(() -> global stfn_calls += 1),
    0.5,
)

ssfn_calls = 0
ssfncb = GenericCallbacks.EveryXSimulationSteps(
    AtInit(() -> global ssfn_calls += 1),
    5,
)

callbacks = ((wtcb, stcb, sscb), fncb, (wtfncb, stfncb, ssfncb))

@testset "GenericCallbacks" begin
    GenericCallbacks.init!(callbacks, ps, nothing, nothing, ps.t)

    @test wtcb.callback.initialized
    @test stcb.callback.initialized
    @test sscb.callback.initialized

    @test wtcb.callback.calls == 0
    @test stcb.callback.calls == 0
    @test sscb.callback.calls == 0
    @test fn_calls == 0
    @test wtfn_calls >= 1
    @test stfn_calls == 1
    @test ssfn_calls == 1

    @test GenericCallbacks.call!(callbacks, ps, nothing, nothing, ps.t) in
          (0, nothing)

    @test wtcb.callback.calls >= 0
    @test stcb.callback.calls == 0
    @test sscb.callback.calls == 0
    @test fn_calls == 1
    @test wtfn_calls >= 1
    @test stfn_calls == 1
    @test ssfn_calls == 1

    ps.t += 0.5
    @test GenericCallbacks.call!(callbacks, ps, nothing, nothing, ps.t) in
          (0, nothing)

    @test wtcb.callback.calls >= 0
    @test stcb.callback.calls == 1
    @test sscb.callback.calls == 0
    @test fn_calls == 2
    @test wtfn_calls >= 1
    @test stfn_calls == 2
    @test ssfn_calls == 1

    sleep(max((2.1 + wtcb.lastcbtime_ns / 1e9 - time_ns() / 1e9), 0.0))
    ps.t += 0.5
    @test GenericCallbacks.call!(callbacks, ps, nothing, nothing, ps.t) in
          (0, nothing)

    @test wtcb.callback.calls >= 1
    @test stcb.callback.calls == 2
    @test sscb.callback.calls == 0
    @test fn_calls == 3
    @test wtfn_calls >= 2
    @test stfn_calls == 3
    @test ssfn_calls == 1

    for i in 1:7
        @test GenericCallbacks.call!(callbacks, ps, nothing, nothing, ps.t) in
              (0, nothing)
    end

    @test wtcb.callback.calls >= 1
    @test stcb.callback.calls == 2
    @test sscb.callback.calls == 1
    @test fn_calls == 10
    @test wtfn_calls >= 2
    @test stfn_calls == 3
    @test ssfn_calls == 3

    GenericCallbacks.fini!(callbacks, ps, nothing, nothing, ps.t)
    @test wtcb.callback.finished
    @test stcb.callback.finished
    @test sscb.callback.finished
    @test wtfn_calls >= 2
    @test stfn_calls == 4
    @test ssfn_calls == 3
end


text:
module InteractVanilla

using WebIO: node, Node, instanceof, props, children, Scope, JSString, @js_str, onimport,
             setobservable!, onjs, onmount, WebIO
using Widgets: AbstractWidget, Widget, Widgets, scope
using Observables: on, Observable, AbstractObservable, ObservablePair, observe, Observables
using Dates
using Colors: Colorant, hex
using UUIDs: uuid4

import JSON

include("component.jl")
include("input.jl")

end # module


text:
export SegmentTree
export getindex, setindex!, top, range, equals

mutable struct SegmentTree{T}
    capacity::Int64
    values::Vector{T}
    select
    limit::T

    function SegmentTree{T}(values::Vector{T}, select, limit::T) where {T}
        local st = SegmentTree{T}(length(values), select, limit)

        for i = 1:length(values)
            st.values[st.capacity-1+i] = values[i]
        end

        for i = st.capacity-2:-1:0
            st.values[i+1] = select(st.values[2+(i<<1)], st.values[3+(i<<1)])
        end

        return st
    end

    function SegmentTree{T}(n::Int64, select, limit::T) where {T}
        local st = new()

        st.limit = limit
        st.select = select
        st.capacity = 1 << Int64(ceil(log2(n)))
        st.values = Vector{T}()
        resize!(st.values, st.capacity << 1)

        fill!(st.values, limit)

        return st
    end
end

function equals(st::SegmentTree{T}, tt::SegmentTree{T}) where {T}
    if st.capacity != tt.capacity
        return false
    end

    if st.values != tt.values
        return false
    end

    return true
end

function Base.:setindex!(st::SegmentTree{T}, value::T, index::Int64) where {T}
    @assert index <= st.capacity

    index += st.capacity - 1
    st.values[index] = value

    while index > 1
        local val = st.values[index]
        local neighbor

        if 1 != index % 2
            neighbor = index + 1
            index = index >> 1
        else
            neighbor = index - 1
            index = index >> 1
        end

        value = st.select(val, st.values[neighbor])

        if (st.values[index] != value)
            st.values[index] = value
        else
            break
        end
    end
end

function top(st::SegmentTree{T}) where {T}
    return st.values[1]
end

function Base.:getindex(st::SegmentTree{T}, index::Int64) where {T}
    return range(st, index, index)
end

function Base.:range(st::SegmentTree{T}, left::Int64, right::Int64) where {T}
    @assert left <= right
    @assert left <= st.capacity
    @assert right <= st.capacity

    left += (st.capacity - 2)
    right += (st.capacity - 2)
    local left_value = st.limit
    local right_value = st.limit

    while left < right
        if 0 == left % 2
            left_value = st.select(left_value, st.values[left+1])
        end

        left >>= 1

        if 1 == right % 2
            right_value = st.select(st.values[right+1], right_value)
        end
        right = (right >> 1) - 1
    end

    if left == right
        left_value = st.select(left_value, st.values[left+1])
    end

    return st.select(left_value, right_value)
end


text:
# Computes transit timing variations to linear order
# in eccentricity.  Please cite Agol & Deck (2015) if
# you make use of this in published research.

# module TTVFaster

# # VERSION < v"0.4-dev" && using Docile

# export Planet_plane_hk, compute_ttv!
# export Planet_plane                   # Deprecated, only exported to make the error message work

include("ttv_succinct.jl")

struct Planet_plane
  mass_ratio :: Float64
  period   :: Float64
  trans0   :: Float64
  ecc    :: Float64
  omega    :: Float64
end

struct Planet_plane_hk{T<:Real} # Parameters of a planet in a plane-parallel system
  # Mass ratio of the planet to the star:
  mass_ratio :: T
  # Initial time of transit:
  period   :: T
  trans0   :: T
  # e times cos or sin of longitude of periastron measured from line of sight, in radians:
  ecosw    :: T
  esinw    :: T
end

# """
# # Error message to explain to anyone who tries to use the old version
# """
function compute_ttv!(jmax::Integer,p1::Planet_plane,p2::Planet_plane,time1::Vector,time2::Vector,ttv1::Vector,ttv2::Vector)
  error("The Planet_plane data structure has been deprecated in favor of Planet_plane_hk")
end

# """
# # Computes transit-timing variations to linear order in
# # eccentricity for non-resonant, plane-parallel planets.
# # Input:
# #   jmax:  Maximum j over which to sum the TTV calculation for both planets
# #     p1:  Planet type for inner planet
# #     p2:  Planet type for outer planet
# #  time1:  Transit times for inner planet
# #  time2:  Transit times for outer planet
# #
# # Output:
# #   ttv1: TTVs of the inner planet
# #   ttv2: TTVs of the outer planet
# """

function compute_ttv!(jmax::Integer,p1::Planet_plane_hk{T},p2::Planet_plane_hk{T},time1::Vector{T},time2::Vector{T},ttv1::Vector{T},ttv2::Vector{T}) where T<:Real

  # Compute the semi-major axis ratio of the planets:
  # println(p1.period,p2.period)
  global alpha = (p1.period/p2.period)^(2//3)  #check memory allocation >>>>>>>>>>>>
  #println(alpha, p1.period, p2.period)
  @assert(alpha < 1)
  @assert(alpha > 0)
  # Number of times:
   ntime1 = length(time1)
   ntime2 = length(time2)
  f1=zeros(T,jmax+2,5)  #check memory allocation >>>>>>>>>>>>
  f2=zeros(T,jmax+2,5)  #check memory allocation >>>>>>>>>>>>
  # Compute the coefficients:
  ttv_succinct!(jmax+1,alpha,f1,f2)  # I need to compute coefficients one higher than jmax
  # Compute TTVs for inner planet (equation 33):
  # Compute since of \pomegas:
  e1 = sqrt(p1.esinw*p1.esinw + p1.ecosw*p1.ecosw)
  e2 = sqrt(p2.esinw*p2.esinw + p2.ecosw*p2.ecosw)
  sin1om=p1.esinw/e1
  sin2om=p2.esinw/e2
  cos1om=p1.ecosw/e1
  cos2om=p2.ecosw/e2

  # Compute mean motions:
  n1=2pi/p1.period
  n2=2pi/p2.period
  # Compute initial longitudes:
  lam10=-n1*p1.trans0 + 2*p1.esinw # 2*p1.eccen*sin1om
  lam20=-n2*p2.trans0 + 2*p2.esinw # 2*p2.eccen*sin2om
  @inbounds for i=1:ntime1
  # Compute the longitudes of the planets at times of transit of planet 1 (equation 49):
    lam11 = n1*time1[i]+lam10
    lam21 = n2*time1[i]+lam20
    psi1  = lam11-lam21 # Compute difference in longitudes at times of transit of planet 1
    sinpsi1=sin(psi1)
    cospsi1=cos(psi1)
    sinlam11 = sin(lam11)
    coslam11 = cos(lam11)
    sinlam1om1=sinlam11*cos1om - coslam11*sin1om
    coslam1om1=coslam11*cos1om + sinlam11*sin1om
    sinlam1om2=sinlam11*cos2om - coslam11*sin2om
    coslam1om2=coslam11*cos2om + sinlam11*sin2om
    ttv1[i]=zero(T) #0.0
    sinjm1psi1=zero(T) #0.0
    cosjm1psi1=one(T) #1.0
  # Sum over j:
    for j=1:jmax
      sinjpsi1=sinjm1psi1*cospsi1 + cosjm1psi1*sinpsi1
      cosjpsi1=cosjm1psi1*cospsi1 - sinjm1psi1*sinpsi1
      ttv1[i] += f1[j+1,1]*sinjpsi1
      ttv1[i] += f1[j+1,2]*e1*(sinjpsi1*coslam1om1 - cosjpsi1*sinlam1om1)
      ttv1[i] += f1[j+1,3]*e1*(sinjpsi1*coslam1om1 + cosjpsi1*sinlam1om1)
      ttv1[i] += f1[j  ,4]*e2*(sinjpsi1*coslam1om2 - cosjpsi1*sinlam1om2)
      ttv1[i] += f1[j+2,5]*e2*(sinjpsi1*coslam1om2 + cosjpsi1*sinlam1om2)
      sinjm1psi1=sinjpsi1
      cosjm1psi1=cosjpsi1
    end
  # Multiply by period and mass ratio, and divide by 2*Pi:
    ttv1[i] = ttv1[i]*p1.period*p2.mass_ratio/(2pi)
  end
  # Compute TTVs for outer planet (equation 33):
  @inbounds for i=1:ntime2
  # Compute the longitudes of the planets at times of transit of planet 2:
    lam12 = n1*time2[i]+lam10
    lam22 = n2*time2[i]+lam20
    sinlam22 = sin(lam22)
    coslam22 = cos(lam22)
    psi2  = lam12-lam22 # Compute difference in longitudes at times of transit of planet 2
    sinpsi2=sin(psi2)
    cospsi2=cos(psi2)
    sinlam2om1=sinlam22*cos1om - coslam22*sin1om
    coslam2om1=coslam22*cos1om + sinlam22*sin1om
    sinlam2om2=sinlam22*cos2om - coslam22*sin2om
    coslam2om2=coslam22*cos2om + sinlam22*sin2om
    ttv2[i]=zero(p2.period) #0.0
    sinjm1psi2=zero(p2.period) #0.0
    cosjm1psi2=one(p2.period) #1.0
  # Sum over j:
    for j=1:jmax
      sinjpsi2=sinjm1psi2*cospsi2 + cosjm1psi2*sinpsi2
      cosjpsi2=cosjm1psi2*cospsi2 - sinjm1psi2*sinpsi2
      ttv2[i] += f2[j+1,1]*sinjpsi2
      ttv2[i] += f2[j+1,2]*e2*(sinjpsi2*coslam2om2 - cosjpsi2*sinlam2om2)
      ttv2[i] += f2[j+1,3]*e2*(sinjpsi2*coslam2om2 + cosjpsi2*sinlam2om2)
      ttv2[i] += f2[j+2,4]*e1*(sinjpsi2*coslam2om1 - cosjpsi2*sinlam2om1)
      ttv2[i] += f2[j  ,5]*e1*(sinjpsi2*coslam2om1 + cosjpsi2*sinlam2om1)
      sinjm1psi2=sinjpsi2
      cosjm1psi2=cosjpsi2
    end
  # Multiply by period and mass ratio, and divide by 2*Pi:
    ttv2[i] = ttv2[i]*p2.period*p1.mass_ratio/(2pi)
  end
  # Finished!
  #ttv1,ttv2 are already allocated
  return
end  # compute_ttv!
# end # module

max_stars_repo_name: bmlindor/TTVFaster

max_issues_repo_name: bmlindor/TTVFaster

max_forks_repo_name: bmlindor/TTVFaster


text:
include("../src/MPSReader.jl")

fixed = true

files = Dict{String, Bool}([
    "lpsolve.mps" => fixed,
    "min4m3.mps" => fixed,
    "min4m3i.mps" => fixed,
    "min4m3b.mps" => fixed,
    "min4m3mip.mps" => fixed,
    "min4m3c0.mps" => fixed,
    "max4m3.mps" => !fixed
])
for file in keys(files)
    println("--- "*file*" ---")
    myModel, myVariables = mpstomodel("data/"*file, files[file])
    @info myModel
    println()
end

#= Expected output:
--- lpsolve.mps ---
┌ Info: Min x[1] + 4 x[2] + 9 x[3]
│ Subject to
│  x[1] ≥ 0.0
│  x[2] ≥ -1.0
│  x[3] ≥ 0.0
│  x[1] ≤ 4.0
│  x[2] ≤ 1.0
│  -x[2] + x[3] = 7.0
│  -x[1] - x[2] ≥ -5.0
└  x[1] + x[3] ≥ 10.0

--- min4m3mip.mps ---
┌ Info: Min -2 x[1] - 2 x[2] - x[3]
│ Subject to
│  x[1] binary
│  x[3] binary
│  x[4] integer
│  x[4] ≥ 0.0
│  x[2] ≥ 0.0
│  x[4] ≤ 10.0
│  x[2] ≤ 50.0
│  x[1] - x[2] + x[3] + 2 x[4] = -8.0
│  -x[1] - 2 x[2] - x[3] - x[4] ≥ 5.0
└  x[1] + x[2] + x[3] - x[4] ≥ 2.0

--- min4m3c0.mps ---
┌ Info: Min -2 x[1] - 2 x[2] - x[3] - 5
│ Subject to
│  x[1] ≥ 0.0
│  x[2] ≥ 0.0
│  x[3] ≥ 0.0
│  x[4] ≥ 0.0
│  x[1] - x[2] + x[3] + 2 x[4] = -8.0
│  -x[1] - 2 x[2] - x[3] - x[4] ≥ 5.0
└  x[1] + x[2] + x[3] - x[4] ≥ 2.0

--- min4m3i.mps ---
┌ Info: Min -2 x[1] - 2 x[2] - x[3]
│ Subject to
│  x[1] integer
│  x[2] integer
│  x[3] integer
│  x[4] integer
│  x[1] ≥ 0.0
│  x[2] ≥ 0.0
│  x[3] ≥ 0.0
│  x[4] ≥ 0.0
│  x[1] - x[2] + x[3] + 2 x[4] = -8.0
│  -x[1] - 2 x[2] - x[3] - x[4] ≥ 5.0
└  x[1] + x[2] + x[3] - x[4] ≥ 2.0

--- min4m3b.mps ---
┌ Info: Min -2 x[1] - 2 x[2] - x[3]
│ Subject to
│  x[1] binary
│  x[2] binary
│  x[3] binary
│  x[4] binary
│  x[1] - x[2] + x[3] + 2 x[4] = -8.0
│  -x[1] - 2 x[2] - x[3] - x[4] ≥ 5.0
└  x[1] + x[2] + x[3] - x[4] ≥ 2.0

--- max4m3.mps ---
┌ Info: Max 2 x[1] + 2 x[2] + x[3]
│ Subject to
│  x[1] ≥ 0.0
│  x[2] ≥ 0.0
│  x[3] ≥ 0.0
│  x[4] ≥ 0.0
│  x[1] - x[2] + x[3] + 2 x[4] = -8.0
│  x[1] + 2 x[2] + x[3] + x[4] ≤ -5.0
└  -x[1] - x[2] - x[3] + x[4] ≤ -2.0

--- min4m3.mps ---
┌ Info: Min -2 x[1] - 2 x[2] - x[3]
│ Subject to
│  x[1] ≥ 0.0
│  x[2] ≥ 0.0
│  x[3] ≥ 0.0
│  x[4] ≥ 0.0
│  x[1] - x[2] + x[3] + 2 x[4] = -8.0
│  -x[1] - 2 x[2] - x[3] - x[4] ≥ 5.0
└  x[1] + x[2] + x[3] - x[4] ≥ 2.0
=#

max_stars_repo_name: remi-garcia/MPSReader

max_issues_repo_name: remi-garcia/MPSReader

max_forks_repo_name: remi-garcia/MPSReader


text:
using DLSMethod
using FITSIO
using HDF5
using LinearAlgebra
using Random
using Statistics

include("BayesianSZNet.jl")
using .BayesianSZNet

ε = 0.00005
N_VAL, N_TEST = 50000, 50000

@time begin
    fid = h5open("data/dr12q_superset.hdf5", "r+")
    id = read(fid["id"])
    z_vi = read(fid["z_vi"])
    z_pipe = read(fid["z_pipe"])
    n = size(id, 2)

    fluxes = zeros(3752, n)

    Threads.@threads for i = 1:n
        loglam, flux = get_spectrum("Superset_DR12Q", id[:, i]...)

        # standardise flux
        flux_mean = mean(flux)
        standard_flux = Vector{Float64}((flux .- flux_mean) ./ std(flux, mean=flux_mean))

        X = polynomial_features(loglam, 3)
        n_best, a_best = dls_fit(X, standard_flux, 2, 0.9)

        # padding
        idx = (LOGLAMMIN - ε .< loglam) .& (loglam .< LOGLAMMAX + ε)
        pixels = sum(idx)    # count number of valid pixels
        Δ_loglam = minimum(loglam) - LOGLAMMIN
        offset = Δ_loglam <= 0 ? 1 : round(Int, Δ_loglam * 1e4) + 1
        fluxes[offset:offset + pixels - 1, i] = (standard_flux - X * a_best)[idx]
    end

    write_dataset(fid, "X", fluxes)

    # split into training, validation and test set (sizes almost according to ILSVRC)
    # seed from random.org
    Random.seed!(66)
    rnd_idx = randperm(n)
    n_tr = n - N_VAL - N_TEST
    idx_tr = rnd_idx[1:n_tr]
    idx_va = rnd_idx[n_tr + 1:n_tr + N_VAL]
    idx_te = rnd_idx[n_tr + N_VAL + 1:end]

    for (name, idx) in [("tr", idx_tr), ("va", idx_va), ("te", idx_te)]
        write_dataset(fid, "id_" * name, id[:, idx])
        write_dataset(fid, "X_" * name, fluxes[:, idx])
        write_dataset(fid, "z_vi_" * name, z_vi[idx])
        write_dataset(fid, "z_pipe_" * name, z_pipe[idx])
    end

    close(fid)
end

max_stars_repo_name: podondra/bdl-redshift

max_issues_repo_name: podondra/bdl-redshift

max_forks_repo_name: podondra/bdl-redshift


text:
# --------------------------------------------------------------------------
# ACE.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

A = ACE.Rotations3D.Rot3DCoeffs()
ACE.Rotations3D.basis(A, SVector(2, 2, 1, 1))

using ACE, LinearAlgebra

trans = PolyTransform(2, 1.0)
fcut = PolyCutoff1s(2, 0.1, 2.0)
basis = SHIPBasis(SparseSHIP(5, 8, 1.0), trans, fcut)

length(basis)

A = ACE.Rotations.Rot3DCoeffs(20)

for bo = 2:5
   ctr = 0
   println();
   @show bo
   for νz in basis.NuZ[bo]
      ν = νz.ν
      kk = getfield.(basis.KL[1][ν], :k)
      if norm(kk) != 0
         continue
      end
      ll = getfield.(basis.KL[1][ν], :l)
      U = ACE.Rotations.basis(A, ll)
      c = [ ACE._Bcoeff(ll, mm, A.cg) for mm in ACE._mrange(ll) ]
      def = norm(U * (U' * c) - c) > 1e-12
      if def > 1e-12
         @show ll, def
         ctr += 1
         if ctr > 10
            break
         end
      else
         print("✓")
      end
   end
end

using StaticArrays

for bo = 2:5
   println()
   @show bo
   ctr = 0
   cartrg = CartesianIndices( ntuple(_->0:3, bo) )
   for ill in cartrg
      ll = SVector(ill.I...)
      if issorted(ll) && sum(ll) <= 20
         U =  try
            ACE.Rotations.basis(A, ll)
         catch
            println()
            @show ll
            continue
         end
         c = [ ACE._Bcoeff(ll, mm, A.cg) for mm in ACE._mrange(ll) ]
         def = norm(U * (U' * c) - c)
         if def > 1e-12
            @show ll, def
            ctr += 1
            if ctr > 10
               break
            end
         else
            print("✓")
         end
      end
   end
end


text:
# default parameter values

𝒫ₕₛ = ComponentVector{Any}(
        # Steady state related
        rA = (0.5,),
        # Endogenous propagation
        τ⁻¹ = InverseGamma(8,8),
        κ   = Uniform(0.0,1.0),
        ψ_1 = Gamma(4,1/8),
        ψ_2 = Gamma(4,1/8),
        ρ_R = Uniform(0.5,0.9),
        # Exogenous shock parameters
        ρ_g = Uniform(0.90,0.99),
        ρ_z = Uniform(0.90,0.99),
        σ_R = InverseGamma(4,0.32),
        σ_g = InverseGamma(4,2.0),
        σ_z = InverseGamma(4,0.5),
)

"""
HerbstSchorfheide()

Small scale New Keynesian model from Herbst and Schorfheide,
"Bayesian Estimation of DSGE Models" (2016, Ch 2.1), in
Blanchard-Kahn canonical form.

"""

function HerbstSchorfheide(θ=Distributions.mode.(𝒫ₕₛ); output_indexes=false)

        eqn= (:mr,:exog_z,:exog_g,:y₋₁,:Δy,:π_obs,:euler,:phillips)  # equations (equations with expectations last)
        x  = (:R,:g,:z,:y₋₁,:Δy,:π_obs)                         # state vars
        y  = (:y,:π)                 # jump vars
        ε  = (:ε_z,:ε_g,:ε_R)        # exogenous shocks

        neqn,nx,ny,nε = length.((eqn,x,y,ε))
        ieqn = (; zip(eqn,(1:neqn...,))...)
        ix   = (; zip(x,(1:nx...,))...)
        iy   = (; zip(y,nx .+ (1:ny...,))...)
        iε   = (; zip(ε,(1:nε...,))...)

        β  = 1/(1+θ.rA/400)

        B = zeros(neqn,nx+ny)
        A = zeros(neqn,nx+ny)
        G = zeros(neqn,nε)

        B[ieqn[:euler],ix[:R]] =  θ.τ⁻¹
        B[ieqn[:euler],ix[:z]] = -θ.ρ_z*θ.τ⁻¹
        B[ieqn[:euler],ix[:g]] = -(1-θ.ρ_g)
        B[ieqn[:euler],iy[:y]] = -1
        B[ieqn[:euler],iy[:π]] = -θ.τ⁻¹
        A[ieqn[:euler],iy[:y]] = -1

        B[ieqn[:y₋₁],ix[:y₋₁]] = 1
        A[ieqn[:y₋₁],iy[:y]]   = 1

        B[ieqn[:Δy],ix[:Δy]]  =  1
        A[ieqn[:Δy],ix[:y₋₁]] = -1
        A[ieqn[:Δy],iy[:y]]   =  1

        B[ieqn[:π_obs],ix[:π_obs]] = 1
        A[ieqn[:π_obs],iy[:π]]     = 1

        B[ieqn[:phillips],iy[:π]] = -β
        B[ieqn[:phillips],ix[:g]] =  θ.κ
        A[ieqn[:phillips],iy[:y]] =  θ.κ
        A[ieqn[:phillips],iy[:π]] = -1

        B[ieqn[:mr],ix[:R]]   =  1
        B[ieqn[:mr],ix[:g]]   = (1-θ.ρ_R)*θ.ψ_2
        A[ieqn[:mr],ix[:R]]   =  θ.ρ_R
        A[ieqn[:mr],iy[:y]]   = (1-θ.ρ_R)*θ.ψ_2
        A[ieqn[:mr],iy[:π]]   = (1-θ.ρ_R)*(1+θ.ψ_1)
        G[ieqn[:mr],iε[:ε_R]] = -θ.σ_R

        B[ieqn[:exog_z],ix[:z]]   = 1
        A[ieqn[:exog_z],ix[:z]]   = θ.ρ_z
        G[ieqn[:exog_z],iε[:ε_z]] = θ.σ_z

        B[ieqn[:exog_g],ix[:g]]   = 1
        A[ieqn[:exog_g],ix[:g]]   = θ.ρ_g
        G[ieqn[:exog_g],iε[:ε_g]] = θ.σ_g

        B₁₁=vec(B[1:nx,1:nx])
        B₁₂=vec(B[1:nx,nx+1:neqn])
        B₂₁=vec(B[nx+1:neqn,1:nx])
        B₂₂=vec(B[nx+1:neqn,nx+1:neqn])
        A₁₁=vec(A[1:nx,1:nx])
        A₁₂=vec(A[1:nx,nx+1:neqn])
        A₂₁=vec(A[nx+1:neqn,1:nx])
        A₂₂=vec(A[nx+1:neqn,nx+1:neqn])
        G₁=vec(G[1:nx,:])
        G₂=vec(G[nx+1:neqn,:])

        Γ_vec = [B₁₁;B₁₂;B₂₁;B₂₂;A₁₁;A₁₂;A₂₁;A₂₂;G₁;G₂]

        axes = (Axis(
                B₁₁=ViewAxis(1:nx^2, ShapedAxis((nx, nx), NamedTuple())),
                B₁₂=ViewAxis(nx^2+1:nx^2+nx*ny, ShapedAxis((nx, ny), NamedTuple())),
                B₂₁=ViewAxis(nx^2+nx*ny+1:nx^2+2*nx*ny, ShapedAxis((ny, nx), NamedTuple())),
                B₂₂=ViewAxis(nx^2+2*nx*ny+1:neqn^2, ShapedAxis((ny, ny), NamedTuple())),
                A₁₁=ViewAxis(neqn^2+1:neqn^2+nx^2, ShapedAxis((nx, nx), NamedTuple())),
                A₁₂=ViewAxis(neqn^2+nx^2+1:neqn^2+nx^2+nx*ny, ShapedAxis((nx, ny), NamedTuple())),
                A₂₁=ViewAxis(neqn^2+nx^2+nx*ny+1:neqn^2+nx^2+2*nx*ny, ShapedAxis((ny, nx), NamedTuple())),
                A₂₂=ViewAxis(neqn^2+nx^2+2*nx*ny+1:2*neqn^2, ShapedAxis((ny, ny), NamedTuple())),
                G₁ =ViewAxis(2*neqn^2+1:2*neqn^2+nx*nε, ShapedAxis((nx, nε), NamedTuple())),
                G₂ =ViewAxis(2*neqn^2+nx*nε+1:2*neqn^2+neqn*nε, ShapedAxis((ny, nε), NamedTuple()))
                ),)

        Γ = ComponentArray(Γ_vec,axes)

        if output_indexes
                return Γ,neqn,nx,ny,nε,ieqn,ix,iy,iε
        else
                return Γ
        end
end

function ChainRulesCore.rrule(::typeof(HerbstSchorfheide), θ)

    Γ,neqn,nx,ny,nε,ieqn,ix,iy,iε = HerbstSchorfheide(θ; output_indexes=true)

    function HerbstSchorfheide_pullback(ΔΓ)

        CAΔΓ = ComponentArray(ΔΓ,getaxes(Γ))

        Ḡ = [reshape(CAΔΓ.G₁,nx,nε) ; reshape(CAΔΓ.G₂,ny,nε)]
        B̄ = [reshape(CAΔΓ.B₁₁,nx,nx) reshape(CAΔΓ.B₁₂,nx,ny) ;
             reshape(CAΔΓ.B₂₁,ny,nx) reshape(CAΔΓ.B₂₂,ny,ny)]
        Ā = [reshape(CAΔΓ.A₁₁,nx,nx) reshape(CAΔΓ.A₁₂,nx,ny) ;
             reshape(CAΔΓ.A₂₁,ny,nx) reshape(CAΔΓ.A₂₂,ny,ny)]

        β  = 1/(1+θ.rA/400)

        ∂σ_g =  Ḡ[ieqn[:exog_g],iε[:ε_g]]
        ∂ρ_g =  Ā[ieqn[:exog_g],ix[:g]] + B̄[ieqn[:euler],ix[:g]]
        ∂σ_z =  Ḡ[ieqn[:exog_z],iε[:ε_z]]
        ∂ρ_z =  Ā[ieqn[:exog_z],ix[:z]] - B̄[ieqn[:euler],ix[:z]]*θ.τ⁻¹
        ∂σ_R = -Ḡ[ieqn[:mr],iε[:ε_R]]
        ∂κ   =  B̄[ieqn[:phillips],ix[:g]] + Ā[ieqn[:phillips],iy[:y]]
        ∂β   = -B̄[ieqn[:phillips],iy[:π]]
        ∂rA  = -∂β*(1/400)*β^2
        ∂ρ_R = -(B̄[ieqn[:mr],ix[:g]]+Ā[ieqn[:mr],iy[:y]])*θ.ψ_2 +
                 Ā[ieqn[:mr],ix[:R]] - Ā[ieqn[:mr],iy[:π]]*(1+θ.ψ_1)
        ∂τ⁻¹ =  B̄[ieqn[:euler],ix[:R]] - θ.ρ_z*B̄[ieqn[:euler],ix[:z]] - B̄[ieqn[:euler],iy[:π]]

        ∂ψ_1 =  Ā[ieqn[:mr],iy[:π]]*(1-θ.ρ_R)
        ∂ψ_2 =  (B̄[ieqn[:mr],ix[:g]] + Ā[ieqn[:mr],iy[:y]])*(1-θ.ρ_R)

        ∂θ = ComponentArray(
                [∂rA,∂τ⁻¹,∂κ,∂ψ_1,∂ψ_2,∂ρ_R,∂ρ_g,∂ρ_z,∂σ_R,∂σ_g,∂σ_z],
                getaxes(θ)
             )

        return NO_FIELDS, ∂θ
    end
    return Γ, HerbstSchorfheide_pullback
end


text:
using StationXML, Test
using Dates: DateTime
import EzXML

@testset "I/O" begin
    @testset "Read(string)" begin
        let file = datapath("JSA.xml"), str = String(read(file))
            @test StationXML.read(file) == StationXML.readstring(str)
        end
    end

    @testset "Read" begin
        let badstring = """
            <?xml version="1.0" encoding="UTF8"?>
            <SomeWeirdThing xmlns="http://www.fdsn.org/xml/station/1">
            </SomeWeirdThing>
            """
            @test_throws ArgumentError StationXML.readstring(badstring)
        end
        # GB.JSA from IRIS
        let sxml = gzipped_read("JSA.xml.gz")
            @test typeof(sxml) == FDSNStationXML
            @test sxml.source == sxml.sender == "IRIS-DMC"
            @test sxml.module_name == "IRIS WEB SERVICE: fdsnws-station | version: 1.1.36"
            @test sxml.schema_version == "1.0"
            @test length(sxml.network) == 1
            @test length(sxml.network[1].station) == 1
            @test length(sxml.network[1].station[1].channel) == 6
            @test ismissing(sxml.network[1].historical_code)
            @test sxml.network[1].station[1].site.name == "ST AUBINS, JERSEY"
            @test sxml.network[1].station[1].start_date == DateTime(2007, 09, 06, 0, 0, 0)
            # Channel
            bhe = channels(sxml)[1]
            @test bhe.latitude.value ≈ 49.187801
            @test bhe.longitude.value ≈ -2.171698
            @test bhe.elevation.value ≈ 39.0
            @test bhe.depth.value ≈ 0.0 atol=1e-5
            @test bhe.azimuth.value ≈ 90.0
            @test bhe.dip.value ≈ 0.0 atol=1e-5
            @test bhe.clock_drift.value ≈ 0.0 atol=1e-5
            @test bhe.calibration_units.name == "A"
            @test bhe.calibration_units.description == "Amperes"
            @test bhe.restricted_status.value == "open"
            # Response
            @test bhe.response isa StationXML.Response
            @test !ismissing(bhe.response.instrument_sensitivity)
            @test ismissing(bhe.response.instrument_polynomial)
            @test length(bhe.response.stage) == 6
            @test count(!ismissing, getfield.((bhe.response.stage[1],),
                (:poles_zeros, :coefficients, :response_list, :fir, :polynomial))) == 1
            pz = bhe.response.stage[1].poles_zeros
            @test !ismissing(pz)
            @test pz.input_units.name == "M/S"
            @test pz.output_units.name == "V"
            @test pz.pz_transfer_function_type == StationXML.PZTransferFunction("LAPLACE (RADIANS/SECOND)")
            @test pz.normalization_factor ≈ 2.31323e9
            @test pz.normalization_frequency ≈ 1.0
            @test length(pz.zero) == 4
            @test length(pz.pole) == 7
            @test [p.real.value for p in pz.pole] ≈ [-0.01815, -0.01815, -196.0, -196.0, -732.0, -732.0, -173.0]
            @test [p.imaginary.value for p in pz.pole] ≈ [0.01799, -0.01799, 231.0, -231.0, 1415.0, -1415.0, 0.0]
            cf = bhe.response.stage[2].coefficients
            @test !ismissing(cf)
            @test cf.input_units.name == pz.output_units.name == "V"
            @test cf.output_units.name == "COUNTS"
            @test cf.cf_transfer_function_type == StationXML.CfTransferFunction("DIGITAL")
            @test length(cf.numerator) == length(cf.denominator) == 0
            dec = bhe.response.stage[2].decimation
            @test !ismissing(dec)
            @test dec.input_sample_rate.value ≈ 30_000.0
            @test dec.factor ≈ 1.0
            @test dec.offset ≈ 0.0 atol=1e-5
            @test dec.delay ≈ 0.0 atol=1e-5
            @test dec.correction ≈ 0.0 atol=1e-5
            @test bhe.response.stage[2].stage_gain.value ≈ 400_000.0
            @test bhe.response.stage[2].stage_gain.frequency ≈ 0.0 atol=1e-5
        end
        # NL.HGN from Orfeus
        let sxml = gzipped_read("orfeus_NL_HGN.xml.gz")
            @test sxml.schema_version == "1.0"
            @test sxml.source == "SeisComP3"
            @test sxml.sender == "ODC"
            @test sxml.created == DateTime(2020, 4, 28, 21, 17, 08, 868)
            @test length(sxml.network) == 1
            net = first(sxml.network)
            @test net.code == "NL"
            @test net.restricted_status == StationXML.RestrictedStatus("open")
            @test net.start_date == DateTime(1993)
            @test net.description == "Netherlands Seismic and Acoustic Network"
            @test length(net.station) == 1
            sta = first(net.station)
            @test sta.code == "HGN"
            @test sta.restricted_status == StationXML.RestrictedStatus("open")
            @test sta.start_date == DateTime(2001, 06, 06)
            @test sta.latitude.value == 50.764
            @test sta.latitude.plus_error === missing
            @test sta.latitude.minus_error === missing
            @test sta.latitude.unit === missing
            @test sta.latitude.datum == "WGS84"
            @test sta.longitude.value == 5.9317
            @test sta.longitude.plus_error === missing
            @test sta.longitude.minus_error === missing
            @test sta.longitude.unit === missing
            @test sta.longitude.datum == "WGS84"
            @test length(sta.channel) == 12
            cha = first(sta.channel)
            @test cha.end_date == DateTime(2003, 10, 24)
            @test cha.longitude == sta.longitude
            @test cha.latitude == sta.latitude
            @test cha.elevation == sta.elevation
            @test cha.dip.value == 0
            @test cha.dip.plus_error === missing
            @test cha.dip.minus_error === missing
            @test cha.dip.unit === missing
            @test cha.sample_rate.value == 40
            @test cha.sample_rate.unit === missing
            @test cha.sample_rate_ratio.number_samples == 40
            @test cha.sample_rate_ratio.number_seconds == 1
            @test cha.sensor.resource_id == "Sensor#20140909071419.913703.12"
            @test cha.sensor.type == "STS-1"
            @test cha.sensor.description == "STS-1"
            @test cha.sensor.model == "STS-1"
            @test cha.data_logger.resource_id == "Datalogger#20140909071419.913633.11"
            @test cha.data_logger.description == "HGN.1993.307.BHE"
            resp = cha.response
            @test resp.instrument_sensitivity.value == 801102000
            @test resp.instrument_sensitivity.frequency == 1
            @test resp.instrument_sensitivity.input_units == StationXML.Units("M/S")
            @test resp.instrument_sensitivity.output_units == StationXML.Units("COUNTS")
            @test length(resp.stage) == 2
            stage1 = resp.stage[1]
            @test stage1.number == 1
            @test stage1.poles_zeros !== missing
            @test stage1.stage_gain.value == 801102000
            @test stage1.stage_gain.frequency == 1
            pz = stage1.poles_zeros
            @test pz.name == "HGN.1993.307.HE"
            @test pz.resource_id == "ResponsePAZ#20140909071419.91376.13"
            @test pz.input_units == StationXML.Units("M/S")
            @test pz.output_units == StationXML.Units("V")
            @test pz.pz_transfer_function_type == StationXML.PZTransferFunction("LAPLACE (RADIANS/SECOND)")
            @test pz.normalization_factor == 3.86603e12
            @test pz.normalization_frequency == 1
            @test [p.real.value for p in pz.pole] == [-0.01234, -0.01234, -62.832,
                -39.144, -39.144, -14.012, -14.012, -56.612, -56.612]
            @test [p.imaginary.value for p in pz.pole] == [0.01234, -0.01234, 0,
                49.148, -49.148, 61.25, -61.25, 27.258, -27.258]
            @test [(z.real.value, z.imaginary.value) for z in pz.zero] ==
                [(0.0, 0.0), (0.0, 0.0)]
            @test [p.number for p in pz.pole] == 0:8
            @test [z.number for z in pz.zero] == 9:10
            @test all(x -> x===missing, p.real.plus_error for p in pz.pole)
            stage2 = resp.stage[2]
            @test stage2.number == 2
            @test stage2.coefficients !== missing
            @test stage2.decimation !== missing
            coefs = stage2.coefficients
            @test coefs.input_units == StationXML.Units("V")
            @test coefs.output_units == StationXML.Units("COUNTS")
            @test coefs.cf_transfer_function_type == StationXML.CfTransferFunction("DIGITAL")
            @test isempty(coefs.numerator)
            @test isempty(coefs.denominator)
            dec = stage2.decimation
            @test dec.input_sample_rate == StationXML.Frequency(40)
            @test dec.factor == 1
            @test dec.offset == 0
            @test dec.delay == 0
            @test dec.correction == 0

            # NL.HGN.02.BHZ
            cha = last(sta.channel)
            @test cha.code == "BHZ"
            @test cha.location_code == "02"
            @test cha.start_date == DateTime(2009, 4, 27)
            @test cha.dip == StationXML.Dip(-90)
            @test length(cha.response.stage) == 3
            stage3 = last(cha.response.stage)
            @test stage3.fir !== missing
            fir = stage3.fir
            @test fir.input_units == fir.output_units == StationXML.Units("COUNTS")
            @test fir.symmetry == StationXML.FIRSymmetry("NONE")
            @test length(fir.numerator_coefficient) == 39
            @test fir.numerator_coefficient == [StationXML.NumeratorCoefficient(c, missing)
                for c in [-5.09997e-08, -2.63324e-08, 3.21041e-07, -1.2548e-05,
                          -1.56874e-06, 0.000633072, 0.000467293, -0.00375769,
                          0.00242649, -0.00181103, -0.00255789, 0.00852735,
                          -0.0154895, 0.0199414, -0.0183877, 0.00706362, 0.0161458,
                          -0.0509218, 0.095284, -0.166693, 0.266723, 0.837233,
                          0.0500619, -0.0891913, 0.0853716, -0.0651036, 0.0395074,
                          -0.0158729, -0.00144455, 0.0107821, -0.0131202, 0.0104801,
                          -0.00623193, 0.00152521, 0.000205709, -0.00314123, 0.00102921,
                          0.000330318, 4.18952e-13]]
        end
    end

    @testset "Schema version" begin
        # Version 1
        for version in ("v1.0", "v1.1", "v1.1.2", "v1.2", "v1.99")
            str = """
                <FDSNStationXML xmlns="http://www.fdsn.org/xml/station/1" schemaVersion="$version">
                  <Source>StationXML.jl</Source>
                  <Sender>Test</Sender>
                  <Created>2019-04-01T09:58:21.123456789+01:00</Created>
                </FDSNStationXML>
                """
            xml = EzXML.parsexml(str)
            if version in ("v1.0", "v1.1", "v1.1.2")
                @test StationXML.schema_version_is_okay(xml)
            else
                @test (@test_logs (:warn,
                    "document is StationXML version $(VersionNumber(version)); " *
                    "only v1.1 data will be read") StationXML.schema_version_is_okay(xml))
            end
        end
        # Version 2+
        @test !StationXML.schema_version_is_okay(EzXML.parsexml("""
            <FDSNStationXML schemaVersion="2.0">
              <Source>StationXML.jl</Source>
              <Sender>AN</Sender>
              <Created>3000-01-01T00:00:00.00</Created>
            </FDSNStationXML>
            """))
    end

    @testset "Writing strings" begin
        let sxml = FDSNStationXML(source="StationXML.jl",
                created=DateTime(3000), schema_version=StationXML.DEFAULT_SCHEMA_VERSION), xml=xmldoc(sxml)
            push!(sxml.network, StationXML.Network(code="AN"))
            @test string(xml) == """
                <?xml version="1.0" encoding="UTF-8"?>
                <FDSNStationXML xmlns="http://www.fdsn.org/xml/station/1" schemaVersion="1.1"><Source>StationXML.jl</Source><Created>3000-01-01T00:00:00</Created></FDSNStationXML>
                """
        end
    end

    @testset "Writing" begin
        filenames = ("irisws_AK.xml.gz", "irisws_IU_BHx.xml.gz", "JSA.xml.gz",
            "orfeus_NL_HGN.xml.gz")

        @testset "To file" begin
            let sxml = sxml = FDSNStationXML(source="AN",
                    created=DateTime("2000-01-01"), schema_version=StationXML.DEFAULT_SCHEMA_VERSION)
                push!(sxml.network, StationXML.Network(code="AN"))
                push!(sxml.network[1].comment, StationXML.Comment("A comment"))
                io = IOBuffer()
                write(io, sxml)
                seekstart(io)
                sxml′ = StationXML.read(io)
                @test sxml == sxml′
                mktemp() do tempfile, f
                    write(f, sxml)
                    seekstart(f)
                    sxml″ = StationXML.read(f)
                    @test sxml == sxml″
                end
            end
        end

        @testset "Round trip file $file" for file in filenames
            sxml = gzipped_read(file)
            # Skip test files which are not expected to be written losslessly
            if VersionNumber(sxml.schema_version) < VersionNumber(StationXML.DEFAULT_SCHEMA_VERSION)
                continue
            end
            io = IOBuffer()
            write(io, sxml)
            seekstart(io)
            sxml′ = StationXML.read(io)
            @test sxml == sxml′
        end

        @testset "Warn removed fields" begin
            let sxml = gzipped_read("orfeus_NL_HGN.xml.gz")
                deleteat!(sxml.network[1].station[1].channel, 2:length(sxml.network[1].station[1].channel))
                @test_logs (:warn,
                    "Not writing field storage_format: removed in StationXML v$(StationXML.DEFAULT_SCHEMA_VERSION)"
                    ) write(devnull, sxml, warn=true)
                @test_logs write(devnull, sxml)
            end
        end
    end
end


text:
module CollectiveSpins

export System, Spin, SpinCollection, CavityMode, CavitySpinCollection,
        geometry,
        interaction, field, GreenTensor, OmegaMatrix, GammaMatrix,
        reducedspin, ReducedSpinBasis, reducedspintransition, reducedspinstate,
                reducedsigmap, reducedsigmam, reducedsigmax, reducedsigmay,
                reducedsigmaz, reducedsigmapsigmam,
        collective_modes, Omega_k_chain, Gamma_k_chain, Omega_k_2D, Gamma_k_2D

include("system.jl")
include("timeevolution_base.jl")
include("geometry.jl")
include("interaction.jl")
include("effective_interaction.jl")
include("effective_interaction_simple.jl")
include("effective_interaction_rotated.jl")
include("quantum.jl")
include("reducedspin.jl")
include("independent.jl")
include("meanfield.jl")
include("mpc.jl")
include("collective_modes.jl")
include("field.jl")

using .quantum
using .meanfield
using .mpc
using .reducedspin
using .interaction
using .collective_modes
using .field

end # module


text:
#=
    Variable bounds
        SingleVariable -in- LessThan
        SingleVariable -in- GreaterThan
        SingleVariable -in- EqualTo
        SingleVariable -in- Interval

        SingleVariable -in- ZeroOne
        SingleVariable -in- Integer
        SingleVariable -in- Semiinteger
        SingleVariable -in- Semicontinuous
=#
constrdict(model::LinQuadOptimizer, ::SVCI{LE}) = cmap(model).upper_bound
constrdict(model::LinQuadOptimizer, ::SVCI{GE}) = cmap(model).lower_bound
constrdict(model::LinQuadOptimizer, ::SVCI{EQ}) = cmap(model).fixed_bound
constrdict(model::LinQuadOptimizer, ::SVCI{IV}) = cmap(model).interval_bound

constrdict(model::LinQuadOptimizer, ::SVCI{MOI.ZeroOne}) = cmap(model).binary
constrdict(model::LinQuadOptimizer, ::SVCI{MOI.Integer}) = cmap(model).integer

constrdict(model::LinQuadOptimizer, ::SVCI{MOI.Semicontinuous{Float64}}) = cmap(model).semicontinuous
constrdict(model::LinQuadOptimizer, ::SVCI{MOI.Semiinteger{Float64}}) = cmap(model).semiinteger

function set_variable_bound(model::LinQuadOptimizer, column::Int,
                            bound::Float64, sense)
    change_variable_bounds!(model, [column], [bound], [sense])
end

function set_variable_bound(model::LinQuadOptimizer, v::SinVar, set::LE)
    set_variable_bound(model, get_column(model, v), set.upper,
                      backend_type(model, Val{:Upperbound}()))
end
function set_variable_bound(model::LinQuadOptimizer, v::SinVar, set::GE)
    set_variable_bound(model, get_column(model, v), set.lower,
                      backend_type(model, Val{:Lowerbound}()))
end
function set_variable_bound(model::LinQuadOptimizer, v::SinVar, set::EQ)
    set_variable_bound(model, get_column(model, v), set.value,
                      backend_type(model, Val{:Upperbound}()))
    set_variable_bound(model, get_column(model, v), set.value,
                      backend_type(model, Val{:Lowerbound}()))
end
function set_variable_bound(model::LinQuadOptimizer, v::SinVar, set::IV)
    set_variable_bound(model, get_column(model, v), set.upper,
                      backend_type(model, Val{:Upperbound}()))
    set_variable_bound(model, get_column(model, v), set.lower,
                      backend_type(model, Val{:Lowerbound}()))
end

"""
    has_value(dict::Dict{K, V}, value::V) where {K, V}

Return true if `dict` has a `key=>value` pair with `value`.
"""
function has_value(dict::Dict{K, V}, value::V) where {K, V}
    return value in values(dict)
end

"""
    __check_for_conflicting__(model::LinQuadOptimizer, variable::SinVar, set,
                              conflicting_type::Type{<:AbstractSet})

Throw an error if `variable` is already constrained to be in a set of type
`conflicting_type`.
"""
function __check_for_conflicting__(model::LinQuadOptimizer, variable::SinVar, set,
                                   conflict_type::Type{<:MOI.AbstractSet})
    if has_value(constrdict(model, SVCI{conflict_type}(0)), variable.variable)
        error("Cannot add constraint $(variable)-in-$(set) as it is already " *
              "constrained by a set of type $(conflict_type).")
    end
end

function __check_for_conflicting__(model::LinQuadOptimizer, variable::SinVar, set,
                                   conflict_type::Type{MOI.ZeroOne})
    for (index, lower, upper) in values(constrdict(model, SVCI{conflict_type}(0)))
        if index == variable.variable
            error("Cannot add constraint $(variable)-in-$(set) as it is already " *
                  "constrained by a set of type $(conflict_type).")
        end
    end
end

"""
    __check_for_conflicting__(model::LinQuadOptimizer, variable::SinVar, set,
                              conflicting_types::AbstractSet...)

Throw an error if `variable` is constrained to be in a set whose type is one of
`conflicting_types`.
"""
function __check_for_conflicting__(model::LinQuadOptimizer, variable::SinVar, set,
                                   conflicting_types...)
    for conflict_type in conflicting_types
        __check_for_conflicting__(model, variable, set, conflict_type)
    end
end

function MOI.addconstraint!(model::LinQuadOptimizer, variable::SinVar, set::S) where S <: LinSets
    __assert_supported_constraint__(model, SinVar, S)
    __check_for_conflicting__(model, variable, set,
        S, MOI.Semicontinuous{Float64}, MOI.Semiinteger{Float64}, MOI.ZeroOne)
    set_variable_bound(model, variable, set)
    model.last_constraint_reference += 1
    index = SVCI{S}(model.last_constraint_reference)
    dict = constrdict(model, index)
    dict[index] = variable.variable
    return index
end

function MOI.delete!(model::LinQuadOptimizer, index::SVCI{S}) where S <: LinSets
    __assert_valid__(model, index)
    delete_constraint_name(model, index)
    dict = constrdict(model, index)
    variable_index = dict[index]
    set_variable_bound(model, SinVar(variable_index), IV(-Inf, Inf))
    delete!(dict, index)
end

# constraint set
MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{SVCI{S}}) where S <: LinSets = true
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintSet, index::SVCI{LE})
    MOI.LessThan{Float64}(
        get_variable_upperbound(model, get_column(model, model[index]))
    )
end
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintSet, index::SVCI{GE})
    MOI.GreaterThan{Float64}(
        get_variable_lowerbound(model, get_column(model, model[index]))
    )
end
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintSet, index::SVCI{EQ})
    MOI.EqualTo{Float64}(
        get_variable_lowerbound(model, get_column(model, model[index]))
    )
end
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintSet, index::SVCI{IV})
    column = get_column(model, model[index])
    return MOI.Interval{Float64}(get_variable_lowerbound(model, column),
                                 get_variable_upperbound(model, column))
end

# constraint function
MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintFunction, ::Type{SVCI{S}}) where S <: LinSets = true
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintFunction, index::SVCI{<: LinSets})
    return SinVar(model[index])
end

# modify
MOI.supports(::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{SVCI{S}}) where S<:LinSets = true
function MOI.set!(model::LinQuadOptimizer, ::MOI.ConstraintSet, index::SVCI{S}, newset::S) where S <: LinSets
    set_variable_bound(model, SinVar(model[index]), newset)
end

#=
    Binary constraints

For some reason CPLEX doesn't respect bounds on a binary variable, so we
should store the previous bounds so that if we delete the binary constraint
we can revert to the old bounds

Xpress is worse, once binary, the bounds are changed independently of what the user does
=#
function MOI.addconstraint!(model::LinQuadOptimizer, variable::SinVar, set::MOI.ZeroOne)
    __assert_supported_constraint__(model, SinVar, MOI.ZeroOne)
    __check_for_conflicting__(model, variable, set, MOI.ZeroOne, MOI.Integer,
        MOI.Semicontinuous{Float64}, MOI.Semiinteger{Float64})
    model.last_constraint_reference += 1
    index = SVCI{MOI.ZeroOne}(model.last_constraint_reference)
    column = get_column(model, variable)
    dict = constrdict(model, index)
    dict[index] = (variable.variable, get_variable_lowerbound(model, column),
                    get_variable_upperbound(model, column))
    change_variable_types!(model, [column], [backend_type(model, set)])
    set_variable_bound(model, column, 1.0,
                       backend_type(model, Val{:Upperbound}()))
    set_variable_bound(model, column, 0.0,
                       backend_type(model, Val{:Lowerbound}()))
    make_problem_type_integer(model)
    return index
end

function MOI.delete!(model::LinQuadOptimizer, index::SVCI{MOI.ZeroOne})
    __assert_valid__(model, index)
    delete_constraint_name(model, index)
    dict = constrdict(model, index)
    (variable, lower, upper) = dict[index]
    column = get_column(model, variable)
    change_variable_types!(model, [column],
                           [backend_type(model, Val{:Continuous}())])
    set_variable_bound(model, column, upper,
                       backend_type(model, Val{:Upperbound}()))
    set_variable_bound(model, column, lower,
                       backend_type(model, Val{:Lowerbound}()))
    delete!(dict, index)
    if !has_integer(model)
        make_problem_type_continuous(model)
    end
end

MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{SVCI{MOI.ZeroOne}}) = true
function MOI.get(::LinQuadOptimizer, ::MOI.ConstraintSet, ::SVCI{MOI.ZeroOne})
    return MOI.ZeroOne()
end

MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintFunction, ::Type{SVCI{MOI.ZeroOne}}) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintFunction, index::SVCI{MOI.ZeroOne})
    return SinVar(model[index][1])
end

#=
    Integer constraints
=#

function MOI.addconstraint!(model::LinQuadOptimizer, variable::SinVar, set::MOI.Integer)
    __assert_supported_constraint__(model, SinVar, MOI.Integer)
    __check_for_conflicting__(model, variable, set, MOI.ZeroOne,
        MOI.Semicontinuous{Float64}, MOI.Semiinteger{Float64})
    change_variable_types!(model, [get_column(model, variable)],
                           [backend_type(model, set)])
    model.last_constraint_reference += 1
    index = SVCI{MOI.Integer}(model.last_constraint_reference)
    dict = constrdict(model, index)
    dict[index] = variable.variable
    make_problem_type_integer(model)
    return index
end

function MOI.delete!(model::LinQuadOptimizer, index::SVCI{MOI.Integer})
    __assert_valid__(model, index)
    delete_constraint_name(model, index)
    dict = constrdict(model, index)
    variable = dict[index]
    change_variable_types!(model, [get_column(model, variable)],
                           [backend_type(model, Val{:Continuous}())])
    delete!(dict, index)
    if !has_integer(model)
        make_problem_type_continuous(model)
    end
end

MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{SVCI{MOI.Integer}}) = true
function MOI.get(::LinQuadOptimizer, ::MOI.ConstraintSet, ::SVCI{MOI.Integer})
    return MOI.Integer()
end

MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintFunction, ::Type{SVCI{MOI.Integer}}) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintFunction,
                 index::SVCI{MOI.Integer})
    return SinVar(model[index])
end

#=
    Semicontinuous / Semiinteger constraints
=#
const SEMI_TYPES = Union{MOI.Semicontinuous{Float64}, MOI.Semiinteger{Float64}}
function MOI.addconstraint!(model::LinQuadOptimizer, variable::SinVar, set::S) where S <: SEMI_TYPES
    __assert_supported_constraint__(model, SinVar, S)
    __check_for_conflicting__(model, variable, set, S, MOI.ZeroOne, MOI.Integer)
    if S == MOI.Semicontinuous{Float64}
        __check_for_conflicting__(model, variable, set, MOI.Semiinteger{Float64})
    else
        __check_for_conflicting__(model, variable, set, MOI.Semicontinuous{Float64})
    end
    column = get_column(model, variable)
    change_variable_types!(model, [column], [backend_type(model, set)])
    set_variable_bound(model, column, set.upper,
                       backend_type(model, Val{:Upperbound}()))
    set_variable_bound(model, column, set.lower,
                       backend_type(model, Val{:Lowerbound}()))
    model.last_constraint_reference += 1
    index = SVCI{S}(model.last_constraint_reference)
    dict = constrdict(model, index)
    dict[index] = variable.variable
    make_problem_type_integer(model)
    return index
end

function MOI.delete!(model::LinQuadOptimizer, index::SVCI{<:SEMI_TYPES})
    __assert_valid__(model, index)
    delete_constraint_name(model, index)
    dict = constrdict(model, index)
    column = get_column(model, dict[index])
    change_variable_types!(model, [column], [backend_type(model, Val{:Continuous}())])
    set_variable_bound(model, column, Inf,
                       backend_type(model, Val{:Upperbound}()))
    set_variable_bound(model, column, -Inf,
                       backend_type(model, Val{:Lowerbound}()))
    delete!(dict, index)
    if !has_integer(model)
        make_problem_type_continuous(model)
    end
end

MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{SVCI{S}}) where S <:SEMI_TYPES = true
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintSet, index::SVCI{S}) where S <: SEMI_TYPES
    dict = constrdict(model, index)
    column = get_column(model, dict[index])
    return S(get_variable_lowerbound(model, column),
             get_variable_upperbound(model, column))
end

MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintFunction, ::Type{SVCI{S}}) where S <:SEMI_TYPES = true
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintFunction, index::SVCI{<:SEMI_TYPES})
    return SinVar(model[index])
end


text:
import MPI

# ------------------------------------
# mpi
# ------------------------------------
function set_mpi(cellxmax, cellymax)
    comm = MPI.COMM_WORLD        # MPIプロセスの呼び出し
    rank = MPI.Comm_rank(comm)   # MPIの分割, 各プロセス番号(0~n)
    size = MPI.Comm_size(comm)   # MPIの分割, 各プロセス番号(0~n)

    num_div = 0.0      # 各方向の分割数
    if floor(size^0.5 * 10000) % 10000 != 0 || size <= 0.0
        println("------------------------\n")
        println(" mpi error \n")
        println(" please change parallel number to be an number of squares \n")
        println("------------------------")
        throw(UndefVarError(:x))
    else
        num_div = Int(size^0.5)
    end

    Nx = cellxmax-4 # 仮想セル分を引く
    Ny = cellymax-4 # 仮想セル分を引く

    #=
    並列時のstart i, end i, 分割数cellNxを定義
    =#
    cellNx, ista, iend, cellNy, jsta, jend = start_and_end0(Nx, Ny, num_div, rank)

    return cellNx, ista, iend, cellNy, jsta, jend, num_div, rank, size, comm
end

# ---------------------------------------------
# mpi division
# ---------------------------------------------
function start_and_end0(Nx, Ny, num_div, rank)

    if Nx % num_div != 0
        println("------------------------\n")
        println(" mpi error \n ")
        println(" please change Nx \n ")
        println(" Nx = ")
        println(Nx)
        println(" \n ")
        println(" num_div = ")
        println(num_div)
        println("------------------------")
    elseif Ny % num_div != 0
        println("------------------------\n")
        println(" mpi error \n ")
        println(" please change Ny \n ")
        println(" Ny = ")
        println(Ny)
        println(" \n ")
        println(" num_div = ")
        println(num_div)
        println("------------------------")
    end

    cellNx = div(Nx, num_div)   # プロセスの分割数
    cellNy = div(Ny, num_div)   # プロセスの分割数
    ista = rank * cellNx + 1  # プロセス毎のstart i (0-> 0*n+1 = 1,        1-> 1*n+1 = n+1)
    iend = ista + cellNx - 1  # プロセス毎のend i   (0-> 0*n+1 + n-1 = n , 1-> 1*n+1 + n-1 = 2n)
    jsta = rank * cellNy + 1  # プロセス毎のstart i (0-> 0*n+1 = 1,        1-> 1*n+1 = n+1)
    jend = ista + cellNy - 1  # プロセス毎のend i   (0-> 0*n+1 + n-1 = n , 1-> 1*n+1 + n-1 = 2n)

    return cellNx, ista, iend, cellNy, jsta, jend
end

function set_value_for_mpi(Qbase, cellx_mpi, celly_mpi, volume, cellcenter, wally, dx, dy,
                        Qbase_all, volume_all, cellcenter_all, wally_all, dx_all, dy_all, nval, num_div, rank)

    is = 3
    ie = cellx_mpi-2
    js = 3
    je = celly_mpi-2
    for l in 1:nval
        for i in is:ie
            for j in js:je
                ii, jj = Conversion_to_coordinate_numbers(i, j, cellx_mpi, celly_mpi, num_div, rank)
                Qbase[i,j,l]  = Qbase_all[ii, jj, l]
            end
        end
    end
    for i in is:ie
        for j in js:je
            ii, jj = Conversion_to_coordinate_numbers(i, j, cellx_mpi, celly_mpi, num_div, rank)
            volume[i,j] = volume_all[ii, jj]
            wally[i,j]  = wally_all[ii, jj]
            dx[i,j]     = dx_all[ii, jj]
            dy[i,j]     = dy_all[ii, jj]
        end
    end

    return Qbase, volume, cellcenter, wally, dx, dy
end

function Conversion_to_coordinate_numbers(i, j, cellx_mpi, celly_mpi, num_div, rank)
    #=
    div system

      rank      div(rank, size^0.5)    rem(rank, size^0.5)

      2 | 3            1 | 1                  0 | 1
     -------          -------                -------
      0 | 1            0 | 0                  0 | 1

      y
      ↑
       →x

    =#
    rankx = rem(rank, num_div)
    ranky = div(rank, num_div)

    ii = i + rankx*(cellx_mpi-2)
    jj = j + ranky*(celly_mpi-2)

    return ii, jj
end

function set_vecA_mpi(vecAx_all, vecAy_all, vecAx, vecAy, cellx_mpi, celly_mpi, num_div, rank)
    #=
    vecA system

                    vecAy[2,3]
                |---------------|
                |               |
                |               |
    vecAx[2,2]  |   cell(2,2)   | vecAx[3,2]
                |               |
                |               |
                |---------------|
                    vecAy[2,2]

    =#
    is = 1
    ie = cellx_mpi+1
    js = 1
    je = celly_mpi
    for l in 1:2
        for i in is:ie
            for j in js:je
                ii, jj = Conversion_to_coordinate_numbers_atbound(i, j, cellx_mpi, celly_mpi, num_div, rank)
                vecAx[i,j,l]  = vecAx_all[ii, jj, l]
            end
        end
    end

    is = 1
    ie = cellx_mpi
    js = 1
    je = celly_mpi+1
    for l in 1:2
        for i in is:ie
            for j in js:je
                ii, jj = Conversion_to_coordinate_numbers_atbound(i, j, cellx_mpi, celly_mpi, num_div, rank)
                vecAy[i,j,l]  = vecAy_all[ii, jj, l]
            end
        end
    end

    return vecAx, vecAy
end

function Conversion_to_coordinate_numbers_atbound(i, j, cellx_mpi, celly_mpi, num_div, rank)
    #=
    div system

      rank      div(rank, size^0.5)    rem(rank, size^0.5)

      2 | 3            1 | 1                  0 | 1
     -------          -------                -------
      0 | 1            0 | 0                  0 | 1

    =#
    rankx = rem(rank, num_div)
    ranky = div(rank, num_div)

    ii = i + rankx*(cellx_mpi)
    jj = j + ranky*(celly_mpi)

    return ii, jj
end

function mpi_div_3d(ite, ls, le, js, je, is, ie)
    #= ite: iteration number (ista, iend generated by start_and_end(N,size))
        ls: start number at l
        le: end number at l
        js:
        je:
        is:
        ie:
        div():商
        rem():余り
    =#
    Ni = ie - is + 1 # number of i ite
    Nj = je - js + 1 # number of j ite
    #Nl = le - ls + 1 # number of l ite

    l = ls + div(ite-1, Ni*Nj)

    ite = ite - div(ite-1, Ni*Nj) * Ni*Nj

    j = js + div(ite-1, Ni)

    ite = ite - div(ite-1, Ni) * Ni

    i = is + div(ite-1, 1)

    return i, j, l
end

function mpi_div_2d(ite, js, je, is, ie)
    #= ite: iteration number (ista, iend generated by start_and_end(N,size))
        js:
        je:
        is:
        ie:
        div():商
        rem():余り
    =#
    Ni = ie - is + 1 # number of i ite

    j = js + div(ite-1, Ni)

    ite = ite - div(ite-1, Ni) * Ni

    i = is + div(ite-1, 1)

    return i, j
end


text:
using Test
using TreeTools

println("##### iterators #####")

begin
	nwk = "(((A1:1,A2:1,B:1,C:1)i_2:0.5,(D1:0,D2:0)i_3:1.5)i_1:1,(E:1.5,(F:1.0,(G:0.,H:0.)i_6:0.5)i_5:1)i_4:1)i_r"
	tree = parse_newick_string(nwk)
	tnodes = sort(["A1", "A2", "B", "C", "i_2", "D1", "D2", "i_3", "i_1", "E", "F", "G", "H", "i_6", "i_5", "i_4", "i_r"])
	tleaves = sort(["A1", "A2", "B", "C", "D1", "D2", "E", "F", "G", "H"])

	@testset "1" begin
		@test sort([x.label for x in POT(tree)]) == tnodes
		@test sort([x.label for x in POTleaves(tree)]) == tleaves
	end

	@testset "in" begin
		for n in nodes(tree)
			@test in(n, tree)
			@test in(n.label, tree)
			if !n.isleaf
				@test !in(n, tree; exclude_internals=true)
			end
		end
	end

	@testset "indexing in tree" begin
		for n in nodes(tree)
			@test tree[n.label] == tree.lnodes[n.label]
		end
		@test_throws KeyError tree["some_random_string"]
	end
end


text:
using PlotlyJS
include(Pkg.dir("CompScienceMeshes","examples","plotlyjs_patches.jl"))

using CompScienceMeshes, BEAST
o, x, y, z = euclidianbasis(3)

#D, Δx = 1.0, 0.25
D, Δx = 1.0, 0.35
Γ = meshsphere(D, Δx)
X, Y = raviartthomas(Γ), buffachristiansen(Γ)

#Δt, Nt = 0.1, 600
Δt, Nt = 0.11, 200
T, δ = timebasisshiftedlagrange(Δt, Nt, 2), timebasisdelta(Δt, Nt)

V = X ⊗ T; W = Y ⊗ δ

width, delay, scaling = 8.0, 12.0, 1.0
gaussian = creategaussian(width, delay, scaling)

direction, polarisation = z, x
E = planewave(polarisation, direction, gaussian, 1.0); H = direction × E

@hilbertspace j; @hilbertspace m′
K, I, N = MWDoubleLayerTDIO(1.0, 1.0, 0), Identity(), NCross()
tdmfie = @discretise (0.5*(N⊗I) + 1.0*K)[m′,j] == -1E[m′]   j∈V  m′∈W
xmfie = solve(tdmfie)

t3 = scatter(x=0:Δt:(Nt-1)Δt, y=xmfie[1,:])

Xmfie, Δω, ω0 = fouriertransform(xmfie, Δt, 0.0, 2)
ω = collect(ω0 + (0:Nt-1)*Δω)
_, i1 = findmin(abs(ω-1.0))

ω1 = ω[i1]
um = Xmfie[:,i1] / fouriertransform(gaussian)(ω1)

fcrm, geo = facecurrents(um, X)
t4 = patch(geo, real.(norm.(fcrm)))
PlotlyJS.plot(t4)


text:
"""
    Layer{F, N}

Consists of `weights` and `bias` for linear mapping, and `activation` for nonlinear mapping.
### Fields
 - `weights::Matrix{N}`
 - `bias::Vector{N}`
 - `activation::F`

See also: [`Network`](@ref)
"""
struct Layer{F<:ActivationFunction, N<:Number}
    weights::Matrix{N}
    bias::Vector{N}
    activation::F
end

"""
A Vector of layers.

    Network([layer1, layer2, layer3, ...])

See also: [`Layer`](@ref)
"""
struct Network
    layers::Vector{Layer} # layers includes output layer
end

"""
    n_nodes(L::Layer)

Returns the number of neurons in a layer.
"""
n_nodes(L::Layer) = length(L.bias)


text:
using AlgebraPDF
using Test

@testset "FunctionWithParameters{NamedTuples}" begin
    d0 = FunctionWithParameters(
        (x;p)->p.a+cos(x)*p.b;
        p=(a=2,b=1))
    #
    @test freepars(d0) == (a=2,b=1)
    @test fixedpars(d0) == ∅
    @test pars(d0) == (a=2,b=1)

    @test nfreepars(d0) == 2

    @test func(d0,0) == 3
    @test func(d0,repeat([0],10)) == repeat([3],10)
    @test func(d0, 1:10) == func(d0, collect(1:10))

    d1 = updatepar(d0, :a, 4)
    @test freepars(d1) == (a=4,b=1)
    @test releasepar(d0, :b) == d0
    @test_throws ArgumentError fixpar(d0, :b)
    #
    @test v2p([1,2],d0) == (a=1,b=2)
    @test p2v(d0) == [2,1]
    @test p2v((b=6,a=5,c=33), d0) == [5,6]
end

@testset "FunctionWithParameters{FlaggedNamedTuple}" begin
    g = FunctionWithParameters(
        (x;p)->exp((x-p.μ)^3/(2*p.σ^2));
        p=Ext(μ=1.2, σ=0.1))
    #
    g′ = fixpar(g, :μ)
    @test fixedpars(g′).μ == 1.2
    g′ = fixpar(g, :μ, 1.3)
    @test fixedpars(g′).μ == 1.3
    #
    @test freepars(g′) == (σ=0.1,)
    @test fixedpars(g′) == (μ=1.3,)
    #
    @test g′(1.3; p=(σ=0.2,)) == 1.0 # now works since μ is fixed
    @test g′(1.3; p=(σ=0.2, μ=1.6)) == 1.0 # will use the fixed value of μ=1.2
    g′′ = updatepar(g′, :μ, 1.2) # however, the update fill do
    @test g′′(1.2) == 1.0
    #
    g′′′ = releasepar(g′′, :μ)
    @test g′′′ == g # true
end

@testset "FunctionWithParameters{FlaggedNamedTuple}" begin
    m0 = FunctionWithParameters(
        (x;p)->p.a+cos(x)*p.b;
        p=Ext(a=2,b=1))
    #
    m1 = updatepar(m0, :a, 4)
    @test freepars(m1) == (a=4,b=1)

    m2 = fixpar(m1, :a)
    @test freepars(m2) == (b=1,)
    @test fixedpars(m2) == (a=4,)
    #
    m3 = releasepar(m0, :a)
    @test m3 == m0
    #
    fixpars(m1, (:a,:b)) == fixpars(m1, [:a,:b])
    #
    @test m1 == releasepars(m2, (:a,))
    @test m1 == releasepars(m2, [:a])
    #
    m4 = updatepars(m0, (a=4, b=5))
    @test freepars(m4) == (a=4,b=5)
end

@testset "func on scalars" begin
    scalar_function = 1.1
    @test func(scalar_function, 3.3) == 1.1
    @test pars(scalar_function, rand([true,false])) == ∅
    #
    Function_function = sin
    @test func(Function_function, 0.0) ≈ 0.0
    @test pars(Function_function, rand([true,false])) == ∅
end

@testset "FAbs2" begin
    m0 = FunctionWithParameters(
        (x;p)->p.a+cos(x)*p.b;
        p=Ext(a=2,b=1))
    #
    am0 = abs2(m0)
    @test func(am0,π/2) == 4
    am1 = updatepar(am0, :a, 5)
    @test func(am1,π/2) == 25
    @test freepars(am0) == freepars(m0)
end

@testset "FLog" begin
    m0 = FunctionWithParameters(
        (x;p)->exp(p.b*x);
        p=Ext(b=2,))
    #
    am0 = log(m0)
    @test am0(1) == 2
    am1 = updatepar(am0, :b, 4)
    @test am1(π/2) ≈ 2π
    @test freepars(am0) == freepars(m0)
end

@testset "NegativeLogLikelihood" begin
    d = Normalized(abs2(FBreitWigner((m=1.1, Γ=0.3))), (-3, 5))
    N = 1000
    data = randn(N)
    nll = NegativeLogLikelihood(d, data)
    #
    @test nll == minussum(log(d), data)
    @test nll(0.0) != 0.0
    @test nll(0.0) == nll(10.0)
    @test nll(0.0; p=(m=1, Γ=2)) == nll(0.0, [1,2])
    @test pars(nll) == pars(d)
    #
    g = Normalized(abs2(FBreitWigner((M=1.1, G=0.3))), (-3, 5))
    s = d+g
    nllsum = NegativeLogLikelihood(s, data)
    @test nllsum(0.0; p=pars(nllsum)+(α1=0.3, α2=0.7)) ≈ nll(0.0)
    #
    enll = Extended(NegativeLogLikelihood(s, data))
    @test pars(enll) == pars(s)
    #
    enll_unit = updatepars(enll, (α1=300, α2=700))
    @test enll_unit(0.0) ≈ nll(0.0) - log(N)*N + N
end

@testset "FSum and FProd" begin
    a1 = FunctionWithParameters(
        (x;p)->p.a+cos(x)*p.b; p=Ext(a=2,b=1))
    a2 = FunctionWithParameters(
        (x;p)->p.a+sin(x)*p.b; p=Ext(a=2,b=1))
    #
    sum1 = a1+a2
    sum2 = +(a1,a2;p=(c1=1.0,c2=1.0))
    sum3 = +(a1,a2;p=Ext(d1=1.0, d2=1.0))

    @test func(sum1,2) == func(sum2,2)
    @test func(sum1,2) == func(sum3,2)

    @test freepars(sum1) == (a=2, b=1, α1=1.0, α2=1.0)
    @test freepars(sum2) == (a=2, b=1, c1=1.0, c2=1.0)
    @test freepars(sum3) == (a=2, b=1, d1=1.0, d2=1.0)
    #
    @test freepars(fixpar(sum3,:a)) == (b=1, d1=1.0, d2=1.0)
    #
    sum4 = fixpar(sum3, :a, 3)
    @test fixedpars(sum4) == (a=3,)
    @test fixedpars(sum4[1]) == (a=3,)
    @test fixedpars(sum4[2]) == (a=3,)
    #
    sum5 = fixpar(sum3,:d2, 2)
    @test func(sum5,3.3) == func(a1,3.3) + 2* func(a2,3.3)
    #
    sum6 = a1-a2
    @test sum6(1.1) == a1(1.1)-a2(1.1)
    #
    prd = a1*a2
    @test prd(1.1) == a1(1.1)*a2(1.1)
end

@testset "no-parameters f and no-parameters normalized f" begin
    d0 = FunctionWithParameters((e;p)->e^2+p.a; p=(a=1.0,))
    f0 = noparsf(d0)
    xv = -1:10
    @test func(d0,xv;p=freepars(d0)) ≈ f0(xv)
end

# implementation with NAMES of parameters build into the funciton call
struct BW1{P} <: AbstractFunctionWithParameters
    p::P
end
import AlgebraPDF:func
func(bw::BW1, x::NumberOrTuple; p=pars(bw)) = p.m*p.Γ/(p.m^2-x^2-1im*p.m*p.Γ)

@testset "User-def type with fixed name" begin
    bw = BW1((m=3.1,Γ=0.1))
    @test pars(bw).m == 3.1
    @test pars(bw).Γ == 0.1
    @test bw(1.1) != 0.0
end

# implementation with ORDER of parameters build into the funciton call
struct BW2{P} <: AbstractFunctionWithParameters
    p::P
end
import AlgebraPDF:func
function func(bw::BW2, x::NumberOrTuple; p=pars(bw))
    m,Γ = (getproperty(p,s) for s in keys(bw.p))
    m*Γ/(m^2-x^2-1im*m*Γ)
end

@testset "User-def type with fixed order" begin
    bw = BW2((m1=3.1,Γ1=0.1))
    #
    @test pars(bw).m1 == 3.1
    @test pars(bw).Γ1 == 0.1
    @test real(bw(3.1)) ≈ 0.0
    @test real(bw(3.1; p=(m1=1.1,Γ1=3.3))) != 0
    #
    bw = BW2((m2=3.1,Γ2=0.1))
    @test pars(bw).m2 == 3.1
    @test pars(bw).Γ2 == 0.1
end

@makefuntype SuperF(x;p) = x^2+p.a*x^3
g = SuperF(p=(a=0.5,))

@testset "@maketype macro" begin
    @test keys(pars(g)) == (:a,)
    @test g(1) == 1.5
end


text:
module TestInterface
    using Test
    using CBSOData4

    include("mock_getodata.jl")

    tbl = CBSOData4.ODataTable_long("80072ned")
    @test length(tbl) == 4409
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct GasTG <: TurbineGov
        R::Float64
        T1::Float64
        T2::Float64
        T3::Float64
        AT::Float64
        Kt::Float64
        V_lim::Tuple{Float64, Float64}
        D_turb::Float64
        P_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        states_types::Vector{StateTypes.StateType}
        internal::InfrastructureSystemsInternal
    end

Parameters of Gas Turbine-Governor. GAST in PSSE and GAST_PTI in PowerWorld.

# Arguments
- `R::Float64`: Speed droop parameter, validation range: `(eps(), 0.1)`
- `T1::Float64`: Governor time constant in s, validation range: `(eps(), 0.5)`
- `T2::Float64`: Combustion chamber time constant, validation range: `(eps(), 0.5)`
- `T3::Float64`: Load limit time constant (exhaust gas measurement time), validation range: `(eps(), 5)`
- `AT::Float64`: Ambient temperature load limit, validation range: `(0, 1)`
- `Kt::Float64`: Load limit feedback gain, validation range: `(0, 5)`
- `V_lim::Tuple{Float64, Float64}`: Operational control limits on fuel valve opening (V_min, V_max)
- `D_turb::Float64`: Speed damping coefficient of gas turbine rotor, validation range: `(0, 0.5)`
- `P_ref::Float64`: Reference Load Set-point, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states of the GAST model are:
	x_g1: Fuel valve opening,
	x_g2: Fuel flow,
	x_g3: Exhaust temperature load
- `n_states::Int`: GasTG has 3 states
- `states_types::Vector{StateTypes.StateType}`: GAST has 3 differential states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct GasTG <: TurbineGov
    "Speed droop parameter"
    R::Float64
    "Governor time constant in s"
    T1::Float64
    "Combustion chamber time constant"
    T2::Float64
    "Load limit time constant (exhaust gas measurement time)"
    T3::Float64
    "Ambient temperature load limit"
    AT::Float64
    "Load limit feedback gain"
    Kt::Float64
    "Operational control limits on fuel valve opening (V_min, V_max)"
    V_lim::Tuple{Float64, Float64}
    "Speed damping coefficient of gas turbine rotor"
    D_turb::Float64
    "Reference Load Set-point"
    P_ref::Float64
    ext::Dict{String, Any}
    "The states of the GAST model are:
	x_g1: Fuel valve opening,
	x_g2: Fuel flow,
	x_g3: Exhaust temperature load"
    states::Vector{Symbol}
    "GasTG has 3 states"
    n_states::Int
    "GAST has 3 differential states"
    states_types::Vector{StateTypes.StateType}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function GasTG(R, T1, T2, T3, AT, Kt, V_lim, D_turb, P_ref=1.0, ext=Dict{String, Any}(), )
    GasTG(R, T1, T2, T3, AT, Kt, V_lim, D_turb, P_ref, ext, [:x_g1, :x_g2, :x_g3], 3, [StateTypes.Differential, StateTypes.Differential, StateTypes.Differential], InfrastructureSystemsInternal(), )
end

function GasTG(; R, T1, T2, T3, AT, Kt, V_lim, D_turb, P_ref=1.0, ext=Dict{String, Any}(), states=[:x_g1, :x_g2, :x_g3], n_states=3, states_types=[StateTypes.Differential, StateTypes.Differential, StateTypes.Differential], internal=InfrastructureSystemsInternal(), )
    GasTG(R, T1, T2, T3, AT, Kt, V_lim, D_turb, P_ref, ext, states, n_states, states_types, internal, )
end

# Constructor for demo purposes; non-functional.
function GasTG(::Nothing)
    GasTG(;
        R=0,
        T1=0,
        T2=0,
        T3=0,
        AT=0,
        Kt=0,
        V_lim=(0.0, 0.0),
        D_turb=0,
        P_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`GasTG`](@ref) `R`."""
get_R(value::GasTG) = value.R
"""Get [`GasTG`](@ref) `T1`."""
get_T1(value::GasTG) = value.T1
"""Get [`GasTG`](@ref) `T2`."""
get_T2(value::GasTG) = value.T2
"""Get [`GasTG`](@ref) `T3`."""
get_T3(value::GasTG) = value.T3
"""Get [`GasTG`](@ref) `AT`."""
get_AT(value::GasTG) = value.AT
"""Get [`GasTG`](@ref) `Kt`."""
get_Kt(value::GasTG) = value.Kt
"""Get [`GasTG`](@ref) `V_lim`."""
get_V_lim(value::GasTG) = value.V_lim
"""Get [`GasTG`](@ref) `D_turb`."""
get_D_turb(value::GasTG) = value.D_turb
"""Get [`GasTG`](@ref) `P_ref`."""
get_P_ref(value::GasTG) = value.P_ref
"""Get [`GasTG`](@ref) `ext`."""
get_ext(value::GasTG) = value.ext
"""Get [`GasTG`](@ref) `states`."""
get_states(value::GasTG) = value.states
"""Get [`GasTG`](@ref) `n_states`."""
get_n_states(value::GasTG) = value.n_states
"""Get [`GasTG`](@ref) `states_types`."""
get_states_types(value::GasTG) = value.states_types
"""Get [`GasTG`](@ref) `internal`."""
get_internal(value::GasTG) = value.internal

"""Set [`GasTG`](@ref) `R`."""
set_R!(value::GasTG, val) = value.R = val
"""Set [`GasTG`](@ref) `T1`."""
set_T1!(value::GasTG, val) = value.T1 = val
"""Set [`GasTG`](@ref) `T2`."""
set_T2!(value::GasTG, val) = value.T2 = val
"""Set [`GasTG`](@ref) `T3`."""
set_T3!(value::GasTG, val) = value.T3 = val
"""Set [`GasTG`](@ref) `AT`."""
set_AT!(value::GasTG, val) = value.AT = val
"""Set [`GasTG`](@ref) `Kt`."""
set_Kt!(value::GasTG, val) = value.Kt = val
"""Set [`GasTG`](@ref) `V_lim`."""
set_V_lim!(value::GasTG, val) = value.V_lim = val
"""Set [`GasTG`](@ref) `D_turb`."""
set_D_turb!(value::GasTG, val) = value.D_turb = val
"""Set [`GasTG`](@ref) `P_ref`."""
set_P_ref!(value::GasTG, val) = value.P_ref = val
"""Set [`GasTG`](@ref) `ext`."""
set_ext!(value::GasTG, val) = value.ext = val
"""Set [`GasTG`](@ref) `states_types`."""
set_states_types!(value::GasTG, val) = value.states_types = val


text:
##
# Check that patches between BinaryBuilder and Julia are consistent.
# julia check_patches.jl .../path/to/julia/patchdir
##

@assert length(ARGS) == 1

DIR = joinpath(@__DIR__, "patches", "llvm_patches")
patches = readdir(DIR)

files = Dict(map(x -> x => joinpath(DIR, x), patches))
orig_files = Dict(map(x -> x => joinpath(ARGS[1], x[6:end]), patches))

for patch in patches
    orig = orig_files[patch]
    file = files[patch]

    if !isfile(orig)
        warn("$patch is not present in $(ARGS[1])")
        continue
    end
    try
        run(`diff $orig $file`)
    catch
        @warn("$patch differs")
    end
end

patches = map(p->string(p[6:end]), patches)

for patch in readdir(ARGS[1])
    if startswith(patch, "llvm") && !(patch in patches)
        warn("$patch not present for LLVMBuilder")
    end
end

max_stars_repo_name: vchuravy/LLVMBuilder

max_issues_repo_name: vchuravy/LLVMBuilder

max_forks_repo_name: vchuravy/LLVMBuilder