text:
### LX_PAT
@testset "latex" begin
    for s in (
            raw"\com",
            raw"  \com ",
        )
        m = match(F.LX_NAME_PAT, s)
        @test m.captures[1] == raw"\com"
    end
    for s in (
            raw"\cöm",
            raw"  \cöm ",
        )
        m = match(F.LX_NAME_PAT, s)
        @test m.captures[1] == raw"\cöm"
    end
    m = match(F.LX_NAME_PAT, raw"\com*")
    @test m.captures[1] == raw"\com*"
    for s in (
            raw"[2]",
            raw" [2] ",
            raw" [ 2]",
            raw" [ 2 ] "
        )
        m = match(F.LX_NARG_PAT, s)
        @test m.captures[2] == "2"
    end
end

### Assign pat
@testset "assign" begin
    for s in (
            raw"cöm_012 = blah",
        )
        m = match(F.ASSIGN_PAT, s)
        @test m.captures[1] == raw"cöm_012"
        @test strip(m.captures[2]) == raw"blah"
    end
end

### LINK
@testset "esclink" begin
    s = Markdown.htmlesc("[hello]")
    a,b,c = match(F.ESC_LINK_PAT, s).captures
    @test isnothing(a)
    @test b == "hello"
    @test isnothing(c)
    s = Markdown.htmlesc("[hello][]")
    a,b,c = match(F.ESC_LINK_PAT, s).captures
    @test isnothing(a)
    @test b == "hello"
    @test isempty(c)
    s = Markdown.htmlesc("[hello][id]")
    a,b,c = match(F.ESC_LINK_PAT, s).captures
    @test isnothing(a)
    @test b == "hello"
    @test c == "id"
    s = Markdown.htmlesc("![hello][id]")
    a,b,c = match(F.ESC_LINK_PAT, s).captures
    @test a == Markdown.htmlesc("!")
    @test b == "hello"
    @test c == "id"
    s = Markdown.htmlesc("[hello]:")
    @test isnothing(match(F.ESC_LINK_PAT, s))
end

@testset "fndef" begin
    for s in (
            raw"[^örα_01]:",
        )
        m = match(F.FN_DEF_PAT, s)
        @test !isnothing(m)
    end
end

