text:
using Test
using Clustering
using Distances
using Random
using LinearAlgebra

import Distances.pairwise!

# custom distance metric
struct MySqEuclidean <: SemiMetric end

# redefinition of Distances.pairwise! for MySqEuclidean type
function pairwise!(r::AbstractMatrix, dist::MySqEuclidean,
                   a::AbstractMatrix, b::AbstractMatrix; dims::Integer=2)
    dims == 2 || throw(ArgumentError("only dims=2 supported for MySqEuclidean distance"))
    mul!(r, transpose(a), b)
    sa2 = sum(abs2, a, dims=1)
    sb2 = sum(abs2, b, dims=1)
    @inbounds r .= sa2' .+ sb2 .- 2r
end

@testset "kmeans() (k-means)" begin

Random.seed!(34568)

m = 3
n = 1000
k = 10

x = rand(m, n)
xt = copy(transpose(x))

equal_kmresults(km1::KmeansResult, km2::KmeansResult) =
    all(getfield(km1, η) == getfield(km2, η) for η ∈ fieldnames(KmeansResult))

@testset "non-weighted" begin
    Random.seed!(34568)
    r = kmeans(x, k; maxiter=50)
    @test isa(r, KmeansResult{Matrix{Float64}, Float64, Int})
    @test size(r.centers) == (m, k)
    @test length(r.assignments) == n
    @test all(a -> 1 <= a <= k, r.assignments)
    @test length(r.costs) == n
    @test length(r.counts) == k
    @test sum(r.counts) == n
    @test r.cweights == map(Float64, r.counts)
    @test sum(r.costs) ≈ r.totalcost

    Random.seed!(34568)
    r_t = kmeans(xt', k; maxiter=50)
    @test equal_kmresults(r, r_t)
end

@testset "non-weighted (float32)" begin
    Random.seed!(34568)
    x32 = map(Float32, x)
    x32t = copy(x32')
    r = kmeans(x32, k; maxiter=50)
    @test isa(r, KmeansResult{Matrix{Float32}, Float32, Int})
    @test size(r.centers) == (m, k)
    @test length(r.assignments) == n
    @test all(a -> 1 <= a <= k, r.assignments)
    @test length(r.costs) == n
    @test length(r.counts) == k
    @test sum(r.counts) == n
    @test r.cweights == map(Float64, r.counts)
    @test sum(r.costs) ≈ r.totalcost

    Random.seed!(34568)
    r_t = kmeans(x32t', k; maxiter=50)
    @test equal_kmresults(r, r_t)
end

@testset "weighted" begin
    w = rand(n)
    Random.seed!(34568)
    r = kmeans(x, k; maxiter=50, weights=w)
    @test isa(r, KmeansResult{Matrix{Float64}, Float64, Float64})
    @test size(r.centers) == (m, k)
    @test length(r.assignments) == n
    @test all(a -> 1 <= a <= k, r.assignments)
    @test length(r.costs) == n
    @test length(r.counts) == k
    @test sum(r.counts) == n

    cw = zeros(k)
    for i = 1:n
        cw[r.assignments[i]] += w[i]
    end
    @test r.cweights ≈ cw
    @test dot(r.costs, w) ≈ r.totalcost

    Random.seed!(34568)
    r_t = kmeans(xt', k; maxiter=50, weights=w)
    @test equal_kmresults(r, r_t)
end

@testset "custom distance" begin
    Random.seed!(34568)
    r = kmeans(x, k; maxiter=50, init=:kmcen, distance=MySqEuclidean())
    r2 = kmeans(x, k; maxiter=50, init=:kmcen)
    @test isa(r, KmeansResult{Matrix{Float64}, Float64, Int})
    @test size(r.centers) == (m, k)
    @test length(r.assignments) == n
    @test all(a -> 1 <= a <= k, r.assignments)
    @test length(r.costs) == n
    @test length(r.counts) == k
    @test sum(r.counts) == n
    @test r.cweights == map(Float64, r.counts)
    @test sum(r.costs) ≈ r.totalcost
    @test equal_kmresults(r, r2)

    Random.seed!(34568)
    r_t = kmeans(xt', k; maxiter=50, init=:kmcen, distance=MySqEuclidean())
    @test equal_kmresults(r, r_t)
end

@testset "Integer data" begin
    x = rand(Int16, m, n)
    Random.seed!(654)
    r = kmeans(x, k; maxiter=50)

    @test isa(r, KmeansResult{Matrix{Float64}, Float64, Int})
end

@testset "kmeans! data types" begin
    Random.seed!(1101)
    for TX in (Int, Float32, Float64)
        for TC in (Float32, Float64)
            for TW in (Nothing, Int, Float32, Float64)
                x = rand(TX, m, n)
                c = rand(TC, m, k)
                if TW == Nothing
                    r = kmeans!(x, c; maxiter=1)
                    @test isa(r, KmeansResult{Matrix{TC},<:Real,Int})
                else
                    w = rand(TW, n)
                    r = kmeans!(x, c; weights=w, maxiter=1)
                    @test isa(r, KmeansResult{Matrix{TC},<:Real,TW})
                end
            end
        end
    end
end
end


text:
# This is part of Quiz.jl. License is MIT.

"Returns the engine metadata of the engine type argument."
function meta end

export meta

ext: jl


text:
module Nexus

using SimpleTraits
using Lazy

import DispatchOps
import Base:
    show, eltype, Tuple, iterate, length, reverse,
    getindex, rand

include("interface.jl")

include("./simplenexus/arc.jl")
include("./simplenexus/digraph.jl")
include("./simplenexus/graph.jl")
include("./simplenexus/arciter.jl")

include("./metanexus/extension.jl")
include("./metanexus/metadigraph.jl")
include("./metanexus/metagraph.jl")

export
    Graph, Digraph, Arc, ArcIter,
    add_arc!, rem_arc!, add_node!, rem_node!,
    nn, na, nodes, arcs, fadj, badj, has_arc, has_node,
    src, tgt, key,

    MetaGraph, MetaDigraph, props, get_prop, has_prop,
    set_props!, set_prop!, rem_prop!, clear_props!,
    filter_arcs, filter_nodes

end


text:
# Combining invalidation and snoopi_deep data

struct StaleTree
    method::Method
    reason::Symbol   # :inserting or :deleting
    mt_backedges::Vector{Tuple{Any,Union{InstanceNode,MethodInstance},Vector{InferenceTimingNode}}}   # sig=>root
    backedges::Vector{Tuple{InstanceNode,Vector{InferenceTimingNode}}}
end
StaleTree(method::Method, reason) = StaleTree(method, reason, staletree_storage()...)

staletree_storage() = (
    Tuple{Any,Union{InstanceNode,MethodInstance},Vector{InferenceTimingNode}}[],
    Tuple{InstanceNode,Vector{InferenceTimingNode}}[])

function Base.show(io::IO, tree::StaleTree)
    iscompact = get(io, :compact, false)::Bool

    print(io, tree.reason, " ")
    printstyled(io, tree.method, color = :light_magenta)
    println(io, " invalidated:")
    indent = iscompact ? "" : "   "

    if !isempty(tree.mt_backedges)
        print(io, indent, "mt_backedges: ")
        showlist(io, tree.mt_backedges, length(indent)+length("mt_backedges")+2)
    end
    if !isempty(tree.backedges)
        print(io, indent, "backedges: ")
        showlist(io, tree.backedges, length(indent)+length("backedges")+2)
    end
    iscompact && print(io, ';')
end

function printdata(io, tnodes::AbstractVector{InferenceTimingNode})
    if length(tnodes) == 1
        print(io, tnodes[1])
    else
        print(io, sum(inclusive, tnodes), " inclusive time for $(length(tnodes)) nodes")
    end
end

"""
    staletrees = precompile_blockers(invalidations, tinf::InferenceTimingNode)

Select just those invalidations that contribute to "stale nodes" in `tinf`, and link them together.
This can allow one to identify specific blockers of precompilation for particular MethodInstances.

# Example

```julia
using SnoopCompileCore
invalidations = @snoopr using PkgA, PkgB;
using SnoopCompile
trees = invalidation_trees(invalidations)
tinf = @snoopi_deep begin
    some_workload()
end
staletrees = precompile_blockers(trees, tinf)
```
"""
function precompile_blockers(trees::Vector{MethodInvalidations}, tinf::InferenceTimingNode; kwargs...)
    sig2node = nodedict!(IdDict{Type,InferenceTimingNode}(), tinf)
    snodes = stalenodes(tinf; kwargs...)
    mi2stalenode = Dict(MethodInstance(node) => i for (i, node) in enumerate(snodes))
    # Prepare "thinned trees" focusing just on those invalidations that blocked precompilation
    staletrees = StaleTree[]
    for tree in trees
        mt_backedges, backedges = staletree_storage()
        for (sig, root) in tree.mt_backedges
            triggers = Set{MethodInstance}()
            thinned = filterbranch(node -> haskey(mi2stalenode, convert(MethodInstance, node)), root, triggers)
            if thinned !== nothing
                push!(mt_backedges, (sig, thinned, [snodes[mi2stalenode[mi]] for mi in triggers]))
            end
        end
        for root in tree.backedges
            triggers = Set{MethodInstance}()
            thinned = filterbranch(node -> haskey(mi2stalenode, convert(MethodInstance, node)), root, triggers)
            if thinned !== nothing
                push!(backedges, (thinned, [snodes[mi2stalenode[mi]] for mi in triggers]))
            end
        end
        if !isempty(mt_backedges) || !isempty(backedges)
            sort!(mt_backedges; by=suminclusive)
            sort!(backedges; by=suminclusive)
            push!(staletrees, StaleTree(tree.method, tree.reason, mt_backedges, backedges))
        end
    end
    sort!(staletrees; by=inclusive)
    return staletrees
end

precompile_blockers(invalidations, tinf::InferenceTimingNode; kwargs...) =
    precompile_blockers(invalidation_trees(invalidations)::Vector{MethodInvalidations}, tinf; kwargs...)

function nodedict!(d, tinf::InferenceTimingNode)
    for child in tinf.children
        sig = MethodInstance(child).specTypes
        oldinf = get(d, sig, nothing)
        if oldinf === nothing || inclusive(child) > inclusive(oldinf)
            d[sig] = child
        end
        nodedict!(d, child)
    end
    return d
end

suminclusive(t::Tuple) = sum(inclusive, last(t))
SnoopCompileCore.inclusive(tree::StaleTree) =
    sum(suminclusive, tree.mt_backedges; init=0.0) +
    sum(suminclusive, tree.backedges; init=0.0)


text:
using DataDeps

register(DataDep(
  "parkinsons-telemonitoring-total",
  "https://archive.ics.uci.edu/ml/datasets/Parkinsons+Telemonitoring",
  "https://archive.ics.uci.edu/ml/machine-learning-databases/parkinsons/telemonitoring/parkinsons_updrs.data",
  "f2c7d5025dec4e92e7feae367a5f7ccf58789a10ac6b54bdf15976c599f9dd39",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=6,
        feature_indices=7:22,
        skipto=2,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
# This file is part of IntegerTriangles.jl.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module TrianglesTables

using PrettyTables, TrianglesExplorer
using TrianglesBase, TrianglesExamples, TrianglesTraitCard, TrianglesUtils

export PrettyTraits, AllPrettyTraits, UpdateDocs

"""
Pretty printing of triangles trait cards.

| A-number | Triangle   | Form | Function  | Sequence                                    |
|----------|------------|------|-----------|---------------------------------------------|
| A000302  | Binomial   | Std  | PolyVal3  | 1, 4, 16, 64, 256, 1024, 4096, 16384        |
| A001333  | SchroederB | Inv  | AltSum    | 1, -1, 3, -7, 17, -41, 99, -239             |
| A006012  | SchroederL | Inv  | AltSum    | 1, -2, 6, -20, 68, -232, 792, -2704         |
| A026302  | Motzkin    | Rev  | Central   | 1, 2, 9, 44, 230, 1242, 6853, 38376         |
| A103194  | Laguerre   | Std  | TransNat0 | 0, 1, 6, 39, 292, 2505, 24306, 263431       |
| A111884  | Lah        | Std  | TransAlts | 1, -1, -1, -1, 1, 19, 151, 1091             |
| nothing  | Laguerre   | Rev  | TransNat1 | 1, 3, 15, 97, 753, 6771, 68983, 783945      |
"""
const ModuleTrianglesTables = ""

const CSS = ["<style> table, td, th {border-collapse: collapse; font-family: sans-serif; color: blue;}",
 "td, th {border-bottom: 0; padding: 4px}",
 "tr:nth-child(odd) {background: #eee;}", "tr:nth-child(even) {background: #fff;}",
 "tr.header {background: orange !important; color: white; font-weight: bold;}",
 "tr.subheader {background: lightgray !important; color: black;}",
 "tr.headerLastRow {border-bottom: 2px solid black;}",
 "th.rowNumber, td.rowNumber {text-align: right;} </style><body>"]

function PrettyTraits(T, name)

    header = [ "Trait", "ANumber", "Sequence" ]
    mat = String[name "" ""]

    for tra in TRAITS
        s = Explore(T, tra)
        s == [] && continue
        mat = [mat; reshape(s, 1, 3)]
    end

    path = datapath(name * ".md")
    rm(path; force=true)
    open(path, "w") do io
        println(io, "```@raw html")
        for ln in CSS
            println(io, ln)
        end
        pretty_table(io, mat, header, backend = :html, standalone= false, alignment=[:l,:l,:l])
        println(io, "```")
    end
end

function AllPrettyTraits()
    for T in Triangles
        println("-- ", T[1])
        PrettyTraits(Triangles[T[1]](32), T[1])
    end
end

function AllCSVtoTable()
    header = ["ANumber", "Triangle", "Type", "Trait", "Sequence"]

    mat = String["" "" "" "" ""]

    inpath = datapath("BIGLIST.csv")
    open(inpath, "r") do csv
        for L ∈ eachline(csv, keep=false)
            mat = [mat; reshape(split(L, ","), 1, 5)]
        end
    end

    outpath = datapath("BIGLISTTEST.md")
    rm(outpath; force=true)
    open(outpath, "w") do html
        println(html, "```@raw html")
        for ln in CSS
            println(html, ln)
        end
        pretty_table(html, mat, header, backend = :html, standalone= false,
                     alignment=[:l,:l,:l,:l,:l])
        println(html, "```")
    end
end

function CSVtoTable(name)
    header = ["Trait", "ANumber", "Sequence"]

    # Depending on the number of csv data:
    #mat = String["" "" "" "" ""]
    #mat = [mat; reshape(split(L, ","), 1, 5)]
    #alignment=[:l,:l,:l,:l,:l])

    p = splitext(name)[begin]
    p = splitpath(p)[end]

    mat = String[p "" ""]
    firstline = true

    inpath = profilespath(name)
    open(inpath, "r") do csv
        for L ∈ eachline(csv, keep=false)
            firstline && (firstline = false; continue)
            s = split(L, ",")
            s[2] = MakeSeqUri(s[2], s[3])
            mat = [mat; reshape(s, 1, 3)]
        end
    end

    outpath = docsrcpath(p * ".md")
    rm(outpath; force=true)
    open(outpath, "w") do html
        println(html, "```@raw html")
        for ln in CSS
            println(html, ln)
        end
        pretty_table(html, mat, header, backend = :html, standalone= false, alignment=[:l,:l,:l])
        println(html, "```")
    end
end

function UpdateDocs()
    csvfiles = csv_files()
    for filename ∈ csvfiles
        fulname = profilespath(filename)
        CSVtoTable(fulname)
    end
end

# START-TEST-########################################################

function test()
    header = ["A-number", "Triangle", "Form", "Function", "Sequence"]
    mat = [  # Type Matrix{String}, or Array{String,2}
    "A000302" "Binomial" "Std" "PolyVal3" "1, 4, 16, 64, 256, 1024, 4096, 16384" ;
    ]
    S = String["A111884" "Lah" "Std" "TransAlts" "1, -1, -1, -1, 1, 19, 151, 1091"]
    mat = [mat; S]
    L = String["A111884", "Lah", "Std", "TransAlts", "1, -1, -1, -1, 1, 19, 151, 1091"]
    mat = [mat; reshape(L, 1, 5)]
    L = Explore("Laguerre", "Rev", "TransNat1")
    mat = [mat; reshape(L, 1, 5)]
    pretty_table(mat, header, alignment=[:l,:l,:l,:l,:l])
    mat
end

function demo()
    PrettyTraits(LahTriangle(32), "Lah")
    PrettyTraits(A046802Triangle(32), "A046802")
end

function perf()
   #AllPrettyTraits()
   #AllCSVtoTable()
   #UpdateDocs()
end

function main()
   test()
   #demo()
   #perf()
end

main()

end # module


text:
using Test
using AtomicLevels
# ^^^ -- to make it possible to run the test file separately

module ATSPParser
    using Test
    using AtomicLevels
    import AtomicLevels: CSF, csfs
    include("atsp/csfparser.jl")

    compare_with_atsp(f, csfout) = @testset "ATSP CSFs: $(csfout)" begin
        atsp_csfs = parse_csf(joinpath(@__DIR__, "atsp", csfout))
        atlev_csfs = f()
        @test length(atlev_csfs) == length(atsp_csfs)
        for atlev_csf in atlev_csfs
            @test atlev_csf in atsp_csfs
        end
    end

    export compare_with_atsp
end

using .ATSPParser

@testset "CSFs" begin
    @testset "Construction" begin
        csf = CSF(rc"1s2 2p- 2p", [IntermediateTerm(0,Seniority(0)), IntermediateTerm(1//2,Seniority(1)), IntermediateTerm(3//2,Seniority(1))], [0, 1//2, 2])
        @test csf isa CSF{RelativisticOrbital{Int},HalfInt}
        @test csf isa RelativisticCSF
        @test csf == csf
        @test csf != CSF(rc"1s2 2p- 2p", [IntermediateTerm(0,Seniority(0)), IntermediateTerm(1//2,Seniority(1)), IntermediateTerm(3//2,Seniority(1))], [0, 1//2, 1])
        @test num_electrons(csf) == 4

        @test CSF(c"1s2 2p", [IntermediateTerm(T"1S",Seniority(0)), IntermediateTerm(T"2Po", Seniority(1))], [T"1S", T"2Po"]) isa NonRelativisticCSF
    end

    @testset "CSF list generation" begin
        @testset "ls coupling" begin
            csf_1s2 = CSF(c"1s2", [IntermediateTerm(T"1S",Seniority(0))],[T"1S"])
            @test csfs(c"1s2") == [csf_1s2]
            @test string(csf_1s2) == "1s²(₀¹S|¹S)+"

            csfs_1s_kp = [CSF(c"1s kp",
                              [IntermediateTerm(T"2S",Seniority(1)),IntermediateTerm(T"2Po",Seniority(1))],
                              [T"2S", T"1Po"]),
                          CSF(c"1s kp",
                              [IntermediateTerm(T"2S",Seniority(1)),IntermediateTerm(T"2Po",Seniority(1))],
                              [T"2S", T"3Po"])]
            @test csfs(c"1s kp") == csfs_1s_kp

            #= zgenconf input
              OI
                1s  2s
                 2   4  -1   2   0   2   n_orbitals,  n_electrons, parity, n_ref, k_min, k_max
                2p  3d
                 0   0
                 6  10
            =#
            compare_with_atsp("1s2c_2s2c_2p1_3_3d3_1.csfs") do
                csfs(c"1s2c 2s2c" ⊗ [c"2p1 3d3", c"2p3 3d1"])
            end
        end
        @testset "jj coupling" begin
            # These are not real tests, they only make sure that the
            # generation does not crash.
            sort(csfs([rc"3s 3p2",rc"3s 3p- 3p"]))
            csfs(rc"[Kr] 5s2 5p-2 5p3 6s")
            csfs(rc"1s kp")

            @test string(CSF(rc"[Kr] 5s2 5p- 5p4 kd",
                             [IntermediateTerm(0//1, Seniority(0)), IntermediateTerm(1//2, Seniority(1)), IntermediateTerm(0//1, Seniority(0)), IntermediateTerm(5//2, Seniority(1))],
                             [0//1, 1//2, 1//2, 3//1])) ==
                "[Kr]ᶜ 5s²(₀0|0) 5p-(₁1/2|1/2) 5p⁴(₀0|1/2) kd(₁5/2|3)-"
        end
    end
end


text:
using QuantumOpticsControlExtensions
using Documenter

DocMeta.setdocmeta!(QuantumOpticsControlExtensions, :DocTestSetup, :(using QuantumOpticsControlExtensions); recursive=true)

makedocs(;
    modules=[QuantumOpticsControlExtensions],
    authors="Michael Goerz <mail@michaelgoerz.net> and contributors",
    repo="https://github.com/JuliaQuantumControl/QuantumOpticsControlExtensions.jl/blob/{commit}{path}#{line}",
    sitename="QuantumOpticsControlExtensions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://juliaquantumcontrol.github.io/QuantumOpticsControlExtensions.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaQuantumControl/QuantumOpticsControlExtensions.jl",
)


text:
@testset "Testing operations and Functions" begin

    layout(3,1,"Conformal")
    u=1.0*id+2.0*e1-e0
    v=2.0*id+1.0*e1+e∞
    # sum
    @test u+v == 3.0*id + 3.0*e1 + 1.5*ep + 0.5*en

    # difference
    @test u-v ==-1.0*id + 1.0*e1 - 0.5*ep - 1.5*en

    # scalar
    @test 1.5*v ==3.0*id + 1.5*e1 + 1.5*ep + 1.5*en
    @test scalar(1.5,v) == 3.0*id + 1.5*e1 + 1.5*ep + 1.5*en

    # reverse
    @test reverse(v) == 2.0*id + 1.0*e1 + 1.0*ep + 1.0*en

    # geometric
    @test u∘v == 5.0*id + 5.0*e1 + 2.0*ep + 1.5*e1p + 2.5*e1n + 1.0*epn
    @test geometric(u,v) == 5.0*id + 5.0*e1 + 2.0*ep + 1.5*e1p + 2.5*e1n + 1.0*epn

    # inner
    @test u⋅v == 3.0*id
    @test inner(u,v) == 3.0*id

    # outer
    @test u^v == 2.0*id + 5.0*e1 + 2.0*ep + 1.5*e1p + 2.5*e1n + 1.0*epn
    @test outer(u,v) == 2.0*id + 5.0*e1 + 2.0*ep + 1.5*e1p + 2.5*e1n + 1.0*epn

    # gradeprojection
    @test gradeprojection(u,1) == 2.0*e1 + 0.5*ep - 0.5*en
    @test gradeprojection(v,1) == 1.0*e1 + 1.0*ep + 1.0*en

    # reverse
    @test reverse(u) == 1.0*id + 2.0*e1 + 0.5*ep - 0.5*en
    @test reverse(v) == 2.0*id + 1.0*e1 + 1.0*ep + 1.0*en

    # dual
    @test dual(u) ==  0.5*e12p - 0.5*e12n - 2.0*e2pn - 1.0*e12pn
    @test dual(v) == - 1.0*e12p - 1.0*e12n - 1.0*e2pn - 2.0*e12pn

    # inverse
    @test inverse(u)∘u ≈ 1.0*id
    @test inverse(v)∘v ≈ 1.0*id

end


text:
module MOITests

using BARON
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

const optimizer = MOIU.CachingOptimizer(BARON.Model{Float64}(), BARON.Optimizer(PrLevel=0));

# TODO: test infeasibility certificates, duals.

@testset "MOI Continuous Linear" begin
    config = MOIT.TestConfig(atol=1e-5, rtol=1e-4, infeas_certificates=false, duals=false)
    excluded = String[
        "linear7", # vector constraints
        "linear8b", # certificate provided in this case (result count is 1)
        "linear8c", # should be unbounded below, returns "Preprocessing found feasible solution with value -.200000000000E+052"
        "linear15", # vector constraints
        "partial_start" # TODO
    ]
    # MOIT.partial_start_test(optimizer, config)
    MOIT.contlineartest(optimizer, config, excluded)
    MOIT.linear8btest(optimizer, MOIT.TestConfig(atol=1e-5, rtol=1e-4, infeas_certificates=true, duals=false))
end

@testset "MOI Integer Linear" begin
    config = MOIT.TestConfig(atol=1e-5, rtol=1e-4, infeas_certificates=false, duals=false)
    excluded = String[
        "int2" # SOS1
    ]
    MOIT.intlineartest(optimizer, config, excluded)
end

@testset "MOI Continuous Quadratic" begin
    # TODO: rather high tolerances
    config = MOIT.TestConfig(atol=1e-3, rtol=1e-3, infeas_certificates=false, duals=false)
    excluded = String[
        "qcp1", # vector constraints
    ]
    MOIT.contquadratictest(optimizer, config, excluded)
end

# @testset "MOI Nonlinear" begin
#     MOIT.nonlineartest(optimizer, config)
# end

end # module


text:
"""
PolaronMobility.jl - https://github.com/jarvist/PolaronMobility.jl
  Codes by Jarvist Moore Frost, 2017-2022; Bradley A.A. Martin 2020-2022
  PolaronMobility.jl is a Julia package which calculates the temperature-dependent polaron
  mobility for a material. It implements the Feynman variational approach with
  Osaka/Hellwarth free-energy for the polaron system.  """
module PolaronMobility

export Polaron, NewPolaron # Type to hold the data
export frohlichalpha, feynmanvw, F, polaronmobility, savepolaron, plotpolaron
export HellwarthBScheme, HellwarthAScheme
export polaron_memory_function  # Polaron memory functions
export optical_absorption       # Polaron optical absorption
export ϵ_ionic_mode, multi_frohlichalpha, var_params, multi_F, polaron_mobility, polaron_complex_impedence, polaron_complex_conductivity
export frohlichPartial, IRtoDielectric, IRtoalpha, DielectricFromIRmode
export Hellwarth1999mobilityRHS
export make_polaron, save_polaron, load_polaron

##### load in library routines... #####
# stdlib
using LinearAlgebra
using Printf
using JLD
using Tullio, LoopVectorization

#Key numerics used:
# one-dimensional numerical integration in Julia using adaptive Gauss-Kronrod quadrature
import QuadGK.quadgk
# Using the powerful Julia Optim package to optimise the variational parameters
using Optim

# Physical constants
const hbar = const ħ = 1.05457162825e-34;                   # kg m2 / s
const eV = const q = const ElectronVolt = 1.602176487e-19;  # kg m2 / s2
const me = MassElectron = 9.10938188e-31;                   # kg
const Boltzmann = const kB =  1.3806504e-23;                # kg m2 / K s2
const ϵ_0 = 8.854E-12;     # permittivity of free space     # C2N−1m−2,
const amu = 1.660_539_066_60e-27;                           # kg

include("types.jl")            # Polaron types
include("FeynmanTheory.jl")    # Actions + variational functions
include("HellwarthTheory.jl")  # multimode -> equivalent mode.
include("MobilityTheories.jl") # Main polaronmobility function
include("MemoryFunction.jl")   # Memory function X calculation.
include("Susceptibility.jl")   # ImX calculation
include("OedipusRex.jl")       # Optical Absorption
include("MultipleBranches.jl")  # Oct 2019 extension to multiple phonon branches

end # module


text:
"""
    angle_adjacent(start::Position, mid::Position, stop::Position, explementary::Bool=false, mercator::Bool=false)

Find the angle formed by two adjacent segments defined by 3 points. The result will be the (positive clockwise)
angle with origin on the `start-mid` segment, or its explementary angle if required.

# Examples
```jldoctest
julia> p1 = Point([-35, 55])
Point([-35.0, 55.0])

julia> p2 = Point([-34.8, 57.5])
Point([-34.8, 57.5])

julia> p3 = Point([-33.4, 59.1])
Point([-33.4, 59.1])

julia> angle_adjacent(p1, p2, p3, false, false)
202.8279033760424

julia> angle_adjacent(p1, p2, p3, true, false) # explementary angle
157.1720966239576
```
"""
function angle_adjacent(start::Position, mid::Position, stop::Position, explementary::Bool=false, mercator::Bool=false)
    azimuth1 = bearing_to_azimuth((mercator !== true) ? bearing(start, mid, false) : rhumb_bearing(start, mid))
    azimuth2 = bearing_to_azimuth((mercator !== true) ? bearing(stop, mid, false) : rhumb_bearing(stop, mid))

    res = abs(azimuth1 - azimuth2)

    if explementary === true
        return 360 - res
    end

    return res
end

angle_adjacent(start::Point, mid::Point, stop::Point, explementary::Bool=false, mercator::Bool=false) = angle_adjacent(start.coordinates, mid.coordinates, stop.coordinates, explementary, mercator)


text:
using Distributed
include("testenv.jl")
addprocs_with_testenv(2)

using JuliaDB, Test, TextParse, IndexedTables, PooledArrays, Dagger, OnlineStats,
    Statistics, MemPool, Random, Serialization, Dagger, Dates, WeakRefStrings, DataValues

include("test_iteration.jl")
include("test_util.jl")

println("<HEARTBEAT>")

include("test_table.jl")
include("test_query.jl")
include("test_join.jl")
include("test_misc.jl")
include("test_rechunk.jl")
include("test_readwrite.jl")
include("test_onlinestats.jl")
include("test_ml.jl")
include("test_weakrefstrings.jl")


text:
using RobustModels: rho,
                    psi,
                    psider,
                    weight,
                    values,
                    estimator_high_breakdown_point_constant,
                    estimator_high_efficiency_constant,
                    estimator_tau_efficient_constant,
                    isbounded,
                    isconvex,
                    efficiency_tuning_constant,
                    breakdown_point_tuning_constant,
                    tau_efficiency_tuning_constant,
                    nothing  # stopper

@testset "Methods loss functions: $(name)" for name in ("L2", "L1", "Huber", "L1L2", "Fair", "Logcosh", "Arctan", "Cauchy", "Geman", "Welsch", "Tukey", "YohaiZamar")
    typeloss = getproperty(RobustModels, Symbol(name * "Loss"))
    l = typeloss()

    @testset "Methods estimators: $(estimator)" for estimator in (nothing, "M", "S", "MM", "Tau", "GeneralizedQuantile")
        typest = if isnothing(estimator)
            # Check LossFunction methods
            typeloss
        else
            # Check AbstractEstimator methods
            T = getproperty(RobustModels, Symbol(estimator * "Estimator"))
            if estimator in ("S", "MM", "Tau")
                if !in(name, ("Geman", "Welsch", "Tukey", "YohaiZamar"))
                    @test_throws TypeError T{typeloss}
                    continue
                end
            end
            T{typeloss}
        end
        est = typest()
        @test_nowarn println(est)

        if !isnothing(estimator)
            if estimator == "Tau"
#                @test isa(loss(est), Tuple{BoundedLossFunction, BoundedLossFunction})
                @test isa(loss(est), CompositeLossFunction)
                @test typeof(first(loss(est))) == typeloss
                @test typeof(last(loss(est))) == typeloss
            else
                @test typeof(loss(est)) == typeloss
            end
        end

        @testset "Bounded $(estimator) estimators: $(name)" begin
            if name in ("Geman", "Welsch", "Tukey", "YohaiZamar")
                @test isbounded(est)
            else
                @test !isbounded(est)
            end
        end

        @testset "Convex $(estimator) estimators: $(name)" begin
            if name in ("Cauchy", "Geman", "Welsch", "Tukey", "YohaiZamar")
                @test !isconvex(est)
            else
                @test isconvex(est)
            end
        end

        @testset "$(estimator) Estimator values: $(name)" begin
            ρ  = rho(est, 1)
            ψ  = psi(est, 1)
            ψp = psider(est, 1)
            w  = weight(est, 1)

            vals = values(est, 1)
            @test length(vals) == 3
            @test vals[1] ≈ ρ  rtol=1e-6
            @test vals[2] ≈ ψ  rtol=1e-6
            @test vals[3] ≈ w  rtol=1e-6
        end

        # Only for LossFunction
        if isnothing(estimator)
            if !in(name, ("L2", "L1"))
                @testset "Estimator high efficiency: $(name)" begin
                    vopt = estimator_high_efficiency_constant(typest)
                    v = efficiency_tuning_constant(typest; eff=0.95, c0=0.9*vopt)
                    @test isapprox(v, vopt; rtol=1e-3)
                end
            end

            if isbounded(est)
                @testset "Estimator high breakdown point: $(name)" begin
                    vopt = estimator_high_breakdown_point_constant(typest)
                    v = breakdown_point_tuning_constant(typest; bp=0.5, c0=1.1*vopt)
                    @test isapprox(v, vopt; rtol=1e-3)
                end

                @testset "τ-Estimator high efficiency: $(name)" begin
                    vopt = estimator_tau_efficient_constant(typest)
                    v = tau_efficiency_tuning_constant(typest; eff=0.95, c0=1.1*vopt)
                    @test isapprox(v, vopt; rtol=1e-3)
                end
            end
        end

        if typest <: AbstractQuantileEstimator
            @testset "MQuantile: $(name)" begin
                τ = 0.5
                qest1 = GeneralizedQuantileEstimator(l, τ)
                qest2 = GeneralizedQuantileEstimator{typeloss}(τ)
                @test qest1 == qest2
                if name == "L2"
                    qest4 = ExpectileEstimator(τ)
                    @test qest1 == qest4
                elseif name == "L1"
                    qest4 = RobustModels.QuantileEstimator(τ)
                    @test qest1 == qest4
                end

                @testset "Pass through of methods" for fun in (rho, psi, psider, weight)
                    ρ1 = fun(est, 1)
                    ρ2 = fun(qest1, 1)
                    @test ρ1 == ρ2
                end

            end
        end
    end

end


text:
@testset "Mutation operators" begin
    ss = RectSearchSpace([(-1.0, 1.0), (0.0, 100.0), (-5.0, 0.0)])

    @testset "UniformMutation" begin
        gibbs = UniformMutation(ss)

        @test search_space(gibbs) == ss
        # incorrect dimensions
        @test_throws BoundsError BlackBoxOptim.apply(gibbs, 2.0, 0, 1)
        @test_throws BoundsError BlackBoxOptim.apply(gibbs, 2.0, 4, 1)

        for i in 1:numdims(ss)
            dim_range = dimrange(ss, i)
            for _ in 1:NumTestRepetitions
                t = BlackBoxOptim.apply(gibbs, 0.0, i, 1)
                @test t >= dim_range[1]
                @test t <= dim_range[2]
            end
        end
    end

    @testset "MutationClock" begin
        mc = MutationClock(UniformMutation(ss), 0.05)

        mutations_per_param = zeros(numdims(ss))
        NumReps = 2_000
        for i in 1:NumReps
            ref_ind = rand_individual(ss)
            ind = copy(ref_ind)
            BlackBoxOptim.apply!(mc, ind, 1)
            @test in(ind, ss)
            mutations_per_param[ind .!= ref_ind] .+= 1
        end
        mut_frequencies = mutations_per_param ./ NumReps
        for p in 1:numdims(ss)
            @test mut_frequencies[p] > 0.03
            @test mut_frequencies[p] < 0.07 # roughly matches the probability
        end
    end

    @testset "FixedGeneticOperatorsMixture" begin
        mx = FixedGeneticOperatorsMixture(GeneticOperator[NoMutation(), MutationClock(UniformMutation(ss), 0.05)], [0.3, 0.7])

        ref_ind = rand_individual(ss)
        n_params_mutated = 0
        for i in 1:10000
            ind = copy(ref_ind)
            BlackBoxOptim.apply!(mx, ind, 1)
            @test in(ind, ss)
            n_params_mutated += ind != ref_ind
        end
        # the number of parameters changed should roughly match the weight of MutationClock multiplied by its mutation probability
        @test n_params_mutated/numdims(ss) > 200
        @test n_params_mutated/numdims(ss) < 500
    end

    @testset "FAGeneticOperatorsMixture" begin
        mx = FAGeneticOperatorsMixture(GeneticOperator[NoMutation(), MutationClock(UniformMutation(ss), 0.05)], pmin = 0.05, pmax = 20.0)

        ref_ind = rand_individual(ss)
        n_params_mutated = 0
        for i in 1:10000
            ind = copy(ref_ind)
            sel_op, tag = BlackBoxOptim.next(mx)
            BlackBoxOptim.apply!(sel_op, ind, 1)
            @test in(ind, ss)
            is_mutated = ind != ref_ind
            n_params_mutated += is_mutated
            is_improved = rand() < 0.25 && is_mutated
            # fake improvement in 0.25 times mutation clock is applied
            BlackBoxOptim.adjust!(mx, tag, 1, is_improved ? 1.0 : 0.0, 0.0, is_mutated)
        end
        # FA should adjust frequencies to [almost] always apply mutation clock
        # the number of parameters changed should roughly match the MutationClock mutation probability
        @test n_params_mutated/numdims(ss) > 300
        @test n_params_mutated/numdims(ss) < 700
        @test frequencies(mx)[1] < 0.1
        @test frequencies(mx)[2] > 0.9
    end

end


text:
struct MonomerIdealParam <: EoSParam
    Mw::SingleParam{Float64}
end

abstract type MonomerIdealModel <: IdealModel end
@newmodelsimple MonomerIdeal MonomerIdealModel MonomerIdealParam

export MonomerIdeal
function MonomerIdeal(components::Array{String,1}; userlocations::Array{String,1}=String[], verbose=false)
    params = getparams(components, ["properties/molarmass"]; userlocations=userlocations, verbose=verbose)
    Mw = params["Mw"]
    packagedparams = MonomerIdealParam(Mw)
    return MonomerIdeal(packagedparams)
end

function a_ideal(model::MonomerIdealModel, v, T, z)
    x = z/sum(z)
    Mw = model.params.Mw.values
    Λ = @. h/√(k_B*T*Mw/N_A)
    return sum(@. x*log(z*N_A/v*Λ^3))-1.
end


text:
println("run...")

imgfile="../blocks.png"
imgarrayfile="$imgfile.jld"

using Images
img = Images.load(imgfile)
(imgwidth, imgheight) = size(img)

#img = Gray.(img) #grayImage
#imgMatrix = reinterpret.(channelview(img))
#println("write ",typeof(imgMatrix),"->" , size(imgMatrix))

#imga = reinterpret.(vec(channelview(rgbImage)))

#open(f -> print(f, imga), imgarrayfile, "w+")
#writecsv(imgarrayfile, imgMatrix)
#a = readcsv(imgarrayfile)
#println("read ",typeof(a),"->" , size(a))
#a = convert(Array{UInt8},a)
#println("convert ",typeof(a), "->", a[1:100],"\n")

using JLD
#JLD.save(imgarrayfile, "img", imgMatrix)
a = JLD.load(imgarrayfile, "img")
println("read ",typeof(a),"->" , size(a))

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
using ParameterInference
using Test

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


text:
#!/usr/local/bin/julia

# Obligatory worldly salutation
println("o hai world")

#=
Use can use block comments
=#

#=
Maths
=#

# Unary operators
1 + 2 - 4 * 3 / 5 ^ 6 % 7

#=
Bitwise operator
=#

# Not
~true
# And
true & false
# Or
false | true
# Xor
true ⊻ false
# Logical shift
true >>> false
# Arithmetic shift
true >> false
# Logical/arithmetic shift
true << true

#=
Every binary arithmetic and bitwise operator also has an
updating version that assigns the result of the operation
back into its left operand.
=#

# Comparative operatations
1_000 == 1000

# Built-in specials
pi
e
π == pi

# Functions that look like maths
f(x) = x^2 + 2x
f(3)

#=
Strings
=#
typeof('x')
typeof("x")
Int('x')
str = "oh hai wurld"
str[1]
str[4:5]
str[end]

#=
Functions
=#
function f(x, y)
  x + y
end
f(3, 4)

f(x, y) = x + y
f(4, 5)

foo(x::Int) = println("$x is an integer")

@time foo(5)
@time foo(π)
foo(x::String) = println("$x is a string")

#=
Arrays
=#
A = [1, 2, 3, 4, 5]
collect(1:5)
[i for i in 1:5]

letters = collect('a':'z')
letters[1]

1:5 .* 2
length(letters)

[1 2 3; 4 5 6]
[[1, 2, 3] [4, 5, 6]]
Array(1:3)
Array([1 3 4])

eye(5)
[(i,j) for i=1:3 for j=1:i]
reshape(1:16, 4, 4)

eye(Int, 10) .* collect(1:10:100)

#=
Control flow
=#

# Compound expressions
z = begin
  x = 1
  y = 2
  x + y
end

z = (x = 1; y = 2; x + y)

# Conditional expressions
if x < y
    println("x is less than y")
elseif x > y
    println("x is greater than y")
else
    println("x is equal to y")
end

x < y ? "less than" : "not less than"

# While loop
i = 1
while true
  println(i)
  if i >= 5
    break
  end
  i += 1
end

# For loop
for i = 1:10
  if i % 3 != 0
    continue
  end
  println(i)
end

# Try-catch
f(x) =
  try
    sqrt(x)
  catch
    sqrt(complex(x, 0))
end

"""
# Tasks aka Coroutines aka Green Threads

Tasks are a control flow feature that allows computations to be suspended and
resumed in a flexible manner.
"""

"""
# Project Euler

## Problem 9: Special Pythagorean triplet

A Pythagorean triplet is a set of three natural numbers, a < b < c, for which,

```math
a^2 + b^2 = c^2
```

For example, 32 + 42 = 9 + 16 = 25 = 5^2.

There exists exactly one Pythagorean triplet for which a + b + c = 1000.
Find the product abc.
"""

"Generates Pythagorean triplets."
function pythagorean_triplet_producer(ch::Channel)
  m = 2
  while true
    for n = 1:(m-1)
      a = 2*m*n
      b = m^2 - n^2
      c = m^2 + n^2
      put!(ch, [a, b, c])
    end
    m += 1
  end
end

"Solves the 9th Projec Euler programming problem."
function pr009()
  chnl = Channel(pythagorean_triplet_producer)
  while true
    a = take!(chnl)
    if sum(a) == 1_000
      return prod(a)
    end
  end
end

@printf "Time elapsed:\t%4.8f\n" @elapsed pr009()

# EOF

max_stars_repo_name: ellisvalentiner/julia-talk

max_issues_repo_name: ellisvalentiner/julia-talk

max_forks_repo_name: ellisvalentiner/julia-talk


text:
# This file includes the simulation of four node network simulation

import Cairo, Fontconfig
using GraphPlot
using Graphs
using SimpleWeightedGraphs
using LinearAlgebra
using DynamicalNetworks
using DifferentialEquations
using PyPlot
pygui(true)
close("all")

# Helper functions
getwave(x, i, d) = getindex.(x, (i - 1) * d + 1)

# Node settings
nodetype = Lorenz           # Node type
P = diagm([1.0, 0.0, 0])    # State coupling matrix
ηt = 10.0                   # Threshold for maximum Lyapuno exponents

# Topology
n = 5
graph = SimpleWeightedGraph(cycle_graph(n))
add_edge!(graph, 1, 3, 1)
add_edge!(graph, 1, 4, 1)
add_edge!(graph, 2, 5, 5)
gplot(graph, nodelabel = 1:nv(graph), layout = circular_layout)

# Coupling strength for the first symmetry
W1 = laplacian_matrix(graph)
Q1 = [
    1 0 0 0 0
    0 0 0 0 1
    0 0 0 1 0
    0 0 1 0 0
    0 1 0 0 0
]
M1 = eigvecs(Q1)
G1 = inv(M1) * (W1 ./ diag(W1)) * M1
B1 = G1[1:2, 1:2]
D1 = G1[3:5, 3:5]
λt1 = eigvals(B1)
λs1 = eigvals(D1)
λt1[1] > λs1[2] || error("Synchronization is not possible")
ϵr1 = ηt ./ [λt1[1], λs1[2]]

# Coupling strength for the second symmetry
graph = SimpleWeightedGraph(cycle_graph(n))
add_edge!(graph, 1, 3, 1)
add_edge!(graph, 1, 4, 1)
add_edge!(graph, 2, 5, 1)
W2 = laplacian_matrix(graph)
Q2 = [
    1 0 0 0 0
    0 0 0 0 1
    0 0 0 1 0
    0 0 1 0 0
    0 1 0 0 0
]
M2 = eigvecs(Q2)
G2 = inv(M2) * (W2 ./ diag(W2)) * M2
B2 = G2[1:2, 1:2]
D2 = G2[2:5, 2:5]
λt2 = eigvals(B2)
λs2 = eigvals(D2)
# λt2[1] > λs2[2] || error("Synchronization is not possible")
ϵr2 = ηt ./ [λt2[1], λs2[2]]

ϵl = maximum([ϵr1[1], ϵr2[1]])
ϵh = minimum([ϵr2[2], ϵr2[2]])
ϵ = ϵl + abs(ϵh - ϵl) * 0.5

# Callback
nbits = 10
tbit = 100.0
bits = rand(Bool, nbits)
ti = 0.0
tf = nbits * tbit
tswitch = ti:tbit:tf-tbit
srcnode1, dstnode1 = 2, 5     # Control link nodes
high = 5.0          # Control coupling strength
low = 1.0           # Control coupling strength
condition(u, t, integrator) = integrator.t ∈ tswitch
function action!(integrator)
    net, bits = integrator.sol.prob.p
    if popfirst!(bits)
        add_edge!(net.graph, srcnode1, dstnode1, high)
    else
        add_edge!(net.graph, srcnode1, dstnode1, low)
    end
end
callback = DiscreteCallback(condition, action!)

# Network
σ = 0.0  # Noise strength
H = getnoisemat(graph) * σ
net = SDENetwork(nodetype, ϵ, graph, H, P)
sol = solvenet!(net, (ti, tf), (copy(bits),), solverkwargs = (saveat = 0.001,), probkwargs = (callback = callback, tstops = tswitch))
t = sol.t
x = sol.u

# Plots
d = dimension(nodetype)
subplot(221)
plot(t, abs.(getwave(x, 1, d) - getwave(x, 2, d)))
axvline.(tswitch, linestyle = "dashed", color = "red")
xlabel("Time")
ylabel("e12")
title(join(string.(Int.(bits))))
subplot(222)
plot(t, abs.(getwave(x, 2, d) - getwave(x, 5, d)))
axvline.(tswitch, linestyle = "dashed", color = "red")
xlabel("Time")
ylabel("e25")
title(join(string.(Int.(bits))))
subplot(223)
plot(t, getwave(x, 1, d))
xlabel("Time")
ylabel("x1")
subplot(224)
plot(getindex.(x, 1), getindex.(x, 2))
xlabel("x1")
ylabel("x2")
tight_layout()

figure()
plot(λt1, zeros(length(λt1)), marker = "*", label = L"\lambda_t")
plot(λs1, ones(length(λs1)), marker = "o", label = L"\lambda_s")
axvline.([λt1[1], λs1[2]], linestyle = "dashed", color = "red")
axvline(ηt / ϵ, linestyle = "dashed", color = "green")
xlabel(L"\lambda")
title("Eigenvalue distribution")
legend()

figure()
plot(λt2, zeros(length(λt2)), marker = "*", label = L"\lambda_t")
plot(λs2, ones(length(λs2)), marker = "o", label = L"\lambda_s")
axvline.([λt2[1], λs2[2]], linestyle = "dashed", color = "red")
axvline(ηt / ϵ, linestyle = "dashed", color = "green")
xlabel(L"\lambda")
title("Eigenvalue distribution")
legend()


text:
using EnvelopeAnalysis
using Documenter

DocMeta.setdocmeta!(EnvelopeAnalysis, :DocTestSetup, :(using EnvelopeAnalysis); recursive=true)

makedocs(;
    modules=[EnvelopeAnalysis],
    authors="Víctor Manuel Hidalgo",
    repo="https://github.com/VMHidalgo/EnvelopeAnalysis.jl/blob/{commit}{path}#{line}",
    sitename="EnvelopeAnalysis.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://VMHidalgo.github.io/EnvelopeAnalysis.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/VMHidalgo/EnvelopeAnalysis.jl",
)


text:
# Extended Karnaugh Map Representation (EKMR) of matrices

# Matrix Transformation Method for conversion between Traditional Matrix Representation (TMR) and EKMR

struct EKMR{T} where {T}
    M::AbstractArray
end


text:
__precompile__()

module HttpCommon

using Compat
import Compat: UTF8String
import URIParser: URI, unescape

export Headers, Request, Cookie, Response,
       escapeHTML, parsequerystring

# All HTTP status codes, as a Dict of code => description
export STATUS_CODES
include("status.jl")

"""
`Headers` represents the header fields for an HTTP request.
"""
typealias Headers Dict{AbstractString,AbstractString}
headers() = Headers(
    "Server"            => "Julia/$VERSION",
    "Content-Type"      => "text/html; charset=utf-8",
    "Content-Language"  => "en",
    "Date"              => Dates.format(now(Dates.UTC), Dates.RFC1123Format) )

"""
A `Request` represents an HTTP request sent by a client to a server.
It has five fields:

* `method`: an HTTP methods string (e.g. "GET")
* `resource`: the resource requested (e.g. "/hello/world")
* `headers`: see `Headers` above
* `data`: the request data as a vector of bytes
"""
type Request
    method::UTF8String      # HTTP method string (e.g. "GET")
    resource::UTF8String    # Resource requested (e.g. "/hello/world")
    headers::Headers
    data::Vector{UInt8}
    uri::URI
end
Request() = Request("", "", Headers(), UInt8[], URI(""))
Request(method, resource, headers, data) = Request(method, resource, headers, data, URI(""))

Base.show(io::IO, r::Request) = print(io, "Request(", r.uri, ", ",
                                        length(r.headers), " headers, ",
                                        sizeof(r.data), " bytes in body)")

"""
A `Cookie` represents an HTTP cookie. It has three fields:
`name` and `value` are strings, and `attrs` is dictionary
of pairs of strings.
"""
type Cookie
    name::UTF8String
    value::UTF8String
    attrs::Dict{UTF8String, UTF8String}
end
Cookie(name, value) = Cookie(name, value, Dict{UTF8String, UTF8String}())
Base.show(io::IO, c::Cookie) = print(io, "Cookie(", c.name, ", ", c.value,
                                        ", ", length(c.attrs), " attributes)")

"""
A `Response` represents an HTTP response sent to a client by a server.
It has six fields:

* `status`: HTTP status code (see `STATUS_CODES`) [default: `200`]
* `headers`: `Headers` [default: `HttpCommmon.headers()`]
* `cookies`: Dictionary of strings => `Cookie`s
* `data`: the request data as a vector of bytes [default: `UInt8[]`]
* `finished`: `true` if the `Reponse` is valid, meaning that it can be
  converted to an actual HTTP response [default: `false`]
* `requests`: the history of requests that generated the response.
  Can be greater than one if a redirect was involved.

Response has many constructors - use `methods(Response)` for full list.
"""
type Response
    status::Int
    headers::Headers
    cookies::Dict{UTF8String, Cookie}
    data::Vector{UInt8}
    request::Nullable{Request}
    history::Vector{Response}

    # Deprecated fields
    finished::Bool
    requests::Vector{Request}
end
# If a Response is instantiated with all of fields except for `finished`,
# `finished` will default to `false`.
typealias HttpData Union{Vector{UInt8}, AbstractString}
Response(s::Int, h::Headers, d::HttpData) =
  Response(s, h, Dict{UTF8String, Cookie}(), d, Nullable(), Response[], false, Request[])
Response(s::Int, h::Headers)              = Response(s, h, UInt8[])
Response(s::Int, d::HttpData)             = Response(s, headers(), d)
Response(d::HttpData, h::Headers)         = Response(200, h, d)
Response(d::HttpData)                     = Response(d, headers())
Response(s::Int)                          = Response(s, headers(), UInt8[])
Response()                                = Response(200)

function Base.show(io::IO, r::Response)
    print(io, "Response(", r.status, " ", STATUS_CODES[r.status], ", ",
          length(r.headers)," headers, ",
          sizeof(r.data)," bytes in body)")
end

"""
escapeHTML(i::AbstractString)

Returns a string with special HTML characters escaped: &, <, >, ", '
"""
function escapeHTML(i::AbstractString)
    # Refer to http://stackoverflow.com/a/7382028/3822752 for spec. links
    o = replace(i, "&", "&")
    o = replace(o, "\"", """)
    o = replace(o, "'", "'")
    o = replace(o, "<", "<")
    o = replace(o, ">", ">")
    return o
end

"""
parsequerystring(query::AbstractString)

Convert a valid querystring to a Dict:

    q = "foo=bar&baz=%3Ca%20href%3D%27http%3A%2F%2Fwww.hackershool.com%27%3Ehello%20world%21%3C%2Fa%3E"
    parsequerystring(q)
    # Dict{ASCIIString,ASCIIString} with 2 entries:
    #   "baz" => "hello world!"
    #   "foo" => "bar"
"""
function parsequerystring{T<:AbstractString}(query::T)
    q = Dict{AbstractString,AbstractString}()
    length(query) == 0 && return q
    for field in split(query, "&")
        keyval = split(field, "=")
        length(keyval) != 2 && throw(ArgumentError("Field '$field' did not contain an '='."))
        q[unescape(keyval[1])] = unescape(keyval[2])
    end
    q
end

end # module HttpCommon


text:
# ScatteredBlob
immutable ScatteredBlob{NODE<:AbstractNode} <: Mocha.Blob
    nodes::DataStructures.SortedDict{Path,NODE}
end

function Base.show(io::IO, blob::ScatteredBlob)
    n_nodes = length(blob.nodes)
    plural = n_nodes > 1 ? "s" : ""
    print(io, "ScatteredBlob(", n_nodes, " node", plural, ")")
end

function forward!(
        backend,
        blob_out::ScatteredBlob,
        bank::Bank1D,
        blob_in::ScatteredBlob)
    pathdepth(blob_in, bank.behavior.pathkey)
    map(node -> pathdepth(bank.behavior.pathkey, keys(blob_in.nodes)))
    γkey = cons(Literal(:γ, 1), bank.behavior.pathkey)
    for j in bank.behavior.j_range
        for χ in 0:(bank.spec.n_filters_per_octave-1)
            ψ = bank.ψs[1 + θ, 1 + χ, 1 + j]
            for (path_in, node_in) in input.nodes
                path_out = copy(path_in)
                path_out[γkey] = γ
                transform!(blob[path_out], blob_in[path_in], ψ)
            end
        end
    end
end


text:
using Plots
using NumericalIntegration

include(joinpath("..", "src","VlasovSolver","PFC.jl"))
include(joinpath("..", "src","MaxwellSolver","FDTD1D.jl"))
include(joinpath("..", "src","MaxwellSolver","PoissonFourier1D.jl"))

Δx = 0.1*2π
Δt = 0.05*Δx
Δp = 0.1

x_min = -5.0*2π
box_length = 20.0*2π
plasma_thickness = 10.0*2π
plasma_temperature = 0.1
plasma_density = 0.1

laser_amplitude = 1.0
laser_duration = 5*2π

total_time = 2π*10

x = collect(x_min:Δx:(box_length+x_min))
p = collect(-laser_amplitude*4:Δp:laser_amplitude*4)

Nx = length(x)
Np = length(p)

f = plasma_density/sqrt(2π*plasma_temperature) *
    (@. exp(-0.5*(p)^2/plasma_temperature)) *
    (@. 0.5*(tanh(x)-tanh(x-plasma_thickness)))'
    #  (@. (x+0.1)./(x+0.1))'
ni = integrate(p, f)
Ni = integrate(x, ni)

n0 = integrate(p, f)
Ne = integrate(x, n0)

g = copy(f)

advect_x! = PFC.generate_solver(g[1,:], f[1,:]; fₘₐₓ = maximum(f))
advect_p! = PFC.generate_solver(f[:,1], g[:,1]; fₘₐₓ = maximum(f))

em = FDTD1D.YeeMesh1D{Float64}(length(x)-1)

pulse_shape = (
    y = (t,x) -> exp(-((x-t)/laser_duration)^2)*sin(x-t),
    z = (t,x) -> 0.0
)

advance_fields! = FDTD1D.make_advance_fields(em, Δt/Δx, pulse_shape, Δt, Δx, x_min, FDTD1D.PML(0,1.0,Δx,Δt))

t = collect(0.0:Δt:total_time)
n = t * n0'
n[1,:] = n0

ε_e = similar(t)
ε = similar(t)

solve_poisson! = PoissonFourier1D.generate_solver(n0, Δx)

e = similar(x)
solve_poisson!(e, n[1,:]-ni)

ε_e[1] = integrate(x, e.^2)
ε[1] = integrate(x, integrate(p, @. f*p^2)) + ε_e[1]

py = zeros(length(x))
pz = zeros(length(x))

for k in 2:length(t)

    for j = 1:Np
        advect_x!(view(g,j,:), view(f,j,:), p[j]*Δt/Δx)
    end

    n[k,:] = integrate(p, g)
    solve_poisson!(e, n[k,:]-ni)

    py[:] = @. py + em.ey*Δt
    pz[:] = @. pz + em.ez*Δt

    advance_fields!(k*Δt, (y=py.*n[k,:], z=pz.*n[k,:]))

    for i = 1:Nx
        advect_p!(view(f,:,i), view(g,:,i), e[i]*Δt/Δp)
    end

    ε_e[k] = integrate(x, e.^2)
    ε[k] = integrate(x, integrate(p, @. f*p^2)) + ε_e[k]
end

heatmap(x/2π, t/2π, em.ey)


text:
@testset "Brunt-Vaisala" begin

    Γneutral = AvailablePotentialEnergyFramework.g/AvailablePotentialEnergyFramework.Dryair.cp
    Γstable = Γneutral - 10e-3u"K/m"
    Γunstable = Γneutral + 10e-3u"K/m"

    z = 1u"m" .* collect(0:50:10e3)
    tvprofile(Γ,z) = reshape(repeat(300 .- Γ*z,1,2),1,1,length(z),2)
    tvprofile(Γ :: Quantity, z :: Array{ <: Quantity}) = reshape(repeat(300u"K" .- Γ*z,1,2),1,1,length(z),2)

    function plot_N2(Γ,z)
        tv_profile = tvprofile(Γ,z)
        p1 = plot(profile1[1,1,:,1],z,title = "Tv profile, Γ = $Γ")
        p2 = plot(AvailablePotentialEnergyFramework.compute_N2_attempt(tv_profile,z),z, title = "N2")
        plot(p1,p2)
    end

    @test isapprox( 1u"1/s/s" .* zeros(length(z),2),compute_N2(tvprofile(Γneutral,z),z), atol=1e-10u"1/s/s")

    @test all(compute_N2(tvprofile(Γunstable,z),z) .< 0u"1/s/s" )

    @test all(compute_N2(tvprofile(Γstable,z),z) .> 0u"1/s/s" )

    @test isapprox(  zeros(length(z),2),compute_N2(tvprofile(ustrip(Γneutral),ustrip.(z)),ustrip.(z)), atol=1e-10)

    @test all(compute_N2(tvprofile(ustrip(Γunstable),ustrip.(z)),ustrip.(z)) .< 0 )

    @test all(compute_N2(tvprofile(ustrip(Γstable),ustrip.(z)),ustrip.(z)) .> 0 )

end

@testset "Thermodynamics" begin

    @test get_saturation_vapor_pressure(273.15u"K") == 6.112u"hPa"
    @test get_partial_vapor_pressure(0,1000u"hPa") == 0u"hPa"
    @test get_partial_vapor_pressure(1,1000u"hPa") == 1000u"hPa"/(18.016/28.966 + 1.0)
    @test get_mixing_ratio(0u"hPa",1000u"hPa") == 0
    @test get_mixing_ratio(get_partial_vapor_pressure(0.5,1000.0), 1000.0) == 0.5
    @test get_virtual_temperature(300u"K",0,0) == 300u"K"
    @test specific_humidity_to_mixing_ratio(mixing_ratio_to_specific_humidity(0.5)) ≈ 0.5
    @test mixing_ratio_to_specific_humidity(specific_humidity_to_mixing_ratio(0.5)) ≈ 0.5
    @test unit(get_specific_entropy(300u"K",0.2,1000u"hPa"))== u"J/K/kg"
    @test get_potential_temperature(300u"K",1000u"hPa",1000u"hPa") == 300u"K"
    @test get_potential_temperature(300u"K",1010u"hPa",1000u"hPa") < 300u"K"
    @test get_potential_temperature(300u"K",900u"hPa",1000u"hPa") > 300u"K"
    @test get_virtual_temperature(300u"K",0u"g/kg") == 300u"K"
    @test get_virtual_temperature(300u"K",0u"g/kg") == 300u"K"
    @test get_virtual_temperature(300u"K",10u"g/kg") > 300u"K"
    @test 1unit(surface_sensible_heat_flux_to_buoyancy(300u"K", 100u"W/m^2")) == 1unit(g)/u"s"*u"m"
    @test 1unit(surface_latent_heat_flux_to_buoyancy(300u"K", 100u"W/m^2")) == 1unit(g)/u"s"*u"m"
    @test 1unit(get_buoyancy(0.5u"K",10.5u"K")) == 1unit(g)
    @test 1unit(radiative_heating_rate_to_buoyancy(10u"K",5u"K/s")) == 1u"m/s^2/s"
    @test 1unit(get_density_ideal_gas(10.0u"Pa",5.0u"K")) == 1u"kg/m^3"
    # No units

    @test get_saturation_vapor_pressure(273.15) == 6.112
    @test get_partial_vapor_pressure(0,1000) == 0
    @test get_partial_vapor_pressure(1,1000) == 1000/(18.016/28.966 + 1.0)
    @test get_mixing_ratio(0,1000) == 0
    @test get_potential_temperature(300,1000,1000) == 300
    @test get_potential_temperature(300,1010,1000) < 300
    @test get_potential_temperature(300,900,1000) > 300
    @test get_virtual_temperature(300,0) == 300
    @test get_virtual_temperature(300,0) == 300
    @test get_virtual_temperature(300,10) > 300
    @test ustrip(surface_sensible_heat_flux_to_buoyancy(300u"K", 100u"W/m^2")) ==  surface_sensible_heat_flux_to_buoyancy(300, 100)
    @test ustrip(surface_latent_heat_flux_to_buoyancy(300u"K", 100u"W/m^2")) ==  surface_latent_heat_flux_to_buoyancy(300, 100)
    @test ustrip(get_buoyancy(0.5u"K",10.5u"K")) == get_buoyancy(0.5,10.5)
    @test ustrip(radiative_heating_rate_to_buoyancy(10u"K",5u"K/s")) == ustrip(radiative_heating_rate_to_buoyancy(10,5))
    @test ustrip(get_density_ideal_gas(10.0u"Pa",5.0u"K")) == get_density_ideal_gas(10.0,5.0)
    pres = Dataset(joinpath(@__DIR__,"testfiles/thermoprofile.nc")) do ds
        1u"hPa" .* variable(ds, "PRES")[:,:]
    end
    tabs = Dataset(joinpath(@__DIR__,"testfiles/thermoprofile.nc")) do ds
        1u"K" .* variable(ds, "TABS")[:,:]
    end
    qv = 1e-3u"kg/g" .* 1u"g/kg" .* Dataset(joinpath(@__DIR__,"testfiles/thermoprofile.nc")) do ds
        variable(ds, "QV")[:,:] #was originally in g/kg
    end
    @info size(pres) size(qv) size(tabs)
    r = specific_humidity_to_mixing_ratio.(qv)
    timeindex = 1200
    pparcel = pres[1,timeindex]
    tparcel = tabs[1,timeindex]
    rparcel = r[1,timeindex]

    #I will create a similar profile but with a perturbation to see what happens
    tabs_unstable = copy(tabs)
    tabs_unstable[2:40,:] .- 7.0u"K"
    tabs_unstable[41:end,:] .+ 7.0u"K"

    @test_broken unit(get_buoyancy_of_lifted_parcel(tparcel,rparcel,pparcel,tabs[:,timeindex],r[:,timeindex],pres[:,timeindex])[1]) == u"K"
    @test_broken get_buoyancy_of_lifted_parcel(tparcel,rparcel,pparcel,tabs_unstable[:,timeindex],r[:,timeindex],pres[:,timeindex])
end


text:
function watson_res(x, r, p)
    # for watson nex = 13|14|15|16
    for i = 1:29
        ti = Float64(i) / 2.9e1
        r1 = 0.0e0
        r2 = x[1]
        t = 1.0e0
        for j = 2:p
              r1 = r1 + Float64(j-1) * t * x[j]
              t = t * ti
              r2 = r2 + t * x[j]
        end
        r[i] = r1 - r2 * r2 - 1.0e0
    end
    r[30] = x[1]
    r[31] = x[2] - x[1] ^ 2 - 1.0e0
    return r
end

function watson_jac(x, j, p)
    for i = 1:29
        ti = Float64(i) / 2.9e1
        r2 = x[1]
        t = 1.e0
        for k = 2:p
             t = t*ti
             r2 = r2 + t*x[k]
        end
        r2 = -2.e0 * r2
        j[i, 1] = r2
        t = 1.e0
        r2 = ti * r2
        for k = 2:p
            j[i, k] = t * (Float64(k-1) + r2)
            t = t * ti
        end
    end

    for i = 30:31
        for k = 2:p
            j[i, k] = 0.e0
        end
    end
    j[30, 1] = 1.e0
    j[31, 1] = -2.e0*x[1]
    j[31, 2] = 1.e0
    return j
end

const watson6 = let res_init=zeros(31), jac_init=zeros(31,6), x_init=zeros(6)

    res(x, r) = watson_res(x, r, 6)
    jac(x, j) = watson_jac(x, j, 6)

    f(;scale=1, verbose=false, print_steps=false) = testone("watson6", res, jac,
                                                   res_init, jac_init, x_init;
                                                   scale=scale, verbose=verbose,
                                                   print_steps=print_steps)
end

const watson9 = let res_init=zeros(31), jac_init=zeros(31,9), x_init=zeros(9)

    res(x, r) = watson_res(x, r, 9)
    jac(x, j) = watson_jac(x, j, 9)

    f(;scale=1, verbose=false, print_steps=false) = testone("watson9", res, jac,
                                                   res_init, jac_init, x_init;
                                                   scale=scale, verbose=verbose,
                                                   print_steps=print_steps)
end

const watson12 = let res_init=zeros(31), jac_init=zeros(31,12), x_init=zeros(12)

    res(x, r) = watson_res(x, r, 12)
    jac(x, j) = watson_jac(x, j, 12)

    f(;scale=1, verbose=false, print_steps=false) = testone("watson12", res, jac,
                                                   res_init, jac_init, x_init;
                                                   scale=scale, verbose=verbose,
                                                   print_steps=print_steps)
end

# the Fortran driver and test code places a 15 iteration limit on
# watson20 so we will do so as well to match.
const watson20 = let res_init=zeros(31), jac_init=zeros(31,20), x_init=zeros(20)

    res(x, r) = watson_res(x, r, 20)
    jac(x, j) = watson_jac(x, j, 20)

    f(;scale=1, verbose=false, print_steps=false) = testone("watson20", res, jac,
                                                   res_init, jac_init, x_init;
                                                   #mxiter=15, scale=scale, verbose=verbose,
                                                   scale=scale, verbose=verbose,
                                                   print_steps=print_steps)
end


text:
mutable struct BinaryLifeform
    score
    genes
end

"""
    Lifeform(size, onlikelihood, initialscore)

Constructor for a BinaryLifeForm struct. Binary life forms are basically wrappers for a binary vector, which
has a likelihood for being 1(`onlikelihood`). Each life form also has a `score` based on it's "fitness". So
the GA's in this package can be used to minimize or maximize this is an open parameter, but Inf/-Inf is a good
`initialscore`.
"""
Lifeform(size, onlikelihood, initialscore) = BinaryLifeform(initialscore, rbinomial(onlikelihood, size))

"""
    SinglePointCrossOver( L1::BinaryLifeform, L2::BinaryLifeform )

Creates two offspring (new BinaryLifeForms) by mixing the genes from `L1` and `L2` after a random position in the vector.
"""
function SinglePointCrossOver( L1::BinaryLifeform, L2::BinaryLifeform )
    len = length( L1.genes )
    @assert( len .== length( L2.genes ) )
    L = rand( 1 : ( len - 1 ) )
    return ( BinaryLifeform( -Inf, vcat(L1.genes[1:L], L2.genes[ ( L + 1 ) : end]) ),
             BinaryLifeform( -Inf, vcat(L2.genes[1:L], L1.genes[ ( L + 1 ) : end]) )  )
end

"""
    Mutate( L::BinaryLifeform, amount = 0.05 )

Assesses each element in the gene vector (inside of `L`). If a randomly drawn value has a binomial probability
of `amount` the element is mutated.
"""
function Mutate( L::BinaryLifeform, amount = 0.05 )
    len = length(L.genes)
    mutater = rbinomial( amount, len ) .== 1
    L.genes[ mutater ] .= ( 1.0 .- L.genes[ mutater ] )
end


text:
addprocs(3)

@everywhere function approximate_pi(trials::Int64)
  hits = 0
  for i in 1:trials
    if sum(rand(2).^2) < 1
      hits += 1
    end
  end
  hits/trials * 4
end

references = [@spawn approximate_pi(1000000) for i=1:3]
approximation = sum([fetch(ref) for ref in references])/3
@assert 3.1 < approximation < 3.2


text:
export ACEConfig, PositionState, DACEConfig

using ACE, StaticArrays, NamedTupleTools

import Base: *, +, -, zero, rand, randn, show, promote_rule, rtoldefault,
       isapprox, getproperty, real
import LinearAlgebra: norm, promote_leaf_eltypes

abstract type XState{SYMS, TT} <: ACE.AbstractState end

"""
`struct State` the main type for states of input variables (particles).
This type is intended only for storing of information but no arithmetic
should be performed on it. For the latter, we have the DState.
"""
struct State{SYMS, TT} <: XState{SYMS, TT}
   x::NamedTuple{SYMS, TT}

   State{SYMS, TT}(t::NamedTuple{SYMS1, TT1}) where {SYMS, SYMS1, TT, TT1} =
      ( SYMS == SYMS1 ? new{SYMS, TT1}(t)
                      : State{SYMS, TT}( merge( _x(zero(State{SYMS, TT})), t ) ) )

end

"""
`struct DState`: A `State`-like variable but acting like vector with arithmetic
operations defined on it, while `State` acts more like a fixed object that cannot
be manipulated. The main application of `DState` is as a derivative of a
`State`; see also `dstate_type`, ``
"""
struct DState{SYMS, TT} <: XState{SYMS, TT}
   x::NamedTuple{SYMS, TT}

   DState{SYMS, TT}(t::NamedTuple{SYMS1, TT1}) where {SYMS, SYMS1, TT, TT1} =
      ( SYMS == SYMS1 ? new{SYMS, TT1}(t)
                      : DState{SYMS, TT}( merge( _x(zero(State{SYMS, TT})), t ) ) )

end

State(t::NamedTuple{SYMS, TT}) where {SYMS, TT} = State{SYMS, TT}(t)
State{SYMS}(t::NamedTuple{SYMS1, TT}) where {SYMS, SYMS1, TT} = State{SYMS, TT}(t)

State(; kwargs...) = State(NamedTuple(kwargs))
State{SYMS}(; kwargs...) where {SYMS} = State{SYMS}(NamedTuple(kwargs))
State{SYMS, TT}(; kwargs...) where {SYMS, TT} = State{SYMS, TT}(NamedTuple(kwargs))

const CTSTT = Union{AbstractFloat, Complex{<: AbstractFloat},
               SVector{N, <: AbstractFloat},
               SVector{N, <: Complex}} where {N}

dstate_type(X::DState) = typeof(X)

@generated function dstate_type(X::State{SYMS, TT}) where {SYMS, TT}
   syms2 = Symbol[]
   tt2 = DataType[]
   for i = 1:length(SYMS)
      if TT.types[i] <: CTSTT
         push!(syms2, SYMS[i])
         push!(tt2, TT.types[i])
      end
   end
   SYMS2 = tuple(syms2...)
   TT2 = "Tuple{" *
            "$(tuple(tt2...))"[2:end-1] * "}"
   DTX = Meta.parse( "DState{$(SYMS2), $TT2}" )
   quote
      $DTX
   end
end

_mypromrl(T::Type{<: Number}, S::Type{<: Number}) =
      promote_type(T, S)
_mypromrl(T::Type{<: Number}, ::Type{<: SVector{N, P}}) where {N, P} =
      SVector{N, promote_type(T, P)}
_mypromrl(::Type{<: SVector{N, P}}, T::Type{<: Number}) where {N, P} =
      promote_type(T, P)
_mypromrl(T::Type{<: SVector{N, P1}}, ::Type{<: SVector{N, P2}}) where {N, P1, P2} =
      SVector{N, promote_type(P1, P2)}

@generated function dstate_type(x::S, X::ACE.XState{SYMS, TT}
                                ) where {S, SYMS, TT}
   syms2 = Symbol[]
   tt2 = DataType[]
   for i = 1:length(SYMS)
      if TT.types[i] <: CTSTT
         push!(syms2, SYMS[i])
         push!(tt2, _mypromrl(S, TT.types[i]))
      end
   end
   SYMS2 = tuple(syms2...)
   TT2 = "Tuple{" *
            "$(tuple(tt2...))"[2:end-1] * "}"
   DTX = Meta.parse( "DState{$(SYMS2), $TT2}" )
   quote
      $DTX
   end
end

dstate_type(S::Type, X::ACE.State) = dstate_type(zero(S), X)

@generated function _ctssyms(::NamedTuple{SYMS, TT}) where {SYMS, TT}
   syms2 = Symbol[]
   for i = 1:length(SYMS)
      if TT.types[i] <: CTSTT
         push!(syms2, SYMS[i])
      end
   end
   SYMS2 = tuple(syms2...)
   quote
      $SYMS2
   end
end

DState(X::TX) where {TX <: State} =
      (dstate_type(X))( select(_x(X), _ctssyms(_x(X))) )

DState(t::NamedTuple{SYMS, TT}) where {SYMS, TT} = DState{SYMS, TT}(t)

DState(; kwargs...) = DState(NamedTuple(kwargs))

_x(X::XState) = getfield(X, :x)
getproperty(X::XState, sym::Symbol) = getproperty(_x(X), sym)

_myrl(x::Number) = real(x)
_myrl(x::SVector) = real.(x)
Base.real(X::TDX) where {TDX <: DState{SYMS}} where {SYMS} =
      TDX( NamedTuple{SYMS}( ntuple(i -> _myrl(getproperty(X, SYMS[i])), length(SYMS)) ) )

_myim(x::Number) = imag(x)
_myim(x::SVector) = imag.(x)
Base.imag(X::TDX) where {TDX <: DState{SYMS}} where {SYMS} =
      TDX( NamedTuple{SYMS}( ntuple(i -> _myim(getproperty(X, SYMS[i])), length(SYMS)) ) )

_mycplx(x::Number) = complex(x)
_mycplx(x::SVector) = complex.(x)
Base.complex(X::TDX) where {TDX <: DState{SYMS}} where {SYMS} =
      TDX( NamedTuple{SYMS}( ntuple(i -> _mycplx(getproperty(X, SYMS[i])), length(SYMS)) ) )

Base.complex(::Type{TDX}) where {TDX <: DState{SYMS}} where {SYMS} =
      typeof( complex( zero(TDX) ) )

function zero(::Union{TX, Type{TX}}) where {TX <: XState{SYMS, TT}} where {SYMS, TT}
   vals = ntuple(i -> _ace_zero(TT.types[i]), length(SYMS))
   return TX( NamedTuple{SYMS}( vals ) )
end

_ace_zero(args...) = zero(args...)
_ace_zero(::Union{Symbol, Type{Symbol}}) = :O

for f in (:rand, :randn)
   eval( quote
      function $f(::Union{TX, Type{TX}}) where {TX <: XState{SYMS, TT}} where {SYMS, TT}
         vals = ntuple(i -> $f(TT.types[i]), length(SYMS))
         return TX( NamedTuple{SYMS}( vals ) )
      end
   end )
end

const _showdigits = 4
_2str(x) = string(x)
_2str(x::AbstractFloat) = "[$(round(x, digits=_showdigits))]"
_2str(x::Complex) = "[$(round(x, digits=_showdigits))]"
_2str(x::SVector{N, <: AbstractFloat}) where {N} = string(round.(x, digits=_showdigits))
_2str(x::SVector{N, <: Complex}) where {N} = string(round.(x, digits=_showdigits))[11:end]

_showsym(X::State) = ""
_showsym(X::DState) = "′"

show(io::IO, X::XState{SYMS}) where {SYMS} =
      print(io, "{" * prod( "$(sym):$(_2str(getproperty(_x(X), sym))), "
                            for sym in SYMS) * "}" * _showsym(X))

for f in (:+, :-)
   eval( quote
      function $f(X1::TX1, X2::TX2) where {TX1 <: XState{SYMS}, TX2 <: XState{SYMS}} where {SYMS}
         vals = ntuple( i -> $f( getproperty(_x(X1), SYMS[i]),
                                 getproperty(_x(X2), SYMS[i]) ), length(SYMS) )
         return TX1( NamedTuple{SYMS}(vals) )
      end
   end )
end

function -(X::TX) where {TX <: XState{SYMS}} where {SYMS}
      vals = ntuple( i -> -getproperty(_x(X), SYMS[i]) )
      return TX( NamedTuple{SYMS}(vals) )
end

function *(X1::TX, a::Number) where {TX <: XState{SYMS}} where {SYMS}
   vals = ntuple( i -> *( getproperty(_x(X1), SYMS[i]), a ), length(SYMS) )
   return TX( NamedTuple{SYMS}(vals) )
end

function *(a::Number, X1::TX) where {TX <: XState{SYMS}} where {SYMS}
   vals = ntuple( i -> *( getproperty(_x(X1), SYMS[i]), a ), length(SYMS) )
   return TX( NamedTuple{SYMS}(vals) )
end

promote_leaf_eltypes(X::XState{SYMS}) where {SYMS} =
   promote_type( ntuple(i -> promote_leaf_eltypes(getproperty(_x(X), SYMS[i])), length(SYMS))... )

norm(X::XState{SYMS}) where {SYMS} =
      sum( norm( getproperty(_x(X), sym) for sym in SYMS )^2 )

import LinearAlgebra: dot
dot(X1::DState{SYMS}, X2::DState{SYMS}) where {SYMS} =
   sum( dot( getproperty(_x(X1), sym), getproperty(_x(X2), sym) )
        for sym in SYMS )

_contract(X1::DState{SYMS1}, X2::DState{SYMS2}) where {SYMS1, SYMS2} =
   sum( sum( getproperty(_x(X1), sym) .* getproperty(_x(X2), sym) )
              for sym in SYMS1 )

isapprox(X1::TX, X2::TX, args...; kwargs...
         ) where {TX <: XState{SYMS}} where {SYMS} =
   all( isapprox( getproperty(_x(X1), sym), getproperty(_x(X2), sym),
                  args...; kwargs...) for sym in SYMS )

# ----- Implementation of a Position State, as a basic example
PositionState{T} = State{(:rr,), Tuple{SVector{3, T}}}

PositionState(r::AbstractVector{T}) where {T <: AbstractFloat} =
      (@assert length(r) == 3; PositionState{T}(; rr = SVector{3, T}(r)))

promote_rule(::Union{Type{S}, Type{PositionState{S}}},
             ::Type{PositionState{T}}) where {S, T} =
      PositionState{promote_type(S, T)}

# some special functionality for PositionState
*(A::AbstractMatrix, X::TX) where {TX <: PositionState} = TX( (rr = A * X.rr,) )
+(X::TX, u::SVector{3}) where {TX <: PositionState} = TX( (rr = X.rr + u,) )

real(X::PositionState{T}) where {T} =
            PositionState{real(T)}( (rr = real.(X.rr), ) )

# ------------------ Basic Configurations Code

abstract type XACEConfig <: AbstractConfiguration  end

"""
`struct ACEConfig`: The canonical implementation of an `AbstractConfiguration`.
Just wraps a `Vector{<: AbstractState}`
"""
struct ACEConfig{STT} <: XACEConfig
   Xs::Vector{STT}   # list of states
end

struct DACEConfig{STT} <: XACEConfig
   Xs::Vector{STT}   # list of states
end

import Base: *
*(t::Number, dcfg::DACEConfig) = DACEConfig( t * dcfg.Xs )
+(cfg::ACEConfig, dcfg::DACEConfig) = ACEConfig( cfg.Xs + dcfg.Xs )

# --- iterator to go through all states in an abstract configuration assuming
#     that the states are stored in cfg.Xs

Base.iterate(cfg::XACEConfig) =
   length(cfg.Xs) == 0 ? nothing : (cfg.Xs[1], 1)

Base.iterate(cfg::XACEConfig, i::Integer) =
   length(cfg.Xs) == i ? nothing : (cfg.Xs[i+1], i+1)

Base.length(cfg::XACEConfig) = length(cfg.Xs)

Base.eltype(cfg::XACEConfig) = eltype(cfg.Xs)

# ---------------- AD code

# this function makes sure that gradients w.r.t. a State become a DState
function rrule(::typeof(getproperty), X::ACE.XState, sym::Symbol)
   val = getproperty(X, sym)
   return val, w -> ( NoTangent(),
                      dstate_type(w[1], X)( NamedTuple{(sym,)}((w,)) ),
                      NoTangent() )
end


text:
module SurfaceRegistration

using LinearAlgebra: norm, tr, det, ⋅, ×
import HCubature
using Cubature
using PyPlot
using ForwardDiff

include("fourier-basis.jl")
export construct_basis

include("operators.jl")
export l2_squared_norm, l2_squared_distance, l2_inner_product, tangent_vector, unscaled_tangent_vector

include("vectorfields.jl")
export construct_vectorfield, construct_jacobian_matrix
export jacobian_determinant, gradient_jacobian_determinant

include("visual.jl")
export plot_parametric_surface, plot_parametric_wireframe, plot_diffeomorphism, add_point, plot_contourf

include("qmap.jl")
export Qmap, af, Daf, normalize_surface

include("gradient_descent.jl")
export SurfaceRegistrationAlgorithm

include("util.jl")
export check_derivatives

end # module


text:
using MathProgBaseMosek
using MathProgBase
using Base.Test

#const MPB = MathProgBase
include(joinpath(Pkg.dir("MathProgBase","test","linprog.jl")))
include(joinpath(Pkg.dir("MathProgBase","test","mixintprog.jl")))
include(joinpath(Pkg.dir("MathProgBase","test","quadprog.jl")))
include(joinpath(Pkg.dir("MathProgBase","test","nlp.jl")))
include(joinpath(Pkg.dir("MathProgBase","test","conicinterface.jl")))
include(joinpath(Pkg.dir("MathProgBase","test","linproginterface.jl")))

function test_mathprogbase(solver)
    @testset "linprog" begin
        linprogtest(solver)
    end

    @testset "mip" begin
        mixintprogtest(solver)
    end

    @testset "qp" begin
        #quadprogtest(solver)
        #qpdualtest(solver)
        #socptest(solver)
    end

    @testset "nlp" begin
        #nlptest(solver)
        #nlptest_nohessian(solver)
        #convexnlptest(solver)
        #rosenbrocktest(solver)
    end

    @testset "conic" begin
        coniclineartest(solver, duals=true)
        # Test conic fallback for LPs
        #coniclineartest(solver, duals=true)
    end

    @testset "linproginterface" begin
        linprogsolvertest(solver)
        linprogsolvertestextra(solver)
        # Test LP fallback for conics
        linprogsolvertest(solver)
    end

end

test_mathprogbase(MosekSolver(QUIET = true))

if (try Pkg.installed("Convex") != nothing catch false end)
    include("testconvex.jl")
end


text:
###############################################################################
#
# Calculate annual profit
#
###############################################################################
"""
    annual_simulation!(
                node::SPACSimple{FT},
                photo_set::AbstractPhotoModelParaSet{FT},
                weather::DataFrame,
                output::DataFrame
    ) where {FT<:AbstractFloat}

Run annual simulation for a growing season, given
- `node` [`SPACSimple`] type struct
- `photo_set` [`AbstractPhotoModelParaSet`] type struct
- `weather` Weather profile in a growing season
- `output` The predefined output result
"""
function annual_simulation!(
            node::SPACSimple{FT},
            photo_set::AbstractPhotoModelParaSet{FT},
            weather::DataFrame,
            output::DataFrame
) where {FT<:AbstractFloat}
    # 0. unpack required values
    @unpack elevation, gaba, laba, latitude, vtoj = node;

    # 1. update the environmental constants based on the node geographycal info
    ratio            = atmospheric_pressure_ratio(elevation);
    node.envir.p_atm = ratio * P_ATM(FT);
    node.envir.p_O₂  = node.envir.p_atm * FT(0.209);

    # 2. calculate the growing season canopy profit
    gscp   = FT(0);
    for i in eachindex( (weather).Year )
        # 2.1 read and update the hourly data
        day  ::FT = (weather).Day[i]
        hour ::FT = (weather).Hour[i]
        _tair::FT = (weather).Tair[i]
        _dair::FT = (weather).D[i]
        p_co2::FT = (weather).CO2[i] * ratio
        r_all::FT = (weather).Solar[i]
        wind ::FT = (weather).Wind[i]
        rain ::FT = (weather).Rain[i]

        node.envir.t_air = _tair + K_0(FT);
        node.envir.p_sat = saturation_vapor_pressure( node.envir.t_air );
        node.envir.p_a   = p_co2;
        node.envir.vpd   = _dair * 1000;
        node.envir.p_H₂O = node.envir.p_sat - node.envir.vpd;
        node.envir.RH    = node.envir.p_H₂O / node.envir.p_sat;
        node.envir.wind  = wind;

        # 2.2 if day time
        zenith = zenith_angle(latitude, day, hour, FT(0));
        if (r_all>0) & (zenith<=85)
            # 2.2.1 update the leaf partitioning
            big_leaf_partition!(node, zenith, r_all)
            @unpack frac_sh, frac_sl = node.container2L;

            # 2.2.2 optimize flows in each layer
            optimize_flows!(node, photo_set);
            leaf_gas_exchange!(node, photo_set, node.opt_f_sl, node.opt_f_sh);
            flow = node.opt_f_sl + node.opt_f_sh;
            anet = frac_sl * node.container2L.cont_sl.an +
                   frac_sh * node.container2L.cont_sh.an;

            # 2.2.3 update drought history
            pressure_profile!(node.hs, node.p_soil, node.opt_f_sl,
                              node.opt_f_sh, frac_sl);

            # 2.2.4 pass values to DataFrame
            output[i, "LAI_sl"] = (node.container2L).lai_sl
            output[i, "PAR_sl"] = (node.container2L).par_sl
            output[i, "RAD_sl"] = (node.container2L).rad_sl
            output[i, "E_sl"  ] = (node.container2L).cont_sl.e
            output[i, "P_sl"  ] = (node.container2L).cont_sl.p
            output[i, "An_sl" ] = (node.container2L).cont_sl.an
            output[i, "Ag_sl" ] = (node.container2L).cont_sl.ag
            output[i, "C_sl"  ] = (node.container2L).cont_sl.c
            output[i, "G_sl"  ] = (node.container2L).cont_sl.gh
            output[i, "T_sl"  ] = (node.container2L).cont_sl.t

            output[i, "LAI_sh"] = (node.container2L).lai_sh
            output[i, "PAR_sh"] = (node.container2L).par_sh
            output[i, "RAD_sh"] = (node.container2L).rad_sh
            output[i, "E_sh"  ] = (node.container2L).cont_sh.e
            output[i, "P_sh"  ] = (node.container2L).cont_sh.p
            output[i, "An_sh" ] = (node.container2L).cont_sh.an
            output[i, "Ag_sh" ] = (node.container2L).cont_sh.ag
            output[i, "C_sh"  ] = (node.container2L).cont_sh.c
            output[i, "G_sh"  ] = (node.container2L).cont_sh.gh
            output[i, "T_sh"  ] = (node.container2L).cont_sh.t

        # 2.3 if night time
        else
            # 2.3.1 calculate leaf temperature
            tlef = max(200, leaf_temperature(node, r_all, FT(0)));

            # 2.3.2 calculate the gas exchange rates
            node.ps.T = tlef;
            leaf_rd!(photo_set.ReT, node.ps);
            anet = -node.ps.Rd;

            # 2.3.3 update temperature effects and then leaf water potential
            flow = FT(0);
            # TODO morea reasonable functions need to be added
            # vc_temperature_effects!(node.hs.leaf, tlef);
            plef = end_pressure(node.hs, flow);

            # 2.3.4 pass values to DataFrame
            output[i, "LAI_sl"] = FT(0)
            output[i, "PAR_sl"] = FT(0)
            output[i, "RAD_sl"] = FT(0)
            output[i, "E_sl"  ] = FT(0)
            output[i, "P_sl"  ] = plef
            output[i, "An_sl" ] = anet
            output[i, "Ag_sl" ] = FT(0)
            output[i, "C_sl"  ] = p_co2
            output[i, "G_sl"  ] = FT(0)
            output[i, "T_sl"  ] = tlef

            output[i, "LAI_sh"] = node.lai
            output[i, "PAR_sh"] = FT(0)
            output[i, "RAD_sh"] = FT(0)
            output[i, "E_sh"  ] = FT(0)
            output[i, "P_sh"  ] = plef
            output[i, "An_sh" ] = anet
            output[i, "Ag_sh" ] = FT(0)
            output[i, "C_sh"  ] = p_co2
            output[i, "G_sh"  ] = FT(0)
            output[i, "T_sh"  ] = tlef
        end

        # pass more data to DataFrame
        output[i, "T_air" ] = _tair
        output[i, "D_air" ] = _dair
        output[i, "Wind"  ] = wind
        output[i, "Rain"  ] = rain
        output[i, "Ca"    ] = p_co2
        output[i, "SWC"   ] = node.swc
        output[i, "P_soil"] = node.p_soil
        output[i, "H_sun" ] = 90 - zenith
        output[i, "A_net" ] = anet
        output[i, "E_crit"] = node.ec

        # 2.4 update soil moisture by converting flow to Kg per hour
        flow /= KG_H_2_MOL_S(FT);
        rain *= gaba * ρ_H₂O(FT) / 1000;
        soil_moisture!(node, flow-rain);
    end

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
module TestModel

using Distributions
using TimeData
using DataArrays
using DataFrames
using Dates
using Base.Test

import Econometrics

## reload(joinpath(Pkg.dir(), "Econometrics", "src", "Econometrics.jl"))

## create dummy Timenum data with missing values
mod = Econometrics.TlsIID(3.4, -0.8, 1.2)
nSim = 100
da = DataArray(Float64, nSim + 10)
da[11:end] = Econometrics.simulate(mod, nSim)
df = DataFrame(OrigData = da)
dats = Date[Date(2001,1,1)+Dates.Day(ii) for ii=1:nSim+10]
tn = Timenum(df, dats)

## test resimulate with given model
##---------------------------------

dataSim = Econometrics.resimulate(tn, mod)

## test equal missing values
@test all(isna(tn.vals[1]) .== isna(dataSim.vals[1]))

## values should come from similar distribution
mod_fit = Econometrics.estimate(Econometrics.TlsIID,
                                dropna(dataSim.vals[1]))

## compare
[[Econometrics.getParams(mod)...]'
[Econometrics.getParams(mod_fit)...]']

## test resimulate with model fitting
##-----------------------------------

dataSim = Econometrics.resimulate(tn, Econometrics.TlsIID)
dataSim = Econometrics.resimulate(tn, Econometrics.NormIID)

end


text:
n = 15_000
R =  1_000
k =    400
P = n - R

y = randn(n)
x = [ones(n) randn(n, k-1)]

@time b1 = recursive_ols(y, x, R)


text:
# This file is a part of UnsafeArrays.jl, licensed under the MIT License (MIT).

using UnsafeArrays
using Test

using Random

@testset "unsafe_array" begin
    function rand_array(::Type{T}, Val_N::Val{N}) where {T, N}
        sz_max = (8, 7, 5, 4, 5)
        sz = ntuple(i -> sz_max[i], Val_N)
        rand!(Array{T}(undef, sz...))
    end

    function rand_array(::Type{T}, Val_N::Val{N}) where {T<:Integer, N}
        sz_max = (8, 7, 5, 4, 5)
        sz = ntuple(i -> sz_max[i], Val_N)
        rand!(Array{T}(undef, sz...), 0:99)
    end

    function test_A(test_code::Function, ::Type{T}, Val_N::Val{N}) where {T, N}
        A = rand_array(T, Val_N)
        UnsafeArrays.@gc_preserve A test_code(A)
    end

    function test_A_UA(test_code::Function, ::Type{T}, Val_N::Val{N}) where {T, N}
        A = rand_array(T, Val_N)
        UA = UnsafeArray(pointer(A), size(A))
        UnsafeArrays.@gc_preserve A test_code(A, UA)
    end

    @testset "ctors" begin
        test_A(Float64, Val(0)) do A
            ptr = pointer(A)
            sz = size(A)
            @test typeof(@inferred(UnsafeArray{Float64,0}(Val(true), ptr, sz))) == UnsafeArray{Float64,0}
            @test typeof(@inferred(UnsafeArray{Float64,0}(ptr, sz))) == UnsafeArray{Float64,0}
            @test typeof(@inferred(UnsafeArray(ptr, sz))) == UnsafeArray{Float64,0}
            @test_throws MethodError typeof(UnsafeArray{Float64,0}(Val(false), ptr, sz))
        end

        test_A(Int, Val(2)) do A
            ptr = pointer(A)
            sz = size(A)
            @test typeof(@inferred(UnsafeArray{Int,2}(Val(true), ptr, sz))) == UnsafeArray{Int,2}
            @test typeof(@inferred(UnsafeArray{Int,2}(ptr, sz))) == UnsafeArray{Int,2}
            @test typeof(@inferred(UnsafeArray(ptr, sz))) == UnsafeArray{Int,2}
            @test_throws MethodError typeof(UnsafeArray{Int,2}(Val(false), ptr, sz))
        end
    end

    @testset "size, length, index style, etc." begin
        run_test(T::Type, Val_N::Val) = test_A_UA(T, Val_N) do A, UA
            T = eltype(UA)
            @test @inferred(length(UA)) == length(A)
            @test @inferred(IndexStyle(UA)) == IndexLinear()
            @test @inferred(LinearIndices(UA)) == LinearIndices(A)
            @test @inferred(eachindex(UA)) == eachindex(A)

            @test @inferred(Base.unsafe_convert(Ptr{T}, UA)) == Base.unsafe_convert(Ptr{T}, A)
            @test @inferred(pointer(UA)) == pointer(A)

            @test @inferred(Base.iscontiguous(UA)) == true
            @test @inferred(Base.iscontiguous(typeof(UA))) == true
        end

        run_test(Int, Val(0))
        run_test(Float64, Val(2))
    end

    @testset "equality, getindex and setindex!" begin
        test_A_UA(Float32, Val(2)) do A, UA
            test_A_UA(Float32, Val(2)) do B, UB
                @test UA != UB
                A .= B
                @test UA == UB
            end
        end

        test_A_UA(Float32, Val(2)) do A, UA
            B = rand_array(Float32, Val(2))
            rand!(B)
            for i in eachindex(UA, B)
                UA[i] = B[i]
            end
            @test A == B
            @test UA == B
            @test all(i -> UA[i] == B[i], eachindex(UA, B))
            @test all(i -> UA[i] == B[i], CartesianIndices(size(UA)))
        end

        test_A_UA(Int, Val(1)) do A, UA
            @test all(x -> x == 42, @inferred fill!(UA, 42))
        end
    end

    @testset "view" begin
        test_A_UA(Int16, Val(0)) do A, UA
            @test view(UA) === UA
            @test view(UA) == view(A)
        end

        test_A_UA(Int, Val(1)) do A, UA
            @test typeof(@inferred(view(UA, :))) == UnsafeArray{Int,1}
            @test typeof(@inferred(view(UA, 2:4))) == UnsafeArray{Int,1}
            @test typeof(@inferred(view(UA, 2))) == UnsafeArray{Int,0}

            @test view(UA, :) == view(A, :)
            @test view(UA, 2:4) == view(A, 2:4)
            @test view(UA, 2) == view(A, 2)

            @test isbits(view(UA, :)) == true
            @test isbits(view(UA, 2:4)) == true
            @test isbits(view(UA, 2)) == true
        end

        test_A_UA(Float32, Val(3)) do A, UA
            @test typeof(@inferred(view(UA, :))) == UnsafeArray{Float32,1}
            @test typeof(@inferred(view(UA, :, :, :))) == UnsafeArray{Float32,3}
            # TODO: Type inference fails for some reason:
            # @test typeof(@inferred(view(UA, :, 3, :))) <: SubArray
            @test typeof((view(UA, :, 3, :))) <: SubArray
            @test typeof(@inferred(view(UA, :, :, 3))) == UnsafeArray{Float32,2}
            @test typeof(@inferred(view(UA, :, 2:4, 3))) == UnsafeArray{Float32,2}
            @test typeof(@inferred(view(UA, :, 2, 3))) == UnsafeArray{Float32,1}
            @test typeof(@inferred(view(UA, 2, 2, 3))) == UnsafeArray{Float32,0}
            @test_throws BoundsError view(UA, 2, 2:20, 3)

            @test view(UA, :) == view(A, :)
            @test view(UA, :, :, :) == view(A, :, :, :)
            @test view(UA, :, 3, :) == view(A, :, 3, :)
            @test view(UA, :, :, 3) == view(A, :, :, 3)
            @test view(UA, :, 2:4, 3) == view(A, :, 2:4, 3)
            @test view(UA, :, 2, 3) == view(A, :, 2, 3)
            @test view(UA, 2, 2, 3) == view(A, 2, 2, 3)

            @test isbits(view(UA, :)) == true
            @test isbits(view(UA, :, :, :)) == true
            @test isbits(view(UA, :, 3, :)) == true
        end
    end

    @testset "reshape" begin
        test_A_UA(Float32, Val(3)) do A, UA
            @test typeof(@inferred(reshape(UA, size(UA, 1), size(UA, 2) * size(UA, 3)))) == UnsafeArray{Float32, 2}
            @test typeof(@inferred(reshape(UA, size(UA, 1) * size(UA, 2) * size(UA, 3)))) == UnsafeArray{Float32, 1}

            @test reshape(UA, size(UA, 1), size(UA, 2) * size(UA, 3)) == reshape(A, size(A, 1), size(A, 2) * size(A, 3))
            @test reshape(UA, size(UA, 1) * size(UA, 2) * size(UA, 3)) == reshape(A, size(A, 1) * size(A, 2) * size(A, 3))
        end

        test_A_UA(Int, Val(1)) do A, UA
            @test typeof(@inferred(reshape(UA, 2, div(size(UA, 1), 2)))) == UnsafeArray{Int, 2}
            @test reshape(UA, 2, div(size(UA, 1), 2)) == reshape(A, 2, div(size(A, 1), 2))
        end
    end

    @testset "copyto! and conversion" begin
        test_A_UA(Int32, Val(3)) do A, UA
            B = similar(A, Int32)
            @test B === @inferred(copyto!(B, UA))
            @test B == A
        end

        test_A_UA(Int32, Val(3)) do A, UA
            B = similar(A, Int64)
            @test B === @inferred(copyto!(B, UA))
            @test B == A
        end

        test_A_UA(Int32, Val(3)) do A, UA
            B = zeros(Int32, size(A)...)
            @test B === @inferred(copyto!(B, 3, UA, 5, 7))
            C = zeros(Int32, size(A)...)
            copyto!(C, 3, A, 5, 7)
            @test B == C
        end

        test_A_UA(Int32, Val(3)) do A, UA
            B = zeros(Int64, size(A)...)
            @test B === @inferred(copyto!(B, 3, UA, 5, 7))
            C = zeros(Int64, size(A)...)
            copyto!(C, 3, A, 5, 7)
            @test B == C
        end

        test_A_UA(Int32, Val(1)) do A, UA
            B = zeros(Int32, size(A)...)
            @test B === @inferred(copyto!(B, 2, UA, 3, 4))
            C = zeros(Int32, size(A)...)
            copyto!(C, 2, A, 3, 4)
            @test B == C
            @test_throws BoundsError copyto!(B, 3, UA, 5, 7)
        end

        test_A_UA(Int32, Val(1)) do A, UA
            B = zeros(Int64, size(A)...)
            @test B === @inferred(copyto!(B, 2, UA, 3, 4))
            C = zeros(Int32, size(A)...)
            copyto!(C, 2, A, 3, 4)
            @test B == C
            @test_throws BoundsError copyto!(B, 3, UA, 5, 7)
        end

        test_A_UA(Float64, Val(3)) do A, UA
            B = rand(Float64, size(A)...)
            A2 = deepcopy(A)
            copyto!(A2, B)
            @test UA === @inferred(copyto!(UA, B))
            @test UA == A2
        end

        test_A_UA(Float64, Val(3)) do A, UA
            B = rand(Float32, size(A)...)
            A2 = deepcopy(A)
            copyto!(A2, B)
            @test UA === @inferred(copyto!(UA, B))
            @test UA == A2
        end

        test_A_UA(Float64, Val(3)) do A, UA
            B = rand(Float64, size(A)...)
            A2 = deepcopy(A)
            copyto!(A2, 3, B, 5, 7)
            @test UA === @inferred(copyto!(UA, 3, B, 5, 7))
            @test UA == A2
        end

        test_A_UA(Float64, Val(3)) do A, UA
            B = rand(Float32, size(A)...)
            A2 = deepcopy(A)
            copyto!(A2, 3, B, 5, 7)
            @test UA === @inferred(copyto!(UA, 3, B, 5, 7))
            @test UA == A2
        end

        test_A_UA(Float64, Val(1)) do A, UA
            B = rand(Float64, size(A)...)
            A2 = deepcopy(A)
            copyto!(A2, 2, B, 3, 4)
            @test UA === @inferred(copyto!(UA, 2, B, 3, 4))
            @test UA == A2
            @test_throws BoundsError copyto!(UA, 3, B, 5, 7)
        end

        test_A_UA(Float64, Val(1)) do A, UA
            B = rand(Float32, size(A)...)
            A2 = deepcopy(A)
            copyto!(A2, 2, B, 3, 4)
            @test UA === @inferred(copyto!(UA, 2, B, 3, 4))
            @test UA == A2
            @test_throws BoundsError copyto!(UA, 3, B, 5, 7)
       end
    end

    @testset "deepcopy" begin
        test_A_UA(Float32, Val(3)) do A, UA
            @test typeof(@inferred(deepcopy(UA))) == typeof(A)
            @test deepcopy(UA) == A
        end

        test_A_UA(Int16, Val(1)) do A, UA
            @test typeof(@inferred(deepcopy(UA))) == typeof(A)
            @test deepcopy(UA) == A
        end
    end

    # # Disabled, as specialization of Base.unaliascopy is disabled:
    #
    # @testset "unaliascopy" begin
    #     test_A_UA(Float32, Val(3)) do A, UA
    #         @test typeof(@inferred(Base.unaliascopy(UA))) == typeof(A)
    #         @test Base.unaliascopy(UA) == A
    #     end
    #
    #     test_A_UA(Int16, Val(1)) do A, UA
    #         @test typeof(@inferred(Base.unaliascopy(UA))) == typeof(A)
    #         @test Base.unaliascopy(UA) == A
    #     end
    # end

    @testset "conversion" begin
        test_A_UA(Float32, Val(3)) do A, UA
            @test typeof(@inferred(convert(Array, UA))) == typeof(A)
            @test convert(Array, UA) == A
        end

        test_A_UA(Int16, Val(1)) do A, UA
            @test typeof(@inferred(convert(Array, UA))) == typeof(A)
            @test convert(Array, UA) == A
        end
    end
end


text:
struct CtbTree
    label::String
    adj::Vector{CtbTree}
    # prob::Float64
end

isleaf(tree::CtbTree) = length(tree.adj) == 0

function dfstraverse(tree::CtbTree, visitor::Function)
    for chd in tree.adj
        dfstraverse(chd, visitor)
    end
    visitor(tree)
    return nothing
end

function label(s::String) #todo special cases
    t = findfirst(!isletter, s)
    return (t == nothing || t == 1) ? s : s[1:t - 1]
end

function CtbTree(chars::Vector{Char}; l = findfirst(isequal('('), chars), trim = label)
    nchar = length(chars)
    l += 1
    r = l
    while chars[r] != '(' && chars[r] != ')' r += 1 end
    if (chars[r] == ')')
        ss = split(join(chars[l:r - 1]))
        @assert length(ss) == 2
        leaf = CtbTree(String(ss[2]), CtbTree[])
        posn = CtbTree(trim(String(ss[1])), CtbTree[leaf])
        return posn
    else
        @assert chars[r] == '('
        ss = split(join(chars[l:r - 1]))
        # @assert length(ss) == 1 string(l, "  ", chars)
        if length(ss) == 0
            return CtbTree(chars; l = r)
        end
        @assert length(ss) == 1
        fa = CtbTree(trim(String(ss[1])), CtbTree[])
        # setlabel!(obj, cur, ss[1])
        nlb = 1
        while nlb > 0
            if (chars[r] == '(')
                if nlb == 1 l = r end
                nlb += 1
            elseif chars[r] == ')'
                nlb -= 1
                if nlb == 1
                    push!(fa.adj, CtbTree(chars; l = l, trim = trim))
                end
            end
            r += 1
        end
        return fa
    end
end

CtbTree(ct::String; trim = label) = CtbTree(collect(ct); trim = trim)

function ==(ta::CtbTree, tb::CtbTree)
        ta.label == tb.label && length(ta.adj) == length(tb.adj) && all(i -> ta.adj[i] == tb.adj[i], 1:length(ta.adj))
end

function size(tree::CtbTree)
    if isleaf(tree) return (1, 1)
    else
        h = w = 0
        for chd in tree.adj
            ch, cw = size(chd)
            h += ch; w = max(w, cw + 1)
        end
        return (h, w)
    end
end

function display(obj::CtbTree)
    xlim, ylim = size(obj)
    mat = fill("", xlim, ylim)
    ileaf = 1
    colors = Dict{String, Int}()
    edges = Vector{Tuple{Int, Int, Int, Int}}()
    function dfs(cur::CtbTree)
        if !haskey(colors, cur.label)
            colors[cur.label] = (1 + (length(colors) + 1) * 10) % 256
        end
        if isleaf(cur)
            mat[ileaf, ylim] = cur.label
            ileaf += 1
            return (ileaf - 1, ylim)
        else
            cxys = Vector{Tuple{Int, Int}}()
            x = xlim; y = ylim;
            for chd in cur.adj
                cx, cy = dfs(chd)
                x = min(x, cx)
                y = min(y, cy - 1)
                push!(cxys, (cx, cy))
            end
            for cxy in cxys push!(edges, (x, y, cxy[1], cxy[2])) end
            mat[x, y] = cur.label
            return (x, y)
        end
    end
    dfs(obj);
    for (lx, ly, rx, ry) in edges
        for x = lx + 1:rx - 1
            new = old = mat[x, ly]
            if old == "" new = "│" end
            if old == "└" new = "├" end
            mat[x, ly] = new
        end
        if lx < rx mat[rx, ly] = "└" end
        mat[rx, ly + 1:ry - 1] .= "─";
    end
    yw = mapreduce(length, max, mat; dims = 1) .+ 2
    for x = 1:xlim, y = 1:ylim
        s = mat[x, y]; ns = length(s)
        # words
        if y == ylim println(" ", s); continue end
        padl = div(yw[y] - ns, 2)
        if (y == ylim - 1)
            padl = yw[y] - ns
        end
        padr = yw[y] - ns - padl
        padlc = padrc = '─'
        if s == "" || s == "│"
            padlc = padrc = ' '
        elseif s == "└" || s == "├" || x == 1 && y == 1
            padlc = ' '
        end
        ns = padlc ^ padl * s * padrc ^ padr
        for i = 1:padl print(padlc) end
        if haskey(colors, s) printstyled(s; color = colors[s]) else print(s) end
        for i = 1:padr print(padrc) end
    end
    return nothing
end

# ctb = parsectb("data/ctb8.0")

# using DataFrames
# import Base.stat
# function stat(ctb::CtbTreeBank)
#     from = Vector{String}(); to = Vector{Vector{String}}()
#     from2 = Vector{String}(); to2 = Vector{String}()
#     function visitor(tree, cur)
#         nchds = length(tree.adj[cur])
#         if nchds > 0
#             if nchds == 1
#                 chd = tree.adj[cur][1]
#                 if length(tree.adj[chd]) == 0
#                     if tree.label[cur] != "-NONE-"
#                         push!(from2, tree.label[cur])
#                         push!(to2, tree.label[chd])
#                     end
#                 elseif tree.label[chd] != "-NONE-"
#                     push!(from, tree.label[cur])
#                     push!(to, tree.label[tree.adj[cur]])
#                 end
#             else
#                 push!(from, tree.label[cur])
#                 push!(to, tree.label[tree.adj[cur]])
#             end
#         end
#     end
#     for vec in ctb
#         for tree in vec
#             dfstraverse(tree, visitor)
#         end
#     end
#     function f(df::DataFrame)
#         return by(df, [:from, :to], tot = :from => length, sort = true)
#     end
#     inn = f(DataFrame(from = from, to = to))
#     pos = f(DataFrame(from = from2, to = to2))
#     inn, pos
# end

# CtbTreeNode = Union{InnerTreeNode, PosTreeNode, LeafTreeNode}
# CtbTreeNode(label::String, id::String) == InnerTreeNode(label, id, Int[], 0.0)

function cnf(root::CtbTree)::CtbTree
    nchd = length(root.adj)
    if nchd == 0
        return root
    elseif nchd <= 2
        newroot = CtbTree(root.label, CtbTree[])
        for i = 1:nchd push!(newroot.adj, cnf(root.adj[i])) end
        return newroot
    else
        newroot = CtbTree(root.label, CtbTree[])
        newright = CtbTree(join(map(x -> x.label, root.adj[2:end]), "+"), root.adj[2:end])
        push!(newroot.adj, cnf(root.adj[1]))
        push!(newroot.adj, cnf(newright))
        return newroot
    end
end

function decnf(root::CtbTree)
    nodes = decnf2(root)
    return length(nodes) == 1 ? nodes[1] : CtbTree(join(map(x -> x.label, nodes), "+"), nodes)
end

function decnf2(root::CtbTree)::Vector{CtbTree}
    nchd = length(root.adj)
    if nchd == 0 return CtbTree[root] end
    if nchd == 1
        newroot = CtbTree(root.label, decnf2(root.adj[1]))
        return CtbTree[newroot]
    end
    newroot = CtbTree(root.label, append!(decnf2(root.adj[1]), decnf2(root.adj[2])))
    istmp = in('+', root.label)
    return istmp ? newroot.adj : CtbTree[newroot]
end


text:
# This provides the functionality of periodic checkpointing. It uses the
# terminology of Revolve; see Alg. 799 published by Griewank et al.
# A minor extension is the  optional `bundle` parameter that allows to treat as many loop
# iterations in one tape/adjoint sweep. If `bundle` is 1, the default, then the behavior is that of Alg. 799.

"""
    Periodic

Periodic checkpointing scheme.

"""
mutable struct Periodic <: Scheme
    steps::Int
    acp::Int
    period::Int
    verbose::Int
    fstore::Union{Function,Nothing}
    frestore::Union{Function,Nothing}
end

function Periodic(
    steps::Int,
    checkpoints::Int,
    fstore::Union{Function,Nothing} = nothing,
    frestore::Union{Function,Nothing} = nothing;
    anActionInstance::Union{Nothing,Action} = nothing,
    bundle_::Union{Nothing,Int} = nothing,
    verbose::Int = 0
)
    if !isa(anActionInstance, Nothing)
        # same as default init above
        anActionInstance.actionflag = 0
        anActionInstance.iteration  = 0
        anActionInstance.cpNum      = 0
    end
    acp             = checkpoints
    period          = div(steps, checkpoints)

    periodic = Periodic(steps, acp, period, verbose, fstore, frestore)

    forwardcount(periodic)
    return periodic
end

function forwardcount(periodic::Periodic)
    if periodic.acp < 0
        error("Periodic forwardcount: error: checkpoints < 0")
    elseif periodic.steps < 1
        error("Periodic forwardcount: error: steps < 1")
    elseif mod(periodic.steps, periodic.acp) != 0
        error("Periodic forwardcount: error: steps ", periodic.steps, " not divisible by checkpoints ", periodic.acp)
    end
end

function checkpoint_struct(body::Function,
        alg::Periodic,
        model_input::MT,
        shadowmodel::MT
    ) where{MT}
    model = deepcopy(model_input)
    model_final = []
    model_check_outer = Array{MT}(undef, alg.acp)
    model_check_inner = Array{MT}(undef, alg.period)
    check = 0
    for i = 1:alg.acp
        model_check_outer[i] = deepcopy(model)
        for j= (i-1)*alg.period: (i)*alg.period-1
            body(model)
        end
    end
    model_final = deepcopy(model)
    for i = alg.acp:-1:1
        model = deepcopy(model_check_outer[i])
        for j= 1:alg.period
            model_check_inner[j] = deepcopy(model)
            body(model)
        end
        for j= alg.period:-1:1
            model = deepcopy(model_check_inner[j])
            Enzyme.autodiff(body, Duplicated(model,shadowmodel))
        end
    end
    return model_final
end