### HTML ENTITY
@testset "html-ent" begin
    for s in (
        "&tab;",
     	"&newline;",
        "!",
        """,
        "	",
        "'",
        "'",
        ":",
        "{",
        "ĉ",
        "ć",
        "ć"
        )
        @test !isnothing(match(F.HTML_ENT_PAT, s))
    end
end

### Code blocks
@testset "code" begin
    list = (
        "```julia:n some code```", # issue 705
        "```julia:ex some code```",
        "```cpp:πγ some code```",
        "```julia:πγ.jl some code```",
        "```julia:./πγ.jl some code```",
        "```julia:./πγ/λη.jl some code```",
        )
    for s in list
        @test !isnothing(match(F.CODE_3_PAT, s))
    end
    for s in list
        @test !isnothing(match(F.CODE_5_PAT, "``$s``"))
    end
    noklist = (
        "```γulia some code```",
        "```julia:0ex some code```",
        )
    for s in noklist
        @test isnothing(match(F.CODE_3_PAT, s))
    end
end

### HBLOCK REGEXES

@testset "hb-if" begin
    for s in (
            "{{if var1}}",
            "{{if  var1 }}",
            "{{ if  var1 }}",
            "{{ if var1 }}"
        )
        m = match(F.HBLOCK_IF_PAT, s)
        @test m.captures[1] == "var1"
    end
    for s in (
            "{{if var1 var2}}",
            "{{ifvar1}}"
        )
        m = match(F.HBLOCK_IF_PAT, s)
        @test isnothing(m)
    end
end
@testset "hb-else" begin
    for s in (
            "{{else}}",
            "{{ else}}",
            "{{  else   }}",
            "{{ else }}"
        )
        m = match(F.HBLOCK_ELSE_PAT, s)
        @test !isnothing(m)
    end
end
@testset "hb-elseif" begin
    for s in (
            "{{elseif var1}}",
            "{{else if  var1 }}",
            "{{ elseif  var1 }}",
            "{{ elseif var1 }}"
        )
        m = match(F.HBLOCK_ELSEIF_PAT, s)
        @test m.captures[1] == "var1"
    end
    for s in (
            "{{else if var1 var2}}",
            "{{elif var1}}"
        )
        m = match(F.HBLOCK_ELSEIF_PAT, s)
        @test isnothing(m)
    end
end
@testset "hb-end" begin
    for s in (
            "{{end}}",
            "{{ end}}",
            "{{  end   }}",
            "{{ end }}"
        )
        m = match(F.HBLOCK_END_PAT, s)
        @test !isnothing(m)
    end
end
@testset "hb-isdef" begin
    for s in (
            "{{isdef var1}}",
            "{{ isdef  var1 }}",
            "{{ isdef  var1 }}",
            "{{ ifdef  var1 }}",
            "{{ ifdef var1 }}"
        )
        m = match(F.HBLOCK_ISDEF_PAT, s)
        @test m.captures[1] == "var1"
    end
    for s in (
            "{{isdef var1 var2}}",
            "{{is def var1}}",
            "{{if def var1}}"
        )
        m = match(F.HBLOCK_ISDEF_PAT, s)
        @test isnothing(m)
    end
end
@testset "hb-isndef" begin
    for s in (
            "{{isnotdef var1}}",
            "{{ isndef  var1 }}",
            "{{ ifndef  var1 }}",
            "{{ ifnotdef  var1 }}",
            "{{ ifndef var1 }}"
        )
        m = match(F.HBLOCK_ISNOTDEF_PAT, s)
        @test m.captures[1] == "var1"
    end
    for s in (
            "{{isnotdef var1 var2}}",
            "{{isnot def var1}}",
            "{{ifn def var1}}"
        )
        m = match(F.HBLOCK_ISNOTDEF_PAT, s)
        @test isnothing(m)
    end
end
@testset "hb-ispage" begin
    for s in (
            "{{ispage var1 var2}}",
            "{{ ispage var1 var2 }}"
        )
        m = match(F.HBLOCK_ISPAGE_PAT, s)
        @test m.captures[1] == "var1 var2"
    end
end
@testset "hb-isnotpage" begin
    for s in (
            "{{isnotpage var1 var2}}",
            "{{ isnotpage var1 var2 }}"
        )
        m = match(F.HBLOCK_ISNOTPAGE_PAT, s)
        @test m.captures[1] == "var1 var2"
    end
end
@testset "hb-for" begin
    for s in (
            "{{for (v1,v2,v3) in iterate}}",
            "{{for (v1, v2,v3) in iterate}}",
            "{{for ( v1, v2, v3) in iterate}}",
            "{{for ( v1  , v2 , v3 ) in iterate}}",
            "{{ for (v1, v2, v3) in iterate }}"
        )
        m = match(F.HBLOCK_FOR_PAT, s)
        @test isapproxstr(m.captures[1], "(v1, v2, v3)")
    end
    s = "{{for v1 in iterate}}"
    m = match(F.HBLOCK_FOR_PAT, s)
    @test isapproxstr(m.captures[1], "v1")

    # WARNING: NOT RECOMMENDED / NEEDS CARE
    s = "{{for v1,v2 in iterate}}"
    m = match(F.HBLOCK_FOR_PAT, s)
    @test isapproxstr(m.captures[1], "v1,v2")
    s = "{{for (v1,v2 in iterate}}"
    m = match(F.HBLOCK_FOR_PAT, s)
    @test isapproxstr(m.captures[1], "(v1,v2")
    s = "{{for v1,v2) in iterate}}"
    m = match(F.HBLOCK_FOR_PAT, s)
    @test isapproxstr(m.captures[1], "v1,v2)")
end
@testset "hb-fun-0" begin
    for s in (
            "{{toc}}",
            "{{ toc }}",
            "{{ toc }}"
            )
        m = match(F.HBLOCK_FUN_PAT, s)
        @test m.captures[1] == "toc"
        @test m.captures[2] === nothing || isempty(strip(m.captures[2]))
    end
end
@testset "hb-fun-1" begin
    for s in (
            "{{fun p1}}",
            "{{ fun  p1 }}",
            "{{ fun p1 }}"
            )
        m = match(F.HBLOCK_FUN_PAT, s)
        @test m.captures[1] == "fun"
        @test strip(m.captures[2]) == "p1"
    end
end
@testset "hb-fun-2" begin
    for s in (
            "{{fun p1 p2}}",
            "{{ fun  p1 p2 }}"
            )
        m = match(F.HBLOCK_FUN_PAT, s)
        @test m.captures[1] == "fun"
        @test strip(m.captures[2]) == "p1 p2"
    end
end

# ========
# Checkers
# ========
@testset "ch-for" begin
    s = "{{for v in iterate}}"
    m = match(F.HBLOCK_FOR_PAT, s).captures[1]
    @test isnothing(F.check_for_pat(m))
    s = "{{for (v1,v2) in iterate}}"
    m = match(F.HBLOCK_FOR_PAT, s).captures[1]
    @test isnothing(F.check_for_pat(m))
    s = "{{for (v in iterate}}"
    m = match(F.HBLOCK_FOR_PAT, s).captures[1]
    @test_throws F.HTMLBlockError F.check_for_pat(m)
    s = "{{for v1,v2) in iterate}}"
    m = match(F.HBLOCK_FOR_PAT, s).captures[1]
    @test_throws F.HTMLBlockError F.check_for_pat(m)
    s = "{{for v1,v2 in iterate}}"
    m = match(F.HBLOCK_FOR_PAT, s).captures[1]
    @test_throws F.HTMLBlockError F.check_for_pat(m)
end


text:
abstract type BibliographyBlock <: Expanders.ExpanderPipeline end

Selectors.order(::Type{BibliographyBlock}) = 12.0  # Expand bibliography last
Selectors.matcher(::Type{BibliographyBlock}, node, page, doc) = Expanders.iscode(node, r"^@bibliography")

const tex2unicode_chars = Dict(
    'o' => "\u00F8",  # \o 	ø 	latin small letter O with stroke
    'O' => "\u00D8",  # \O 	Ø 	latin capital letter O with stroke
    'l' => "\u0142",  # \l 	ł 	latin small letter L with stroke
    'L' => "\u0141",  # \L 	Ł 	latin capital letter L with stroke
    'i' => "\u0131",  # \i 	ı 	latin small letter dotless I
)

const tex2unicode_replacements = (
    "---" => "—", # em dash needs to go first
    "--"  => "–",

    # replace quoted single letters before the remaining replacements, and do
    # them all at once, as these patterns rely on word boundaries which can
    # change due to the replacements we perform
    r"\\[oOlLi]\b" => c -> tex2unicode_chars[c[2]],

    r"\\`\{(\S{1})\}" => s"\1\u300", # \`{o} 	ò 	grave accent
    r"\\'\{(\S{1})\}" => s"\1\u301", # \'{o} 	ó 	acute accent
    r"\\\^\{(\S{1})\}" => s"\1\u302", # \^{o} 	ô 	circumflex
    r"\\~\{(\S{1})\}" => s"\1\u303", # \~{o} 	õ 	tilde
    r"\\=\{(\S{1})\}" => s"\1\u304", # \={o} 	ō 	macron accent (a bar over the letter)
    r"\\u\{(\S{1})\}" => s"\1\u306",  # \u{o} 	ŏ 	breve over the letter
    r"\\\.\{(\S{1})\}" => s"\1\u307", # \.{o} 	ȯ 	dot over the letter
    r"\\\\\"\{(\S{1})\}" => s"\1\u308", # \"{o} 	ö 	umlaut, trema or dieresis
    r"\\r\{(\S{1})\}" => s"\1\u30A",  # \r{a} 	å 	ring over the letter (for å there is also the special command \aa)
    r"\\H\{(\S{1})\}" => s"\1\u30B",  # \H{o} 	ő 	long Hungarian umlaut (double acute)
    r"\\v\{(\S{1})\}" => s"\1\u30C",  # \v{s} 	š 	caron/háček ("v") over the letter
    r"\\d\{(\S{1})\}" => s"\1\u323",  # \d{u} 	ụ 	dot under the letter
    r"\\c\{(\S{1})\}" => s"\1\u327",  # \c{c} 	ç 	cedilla
    r"\\k\{(\S{1})\}" => s"\1\u328",  # \k{a} 	ą 	ogonek
    r"\\b\{(\S{1})\}" => s"\1\u331",  # \b{b} 	ḇ 	bar under the letter
    r"\\t\{(\S{1})(\S{1})\}" => s"\1\u0361\2",  # \t{oo} 	o͡o 	"tie" (inverted u) over the two letters
    r"\{\}" => s"",  # empty curly braces should not have any effect
    r"\{([[:alnum:]]+)\}" => s"\1",  # {<text>} 	<text> 	bracket stripping after applying all rules

    # TODO:
    # \"{\i} 	ï 	Latin Small Letter I with Diaeresis

    # Sources : https://www.compart.com/en/unicode/U+0131 enter the unicode character into the search box
)

function tex2unicode(s)
    for replacement in tex2unicode_replacements
        s = replace(s, replacement)
    end
    return Unicode.normalize(s)
end

linkify(text, link) = isempty(link) ? text : "$text"

function Selectors.runner(::Type{BibliographyBlock}, x, page, doc)
    @info "Expanding bibliography."
    raw_bib = "
"
    for (id, entry) in doc.plugins[CitationBibliography].bib
        @info "Expanding bibliography entry: $id."

        # Add anchor that citations can link to from anywhere in the docs.
        Anchors.add!(doc.internal.headers, entry, entry.id, page.build)

        authors = xnames(entry) |> tex2unicode
        link = xlink(entry)
        title = xtitle(entry) |> tex2unicode
        published_in = xin(entry) |> tex2unicode

        raw_bib *= """
$id

$authors, $(linkify(title, link)), $published_in

"""
    end
    raw_bib *= "\n
"

    page.mapping[x] = Documents.RawNode(:html, raw_bib)
end

ext: jl


text:
export AbstractMeterProvider, AbstractInstrument, global_meter_provider, Meter

abstract type AbstractMeterProvider end

Base.@kwdef struct DummyMeterProvider <: AbstractMeterProvider end

Base.push!(p::DummyMeterProvider, x) = nothing

const GLOBAL_METER_PROVIDER = Ref(DummyMeterProvider())
global_meter_provider() = GLOBAL_METER_PROVIDER[]
global_meter_provider(p) = GLOBAL_METER_PROVIDER[] = p

abstract type AbstractInstrument{T} end

struct Meter{P<:AbstractMeterProvider}
    name::String
    provider::P
    version::VersionNumber
    schema_url::String
    instruments::Vector{AbstractInstrument}
    instrumentation_info::InstrumentationInfo
    function Meter(
        name::String;
        provider::P = global_meter_provider(),
        version = v"0.0.1-dev",
        schema_url = "",
        instrumentation_info = InstrumentationInfo(),
    ) where {P<:AbstractMeterProvider}
        m = new{P}(
            name,
            provider,
            version,
            schema_url,
            AbstractInstrument[],
            instrumentation_info,
        )
        push!(provider, m)
        m
    end
end


text:
struct ClebschGordan{A,B,C,j₁T,j₂T}
    j₁::A
    m₁::A
    j₂::B
    m₂::B
    J::C
    M::C
end

Base.zero(::Type{ClebschGordan{A,B,C,j₁T,j₂T}}) where {A,B,C,j₁T,j₂T} =
    ClebschGordan{A,B,C,j₁T,j₂T}(zero(A),zero(A),zero(B),zero(B),zero(C),zero(C))

const ClebschGordanℓs{I<:Integer} = ClebschGordan{I,Rational{I},Rational{I},:ℓ,:s}
ClebschGordanℓs(j₁::I, m₁::I, j₂::R, m₂::R, J::R, M::R) where {I<:Integer,R<:Rational{I}} =
    ClebschGordanℓs{I}(j₁, m₁, j₂, m₂, J, M)

Base.convert(::Type{T}, cg::ClebschGordan) where {T<:Real} =
    clebschgordan(cg.j₁,cg.m₁,cg.j₂,cg.m₂,cg.J,cg.M)

# We have chosen the Clebsch–Gordan coefficients to be real.
Base.adjoint(cg::ClebschGordan) = cg

Base.show(io::IO, cg::ClebschGordan) =
    write(io, "⟨$(rs(cg.j₁)),$(rs(cg.j₂));$(rs(cg.m₁)),$(rs(cg.m₂))|$(rs(cg.J)),$(rs(cg.M))⟩")

function Base.show(io::IO, cg::ClebschGordanℓs)
    spin = cg.m₂ > 0 ? "↑" : "↓"
    write(io, "⟨$(spectroscopic_label(cg.j₁));$(cg.m₁),$(spin)|$(rs(cg.J)),$(rs(cg.M))⟩")
end

#=
Relevant equations taken from

- Sakurai, J. J., & Napolitano, J. J. (2010). Modern quantum mechanics
(2nd edition). : Addison-Wesley.
=#

#=

Repeated action [Eq. (3.8.45)] with the ladder operators on the jmⱼ
basis

J±|ℓs;jmⱼ⟩

gives the following recursion relations, Eq. (3.8.49):

√((j∓mⱼ)(j±mⱼ+1))⟨ℓs;mₗ,mₛ|ℓs;j,mⱼ±1⟩
= √((ℓ∓mₗ+1)(ℓ±mₗ))⟨ℓs;mₗ∓1,mₛ|ℓs;j,mⱼ⟩
+ √((s∓mₛ+1)(s±mₛ))⟨ℓs;mₗ,mₛ∓1|ℓs;j,mⱼ⟩.

We know [Eq. (3.8.58)] that

⟨ℓs;±ℓ,±1/2|ℓs;ℓ±1/2,ℓ±1/2⟩ = 1.

=#
function rotate!(blocks::Vector{M}, orbs::RelativisticOrbital...) where {T,M<:AbstractMatrix{T}}
    2length(blocks) == sum(degeneracy.(orbs))-2 ||
        throw(ArgumentError("Invalid block size for $(orbs...)"))
    ℓ = kappa_to_ℓ(first(orbs).κ)
    # We sort by mⱼ and remove the first and last elements since they
    # are pure and trivially unity.
    ℓms = sort(vcat([[(ℓ,m,s) for m ∈ -ℓ:ℓ] for s = -1//2:1//2]...)[2:end-1], by=((ℓms)) -> (ℓms[2],+(ℓms[2:3]...)))
    jmⱼ = sort(vcat([[(j,mⱼ) for mⱼ ∈ -j:j] for j ∈ [convert(Rational, kappa_to_j(o.κ)) for o in orbs]]...), by=last)[2:end-1]
    for (a,(ℓ,m,s)) in enumerate(ℓms)
        bi = cld(a, 2)
        o = 2*(bi-1)
        for (b,(j,mⱼ)) in enumerate(jmⱼ)
            b-o ∉ 1:2 && continue
            blocks[bi][a-o,b-o] = ClebschGordanℓs(ℓ,m,1//2,s,j,mⱼ)
        end
    end
    blocks
end

"""
    jj2lsj([T=Float64, ]orbs...)

Generates the block-diagonal matrix that transforms jj-coupled
configurations to lsj-coupled ones.

The blocks correspond to invariant subspaces, which rotate among
themselves to form new linear combinations.

They are sorted according to n,ℓ and within the blocks, the columns
are sorted according to mⱼ,j (increasing) and the rows according to
s,ℓ,m (increasing).

E.g. the p-block will have the following structure:

```
 ⟨p;-1,↓|3/2,-3/2⟩  │       ⋅                  ⋅             │       ⋅                ⋅           │       ⋅
 ───────────────────┼────────────────────────────────────────┼────────────────────────────────────┼─────────────────
      ⋅             │  ⟨p;0,↓|1/2,-1/2⟩   ⟨p;0,↓|3/2,-1/2⟩   │       ⋅                ⋅           │       ⋅
      ⋅             │  ⟨p;-1,↑|1/2,-1/2⟩  ⟨p;-1,↑|3/2,-1/2⟩  │       ⋅                ⋅           │       ⋅
 ───────────────────┼────────────────────────────────────────┼────────────────────────────────────┼─────────────────
      ⋅             │       ⋅                  ⋅             │  ⟨p;1,↓|1/2,1/2⟩  ⟨p;1,↓|3/2,1/2⟩  │       ⋅
      ⋅             │       ⋅                  ⋅             │  ⟨p;0,↑|1/2,1/2⟩  ⟨p;0,↑|3/2,1/2⟩  │       ⋅
 ───────────────────┼────────────────────────────────────────┼────────────────────────────────────┼─────────────────
      ⋅             │       ⋅                  ⋅             │       ⋅                ⋅           │  ⟨p;1,↑|3/2,3/2⟩
```

"""
function jj2lsj(::Type{T}, orbs::RelativisticOrbital...) where T
    nℓs = map(o -> (o.n, kappa_to_ℓ(o.κ)), sort([orbs...]))
    blocks = map(unique(nℓs)) do (n,ℓ)
        i = findall(isequal((n,ℓ)), nℓs)
        subspace = orbs[i]
        mⱼ = map(subspace) do orb
            j = convert(Rational, kappa_to_j(orb.κ))
            -j:j
        end

        jₘₐₓ = maximum([convert(Rational, kappa_to_j(o.κ)) for o in subspace])
        pure = [Matrix{T}(undef,1,1),Matrix{T}(undef,1,1)]
        pure[1][1]=ClebschGordanℓs(ℓ,-ℓ,1//2,-1//2,jₘₐₓ,-jₘₐₓ)
        pure[2][1]=ClebschGordanℓs(ℓ,ℓ,1//2,1//2,jₘₐₓ,jₘₐₓ)

        n = sum(length.(mⱼ))-2
        if n > 0
            nblocks = div(n,2)
            b = [zeros(T,2,2) for i = 1:nblocks]
            rotate!(b, subspace...)
            [pure[1],b...,pure[2]]
        else
            pure
        end
    end |> b -> vcat(b...)
    rows = size.(blocks,1)
    N = sum(rows)
    R = BlockBandedMatrix(Zeros{T}(N,N), (rows,rows), (0,0))
    for (i,b) in enumerate(blocks)
        R[Block(i,i)] .= b
    end
    R
end
jj2lsj(orbs::RelativisticOrbital...) = jj2lsj(Float64, orbs...)

function jj2lsj(sorb::SpinOrbital)
    @unpack orb,mℓ,spin = sorb
    @unpack n,ℓ = orb
    ms = (spin ? 1 : -1)//2
    mj = mℓ + ms
    map(max(abs(mj),ℓ-1//2):(ℓ+1//2)) do j
        (RelativisticOrbital(n,ℓ,j),mj)=>clebschgordan(ℓ,mℓ,1//2,ms,j,mj)
    end
end

export jj2lsj, ClebschGordan, ClebschGordanℓs


text:
module FluidSynth

# Load shared lib
deps = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(deps)
    include(deps)
else
    error("fluidsynth not properly installed. Please run Pkg.build(\"fluidsynth\")")
end

export FLUID_OK, FLUID_FAILED,
    Settings, setstr, setint,
    Synth, sfload, program_select, noteoff, noteon, write_s16_stereo, write_s16, write_float, set_gain, set_sample_rate,
    Player, add, play, stop, join, get_status,
    PlayerStatus, Playing, Ready, Done,
    FileRenderer, process_block,
    AudioDriver

const FLUID_OK = 0
const FLUID_FAILED = -1

mutable struct Settings
    settings_ptr::Ptr{Cvoid}

    function Settings()
    settings = new(C_NULL)
    settings.settings_ptr = ccall((:new_fluid_settings, libfluidsynth), Ptr{Cvoid}, ())
    finalizer(settings) do settings
    ccall((:delete_fluid_settings, libfluidsynth), Cvoid, (Ptr{Cvoid},), settings.settings_ptr)
    end
    settings
    end
end

function setstr(settings::Settings, name::String, str::String)
    ret = ccall((:fluid_settings_setstr, libfluidsynth), Cint, (Ptr{Cvoid}, Cstring, Cstring), settings.settings_ptr, name, str)
    if ret == FLUID_FAILED
    throw(ErrorException(""))
    end
    ret
end

function setint(settings::Settings, name::String, val::Int32)
    ret = ccall((:fluid_settings_setint, libfluidsynth), Cint, (Ptr{Cvoid}, Cstring, Cint), settings.settings_ptr, name, val)
    if ret == FLUID_FAILED
    throw(ErrorException("fail setint"))
    end
    ret
end

mutable struct Synth
    synth_ptr::Ptr{Cvoid}
    settings::Settings

    function Synth(settings::Settings)
    synth = new(C_NULL, settings)

    synth.synth_ptr = ccall((:new_fluid_synth, libfluidsynth), Ptr{Cvoid}, (Ptr{Cvoid},), settings.settings_ptr)
    finalizer(synth) do synth
    ccall((:delete_fluid_synth, libfluidsynth), Cvoid, (Ptr{Cvoid},), synth.synth_ptr)
    end
    synth
    end
end

# fluidsynth interfaces.
"""
Load SoundFont file.
Return SoundFont ID.

# Arguments
- `filename` : path to SoundFont file.
- `reset_presets` : true to re-assign presets for all MIDI channels
"""
function sfload(synth::Synth, filename::AbstractString, reset_presets=false::Bool)
    c_reset_preset = 0
    if reset_presets
    c_reset_preset = 1
    end
    ret = ccall((:fluid_synth_sfload, libfluidsynth), Cint, (Ptr{Cvoid}, Cstring, Cint), synth.synth_ptr, filename, c_reset_preset);
    if ret == FLUID_FAILED
    throw(ErrorException(""))
    end
    ret
end

"""
Select an instrument on a MIDI channel by SoundFont ID, bank and program numbers.
# Arguments
- `chan` : MIDI channel number
- `sfont_id` : ID of loaded SoundFont
- `bank_num` : MIDI bank number
- `preset_num` : MIDI program number
"""
function program_select(synth::Synth, chan::Int32, sfont_id::Int32, bank_num::Int32, preset_num::Int32)
    ret = ccall((:fluid_synth_program_select, libfluidsynth), Cint, (Ptr{Cvoid}, Cint, Cint, Cint, Cint), synth.synth_ptr, chan, sfont_id, bank_num, preset_num)
    if ret == FLUID_FAILED
    throw(ErrorException(""))
    end
end

"""
Send a note-on event.
# Arguments
- `chan` : MIDI channel number(0-127)
- `key` : MIDI note number(0-127)
- `vel` : MIDI velocity(0-127, 0=noteoff)
"""
function noteon(synth::Synth, chan::Int32, key::Int32, vel::Int32)
    ret = ccall((:fluid_synth_noteon, libfluidsynth), Cint, (Ptr{Cvoid}, Cint, Cint, Cint), synth.synth_ptr, chan, key, vel)
    if ret == FLUID_FAILED
    throw(ErrorException(""))
    end
end

"""
Send a note-off event.
# Arguments
- `chan` : MIDI channel number(0-127)
- `key` : MIDI note number(0-127)
"""
function noteoff(synth::Synth, chan::Int32, key::Int32)
    ret = ccall((:fluid_synth_noteoff, libfluidsynth), Cint, (Ptr{Cvoid}, Cint, Cint), synth.synth_ptr, chan, key)
    if ret == FLUID_FAILED
    throw(ErrorException("no matched noteon event"))
    end
end

"""
Set the MIDI pitch bend controller value on a MIDI channel.
# Arguments
- `chan` : MIDI channel number(0-127)
- `val` : MIDI pitch bend value(0-16383 with 8192 being center)
"""
function pitch_bend(synth::Synth, chan::Int32, val::Int32)
    ret = ccall((:fluid_synth_pitch_bend, libfluidsynth), Cint, (Ptr{Cvoid}, Cint, Cint), synth.synth_ptr, chan, val)
    if ret == FLUID_FAILED
    throw(ErrorException(""))
    end
end

"""
Send a MIDI controller event on a MIDI channel.
# Arguments
- `chan` : MIDI channel number(0-127)
- `num` : MIDI controller number(0-127)
- `val` : MIDI controller value(0-127)
"""
function cc(synth::Synth, chan::Int32, num::Int32, val::Int32)
    ret = ccall((:fluid_synth_cc, libfluidsynth), Cint, (Ptr{Cvoid}, Cint, Cint, Cint), synth.synth_ptr, chan, num, val)
    if ret == FLUID_FAILED
    throw(ErrorException(""))
    end
end

function set_gain(synth::Synth, gain::Float32)
    ccall((:fluid_synth_set_gain, libfluidsynth), Cvoid, (Ptr{Cvoid}, Cfloat), synth.synth_ptr, gain)
end

function set_sample_rate(synth::Synth, sample_rate::Float32)
    ccall((:fluid_synth_set_gain, libfluidsynth), Cvoid, (Ptr{Cvoid}, Cfloat), synth.synth_ptr, sample_rate)
end

"""
Synthesize a block of 16 bit audio samples to audio buffers.
# Arguments
- `len` : Count of audio frames to synthesize
- `lout::Vector{Int16}` : Vector of signed 16 bit to store left channel of audio
- `loff` : Offset index in 'lout' for first sample
- `lincr` : Increment between samples stored to 'lout'.
- `rout::Vector{Int16}` : Vector of signed 16 bit to store right channel of audio
- `roff` : Offset index in 'rout' for first sample
- `rincr` : Increment between samples stored to 'rout'.
"""
function write_s16(synth::Synth, len::Int32, lout::Vector{Int16}, loff::Int32, lincr::Int32, rout::Vector{Int16}, roff::Int32, rincr::Int32)
    ret = ccall((:fluid_synth_write_s16, libfluidsynth),
    Cint,
    (Ptr{Cvoid}, Cint, Ptr{Cvoid}, Cint, Cint, Ptr{Cvoid}, Cint, Cint),
    synth.synth_ptr, len, lout, loff, lincr, rout, roff, rincr)
    if ret == FLUID_FAILED
    throw(ErrorException(""))
    end
end

function write_float(synth::Synth, len::Int32, lout::Vector{Float32}, loff::Int32, lincr::Int32, rout::Vector{Float32}, roff::Int32, rincr::Int32)
    ret = ccall((:fluid_synth_write_float, libfluidsynth),
    Cint,
    (Ptr{Cvoid}, Cint, Ptr{Cvoid}, Cint, Cint, Ptr{Cvoid}, Cint, Cint),
    synth.synth_ptr, len, lout, loff, lincr, rout, roff, rincr)
    if ret == FLUID_FAILED
    throw(ErrorException(""))
    end
end

# player
mutable struct Player
    player_ptr::Ptr{Cvoid}
    synth::Synth

    function Player(synth::Synth)
 		player = new(C_NULL, synth)
    player.player_ptr = ccall((:new_fluid_player, libfluidsynth), Ptr{Cvoid}, (Ptr{Cvoid},), synth.synth_ptr)
    finalizer(player) do player
    ccall((:delete_fluid_player, libfluidsynth), Cvoid, (Ptr{Cvoid},), player.player_ptr)
    end
    player
    end
end

"""
Add a MIDI file to a player queue.
"""
function add(player::Player, midifile::String)
    ret = ccall((:fluid_player_add, libfluidsynth),
    Cint,
    (Ptr{Cvoid}, Cstring),
    player.player_ptr, midifile)
    if ret == FLUID_FAILED
    throw(ErrorException(""))
    end
end

function play(player::Player)
    ret = ccall((:fluid_player_play, libfluidsynth),
    Cint,
    (Ptr{Cvoid}, ),
    player.player_ptr)
    if ret == FLUID_FAILED
    throw(ErrorException(""))
    end
end

function stop(player::Player)
    ccall((:fluid_player_stop, libfluidsynth),
    Cint,
    (Ptr{Cvoid}, ),
    player.player_ptr)
    # this metod always reutrn ok.
end

function join(player::Player)
    ret = ccall((:fluid_player_join, libfluidsynth),
    Cint,
    (Ptr{Cvoid}, ),
    player.player_ptr)
    # this metod always reutrn ok.
end

@enum PlayerStatus begin
    Ready = 0
    Playing = 1
    Done = 2
end

function get_status(player::Player)
    ret = ccall((:fluid_player_get_status, libfluidsynth), Cint, (Ptr{Cvoid},), player.player_ptr)
    if ret == Int(Ready::PlayerStatus)
    return Ready::PlayerStatus
    elseif ret == Int(Playing::PlayerStatus)
    return Playing::PlayerStatus
    else
    return Done::PlayerStatus
    end

end

# file renderer
mutable struct FileRenderer

    ptr::Ptr{Cvoid}
    synth::Synth

    function FileRenderer(synth::Synth)
    file_renderer = new(C_NULL, synth)

    file_renderer.ptr = ccall((:new_fluid_file_renderer, libfluidsynth), Ptr{Cvoid}, (Ptr{Cvoid},), synth.synth_ptr)
    finalizer(file_renderer) do file_renderer
    ccall((:delete_fluid_file_renderer, libfluidsynth), Cvoid, (Ptr{Cvoid},), file_renderer.ptr)
    end
    file_renderer
    end
end

function process_block(file_renderer::FileRenderer)
    ccall((:fluid_file_renderer_process_block, libfluidsynth), Cint, (Ptr{Cvoid},), file_renderer.ptr)
end

# file renderer
mutable struct AudioDriver

    ptr::Ptr{Cvoid}
    settings::Settings
    synth::Synth

    function AudioDriver(settings::Settings, synth::Synth)
    audio_driver = new(C_NULL, settings, synth)

    audio_driver.ptr = ccall((:new_fluid_audio_driver, libfluidsynth), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}), settings.settings_ptr, synth.synth_ptr)
    finalizer(audio_driver) do audio_driver
        ccall((:delete_fluid_audio_driver, libfluidsynth), Cvoid, (Ptr{Cvoid},), audio_driver.ptr)
    end
    audio_driver
    end
end

# support methods
"""
Synthesize a 16 bit stereo audio samples.
Retuen synthesized vector. vector size is 2 * len(Because of stereo).
# Arguments
- `len::Int32` : Count of audio frames to synthesize.
"""
function write_s16_stereo(synth::Synth, len::Int32)
    buf = Vector{Int16}(undef, 2 * len)
    write_s16(synth, len, buf, Int32(0) , Int32(2), buf, Int32(1), Int32(2))
    return buf
end

end # module


text:
using MOCNeutronTransport
@testset "VTK" begin
    @testset "Triangles" begin
        # Three triangles
        #      4---------5
        #     / \       / \
        #    /   \     /   \
        #   /     \   /     \
        #  /       \ /       \
        # 1---------2---------3
        filepath = "./mesh/mesh_files/three_triangles.vtk"
        ref_points = [Point_2D(0.0), Point_2D(2.0), Point_2D(4.0), Point_2D(1.0, 1.0), Point_2D(3.0, 1.0)]
        ref_faces = [
                        SVector(5, 1, 2, 4),
                        SVector(5, 2, 5, 4),
                        SVector(5, 2, 3, 5)
                    ]
        ref_edges = [SVector(1, 2), SVector(1, 4), SVector(2, 3),
                     SVector(2, 4), SVector(2, 5), SVector(3, 5), SVector(4, 5)]

        # read_vtk
        mesh = read_vtk_2d(filepath)
        @test mesh.points == ref_points
        @test mesh.faces == ref_faces
        @test mesh.name == "three_triangles"

        # write_vtk
        write_vtk_2d("write_three_triangle.vtk", mesh)
        ref_file = open(filepath, "r")
        test_file = open("./write_three_triangle.vtk", "r")
        while !eof(ref_file)
            ref_line = readline(ref_file)
            test_line = readline(test_file)
            @test ref_line == test_line
        end
        close(ref_file)
        close(test_file)
        rm("./write_three_triangle.vtk")
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
#######################################################################
# Functions for processing information about packages
# and sorting them by the number of stars
#######################################################################

#--------------------------------------------------
# Data
#--------------------------------------------------

# Expected format of JSON data: `{packages: [...]}`
const PKGS_KEY = "packages"

# Expected format of a single package:
# `{name: ..., metadata: {repo:..., starcount: ...}}`
const NAME_KEY          = "name"
const METADATA_KEY      = "metadata"
const REPO_KEY          = "repo"
const STARTCOUNT_KEY    = "starcount"
const VERSIONS_KEY      = "versions"

# Placeholders if data is not available
const NAME_NA           = "<NA-name>"
const REPO_NA           = "<NA-repo>"
const VERSION_NA        = "<NA-version>"

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Extracting list of packages
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# (Vector{PackageInfo}, Bool, Vector{String}) → Vector{String}
# Returns the list of package repos or package names depending on `showName`,
# sorted from the most to least starred packages,
# excluding packages with repos listed in `excludedRepos`
getSortedPkgs(
    pkgs :: Vector, showName :: Bool = false;
    excludedRepos :: Vector{String} = String[], includeVersion :: Bool = false
) :: Vector{String} = getSortedPkgsInfo(
    pkgs, showName ?
        (includeVersion ? getNameAndVersion : getName) :
        (includeVersion ? getRepoAndVersion : getRepo);
    excludedRepos=excludedRepos
)

# (Vector{PackageInfo}, PackageInfo→Any, Vector{String}) → Vector{Any}
# Returns the list of package information
# sorted from the most to least starred packages,
# excluding packages with repos listed in `excludedRepos`
getSortedPkgsInfo(
    pkgs :: Vector, getInfo :: Function;
    excludedRepos :: Vector{String} = String[]
) :: Vector = begin
    if !isempty(excludedRepos)
        @status "Cleaning packages..."
        pkgs = filter(pkg -> !in(getRepo(pkg), excludedRepos), pkgs)
    end
    @status "Sorting packages..."
    sorted = sortStarred(pkgs)
    @status "Extracting information..."
    unique(map(getInfo, sorted))
end

# (Vector{PackageInfo}, String) → Vector{PackageInfo}
# Sorts packages from the most starred to least starred
sortStarred(pkgs :: Vector) = sort(pkgs, by=getStarCount, rev=true)

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Utilitites for working with packages info
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# Retrieves the list of packages from file `pkgsInfoFilePath`
# ASSUMPTION: the file has the expected JSON format
loadPkgsData(pkgsInfoFilePath :: String) :: Vector = begin
    data = JSON.parsefile(pkgsInfoFilePath)
    data[PKGS_KEY]
end

# (PackageInfo, String) → Any
# Retrieves the value of `key` from `pkgInfo` metadata
# or returns `default` value
# ASSUMPTION: `pkgInfo` contains metadata
getMetaDataValue(pkgInfo :: Dict, key :: String, default :: Any) :: Any = begin
    metaData = pkgInfo[METADATA_KEY]
    val = get(metaData, key, default)
    # in case the value of the element is `null` in JSON
    val !== nothing ? val : default
end

# (PackageInfo, String) → String
getName(pkgInfo :: Dict) :: String = get(pkgInfo, NAME_KEY, NAME_NA)

# (PackageInfo, String) → String
getRepo(pkgInfo :: Dict) :: String =
    getMetaDataValue(pkgInfo, REPO_KEY, REPO_NA)

# (PackageInfo, String) → Int
getStarCount(pkgInfo :: Dict) :: Int =
    getMetaDataValue(pkgInfo, STARTCOUNT_KEY, 0)

# (PackageInfo, String) → String
# List of versions is sorted from the latest version
getLatestVersion(pkgInfo :: Dict) :: String = begin
    versions = getMetaDataValue(pkgInfo, VERSIONS_KEY, [])
    isempty(versions) ? VERSION_NA : versions[1]
end

# (PackageInfo, String) → String
getNameAndVersion(pkgInfo :: Dict) :: String =
    "$(getName(pkgInfo)), $(getLatestVersion(pkgInfo))"

# (PackageInfo, String) → String
getRepoAndVersion(pkgInfo :: Dict) :: String =
    "$(getRepo(pkgInfo)), $(getLatestVersion(pkgInfo))"


text:
module TestPackage0

include("precompile_includer.jl")
end


text:
#=
Please see
https://computationalmindset.com/en/neural-networks/ordinary-differential-equation-solvers.html#ode1
for details
=#

#=Original problem
x′ + x = sin t + 3 cos 2t
x(0) = 0
=#

#=Esplicit form of equation
x′ = sin t  + 3 cos 2t − x
=#

#=Analytical solution
x(t) = 1/2 sin t - 1/2 cos t + 3/5 cos 2t + 6/5 sin 2t - 1/10 e^-t
=#

using DifferentialEquations
using Plots

ode_fn(x,p,t) = sin(t) + 3.0 * cos(2.0 * t) - x

an_sol(t) = (1.0/2.0) * sin(t) - (1.0/2.0) * cos(t) +
            (3.0/5.0) * cos(2.0*t) + (6.0/5.0) * sin(2.0*t) -
            (1.0/10.0) * exp(-t)

t_begin=0.0
t_end=10.0
tspan = (t_begin,t_end)
x_init=0.0

prob = ODEProblem(ode_fn, x_init, tspan)
num_sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

plot(num_sol.t, an_sol.(num_sol.t),
    linewidth=2, ls=:dash,
    title="ODE 1st order IVP solved by D.E. package",
    xaxis="t", yaxis="x",
    label="analytical",
    legend=true)
plot!(num_sol,
    linewidth=1,
    label="numerical")

max_stars_repo_name: ettoremessina/differential-equations

max_issues_repo_name: ettoremessina/differential-equations

max_forks_repo_name: ettoremessina/differential-equations


text:
@testset "Mass matrices and inverses" begin
    R = FiniteDifferences(20,0.2)
    @testset "Mass matrix" begin
        @test R'R == step(R)*I
    end
    @testset "Inverses" begin
        R⁻¹ = pinv(R)
        @test R⁻¹*R == I
        @test R*R⁻¹ == I

        cu = rand(size(R,2))
        cv = rand(size(R,2))
        cuv = [cu cv]

        u = R*cu
        v = R*cv
        uv = R*cuv

        @test R⁻¹*u === cu
        @test R⁻¹*v === cv
        @test R⁻¹*uv === cuv

        ut = u'
        # Fails with: ERROR: MethodError: no method matching axes(::UniformScaling{Float64}, ::Int64)
        # @test ut*R⁻¹' === ut.args[1]
    end
end

@testset "Scalar operators" begin
    B = FiniteDifferences(5,1.0)
    x = axes(B,1)
    x² = B'QuasiDiagonal(x.^2)*B
    @test x² isa Diagonal

    v = ones(5)
    @test x²*v == locs(B).^2

    y = Inclusion(0.0..5.0)
    @test_throws DimensionMismatch B'QuasiDiagonal(y.^2)*B

    @testset "Restricted bases" begin
        N = 10
        ρ = 1.0

        R = FiniteDifferences(N, ρ)
        r = axes(R, 1)

        sel = 3:6
        sel2 = 8:10

        R̃ = R[:, sel]
        R′ = R[:, sel2]

        x = QuasiDiagonal(r)

        A = R'x*R
        @test A isa Diagonal

        a = R̃'x*R
        @test a isa BandedMatrix
        @test a == Matrix(A)[sel,:]
        @test bandwidths(a) == (-2,2)

        a = R'x*R̃
        @test a isa BandedMatrix
        @test a == Matrix(A)[:,sel]
        @test bandwidths(a) == (2,-2)

        a = R̃'x*R̃
        @test a isa BandedMatrix
        @test isdiag(a)
        @test a == Matrix(A)[sel,sel]

        # Non-overlapping restrictions
        a = R̃'x*R′
        @test a isa BandedMatrix
        @test size(a) == (length(sel), length(sel2))
        @test iszero(a)
        @test bandwidths(a) == (3,-5)
    end
end


text:
function solve!(solver::AbstractNonlinearSolver; verbose=false)

    @unpack P, prob, linalg = solver
    init_values!(P, prob)

    for n = 2:length(prob.pviews)
        @unpack Δt = @inbounds prob.pviews[n-1]
        prepare_step!(prob, n, P)
        perform_step!(solver, Δt[], n)

        if verbose
            @printf "n: %d, " n
            @printf "iter: %d, " solver.iter
            @printf "norm_r: %.5e, " solver.norm_r
            @printf "dP_rel: %.5e\n" solver.dP_rel
        end

        accept_step!(prob, n, P)
        solver.success || break
    end
    solver.success || error("Newton solver did not converge")

    return solver
end

function init_values!(P, prob::NonlinearProblem)

    @unpack Pi = @inbounds prob.pviews[1]
    @unpack Vwprev, Voprev, Vwi, Voi = prob.cache

    # Initialization at zero time step
    update_cache!(prob, 1)
    copyto!(P, Pi)
    copyto!(Vwprev, Vwi)
    copyto!(Voprev, Voi)

    return P, prob
end

function prepare_step!(prob::NonlinearProblem, n, P)

    params = @inbounds prob.pviews[n]
    @unpack Pmin, Pmax, Qliq, Qinj, Qliq_h, Qinj_h, λ = params
    @unpack Qsum = prob.cache

    update_cache!(prob, n)
    @turbo for i = 1:length(P)
        # Correction of rates taking into account the upper and lower bounds of pressure
        Qliq[i] = (P[i] ≥ Pmin[i]) * Qliq_h[i]
        Qinj[i] = (P[i] ≤ Pmax[i]) * Qinj_h[i]
        Qsum[i] = Qliq[i] - λ[i] * Qinj[i]
    end

    return prob
end

function accept_step!(prob::NonlinearProblem, n, P)

    params = @inbounds prob.pviews[n]
    @unpack jac_next, Pcalc, Δt = params
    @unpack Vwprev, Voprev, Vw, Vo, cwf, cof = prob.cache

    # Saving data at the current time step
    copyto!(Pcalc, P)
    copyto!(Vwprev, Vw)
    copyto!(Voprev, Vo)

    # The diagonal of the jacobian with respect to reservoir pressure on the previous time step
    @turbo for i = 1:length(P)
        jac_next[i] = -(Vw[i] * cwf[i] + Vo[i] * cof[i]) / Δt[]
    end

    return prob
end

function update_cache!(prob::NonlinearProblem{T}, n) where {T}

    params = @inbounds prob.pviews[n]

    # Initial porous volumes of fluids
    if params.Vupd[]
        @unpack Vpi, Swi, Bwi, Boi = params
        @unpack Vwi, Voi = prob.cache
        @turbo for i = 1:length(Vpi)
            Vwi[i] = Swi[i] * Vpi[i] / Bwi[i]
            Voi[i] = (one(T) - Swi[i]) * Vpi[i] / Boi[i]
        end
    end

    # Elements of the jacobian that do not depend on reservoir pressure
    if params.Tupd[]
        @unpack Tconn = params
        @unpack CTC, Cbuf = prob.cache
        # copyto!(Cbuf, prob.C)
        # lmul!(Diagonal(Tconn), Cbuf)
        # mul!(CTC, prob.C', Cbuf)
        # FIXED: For a sparse matrix 'C', the construction below is faster
        CTC .= prob.C' * Diagonal(Tconn) * prob.C
    end

    # Total compressibility of the reservoir-fluid system
    if params.cupd[]
        @unpack cw, co, cf = params
        @unpack cwf, cof = prob.cache
        @turbo for i = 1:length(cf)
            cwf[i] = cw[i] + cf[i]
            cof[i] = co[i] + cf[i]
        end
    end

    return prob
end


text:
# This file created the plots in the overleaf file. I only changedline 38, n_epochs.
using Flux, DiffEqFlux, OrdinaryDiffEq, DiffEqParamEstim, Plots, Optim, Dates
####################################################### Observation ###############################################################
# Start conditions for the two species in the system
u0 = Float32[2.; 0.]
# Number of evaluations of the neural ODE. It relates to the numbers of layers of the neural net (depth of network).
datasize = 30
# Time span in which of evaluation will be and actual timepoints of evaluations
tspan = (0.0f0, 1.5f0)
t = range(tspan[1], tspan[2], length = datasize)
# The true ODE (with the true parameters) which the neural net should learn
function trueODEfunc(du, u, p, t)
  true_A = [-0.1 2.0; -2.0 -0.1]
  du .= ((u.^3)'true_A)'
end
# Construction of the ODEProblem and solving the ODEProblem with Tsit5 solver
prob = ODEProblem(trueODEfunc, u0, tspan)
ode_data = Array(solve(prob,Tsit5(),saveat=t))
################################################### Building a neural ODE ##########################################################
# Derivative is modeled by a neural net. Chain concatinates the functions ode function and two dense layers.
dudt = Chain(x -> x.^3,
       Dense(2,50,tanh),
       Dense(50,2))
# Parameters of the model which are to be learnt. They are: W1 (2x50), b1 (50), W2 (50x2), b2 (2)
ps = Flux.params(dudt)
# Getting loss function from two stage collocation function
function node_two_stage_function(model, x, tspan, saveat, ode_data,
            args...; kwargs...)
  dudt_(du,u,p,t) = du .= model(u)
  prob_fly = ODEProblem(dudt_,x,tspan)
  two_stage_method(prob_fly, saveat, ode_data)
end
loss_n_ode = node_two_stage_function(dudt, u0, tspan, t, ode_data, Tsit5(), reltol=1e-7, abstol=1e-9)
#  loss function
two_stage_loss_fct()=loss_n_ode.cost_function(ps)
# Defining anonymous function for the neural ODE with the model. in: u0, out: solution with current params.
n_ode = x->neural_ode(dudt, x, tspan, Tsit5(), saveat=t, reltol=1e-7, abstol=1e-9)
n_epochs = 5000 # vary between [500, 1000, 2500, 5000]
data1 = Iterators.repeated((), n_epochs)
opt1 = ADAM(0.1)
# Callback function to observe two stage training.
mutable struct saver
    losses::Array{Float64,1}
    times::Array{Dates.Time,1}
    count_epochs::Int128
end
function saver(n_epochs)
    losses = zeros(n_epochs)
    times = fill(Dates.Time(Dates.now()),n_epochs)
    count_epochs = 0
    return saver(losses,times,count_epochs)
end
function update_saver(saver, loss_i, time_i)
    epoch_i = saver.count_epochs
    saver.losses[epoch_i] = loss_i
    saver.times[epoch_i] = time_i
end

sa = saver(n_epochs)

cb1 = function ()
    sa.count_epochs = sa.count_epochs +  1
    update_saver(sa, Tracker.data(two_stage_loss_fct()),Dates.Time(Dates.now()))
    #println("\"",Tracker.data(two_stage_loss_fct()),"\" \"",Dates.Time(Dates.now()),"\";")
end
#two stage training call
a = @time Flux.train!(two_stage_loss_fct, ps, data1, opt1, cb = cb1)
# Call n_ode to get first prediction and to show startpoint for training.

pred = n_ode(u0)
print(sum(abs2,ode_data .- pred))

scatter(t, ode_data[1,:], label="data")
scatter!(t, Flux.data(pred[1,:]), label="prediction")
scatter!(t, ode_data[2,:], label="data")
scatter!(t, Flux.data(pred[2,:]), label="prediction")

#savefig("sogood.png")
#print("hi")
#create outside references and enclose them if that's what you need

test_u0 = [2.5 , 1.]
pred_t = n_ode(test_u0)

prob_t = ODEProblem(trueODEfunc, test_u0, tspan)
ode_data_t = Array(solve(prob_t,Tsit5(),saveat=t))

scatter(t, ode_data_t[1,:], label="data")
plot!(t, Flux.data(pred_t[1,:]), label="prediction")
scatter!(t, ode_data_t[2,:], label="data")
plot!(t, Flux.data(pred_t[2,:]), label="prediction")


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ eb3b3be0-67d1-11eb-144b-2f9581950503
md"""
#### Math can be easy!
#### When you come on someting new
#### $~~~~$just think in simple steps.
"""

# ╔═╡ 66294030-67d3-11eb-25b1-8d25703546db
md"""
#### sqrt()
##### ??? What is this?
####  Pluto 🛡 to the rescue!
##### Doing Things have a name.$~~$ This name is $~~ sqrt$
##### The $~~$()$~~$ is a place for something to worked on.
##### Doing Things can teach a new name to do the same work.
##### We may like one,$~$ but both will do the same job.
"""

# ╔═╡ f4835c0a-67db-11eb-134c-85092b5a1ced
md"""
##### Looking for a short
##### cut across a square 🔲?
##### Walking on the edges is longer.
##### But corner to corner is shortest.
##### We call the shortest, the root 🥕.
"""

# ╔═╡ 1976e31c-67df-11eb-2d4a-1d8e612a99a9

md"""
##### sqrt() is a Doing Thing
##### and is called a function.
##### It can tells how far it
##### is across a 1 by 1 🔲
##### Instead of going on the sides,
##### place 1 + 1 inside of ()
"""

# ╔═╡ fca3b432-86f8-11eb-3086-b3adfccd0cda
sqrt(1+1)

# ╔═╡ 84554b8a-67ce-11eb-2e10-7550cbde4485
md"""
##### Hard to think that $~$sqrt$~$ means square root?
##### I bet you will never forget
## 🔲🥕
##### Let's use $~$=$~$ to make a new Doing Thing.
"""

# ╔═╡ f8c8c828-67d0-11eb-16b6-f384a4a4a5df
🔲🥕 = sqrt

# ╔═╡ 279bdd4e-67e2-11eb-1adb-b98afc7e812e
md"""
##### Yes, it works!
##### Tomorrow, will you remember?
### 🔲🥕  $~~~~or~~~~~$   sqrt
"""

# ╔═╡ 1c3a7e52-67d1-11eb-0c5a-a5b36b4fc989
🔲🥕(1+1)

# ╔═╡ 7dab943e-86f9-11eb-284d-4fb84d4c9b63
md"""
##### Wow!$~~$ The long way is 2.
##### But the short way is not even 1 and 1/2.
"""

# ╔═╡ Cell order:
# ╠═eb3b3be0-67d1-11eb-144b-2f9581950503
# ╠═66294030-67d3-11eb-25b1-8d25703546db
# ╠═f4835c0a-67db-11eb-134c-85092b5a1ced
# ╠═1976e31c-67df-11eb-2d4a-1d8e612a99a9
# ╠═fca3b432-86f8-11eb-3086-b3adfccd0cda
# ╠═84554b8a-67ce-11eb-2e10-7550cbde4485
# ╠═f8c8c828-67d0-11eb-16b6-f384a4a4a5df
# ╠═279bdd4e-67e2-11eb-1adb-b98afc7e812e
# ╠═1c3a7e52-67d1-11eb-0c5a-a5b36b4fc989
# ╠═7dab943e-86f9-11eb-284d-4fb84d4c9b63


text:
type fileHeader
	jobid  :: Int32
	linnum :: Int32
	renum  :: Int32
	ntrpe  :: Int16
	natrpe :: Int16
	dt     :: Int16
	dtfr   :: Int16
	ns     :: Int16
	nsfr   :: Int16
	fmtc   :: Int16
	expf   :: Int16
	trsc   :: Int16
	vsumc  :: Int16
	sfs    :: Int16
	sfe    :: Int16
	slen   :: Int16
	styp   :: Int16
	tnumsc :: Int16
	stalens:: Int16
	stalene:: Int16
	tyta   :: Int16
	corr   :: Int16
	rgc    :: Int16
	arm    :: Int16
	unit   :: Int16
	pol    :: Int16
	vpol   :: Int16
	fvn    :: Int16
	fltf   :: Int16
	netfh  :: Int16
end

fileHeader_count = Dict{AbstractString, Int32}()
fileHeader_count["jobid"  ] = 3200
fileHeader_count["linnum" ] = 3204
fileHeader_count["renum"  ] = 3208
fileHeader_count["ntrpe"  ] = 3212
fileHeader_count["natrpe" ] = 3214
fileHeader_count["dt"     ] = 3216
fileHeader_count["dtfr"   ] = 3218
fileHeader_count["ns"     ] = 3220
fileHeader_count["nsfr"   ] = 3222
fileHeader_count["fmtc"   ] = 3224
fileHeader_count["expf"   ] = 3226
fileHeader_count["trsc"   ] = 3228
fileHeader_count["vsumc"  ] = 3230
fileHeader_count["sfs"    ] = 3232
fileHeader_count["sfe"    ] = 3234
fileHeader_count["slen"   ] = 3236
fileHeader_count["styp"   ] = 3238
fileHeader_count["tnumsc" ] = 3240
fileHeader_count["stalens"] = 3242
fileHeader_count["stalene"] = 3244
fileHeader_count["tyta"   ] = 3246
fileHeader_count["corr"   ] = 3248
fileHeader_count["rgc"    ] = 3250
fileHeader_count["arm"    ] = 3252
fileHeader_count["unit"   ] = 3254
fileHeader_count["pol"    ] = 3256
fileHeader_count["vpol"   ] = 3258
fileHeader_count["fvn"    ] = 3500
fileHeader_count["fltf"   ] = 3502
fileHeader_count["netfh"  ] = 3504

function InitFileHeader()
	fh = fileHeader(0,0,0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,0,0,)
	return fh
end

function GrabFileHeader(stream)

	position = fileHeader_count["jobid"]
	seek(stream,position)
	fh = InitFileHeader()
	fh.jobid   = read(stream,typeof(fh.jobid))
	fh.linnum  = read(stream,typeof(fh.linnum))
	fh.renum   = read(stream,typeof(fh.renum))
	fh.ntrpe   = read(stream,typeof(fh.ntrpe))
	fh.natrpe  = read(stream,typeof(fh.natrpe))
	fh.dt      = read(stream,typeof(fh.dt))
	fh.dtfr    = read(stream,typeof(fh.dtfr))
	fh.ns      = read(stream,typeof(fh.ns))
	fh.nsfr    = read(stream,typeof(fh.nsfr))
	fh.fmtc    = read(stream,typeof(fh.fmtc))
	fh.expf    = read(stream,typeof(fh.expf))
	fh.trsc    = read(stream,typeof(fh.trsc))
	fh.vsumc   = read(stream,typeof(fh.vsumc))
	fh.sfs     = read(stream,typeof(fh.sfs))
	fh.sfe     = read(stream,typeof(fh.sfe))
	fh.slen    = read(stream,typeof(fh.slen))
	fh.styp    = read(stream,typeof(fh.styp))
	fh.tnumsc  = read(stream,typeof(fh.tnumsc))
	fh.stalens = read(stream,typeof(fh.stalens))
	fh.stalene = read(stream,typeof(fh.stalene))
	fh.tyta    = read(stream,typeof(fh.tyta))
	fh.corr    = read(stream,typeof(fh.corr))
	fh.rgc     = read(stream,typeof(fh.rgc))
	fh.arm     = read(stream,typeof(fh.arm))
	fh.unit    = read(stream,typeof(fh.unit))
	fh.pol     = read(stream,typeof(fh.pol))
	fh.vpol    = read(stream,typeof(fh.vpol))
	fh.fvn     = read(stream,typeof(fh.fvn))
	fh.fltf    = read(stream,typeof(fh.fltf))
	fh.netfh   = read(stream,typeof(fh.netfh))

	return fh
end

function PutFileHeader(stream,fh)
	position = 3600
	seek(stream,position)
	write(stream,fh.jobid)
	write(stream,fh.linnum)
	write(stream,fh.renum)
	write(stream,fh.ntrpe)
	write(stream,fh.natrpe)
	write(stream,fh.dt)
	write(stream,fh.dtfr)
	write(stream,fh.ns)
	write(stream,fh.nsfr)
	write(stream,fh.fmtc)
	write(stream,fh.expf)
	write(stream,fh.trsc)
	write(stream,fh.vsumc)
	write(stream,fh.sfs)
	write(stream,fh.sfe)
	write(stream,fh.slen)
	write(stream,fh.styp)
	write(stream,fh.tnumsc)
	write(stream,fh.stalens)
	write(stream,fh.stalene)
	write(stream,fh.tyta)
	write(stream,fh.corr)
	write(stream,fh.rgc)
	write(stream,fh.arm)
	write(stream,fh.unit)
	write(stream,fh.pol)
	write(stream,fh.vpol)
	write(stream,fh.fvn)
	write(stream,fh.fltf)
	write(stream,fh.netfh)

end

type SegyHeader
	tracl::Int32
	tracr::Int32
	fldr::Int32
	tracf::Int32
	ep::Int32
	cdp::Int32
	cdpt::Int32
	trid::Int16
	nva::Int16
	nhs::Int16
	duse::Int16
	offset::Int32
	gelev::Int32
	selev::Int32
	sdepth::Int32
	gdel::Int32
	sdel::Int32
	swdep::Int32
	gwdep::Int32
	scalel::Int16
	scalco::Int16
	sx::Int32
	sy::Int32
	gx::Int32
	gy::Int32
	counit::Int16
	wevel::Int16
	swevel::Int16
	sut::Int16
	gut::Int16
	sstat::Int16
	gstat::Int16
	tstat::Int16
	laga::Int16
	lagb::Int16
	delrt::Int16
	muts::Int16
	mute::Int16
	ns::Int16
	dt::Int16
	gain::Int16
	igc::Int16
	igi::Int16
	corr::Int16
	sfs::Int16
	sfe::Int16
	slen::Int16
	styp::Int16
	stas::Int16
	stae::Int16
	tatyp::Int16
	afilf::Int16
	afils::Int16
	nofilf::Int16
	nofils::Int16
	lcf::Int16
	hcf::Int16
	lcs::Int16
	hcs::Int16
	year::Int16
	day::Int16
	hour::Int16
	minute::Int16
	sec::Int16
	timbas::Int16
	trwf::Int16
	grnors::Int16
	grnofr::Int16
	grnlof::Int16
	gaps::Int16
	otrav::Int16
	d1::Float32
	f1::Float32
	d2::Float32
	f2::Float32
	ungpow::Float32
	unscale::Float32
	ntr::Int32
	mark::Int16
	unass::Int16
end

segy_count = Dict{AbstractString,Int32}()
segy_count["tracl"]  = 0
segy_count["tracr"]  = 4
segy_count["fldr"]   = 8
segy_count["tracf"]  = 12
segy_count["ep"]     = 16
segy_count["cdp"]    = 20
segy_count["cdpt"]   = 24
segy_count["trid"]   = 28
segy_count["nva"]    = 30
segy_count["nhs"]    = 32
segy_count["duse"]   = 34
segy_count["offset"] = 36
segy_count["gelev"]  = 40
segy_count["selev"]  = 44
segy_count["sdepth"] = 48
segy_count["gdel"]   = 52
segy_count["sdel"]   = 56
segy_count["swdep"]  = 60
segy_count["gwdep"]  = 64
segy_count["scalel"] = 68
segy_count["scalco"] = 70
segy_count["sx"]     = 72
segy_count["sy"]     = 76
segy_count["gx"]     = 80
segy_count["gy"]     = 84
segy_count["counit"] = 88
segy_count["wevel"]  = 90
segy_count["swevel"] = 92
segy_count["sut"]    = 94
segy_count["gut"]    = 96
segy_count["sstat"]  = 98
segy_count["gstat"]  = 100
segy_count["tstat"]  = 102
segy_count["laga"]   = 104
segy_count["lagb"]   = 106
segy_count["delrt"]  = 108
segy_count["muts"]   = 110
segy_count["mute"]   = 112
segy_count["ns"]     = 114
segy_count["dt"]     = 116
segy_count["gain"]   = 118
segy_count["igc"]    = 120
segy_count["igi"]    = 122
segy_count["corr"]   = 124
segy_count["sfs"]    = 126
segy_count["sfe"]    = 128
segy_count["slen"]   = 130
segy_count["styp"]   = 132
segy_count["stas"]   = 134
segy_count["stae"]   = 136
segy_count["tatyp"]  = 138
segy_count["afilf"]  = 140
segy_count["afils"]  = 142
segy_count["nofilf"] = 144
segy_count["nofils"] = 146
segy_count["lcf"]    = 148
segy_count["hcf"]    = 150
segy_count["lcs"]    = 152
segy_count["hcs"]    = 154
segy_count["year"]   = 156
segy_count["day"]    = 158
segy_count["hour"]   = 160
segy_count["minute"] = 162
segy_count["sec"]    = 164
segy_count["timbas"] = 166
segy_count["trwf"]   = 168
segy_count["grnors"] = 170
segy_count["grnofr"] = 172
segy_count["grnlof"] = 174
segy_count["gaps"]   = 176
segy_count["otrav"]  = 178
segy_count["d1"]     = 180
segy_count["f1"]     = 184
segy_count["d2"]     = 188
segy_count["f2"]     = 192
segy_count["ungpow"] = 196
segy_count["unscale"]= 200
segy_count["ntr"]    = 204
segy_count["mark"]   = 208
segy_count["unass"]  = 210
segy_count["trace"]  = 240

function InitSegyHeader()
	h = SegyHeader(0,0,0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,0,0)
	return h
end

function GrabSegyHeader(stream,swap_bytes,nt,file_header_size,j)
	position = file_header_size + (240+nt*4)*(j-1)
	seek(stream,position)
	h = InitSegyHeader()
	if (swap_bytes == false)
		h.tracl  = read(stream,typeof(h.tracl))
		h.tracr  = read(stream,typeof(h.tracr))
		h.fldr   = read(stream,typeof(h.fldr))
		h.tracf  = read(stream,typeof(h.tracf))
		h.ep     = read(stream,typeof(h.ep))
		h.cdp    = read(stream,typeof(h.cdp))
		h.cdpt   = read(stream,typeof(h.cdpt))
		h.trid   = read(stream,typeof(h.trid))
		h.nva    = read(stream,typeof(h.nva))
		h.nhs    = read(stream,typeof(h.nhs))
		h.duse   = read(stream,typeof(h.duse))
		h.offset = read(stream,typeof(h.offset))
		h.gelev  = read(stream,typeof(h.gelev))
		h.selev  = read(stream,typeof(h.selev))
		h.sdepth = read(stream,typeof(h.sdepth))
		h.gdel   = read(stream,typeof(h.gdel))
		h.sdel   = read(stream,typeof(h.sdel))
		h.swdep  = read(stream,typeof(h.swdep))
		h.gwdep  = read(stream,typeof(h.gwdep))
		h.scalel = read(stream,typeof(h.scalel))
		h.scalco = read(stream,typeof(h.scalco))
		h.sx     = read(stream,typeof(h.sx))
		h.sy     = read(stream,typeof(h.sy))
		h.gx     = read(stream,typeof(h.gx))
		h.gy     = read(stream,typeof(h.gy))
		h.counit = read(stream,typeof(h.counit))
		h.wevel  = read(stream,typeof(h.wevel))
		h.swevel = read(stream,typeof(h.swevel))
		h.sut    = read(stream,typeof(h.sut))
		h.gut    = read(stream,typeof(h.gut))
		h.sstat  = read(stream,typeof(h.sstat))
		h.gstat  = read(stream,typeof(h.gstat))
		h.tstat  = read(stream,typeof(h.tstat))
		h.laga   = read(stream,typeof(h.laga))
		h.lagb   = read(stream,typeof(h.lagb))
		h.delrt  = read(stream,typeof(h.delrt))
		h.muts   = read(stream,typeof(h.muts))
		h.mute   = read(stream,typeof(h.mute))
		h.ns     = read(stream,typeof(h.ns))
		h.dt     = read(stream,typeof(h.dt))
		h.gain   = read(stream,typeof(h.gain))
		h.igc    = read(stream,typeof(h.igc))
		h.igi    = read(stream,typeof(h.igi))
		h.corr   = read(stream,typeof(h.corr))
		h.sfs    = read(stream,typeof(h.sfs))
		h.sfe    = read(stream,typeof(h.sfe))
		h.slen   = read(stream,typeof(h.slen))
		h.styp   = read(stream,typeof(h.styp))
		h.stas   = read(stream,typeof(h.stas))
		h.stae   = read(stream,typeof(h.stae))
		h.tatyp  = read(stream,typeof(h.tatyp))
		h.afilf  = read(stream,typeof(h.afilf))
		h.afils  = read(stream,typeof(h.afils))
		h.nofilf = read(stream,typeof(h.nofilf))
		h.nofils = read(stream,typeof(h.nofils))
		h.lcf    = read(stream,typeof(h.lcf))
		h.hcf    = read(stream,typeof(h.hcf))
		h.lcs    = read(stream,typeof(h.lcs))
		h.hcs    = read(stream,typeof(h.hcs))
		h.year   = read(stream,typeof(h.year))
		h.day    = read(stream,typeof(h.day))
		h.hour   = read(stream,typeof(h.hour))
		h.minute = read(stream,typeof(h.minute))
		h.sec    = read(stream,typeof(h.sec))
		h.timbas = read(stream,typeof(h.timbas))
		h.trwf   = read(stream,typeof(h.trwf))
		h.grnors = read(stream,typeof(h.grnors))
		h.grnofr = read(stream,typeof(h.grnofr))
		h.grnlof = read(stream,typeof(h.grnlof))
		h.gaps   = read(stream,typeof(h.gaps))
		h.otrav  = read(stream,typeof(h.otrav))
		h.d1     = read(stream,typeof(h.d1))
		h.f1     = read(stream,typeof(h.f1))
		h.d2     = read(stream,typeof(h.d2))
		h.f2     = read(stream,typeof(h.f2))
		h.ungpow = read(stream,typeof(h.ungpow))
		h.unscale= read(stream,typeof(h.unscale))
		h.ntr    = read(stream,typeof(h.ntr))
		h.mark   = read(stream,typeof(h.mark))
		h.unass  = read(stream,typeof(h.unass))
	else
		h.tracl  = bswap(read(stream,typeof(h.tracl)))
		h.tracr  = bswap(read(stream,typeof(h.tracr)))
		h.fldr   = bswap(read(stream,typeof(h.fldr)))
		h.tracf  = bswap(read(stream,typeof(h.tracf)))
		h.ep     = bswap(read(stream,typeof(h.ep)))
		h.cdp    = bswap(read(stream,typeof(h.cdp)))
		h.cdpt   = bswap(read(stream,typeof(h.cdpt)))
		h.trid   = bswap(read(stream,typeof(h.trid)))
		h.nva    = bswap(read(stream,typeof(h.nva)))
		h.nhs    = bswap(read(stream,typeof(h.nhs)))
		h.duse   = bswap(read(stream,typeof(h.duse)))
		h.offset = bswap(read(stream,typeof(h.offset)))
		h.gelev  = bswap(read(stream,typeof(h.gelev)))
		h.selev  = bswap(read(stream,typeof(h.selev)))
		h.sdepth = bswap(read(stream,typeof(h.sdepth)))
		h.gdel   = bswap(read(stream,typeof(h.gdel)))
		h.sdel   = bswap(read(stream,typeof(h.sdel)))
		h.swdep  = bswap(read(stream,typeof(h.swdep)))
		h.gwdep  = bswap(read(stream,typeof(h.gwdep)))
		h.scalel = bswap(read(stream,typeof(h.scalel)))
		h.scalco = bswap(read(stream,typeof(h.scalco)))
		h.sx     = bswap(read(stream,typeof(h.sx)))
		h.sy     = bswap(read(stream,typeof(h.sy)))
		h.gx     = bswap(read(stream,typeof(h.gx)))
		h.gy     = bswap(read(stream,typeof(h.gy)))
		h.counit = bswap(read(stream,typeof(h.counit)))
		h.wevel  = bswap(read(stream,typeof(h.wevel)))
		h.swevel = bswap(read(stream,typeof(h.swevel)))
		h.sut    = bswap(read(stream,typeof(h.sut)))
		h.gut    = bswap(read(stream,typeof(h.gut)))
		h.sstat  = bswap(read(stream,typeof(h.sstat)))
		h.gstat  = bswap(read(stream,typeof(h.gstat)))
		h.tstat  = bswap(read(stream,typeof(h.tstat)))
		h.laga   = bswap(read(stream,typeof(h.laga)))
		h.lagb   = bswap(read(stream,typeof(h.lagb)))
		h.delrt  = bswap(read(stream,typeof(h.delrt)))
		h.muts   = bswap(read(stream,typeof(h.muts)))
		h.mute   = bswap(read(stream,typeof(h.mute)))
		h.ns     = bswap(read(stream,typeof(h.ns)))
		h.dt     = bswap(read(stream,typeof(h.dt)))
		h.gain   = bswap(read(stream,typeof(h.gain)))
		h.igc    = bswap(read(stream,typeof(h.igc)))
		h.igi    = bswap(read(stream,typeof(h.igi)))
		h.corr   = bswap(read(stream,typeof(h.corr)))
		h.sfs    = bswap(read(stream,typeof(h.sfs)))
		h.sfe    = bswap(read(stream,typeof(h.sfe)))
		h.slen   = bswap(read(stream,typeof(h.slen)))
		h.styp   = bswap(read(stream,typeof(h.styp)))
		h.stas   = bswap(read(stream,typeof(h.stas)))
		h.stae   = bswap(read(stream,typeof(h.stae)))
		h.tatyp  = bswap(read(stream,typeof(h.tatyp)))
		h.afilf  = bswap(read(stream,typeof(h.afilf)))
		h.afils  = bswap(read(stream,typeof(h.afils)))
		h.nofilf = bswap(read(stream,typeof(h.nofilf)))
		h.nofils = bswap(read(stream,typeof(h.nofils)))
		h.lcf    = bswap(read(stream,typeof(h.lcf)))
		h.hcf    = bswap(read(stream,typeof(h.hcf)))
		h.lcs    = bswap(read(stream,typeof(h.lcs)))
		h.hcs    = bswap(read(stream,typeof(h.hcs)))
		h.year   = bswap(read(stream,typeof(h.year)))
		h.day    = bswap(read(stream,typeof(h.day)))
		h.hour   = bswap(read(stream,typeof(h.hour)))
		h.minute = bswap(read(stream,typeof(h.minute)))
		h.sec    = bswap(read(stream,typeof(h.sec)))
		h.timbas = bswap(read(stream,typeof(h.timbas)))
		h.trwf   = bswap(read(stream,typeof(h.trwf)))
		h.grnors = bswap(read(stream,typeof(h.grnors)))
		h.grnofr = bswap(read(stream,typeof(h.grnofr)))
		h.grnlof = bswap(read(stream,typeof(h.grnlof)))
		h.gaps   = bswap(read(stream,typeof(h.gaps)))
		h.otrav  = bswap(read(stream,typeof(h.otrav)))
		h.d1     = bswap(read(stream,typeof(h.d1)))
		h.f1     = bswap(read(stream,typeof(h.f1)))
		h.d2     = bswap(read(stream,typeof(h.d2)))
		h.f2     = bswap(read(stream,typeof(h.f2)))
		h.ungpow = bswap(read(stream,typeof(h.ungpow)))
		h.unscale= bswap(read(stream,typeof(h.unscale)))
		h.ntr    = bswap(read(stream,typeof(h.ntr)))
		h.mark   = bswap(read(stream,typeof(h.mark)))
		h.unass  = bswap(read(stream,typeof(h.unass)))
	end

	return h
end

function PutSegyHeader(stream,h,nt,file_header_size,j)
	position = file_header_size + (240+nt*4)*(j-1)
	seek(stream,position)
	write(stream,h.tracl)
	write(stream,h.tracr)
	write(stream,h.fldr)
	write(stream,h.tracf)
	write(stream,h.ep)
	write(stream,h.cdp)
	write(stream,h.cdpt)
	write(stream,h.trid)
	write(stream,h.nva)
	write(stream,h.nhs)
	write(stream,h.duse)
	write(stream,h.offset)
	write(stream,h.gelev)
	write(stream,h.selev)
	write(stream,h.sdepth)
	write(stream,h.gdel)
	write(stream,h.sdel)
	write(stream,h.swdep)
	write(stream,h.gwdep)
	write(stream,h.scalel)
	write(stream,h.scalco)
	write(stream,h.sx)
	write(stream,h.sy)
	write(stream,h.gx)
	write(stream,h.gy)
	write(stream,h.counit)
	write(stream,h.wevel)
	write(stream,h.swevel)
	write(stream,h.sut)
	write(stream,h.gut)
	write(stream,h.sstat)
	write(stream,h.gstat)
	write(stream,h.tstat)
	write(stream,h.laga)
	write(stream,h.lagb)
	write(stream,h.delrt)
	write(stream,h.muts)
	write(stream,h.mute)
	write(stream,h.ns)
	write(stream,h.dt)
	write(stream,h.gain)
	write(stream,h.igc)
	write(stream,h.igi)
	write(stream,h.corr)
	write(stream,h.sfs)
	write(stream,h.sfe)
	write(stream,h.slen)
	write(stream,h.styp)
	write(stream,h.stas)
	write(stream,h.stae)
	write(stream,h.tatyp)
	write(stream,h.afilf)
	write(stream,h.afils)
	write(stream,h.nofilf)
	write(stream,h.nofils)
	write(stream,h.lcf)
	write(stream,h.hcf)
	write(stream,h.lcs)
	write(stream,h.hcs)
	write(stream,h.year)
	write(stream,h.day)
	write(stream,h.hour)
	write(stream,h.minute)
	write(stream,h.sec)
	write(stream,h.timbas)
	write(stream,h.trwf)
	write(stream,h.grnors)
	write(stream,h.grnofr)
	write(stream,h.grnlof)
	write(stream,h.gaps)
	write(stream,h.otrav)
	write(stream,h.d1)
	write(stream,h.f1)
	write(stream,h.d2)
	write(stream,h.f2)
	write(stream,h.ungpow)
	write(stream,h.unscale)
	write(stream,h.ntr)
	write(stream,h.mark)
	write(stream,h.unass)
end

function MapHeaders(h_in,j,map_type)

	if map_type=="SegyToSeis"
		#scalco = abs(convert(Float32,h_in[1].scalco)) < 10 ? convert(Float32,h_in[1].scalco) : sign(convert(Float32,h_in[1].scalco))*log10(abs(convert(Float32,h_in[1].scalco)))
		scalco = sign(convert(Float32,h_in[1].scalco))*log10(abs(convert(Float32,h_in[1].scalco)))

		#scalel = abs(convert(Float32,h_in[1].scalel)) < 10 ? convert(Float32,h_in[1].scalel) : sign(convert(Float32,h_in[1].scalco))*log10(abs(convert(Float32,h_in[1].scalco)))
		scalel = sign(convert(Float32,h_in[1].scalco))*log10(abs(convert(Float32,h_in[1].scalco)))
		h_out = InitSeisHeader()
		h_out.tracenum = convert(typeof(h_out.tracenum),j)
		h_out.o1 = convert(typeof(h_out.o1),0)
		h_out.n1 = convert(typeof(h_out.n1),h_in[1].ns)
		h_out.d1 = convert(typeof(h_out.d1),h_in[1].dt/1000000)
		h_out.sx = scalco >= 0 ? convert(typeof(h_out.sx),h_in[1].sx)*10^scalco : convert(typeof(h_out.sx),h_in[1].sx)/(10^abs(scalco))
		h_out.sy = scalco >= 0 ? convert(typeof(h_out.sy),h_in[1].sy)*10^scalco : convert(typeof(h_out.sx),h_in[1].sy)/(10^abs(scalco))
		h_out.gx = scalco >= 0 ? convert(typeof(h_out.gx),h_in[1].gx)*10^scalco : convert(typeof(h_out.gx),h_in[1].gx)/(10^abs(scalco))
		h_out.gy = scalco >= 0 ? convert(typeof(h_out.gy),h_in[1].gy)*10^scalco : convert(typeof(h_out.gy),h_in[1].gy)/(10^abs(scalco))
		h_out.h = convert(typeof(h_out.h),h_in[1].offset)
		h_out.isx = convert(typeof(h_out.isx),h_in[1].ep)
		h_out.imx = convert(typeof(h_out.imx),h_in[1].cdp)
		h_out.selev = scalel >= 0 ? convert(typeof(h_out.selev),h_in[1].selev)*(10^scalel) : convert(typeof(h_out.selev),h_in[1].selev)/(10^abs(scalel))
		h_out.gelev = scalel >= 0 ? convert(typeof(h_out.gelev),h_in[1].gelev)*(10^scalel) : convert(typeof(h_out.gelev),h_in[1].gelev)/(10^abs(scalel))
		h_out.trid = convert(typeof(h_out.trid),h_in[1].trid)
	else
		h_out = InitSegyHeader()
		h_out.tracl = convert(typeof(h_out.tracl),j)
		h_out.tracr = convert(typeof(h_out.tracr),j)
		h_out.scalel = convert(typeof(h_out.scalel),-3)
		h_out.scalco = convert(typeof(h_out.scalco),-3)
		h_out.counit = convert(typeof(h_out.counit),1)
		h_out.gain = convert(typeof(h_out.gain),1)
		h_out.corr = convert(typeof(h_out.corr),1)
		h_out.ns = convert(typeof(h_out.ns),Int(h_in[1].n1))
		h_out.dt = convert(typeof(h_out.dt),round(Int,h_in[1].d1*1000000))
		h_out.sx = convert(typeof(h_out.sx),round(Int,h_in[1].sx*1000))
		h_out.sy = convert(typeof(h_out.sy),round(Int,h_in[1].sy*1000))
		h_out.gx = convert(typeof(h_out.gx),round(Int,h_in[1].gx*1000))
		h_out.gy = convert(typeof(h_out.gy),round(Int,h_in[1].gy*1000))
		h_out.offset = convert(typeof(h_out.offset),round(Int,h_in[1].h))
		h_out.ep = convert(typeof(h_out.ep),Int(h_in[1].isx))
		h_out.cdp = convert(typeof(h_out.cdp),Int(h_in[1].imx))
		h_out.selev = convert(typeof(h_out.selev),round(Int,h_in[1].selev*1000))
		h_out.gelev = convert(typeof(h_out.gelev),round(Int,h_in[1].gelev*1000))
		h_out.trid = convert(typeof(h_out.trid),h_in[1].trid)
	end

	return h_out
end

import Base.convert

primitive type IBMFloat32 32 end

ieeeOfPieces(fr::UInt32, exp::Int32, sgn::UInt32) = reinterpret(Float32, convert(UInt32,fr >>> 9) | convert(UInt32,exp << 23) | sgn) :: Float32
import Base.convert

function convert(::Type{Float32}, ibm::IBMFloat32)
  local fr::UInt32 = ntoh(reinterpret(UInt32, ibm))
  local sgn::UInt32 = fr & 0x80000000 # save sign
  fr <<= 1 # shift sign out
  local exp::Int32 = convert(Int32,fr >>> 25) # save exponent
  fr <<= 7 # shift exponent out

  if (fr == convert(UInt32,0))
    zero(Float32)
  else
    # normalize the signficand
    local norm::UInt32 = leading_zeros(fr)
    fr <<= norm
    exp = (exp << 2) - 130 - norm

    # exp <= 0 --> ieee(0,0,sgn)
    # exp >= 255 --> ieee(0,255,sgn)
    # else -> ieee(fr<<1, exp, sgn)
    local clexp::Int32 = exp & convert(Int32,0xFF)
    ieeeOfPieces(clexp == exp ? fr << 1 : convert(UInt32,0), clexp, sgn)
  end
end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
@testset "closest" begin
    x = 1:3:19
    @test closest(x, 3) == 2
    @test closest(x, 4) == 2
    @test closest(x, 5) == 2

    x = 19:-3:1
    @test closest(x, 3) == 6
    @test closest(x, 4) == 6
    @test closest(x, 5) == 6

    x = [4, 1, 7, 10, 13, 16, 19]
    @test closest(x, 3) == 1
    @test closest(x, 4) == 1
    @test closest(x, 5) == 1

    x = Axis([1.5, 2.0, 2.5, 3.0])
    @test x[closest(1.6)] == 1
    @test x[closest(1.9)] == 2
    @test x[closest(2.0)] == 2
    @test checkindex(Bool, x, closest(1.6))
end


text:
const WHITELIST = [' ','-'] # Characters which may be repeated

# This function operates in essentially constant time on real-ish words: 90 ns
# Most of this time (~70 ns) is spent initialising the Set.
function isisogram(s::AbstractString)
    seen = Set{Char}()
    for C in s
        c = lowercase(C) # Most words aren't isograms - this ensures we don't spend time lowercasing the whole string when we probably don't need to
        c ∈ seen && return false
        !(c in WHITELIST) && push!(seen,c)
    end
    true
end

# Lines below are interesting but require the StatsBase.jl and StaticArrays.jl packages

# # cmcaine patented technique
# # (ça ne marche pas pour les mots étrangers)
# #
# # ~5 times faster than above method for short real-ish words.
# # NB: all speedup comes from StaticArrays - a normal array makes it roughly the same speed as the function above.
# using StaticArrays
# function isisogram2(input)
#     counts = @MVector zeros(Int, 26)
#     @inbounds for c in codeunits(input)
#         if codepoint('a') <= c <= codepoint('z')
#             counts[c - codepoint('a')] += 1
#         elseif codepoint('A') <= c <= codepoint('Z')
#             counts[c - codepoint('A')] += 1
#         end
#         # Early return:
#         # This line makes the function operate in approx constant time for real-ish words, ~60 ns
#         # It is more performant without it for words < 50 characters.
#         # Longer words are unlikely to be pangrams anyway given that the alphabet has 26 characters.
# #       any(counts .> 1) && return false
#     end
#     all(counts .< 2)
# end
#
# # For testing
# import StatsBase
# using BenchmarkTools
# const ALPHABET = 'a':'z' |> collect
# const WEIGHTS = [ # from https://www.wikipedia.org/wiki/Letter_frequency
#  8.167e-2,
#  1.492e-2,
#  2.202e-2,
#  4.253e-2,
#  12.702e-2,
#  2.228e-2,
#  2.015e-2,
#  6.094e-2,
#  6.966e-2,
#  0.153e-2,
#  1.292e-2,
#  4.025e-2,
#  2.406e-2,
#  6.749e-2,
#  7.507e-2,
#  1.929e-2,
#  0.095e-2,
#  5.987e-2,
#  6.327e-2,
#  9.356e-2,
#  2.758e-2,
#  0.978e-2,
#  2.560e-2,
#  0.150e-2,
#  1.994e-2,
#  0.077e-2,
# ]
# randword(n) = [ALPHABET[StatsBase.sample(StatsBase.ProbabilityWeights(WEIGHTS))] for i in 1:n] |> join
# @benchmark isisogram(w) setup=(w=randword(10))
# @benchmark isisogram2(w) setup=(w=randword(10))


text:
using ConfigLH, Test

@testset "projectLH" begin
	println("-----")
	println("Testing projectLH")
    pName = "test";
    @test isdir(ConfigLH.project_dir(2018));
    proj = ConfigLH.Project(pName, pwd(), [pwd()]);
    show(proj);

    proj2 = ConfigLH.get_project(pName);
    @test proj2.name == pName;
    project_start(pName)
end

# ------------


text:
radiation_type = typeof(0.1W*m^-2)
snowdepth_type = typeof(0.1m)
airtemperature_type = typeof(0.1K)
relhumidity_type = typeof(0.001)
windspeed_type = typeof(0.1m*s^-1)
soiltemperature_type = typeof(0.1K)
soilwaterpotential_type = typeof(0.0kPa)
soilwatercontent_type = typeof(0.001)

abstract type AbstractMicroclimate end
abstract type AbstractLayeredMicroclimate{S,I,A,R} <: AbstractMicroclimate end
abstract type AbstractMicroclimControl <: AbstractMicroclimate end
abstract type AbstractMicroclimGrid{S,I,A,R} <: AbstractLayeredMicroclimate{S,I,A,R} end
abstract type AbstractMicroclimPoint{S,I,A,R} <: AbstractLayeredMicroclimate{S,I,A,R} end

@udefault_kw @units @limits struct MicroclimControl{RA,SN,AT,RH,WS,ST,WP,WC} <: AbstractMicroclimControl
    radiation::RA          | 200.0  | W*m^-2 | (500.0, 1200.0)
    snowdepth::SN          | 0.0    | m      | (0.0, 10.0)
    airtemperature::AT     | 290.0  | K      | (260.0, 370.0)
    relhumidity::RH        | 0.7    | _      | (0.0, 1.0)
    windspeed::WS          | 0.1    | m*s^-1 | (0.0, 20.0)
    soiltemperature::ST    | 290.0  | K      | (260.0, 370.0)
    soilwaterpotential::WP | -100.0 | kPa    | (-0.0, -10000.0)
    soilwatercontent::WC   | 0.3    | _      | (0.0, 1.0)
end

struct MicroclimGrid{S,I,A,R,RA,SN,AT,RH,WS,ST,WP,WC} <: AbstractMicroclimGrid{S,I,A,R}
    radiation::RA
    snowdepth::SN
    airtemperature::AT
    relhumidity::RH
    windspeed::WS
    soiltemperature::ST
    soilwaterpotential::WP
    soilwatercontent::WC
end

MicroclimGrid{S,I,A,R}(radiation::RA, snowdepth::SN, airtemperature::AT, relhumidity::RH,
                        windspeed::WS, soiltemperature::ST, soilwaterpotential::WP,
                        soilwatercontent::WC) where {S,I,A,R,RA,SN,AT,RH,WS,ST,WP,WC} =
    MicroclimGrid{S,I,A,R,RA,SN,AT,RH,WS,ST,WP,WC}(radiation, snowdepth, airtemperature,
        relhumidity, windspeed, soiltemperature, soilwaterpotential, soilwatercontent)

struct MicroclimPoint{S,I,A,R,RA,SN,AT,RH,WS,ST,WP,WC} <: AbstractMicroclimPoint{S,I,A,R}
    radiation::RA
    snowdepth::SN
    airtemperature::AT
    relhumidity::RH
    windspeed::WS
    soiltemperature::ST
    soilwaterpotential::WP
    soilwatercontent::WC
end

MicroclimPoint{S,I,A,R}(radiation::RA, snowdepth::SN, airtemperature::AT, relhumidity::RH,
                        windspeed::WS, soiltemperature::ST, soilwaterpotential::WP,
                        soilwatercontent::WC) where {S,I,A,R,RA,SN,AT,RH,WS,ST,WP,WC} =
    MicroclimPoint{S,I,A,R,RA,SN,AT,RH,WS,ST,WP,WC}(radiation, snowdepth, airtemperature,
        relhumidity, windspeed, soiltemperature, soilwaterpotential, soilwatercontent)

# From external data without units
MicroclimPoint{S,I,A,R}(radiation::Vector{<:AbstractFloat}, snowdepth::Vector{<:AbstractFloat},
                      airtemperature::Matrix{<:AbstractFloat}, relhumidity::Matrix{<:AbstractFloat},
                      windspeed::Matrix{<:AbstractFloat}, soiltemperature::Matrix{<:AbstractFloat},
                      soilwaterpotential::Matrix{<:AbstractFloat}, soilwatercontent::Matrix{<:AbstractFloat}
                     ) where {S,I,A,R} = begin
    rad = missing_otherwise(to_radiation, radiation)
    snow = missing_otherwise(to_snowdepth, snowdepth)
    airt = missing_otherwise(to_airtemperature, airtemperature)
    relh = missing_otherwise(to_relhumidity, relhumidity)
    wind = missing_otherwise(to_windspeed, windspeed)
    soilt = missing_otherwise(to_soiltemperature, soiltemperature)
    soilwp = missing_otherwise(to_soilwaterpotential, soilwaterpotential)
    soilwc = missing_otherwise(to_soilwatercontent, soilwatercontent)

    args = (rad, snow, airt, relh, wind, soilt, soilwp, soilwc)
    if any(ismissing.(args))
        missing
    else
        MicroclimPoint{S,I,A,R,typeof.(args)...}(args...)
    end
end

# From microclim dataset in tuples and a cartesian index
# Requires conversion to float from Int
MicroclimPoint{S,I,A,R}(radiation::Vector, snowdepth::Vector, airtemperature::Tuple,
                        relhumidity::Tuple, windspeed::Tuple, soiltemperature::Tuple,
                        soilwaterpotential::Tuple, soilwatercontent::Tuple, i::CartesianIndex) where {S,I,A,R} = begin
    rad = missing_otherwise_cat(x -> to_radiation(0.1x), radiation, i)
    snow = missing_otherwise_cat(x -> to_snowdepth(0.001x), snowdepth, i)
    airt = missing_otherwise_cat(x -> to_airtemperature(0.1x), airtemperature, i)
    relh = missing_otherwise_cat(x -> to_relhumidity(0.001x), relhumidity, i)
    wind = missing_otherwise_cat(x -> to_windspeed(0.1x), windspeed, i)
    soilt = missing_otherwise_cat(x -> to_soiltemperature(0.1x), soiltemperature, i)
    soilwp = missing_otherwise_cat(x -> to_soilwaterpotential(0.1x), soilwaterpotential, i)
    soilwc = missing_otherwise_cat(x -> to_soilwatercontent(0.001x), soilwatercontent, i)

    args = (rad, snow, airt, relh, wind, soilt, soilwp, soilwc)
    if any(ismissing.(args))
        missing
    else
        MicroclimPoint{S,I,A,R,typeof.(args)...}(args...)
    end
end

# From MiriclimGrid and a cartesian index
MicroclimPoint(env::AbstractMicroclimGrid{S,I,A,R}, i::CartesianIndex) where {S,I,A,R} =
    MicroclimPoint{S,I,A,R}(radiation(env), snowdepth(env), airtemperature(env),
                   relhumidity(env), windspeed(env), soiltemperature(env),
                   soilwaterpotential(env), soilwatercontent(env), i)

missing_otherwise(f, d) =
    (length(d) == 0 || any(ismissing.(d))) ? missing : f.(d)

missing_otherwise_cat(f, ds::Tuple, i) = begin
    maybedata = map(d -> missing_otherwise_cat(f, d, i), ds)
    if any(map(ismissing, maybedata))
        missing
    else
        hcat(maybedata...)
    end
end
missing_otherwise_cat(f, d, i) =
    if length(d) == 0
        nothing
    else
        p = catpoint(d, i)
        any(ismissing.(p)) ? missing : f.(p)
    end

@inline layer_sizes(env::MicroclimPoint{S,I,NL,R}) where {S,I,NL,R} = layer_sizes(zero(NL), I, NL)
@inline layer_bounds(env::MicroclimPoint{S,I,NL,R}) where {S,I,NL,R} = layer_bounds(I, NL)

abstract type AbstractMicroclimInstant end

struct MicroclimInstant{M,T,I} <: AbstractMicroclimInstant
    microclimate::M
    t::T
    interp::I
end
MicroclimInstant(microclim, t, height::Number) =
    MicroclimInstant(microclim, t, LinearLayerInterpolators(microclim, height))

@inline layermax(layers, x::MicroclimInstant) = layermax(layers, x.t, x.interp.increment)

struct MicroclimateLayers{L<:Tuple}
    layers::L
end

# unit conversions
to_radiation(x) = x * W*m^-2
to_snowdepth(x) = x * m
to_airtemperature(x) = x * °C |> K
to_relhumidity(x) = x
to_windspeed(x) = x * m*s^-1
to_soiltemperature(x) = x * °C |> K
to_soilwaterpotential(x) = x * kPa
to_soilwatercontent(x) = x

# Get layer specifications from type system
get_shade(layers::AbstractLayeredMicroclimate{S,I,A,R}) where {S,I,A,R} = S
get_increments(layers::AbstractLayeredMicroclimate{S,I,A,R}) where {S,I,A,R} = I
get_nextlayer(layers::AbstractLayeredMicroclimate{S,I,A,R}) where {S,I,A,R} = A
get_range(layers::AbstractLayeredMicroclimate{S,I,A,R}) where {S,I,A,R} = R

@inline catpoint(layers::Tuple{}, i::CartesianIndex) = Float64[]
@inline catpoint(layers::Tuple, i::CartesianIndex) =
    cat([[cat((y[i, :] for y in l)..., dims=1)...] for l in layers]..., dims=2)
@inline catpoint(layers::Vector, i::CartesianIndex) =
    cat((getindex(y, i, :) for y in layers)..., dims=1)

# Accessors for single layer data
for name in (:radiation, :snowdepth)
    eval(quote
        $name(env::AbstractMicroclimate) = env.$name
        $name(env::AbstractMicroclimPoint, t, interp) = lin_interp($name(env), t)
        $name(env::AbstractMicroclimInstant) = lin_interp($name(env.microclimate), env.t)
        $name(env::AbstractMicroclimControl, x, args...) = $name(env)
    end)
end

# Accessor functions for range and increment data
for (name, mode) in ((:airtemperature, :range), (:relhumidity, :range), (:windspeed, :range),
                     (:soiltemperature, :increment), (:soilwaterpotential, :increment),
                     (:soilwatercontent, :increment))
    eval(quote
        $name(env::AbstractMicroclimate) = env.$name
        $name(env::AbstractMicroclimPoint, t, interp::LinearLayerInterpolators) =
            $name(env::AbstractMicroclimPoint, t, interp.$mode)
        $name(env::AbstractMicroclimPoint, t, height::Number) =
            $name(env, t, $(Symbol(mode, :_interpolator))(height))
        $name(env::AbstractMicroclimInstant) =
            $name(env.microclimate, env.t, env.interp)
        $name(env::AbstractMicroclimInstant, interp::AbstractLayerInterpolator) =
            $name(env.microclimate, env.t, interp)
        $name(env::AbstractMicroclimControl, x, args...) = $name(env)
        $name(env::AbstractMicroclimPoint, t, interp::LinearLayerInterpolator) =
            interp_layer($name(env), t, interp)
    end)
end

# Accessors for weighted mean of 8 layer data
for name in (:soiltemperature, :soilwatercontent, :soilwaterpotential)
    eval(quote
        $(Symbol(:mean_, name))(env::AbstractMicroclimPoint, t, x) =
            weightedmean(env, $name(env), t, x)
        $(Symbol(:mean_, name))(env::AbstractMicroclimInstant) =
            weightedmean(env, $name(env.microclimate), env.t, env.interp)
        $(Symbol(:mean_, name))(env::AbstractMicroclimInstant, t::Number) =
            weightedmean(env, $name(env.microclimate), t, env.interp)
        $(Symbol(:mean_, name))(env::AbstractMicroclimInstant, interp::AbstractLayerInterpolator) =
            weightedmean(env, $name(env.microclimate), env.t, interp)
        # No mean for control, just return the control value
        $(Symbol(:mean_, name))(env::AbstractMicroclimControl, args...) = $name(env)
    end)
end


text:
using Pardiso
using Dates

function dolinearsolve(A::SparseMatrixCSC, b::Array, matrixsym::FDFDMatSymmetry)
    time0 = now()

    @info "Solving linear system"
    @info "Unknowns: $(length(b))"

    if haskey(ENV, "FDFD_SOLVER") && lowercase(ENV["FDFD_SOLVER"]) == "pardiso"
        @info "Solver: Pardiso"
        if haskey(ENV, "OMP_NUM_THREADS")
            @info "Threads: $(ENV["OMP_NUM_THREADS"])"
        else
            @warn "OMP_NUM_THREADS was not set. Pardiso may not be using the expected number of cores."
        end

        # Prefer the MKL Pardiso solver over the vanilla Pardiso solver
        ps = nothing
        try
            ps = MKLPardisoSolver()
        catch e
            if e isa LoadError
                ps = PardisoSolver()
            else
                error("Error while loading the Pardiso solver.")
            end
        end
        set_matrixtype!(ps, Pardiso.COMPLEX_NONSYM)
        x = Pardiso.solve(ps, A, b);
        set_phase!(ps, Pardiso.RELEASE_ALL)
        pardiso(ps)
    else
        @info "Solver: Julia"
        x = lu(A)\b
    end

    @info "Solve time: $(now()-time0)"

    return x
end


text:
#=
        project1.jl -- This is where the magic happens!

    All of your code must either live in this file, or be `include`d here.
=#

#=
    If you want to use packages, please do so up here.
    Note that you may use any packages in the julia standard library
    (i.e. ones that ship with the julia language) as well as Statistics
    (since we use it in the backend already anyway)
=#
# import Pkg;
# Pkg.add("Plots")
#import Pkg
#Pkg.add("Plots")

#Pkg.add("PyPlot")
#Pkg.add("GR")

#didn't end up adding yet
#Pkg.add("UnicodePlots")
#Pkg.add("PlotlyJS")

# using Plots
# gr()

# Example:
using LinearAlgebra

#=
    If you're going to include files, please do so up here. Note that they
    must be saved in project1_jl and you must use the relative path
    (not the absolute path) of the file in the include statement.

    [Good]  include("somefile.jl")
    [Bad]   include("/pathto/project1_jl/somefile.jl")
=#

# Example
# include("myfile.jl")

"""
    optimize(f, g, x0, n, prob)

Arguments:
    - `f`: Function to be optimized
    - `g`: Gradient function for `f`
    - `x0`: (Vector) Initial position to start from
    - `n`: (Int) Number of evaluations allowed. Remember `g` costs twice of `f`
    - `prob`: (String) Name of the problem. So you can use a different strategy for each problem. E.g. "simple1", "secret2", etc.

Returns:
    - The location of the minimum
"""
function optimize(f, g, x0, n, prob)
#     Basing first attempt on Nesterov's momentum gradient descent as described in Alg. 5.4 in the course textbook
 	v = vec(zeros(length(x0),1)); #set momentum to zero
	i = 1;
#     x_history = x0; #x is
    #x version 2
    x_history = Array{Array{Float64,1}}(undef, 1);
 	x_history[1] = x0;
    if prob == "simple1"
        alpha = 0.4; #learning rate
        beta = 0.3;#momentum decay, 0 = pure gradient method
    else
        alpha = 0.9; #learning rate
        beta = 0.8;#momentum decay, 0 = pure gradient method
    end
    #version 1
#     while count(f, g) < n
#         xi_vec = vec(x_history[:,i]);
# #         if i == 1
# #             xi_vec = x_history;
# #         else
# #             xi_vec = vec(pop(x_history));
# #         end

#         tmp_vec = xi_vec + beta*v;
#         v = beta*v - (alpha^i)*g(tmp_vec); #adjusted Nestov's alg. to incorporate decreasing learning rate
#         normThresh = 1;
#         if norm(v) > normThresh
#             v = v / norm(v)
#         end
#         xpi_vec = xi_vec + v;
#         x_history = hcat(x_history, xpi_vec);
# #         x_history = push!(x_history, xpi_vec);
# #         println(typeof(x_history))
#         i = i + 1;
#     end
#         #end loop
# #     println(x_history)
# #     println(v)
#     return vec(x_history[:,i-1])
# #     return pop(x_history)

    #version 2
    while count(f, g) < n
        xi_vec = vec(x_history[i]);
#         if i == 1
#             xi_vec = x_history;
#         else
#             xi_vec = pop(x_history); #x_history should be an array of arrays, non-muatating pop
#         end

        tmp_vec = xi_vec + beta*v;
        v = beta*v - (alpha^i)*g(tmp_vec); #adjusted Nestov's alg. to incorporate decreasing learning rate
        normThresh = 1;
        if norm(v) > normThresh
            v = v / norm(v)
        end
        xpi_vec = xi_vec + v;
#         x_history = hcat(x_history, xpi_vec);
         x_history = push!(x_history, xpi_vec);
#         println(typeof(x_history))
        i = i + 1;
    end
#     println(collect(x_history))
        #end loop
#     println(x_history)
#     println(v)
#     return vec(x_history[:,i-1])
#     plot!, [0 0 0], [0 10 50], st = :contour)
    return vec(x_history[i-1])
end

#function to return full history of optimization path
function optimize_hist(f, g, x0, n, prob)
#     Basing first attempt on Nesterov's momentum gradient descent as described in Alg. 5.4 in the course textbook
 	v = vec(zeros(length(x0),1)); #set momentum to zero
	i = 1;
#     x_history = x0; #x is
    #x version 2
    x_history = Array{Array{Float64,1}}(undef, 1);
 	x_history[1] = x0;
    if prob == "simple1"
        alpha = 0.4; #learning rate
        beta = 0.3;#momentum decay, 0 = pure gradient method
    else
        alpha = 0.9; #learning rate
        beta = 0.8;#momentum decay, 0 = pure gradient method
    end
    #version 1
#     while count(f, g) < n
#         xi_vec = vec(x_history[:,i]);
# #         if i == 1
# #             xi_vec = x_history;
# #         else
# #             xi_vec = vec(pop(x_history));
# #         end

#         tmp_vec = xi_vec + beta*v;
#         v = beta*v - (alpha^i)*g(tmp_vec); #adjusted Nestov's alg. to incorporate decreasing learning rate
#         normThresh = 1;
#         if norm(v) > normThresh
#             v = v / norm(v)
#         end
#         xpi_vec = xi_vec + v;
#         x_history = hcat(x_history, xpi_vec);
# #         x_history = push!(x_history, xpi_vec);
# #         println(typeof(x_history))
#         i = i + 1;
#     end
#         #end loop
# #     println(x_history)
# #     println(v)
#     return vec(x_history[:,i-1])
# #     return pop(x_history)

    #version 2
    while count(f, g) < n
        xi_vec = vec(x_history[i]);
#         if i == 1
#             xi_vec = x_history;
#         else
#             xi_vec = pop(x_history); #x_history should be an array of arrays, non-muatating pop
#         end

        tmp_vec = xi_vec + beta*v;
        v = beta*v - (alpha^i)*g(tmp_vec); #adjusted Nestov's alg. to incorporate decreasing learning rate
        normThresh = 1;
        if norm(v) > normThresh
            v = v / norm(v)
        end
        xpi_vec = xi_vec + v;
#         x_history = hcat(x_history, xpi_vec);
         x_history = push!(x_history, xpi_vec);
#         println(typeof(x_history))
        i = i + 1;
#           println("test")
    end

#     println(collect(x_history))
        #end loop
#     println(x_history)
#     println(v)
#     return vec(x_history[:,i-1])
#     plot!, [0 0 0], [0 10 50], st = :contour)
    return x_history
end


text:
## Our rule for chosing γ from μ

function γ_from_μ(μ::R, settings::Settings) where {R <: Real}
    if settings.γ_updt_rule == :safe
        return abs(cbrt(settings.μ_min)) # cbrt also works quite well
    elseif settings.γ_updt_rule == :adaptive
        γ_adaptive = (10^-3)*abs(sqrt(μ))
        return γ_adaptive # adaptive actually does not work very well
    elseif settings.γ_updt_rule == :adaptiveplus
        log_γ_adaptiveplus = ((3/(log10(settings.μ_max)+10))*(log10(μ)+10))-6
        γ_adaptiveplus = 10^log_γ_adaptiveplus
        return γ_adaptiveplus
    else
        @error "either select :adaptive or :safe or :adaptiveplus as your γ selection rule"
        return nothing
    end
end

# # Projection onto the sparse set
#  First thing to do is, projecting onto the sparse set. Projection on the sparse set can be computed as follows.

function Π_exact(C::SparseSet, x::AbstractVector{R}) where {R <: Real}
    ## this function takes a vector x and returns a projection of x on the set
    ## card(⋅) <= k, ||⋅||_∞ <= M
    ## Extract the information from teh SparseSet
    M = C.M
    k = C.r

    n = length(x)
    y = zero(x)

    perm = sortperm(abs.(x), rev=true)
    ## sortperm(abs.(x), rev=true) returns a sorted array index set from larger to smaller components of absolute values of x
    for i in 1:n
        if i in perm[1:k]
            if x[i]>M
                y[i] = M
            elseif x[i] < -M
                y[i] = -M
            else
                y[i] = x[i]
            end
        end
    end

    return y
end

# Finally, the projection onto the low rank set.
# # Pseudocode for `Π_exact`
# ---

# ![Algorithm](/literate_assets/low_rank_projection.png)

function Π_exact(C::RankSet, x::AbstractMatrix{R}; iterative::Bool = true) where {R <: Real}

    ## extract the values
    y = similar(x)
    M = C.M # bound on the matrix elements
    r = C.r # rank constraint rank(x) ≦ r
    if iterative == true
        U, S, V = tsvd(x, r)
        for i in 1:r
            if S[i] > M
                S[i] = M
            end
        end
        SVt = S.* V' # this line essentially does the same as Diagonal(X)*V' but more efficiently
        mul!(y, U, SVt)
        return y
    else
        F = svd(x)
        σ_x = F.S
        for i in 1:r
            if σ_x[i] > M
                σ_x[i] = M
            end
        end
        y = F.U[:,1:r]*Diagonal(σ_x[1:r])*V[:,1:r]'
        return y
    end # if else
end

# Projection onto the unit hypercube.

function Π_exact(C::UnitHyperCube, x::AbstractVector{R}) where {R <: Real}
    ## this function takes a vector x and returns a projection of x on the set {0,1}^n
    n = C.n
    y = zero(x) # creates a 0 vector with the same dimension as x
    for i=1:n
        if x[i] > 0.5
            y[i]=1
        end
    end
    return y
end

# update the state

function update_state!(state::State, init_info::InitInfo, problem::Problem, settings::Settings)

    # extract information from the state and init info
    x = state.x
    y = state.y
    z = init_info.z0
    β = settings.β
    γ = init_info.γ
    μ = init_info.μ

    # create best point variables
    best_point_x = x
    best_point_y = y
    best_point_z = z
    best_fxd_pnt_gap = settings.big_M
    best_fsblt_gap = settings.big_M

    i = 1
    while i <= settings.n_iter_max
        x, y, z = inner_iteration(x, y, z, β, γ, μ, problem)
        crnt_fxd_pnt_gap = norm(x-y, Inf)
        y_fsbl = Π_exact(problem.C, y)
        crnt_fsblt_gap = norm(y-y_fsbl, Inf)

        # update the best points so far if we have lower objective value
        if crnt_fxd_pnt_gap <= best_fxd_pnt_gap
            best_point_x = x
            best_point_y = y
            best_point_z = z
            best_fxd_pnt_gap = crnt_fxd_pnt_gap
            best_fsblt_gap = crnt_fsblt_gap
        end # if

        # display important information
        if settings.verbose == true
            if mod(i, settings.freq) == 0
                @info "μ = $μ | iteration = $i | log fixed point gap = $(log10(crnt_fxd_pnt_gap)) | log feasibility gap = $(log10(crnt_fsblt_gap))"
            end # if
        end # if

        # termination condition check
        if crnt_fxd_pnt_gap <= settings.tol
            if settings.verbose == true
                @info "μ = $μ | log fixed point gap reached $(log10(settings.tol))"
            end
            break
        end #if

        # increment the count by 1
        i = i+1

    end # while

    ## update the state
    state.x = best_point_x
    state.y = best_point_y
    state.z = best_point_z
    state.i = i
    state.fxd_pnt_gap = best_fxd_pnt_gap
    state.fsblt_gap = best_fsblt_gap
    state.μ = μ
    state.γ = γ

    # Print information about the best point
    if settings.verbose == true
        @info "information about the best state found for μ = $(state.μ)"
        @info "μ = $(state.μ) | log fixed point gap = $(log10(state.fxd_pnt_gap)) | log feasibility gap = $(log10(state.fsblt_gap)) | inner iterations = $(state.i)"
    end

    return state
end

## the inner iteration function
# The inner iteration function has the following pseudocde
# # Pseudocode for `Π_exact`
# ---

# ![Algorithm](/literate_assets/inner_iteration.png)

function inner_iteration(x::A, y::A, z::A, β::R, γ::R, μ::R, problem::Problem) where {A <: AbstractVecOrMat{<:Real}, R <: Real}
    # old iteration: x = prox_NExOS(problem.f, β, γ, z)
    x = prox_NExOS(problem.f, γ, z)
    y = Π_NExOS(problem.C, β, γ, μ, 2*x - z)
    z = z + y - x
    return x, y, z
end

## proximal function on f: this is the old function to be removed later
# function prox_NExOS(f::ProximableFunction, β::R, γ::R, x::A) where {R <: Real, A <: AbstractVecOrMat{<:Real}}
#     prox_param = γ/((β*γ)+1)
#     scaled_input = (1/((β*γ)+1)).*x
#     prox_point = similar(x) # allocate prox_NExOS output
#     prox!(prox_point, f, scaled_input, prox_param) # this is the prox! function from the ProximalOperators package
#     return prox_point
# end

## New proximal function, no scaling required
function prox_NExOS(f::ProximableFunction, γ::R, x::A) where {R <: Real, A <: AbstractVecOrMat{<:Real}}
    prox_point = similar(x) # allocate prox_NExOS output
    prox!(prox_point, f, x, γ) # this is the prox! function from the ProximalOperators package
    return prox_point
end

## NExOS projection function
# This implements the following pseudocode
# ![relaxed_projection_algorithm](/literate_assets/relaxed_projection.png)

function Π_NExOS(C::ProxRegularSet, β::R, γ::R, μ::R, x::A) where {R <: Real, A <: AbstractVecOrMat{<:Real}}
    dnm = γ + (μ*((β*γ)+1))
    scaled_input = (1/((β*γ)+1)).*x
    pi_x_scaled_input = Π_exact(C, scaled_input)
    t1 = (μ/dnm).*x
    t2 = (γ/dnm)*pi_x_scaled_input
    return t1+t2
end

# Now we write the update init information function.

function update_init_info!(state::State, init_info::InitInfo, problem::Problem, settings::Settings)

    init_info.μ = init_info.μ*settings.μ_mult_fact
    init_info.γ = γ_from_μ(init_info.μ, settings)
    x_μ = Π_NExOS(problem.C, settings.β, init_info.γ, init_info.μ, state.x)
    u1 = similar(x_μ) # allocate memory
    gradient!(u1, problem.f, x_μ)
    u = u1 + settings.β.*x_μ # this is a gradient of the function f + (β/2)*||⋅||^2 at x_μ
    # we are using the function gradient from proximaloperators.jl
    init_info.z0 = x_μ + init_info.γ*u

    return init_info

end

# Experiment with the initial information

function update_init_info_experimental!(state::State, init_info::InitInfo, problem::Problem, settings::Settings)

    init_info.μ = init_info.μ*settings.μ_mult_fact
    init_info.γ = γ_from_μ(init_info.μ, settings)
    init_info.z0 = state.z

    return init_info

end

## Extending the proximal opertor for least sqaures penalty over matrices

function ProximalOperators.prox!(Y::Array{R,2}, f::ProximalOperators.LeastSquaresIterative, X::Array{R,2}, gamma::R=R(1)) where {R <: Real}
    m, n = size(Y)
    y = similar(vec(Y)) # create storaget to store the proximal operator
    x = vec(X)
    prox!(y, f, x, gamma)
    # need to convert y back to Y
    Y[:,:] = reshape(y,m,n) # Y[:,:] is used so that the contents are changed
    return f(y)
end

function ProximalOperators.prox!(Y::Array{R,2}, f::ProximalOperators.LeastSquaresDirect, X::Array{R,2}, gamma::R=R(1)) where {R <: Real}
    m, n = size(Y)
    y = similar(vec(Y)) # create storaget to store the proximal operator
    x = vec(X)
    prox!(y, f, x, gamma)
    # need to convert y back to Y
    Y[:,:] = reshape(y,m,n) # Y[:,:] is used so that the contents are changed
    return f(y)
end

# Computing gradient for the leastsquares over matrix function direct case
function ProximalOperators.gradient!(Y::Array{R,2}, f::ProximalOperators.LeastSquaresDirect, X::Array{R,2}) where {R<: Real}
    m, n = size(Y)
    y = similar(vec(Y)) # create storaget to store the proximal operator
    x = vec(X)
    gradient!(y, f, x)
    Y[:,:] = reshape(y,m,n) # Y[:,:] is used so that the contents are changed
    return f(y)
end

# Computing gradient for the leastsquares over matrix function direct case
function ProximalOperators.gradient!(Y::Array{R,2}, f::ProximalOperators.LeastSquaresIterative, X::Array{R,2}) where {R<: Real}
    m, n = size(Y)
    y = similar(vec(Y)) # create storaget to store the proximal operator
    x = vec(X)
    gradient!(y, f, x)
    Y[:,:] = reshape(y,m,n) # Y[:,:] is used so that the contents are changed
    return f(y)
end


text:
using Raylib
using Test

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


text:
using NeuralQuantum, Random, Plots

Random.seed!(Random.GLOBAL_RNG, 1234)

# Parameters of the transverse field ising model
N = 20
h = 1.0
J = 1.0

# Constructs the Hilbert space for N 1//2 spins.
hilb = HomogeneousSpin(N, 1//2)

# Builds the hamiltonian
H = LocalOperator(hilb)
for i=1:N
    global H  -= h * sigmax(hilb, i)
    global H  += J * sigmaz(hilb, i) * sigmaz(hilb, mod(i, N)+1)
end

# Build the FFNN with one dense layer and a final sum layer
# instead of using sum, we use sum_autobatch, so that we can correctly
# batch calls to the FFNN.
ch = Chain(Dense(ComplexF64, N, 2*N, af_logcosh), sum_autobatch)

# Wrap the chain (FFNN) into an object correctly describing a Pure State
net = PureStateAnsatz(ch, N)

# Randomize all parameters
init_random_pars!(net, sigma=0.01)

sampl = MetropolisSampler(LocalRule(), 75, N, burn=100)
algo  = SR(ϵ=(0.1), algorithm=sr_cg, precision=1e-3)

is = BatchedSampler(net, sampl, H, algo; batch_sz=16)

optimizer = Optimisers.Descent(0.05)

Evalues = Float64[];
Eerr = Float64[];
for i=1:300
    ldata, prec = sample!(is)
    ob = compute_observables(is)

    println("$i - $ldata")

    push!(Evalues, real(ldata.mean))
    push!(Eerr, ldata.error)
    grad = precondition!(prec, algo, i)
    Optimisers.update!(optimizer, net, grad)
end

plot(Evalues, yerr=Eerr)

# N=20 (thanks netket)
exact = -1.274549484318e00 * 20
hline!([exact, exact])