text:
####
# Moste complex code , can probebly be improved
####

abstract type AbstractDesign{T} end

""""
    SingleDesign{T} <: AbstractDesign{T}

A Single design is a structure that contains one special define type of construct.
They can be combined to form more complex design space.
"""
abstract type SingleDesign{T} <: AbstractDesign{T} end

function Base.show(io::IO,m::SingleDesign)
        print(io,"Used modules : ")
        print(io,m.mod ,"\n")
        print(io,"allowed length : ", m.len, "\n",
               "constraints : " , m.con ,"\n")
        print(io,"\t \t designspace \n")
        print(io,m.space)
end

"""
    OrderedDesign

Ordered space is where the possition in the construct has a functional role, the
construct [a,b,c] is different compeared to [c,b,a].
In the General case every module can be used infinite number of times.

Fields:
    - `mod::AbstractMod`: The used modules in this design
    - `len::Int` : The allowed length of the construct
    - `con::AbstractConstraints`: The constraints of the design space
    - `space::AbstractSpace`: Design space containing all constructs, if
                    possible, these constructs are not explicitly generated.
"""
struct OrderedDesign{Tm <: AbstractMod, Tcon <: AbstractConstraints,Tspace <: AbstractSpace{T} where T} <: SingleDesign{Tspace}
    mod::Tm
    len::Int
    con::Tcon
    space:: Tspace
end

#FIXME  makeorderedspace ? make_ordered_space
"""
    makeorderedspace(mod::GroupMod, len::Int)

The function generates all ordered constructs with given length and modules. The
modules can be infinity redrawn. The constructs are not explicitly generated.
With the use of the Kronecker product, every construct can be made on the fly
only on the moment when it is required.
"""
makeorderedspace(mod::GroupMod, len::Int) = reshape(mod.m, length(mod.m),1) |> y -> ⊗(y,len)
#makesinglespace(mod::GroupMod) = reshape(mod.m, length(mod.m),1) |> y -> ⊗(y,ones(Int64,length(mod.m),1))
makesinglespace(mod::GroupMod)  = error("no order design with length 1")
"""
    UnorderedDesign

Unordered space is where the position in the construct has a non-functional role,
the construct `[a,b,c]`` is seen as equal to [c,b,a], so only one of both will be
generated.

Al modules can only be used once.

Fields:
    - `mod::AbstractMod`: The used moduels in this design
    - `len::Int`: The allowed length of the constructs
    - `con::AbstractConstraints`: The constraints of the design space
    - `space::AbstractSpace`: Design space containing all constructs, if possible, these constructs are not explicitly generated.
"""
struct UnorderedDesign{Tm <: AbstractMod,Tcon <: AbstractConstraints,Tspace <: AbstractSpace{T} where T} <: SingleDesign{Tspace}
    mod::Tm
    len::Int
    con::Tcon
    space:: Tspace
end

"""
    makeunorderedspace(mod::GroupMod,len::Int)

The function generates all unordered constructs with given length and modules.
Every module can only be used one single time in a construct.
The constructs aren't explicitly generated.
With the use of the ``Combination`` structure, every construct can be made on
the fly only on the moment it is required.
"""
makeunorderedspace(mod::GroupMod, len::Int) = Combination(mod.m,len)

"""
    constructdesign(mod::GroupMod, len::Int; con::NoConstraint = NoConstraint(nothing) , order::Bool = false)

Returns a correct filled `UnorderedDesign` with given input modules and desired length.
If `order = true`, then an `OrderedDesign` is returned using the same input settings
See [`UnorderedDesign`](@ref) or [`OrderedDesign`](@ref)
"""
function constructdesign(mod::GroupMod, len::Int; order::Bool = false)
    if order == true
        if len == 1
            makesinglespace(mod) |> FullOrderedSpace |> (y -> OrderedDesign(mod,len,NoConstraint(nothing),y))
        else
            makeorderedspace(mod,len) |> FullOrderedSpace |> (y -> OrderedDesign(mod,len,NoConstraint(nothing),y))
        end
    else
        makeunorderedspace(mod,len) |> FullUnorderedSpace |> (y -> UnorderedDesign(mod,len,NoConstraint(nothing),y))
    end
end

#FIXME: this function is too long for a one-liner
# precomuted full space all other field are ignored, currently not used , manby usefull in the future.
function constructdesign(Tspace::ComputedSpace ; mod = GroupMod(nothing),
     len = LenConstraint(nothing), con = NoConstraint(nothing), order= false)
     if order == true
         return  OrderedDesign(mod,len,con,Tspace)
     else
         return UnorderedDesign(mod,len,con,Tspace)
     end
end

"""
    constructdesign(mod::GroupMod, len::Int, con::ConstructConstraints; order = false)

Returns a correct filled `UnorderedDesign` with given input modules and desired
length. The constraints are attributed and evaluated if they are useful for the
given constructed type.

Ordered constraints are only helpful if  `order = true`, then an `OrderedDesign`
is returned.

See [`UnorderedDesign`](@ref) or [`OrderedDesign`](@ref)
"""
function constructdesign(mod::GroupMod, len::Int, con::ConstructConstraints{T} where T; order = false)
    if order == true
        makeorderedspace(mod,len) |> x -> FrameSpace(FullOrderedSpace(x),con) |> (y -> OrderedDesign(mod,len,con,y))
    else
        if isa(con,UnOrderedConstraint) || promotetype(eltype(con),UnOrderedConstraint) == UnOrderedConstraint
            return makeunorderedspace(mod,len) |> x -> FrameSpace(FullUnorderedSpace(x),con) |> (y -> UnorderedDesign(mod,len,con,y))
        else
            error("Ordered contrains can't be used in unorderd design set order to true or remove Constraints")
        end
    end
end

"""
    getspace(Design::SingleDesign; full = false)

If the input design has no Constraints, then an efficient space is generated
This space type allow most `Base` function that you would have on the explicit generated design space.
If constraints are added than no efficient space can be generated.
For most functionalities explicitly, generation is needed.
In some case one can continue using the whole design structure.
If `full = true` the whole space is generated explicitly, which isn't recommended for large designs.
Still for constrained problems is sometimes the best or only option.
"""

function getspace(Design::SingleDesign; full = false)
          if full== false
              if  isa(Design.space, FrameSpace)
                  @info "given space is the closed effienct space that can be made, not all Constraints are used"
              end
              return Design.space
          else
            return  createdspace(Design.space,Design.con)
        end
end

### Create space can be removed in the future because Base.iterate{Framespace} has same functionalities.
#was implented before and it works so kept for now.

"""
    createdspace(space::EffSpace)

Generated explicated all constructed for a space that can be handled efficiently
"""
function createdspace(space::EffSpace,::AbstractConstraints)
    newspace = [construct for construct in space]
    return ComputedSpace(newspace)
end

"""
    createdspace(space::Framespace,con::Constructcostrains)

It generated the space that has one or multiple constraints.
The constructed that aren't allowed base on given constraints are filtered out.
"""
function createdspace(space::FrameSpace,con::ConstructConstraints)
    new_space = eltype(space)[]
    for temp_construct in space.space
        #print(tempconstruct)
            if filterconstraint(temp_construct,con) == false
                 push!(newspace,tempconstruct)
            end
   end
    return ComputedSpace(new_space)
end

"""
    createdspace(space::Comuputedspace)

fall back if space was computed before.
"""
createdspace(space::ComputedSpace,::AbstractConstraints) = space.space

"""
    constraintcheck(lenconstruct::Int,Constraint)

quick check if the Constraint make sense given the length of the constructs.
"""
constraintcheck(len_construct::Int,constraint) = (len_construct >= length(constraint.combination)) && (len_construct >= maximum(constraint.pos))

""""
    filterconstraint(construct::OrderedConstruct, con::OrderedConstraint)

Evaluated if the `construct` is allowed given  `con`.  Returns `true` when the construct isn't allowed.
Because the constraint is ordered, the given modules are only tested on the given position. Only with an exact match, it returns true,if no match is found it returns `false`.
"""
function filterconstraint(construct::OrderedConstruct, con::OrderedConstraint)
    # check if Constraints makes sence to given construct
    if constraintcheck(length(construct),con) == false
        @warn "skiped constraint that didn't make sence"
        return false
    else
         for (pos,element) in  zip(con.pos,con.combination)
             if construct[pos] != element
                 return false
            end
        end
    end
    return true
end

""""
    filterconstraint(construct::AbstractConstruct, con::UnOrderedConstraint)

Evaluated if the `construct` is allowed given  `con`.  Returns `true` when the construct isn't allowed.  If no match is found it returns `false`
Because the constraint is unordered only all modules in the constraint need to be present in construct.
The position is not evaluated.
"""

function filterconstraint(construct::AbstractConstruct, con::UnOrderedConstraint)
    con_count = count_elements(con.combination)
    construct_count =  count_elements(construct.c)
    for (key,n) in con_count
        if haskey(construct_count,key)
            if construct_count[key] < n
                return false
           end
       else
           return false
       end
   end
   return true
end

function count_elements(a::Array)
     n = Dict(key => 0 for key in a)
     for element in a
         n[element] +=1
    end
    return n
end

""""
    filterconstraint(construct::AbstractConstruct, con::ComposedConstructConstraints)

A wrapper to evaluated multiple constraints regarding the same constructs.
The function terminates if one of the constraint matches to the given constructs and returns `true`.
"""

function filterconstraint(construct::AbstractConstruct, con::ComposedConstructConstraints)
    for temp_con in con.constructcon
        if filterconstraint(construct,temp_con)
            return true
        end
    end
    return false
end

"""
    MultiDesign{Tm <: SingleDesign{T} where T} <: AbstractDesign{T}

Designs that are constructed out of different `SingleDesign`.
The underlining Designs are stored  in a Vector.
Currently only used to allow multiple lengths of constructs for the same design rules ( modules and constraints ).
"""

struct MultiDesign{T} <: AbstractDesign{T}
    d::Vector{T}
end

function Base.show(io::IO,d::MultiDesign)
    l = length(d)
   print(io,"Number of single designs: $l \n")
   for i in d
       print(io,"\n", i , "\n")
   end
end

"""
addition of differnt designs
"""
Base.eltype(::MultiDesign{T}) where T = T
Base.length(design::MultiDesign) = length(design.d)
Base.iterate(design::MultiDesign,state = 1) = length(design.d) >= state ? (design[state],state+1) : nothing
Base.getindex(design::MultiDesign,i) = design.d[i]
Base.size(design::MultiDesign) = (length(design), map((i -> length(i.space,nothing)),design) |> sum )

"""
    constructdesign(mod::GroupMod, min::Int, max::Int; order = false)

Returns MulitDesign contianing the design space for ever length between `min` and `max`.
If `order = true`, then an `OrderedDesign` is constructed.
"""

constructdesign(mod::GroupMod, min::Int, max::Int; order = false) = MultiDesign([constructdesign(mod,len,order = order) for len in min:max])

"""
    constructdesign(mod::GroupMod, min::Int, max::Int; order = false)

Returns MulitDesign contianing the design space for ever length in the `lens` array.
If `order = true`, then an `OrderedDesign` is constructed.
"""

constructdesign(mod::GroupMod, lens::Array{Int}; order = false) = MultiDesign([constructdesign(mod,len,order = order) for len in lens])

"""
    constructdesign(mod::GroupMod, min::Int, max::Int; order = false)

Returns MulitDesign contianing the design space for ever length in the `lens` array.
If `order = true`, then an `OrderedDesign` is constructed.
"""

constructdesign(mod::GroupMod, lens::Array{Int}, con::ConstructConstraints{T} where T; order = false) = MultiDesign([constructdesign(mod,len,con,order = order) for len in lens])

"""
    getspace(multipledesign::MultiDesign; full = false )

Simple wrapper to obtain all different space in the MultiDesign object using `getspace(Design::SingleDesign; full = false)` iteratively.
type determined on the space of the first design
"""

function getspace(design::MultiDesign; full = false )

    if full
        all_spaces = collect(design[1].space)
        for i in 2:length(design)
            append!(all_spaces,collect(design[i].space))
        end

        return ComputedSpace(all_spaces)

    else
        all_spaces = AbstractSpace[]
        for d in design
             push!(all_spaces,getspace(d))
        end
        return MultiSpace(all_spaces)
    end

end

# Combining of combinations if not all intermediate-length are desired.

Base.:+(des1::SingleDesign,des2::SingleDesign) where T = MultiDesign([des1,des2])
Base.:+(des1::MultiDesign{T} where T, des2::SingleDesign) = MultiDesign([des1.d...,des2])
Base.:+(des1::SingleDesign,des2::MultiDesign) = +(des2,des1)
Base.:+(des1::MultiDesign{T} where T, des2::MultiDesign{N} where N) = MultiDesign([des1.d...,des2.d...])

max_stars_repo_name: SamDier/BOMoD

max_issues_repo_name: SamDier/BOMoD

max_forks_repo_name: SamDier/BOMoD


text:
function mass(u, q, w, v, x, t, mu, eta)
    m = [1.0, 1.0, 1.0, 1.0, 1.0];
    return m;
end
function flux(u, q, w, v, x, t, mu, eta)
    gam = mu[1];
    gam1 = gam - 1.0;
    Re = mu[2];
    Pr = mu[3];
    Minf = mu[4];
    Re1 = 1/Re;
    M2 = Minf^2;
    c23 = 2.0/3.0;
    fc = 1/(gam1*M2*Re*Pr);
    r = u[1];
    ru = u[2];
    rv = u[3];
    rw = u[4];
    rE = u[5];
    rx = q[1];
    rux = q[2];
    rvx = q[3];
    rwx = q[4];
    rEx = q[5];
    ry = q[6];
    ruy = q[7];
    rvy = q[8];
    rwy = q[9];
    rEy = q[10];
    rz = q[11];
    ruz = q[12];
    rvz = q[13];
    rwz = q[14];
    rEz = q[15];
    r1 = 1/r;
    uv = ru*r1;
    vv = rv*r1;
    wv = rw*r1;
    E = rE*r1;
    ke = 0.5*(uv*uv+vv*vv+wv*wv);
    p = gam1*(rE-r*ke);
    h = E+p*r1;
    fi = [ru, ru*uv+p, rv*uv, rw*uv, ru*h, rv, ru*vv, rv*vv+p, rw*vv, rv*h, rw, ru*wv, rv*wv, rw*wv+p, rw*h];
    ux = (rux - rx*uv)*r1;
    vx = (rvx - rx*vv)*r1;
    wx = (rwx - rx*wv)*r1;
    kex = uv*ux + vv*vx + wv*wx;
    px = gam1*(rEx - rx*ke - r*kex);
    Tx = gam*M2*(px*r - p*rx)*r1^2;
    uy = (ruy - ry*uv)*r1;
    vy = (rvy - ry*vv)*r1;
    wy = (rwy - ry*wv)*r1;
    key = uv*uy + vv*vy + wv*wy;
    py = gam1*(rEy - ry*ke - r*key);
    Ty = gam*M2*(py*r - p*ry)*r1^2;
    uz = (ruz - rz*uv)*r1;
    vz = (rvz - rz*vv)*r1;
    wz = (rwz - rz*wv)*r1;
    kez = uv*uz + vv*vz + wv*wz;
    pz = gam1*(rEz - rz*ke - r*kez);
    Tz = gam*M2*(pz*r - p*rz)*r1^2;
    txx = Re1*c23*(2*ux - vy - wz);
    txy = Re1*(uy + vx);
    txz = Re1*(uz + wx);
    tyy = Re1*c23*(2*vy - ux - wz);
    tyz = Re1*(vz + wy);
    tzz = Re1*c23*(2*wz - ux - vy);
    fv = [0, txx, txy, txz, uv*txx + vv*txy + wv*txz + fc*Tx, 0, txy, tyy, tyz, uv*txy + vv*tyy + wv*tyz + fc*Ty, 0, txz, tyz, tzz, uv*txz + vv*tyz + wv*tzz + fc*Tz];
    f = fi+fv;
    f = reshape(f,(5,3));
    return f;
end
function source(u, q, w, v, x, t, mu, eta)
    s = [0.0, 0.0, 0.0, 0.0, 0.0];
    return s;
end
function ubou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    ub = [0.0, 0.0, 0.0, 0.0, 0.0];
    return ub;
end
function fbou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    f = flux(u, q, w, v, x, t, mu, eta);
    fb = f[:,1]*n[1] + f[:,2]*n[2] + f[:,3]*n[3] + tau[1]*(u-uhat);
    return fb;
end
function initu(x, mu, eta)
    L = 2*pi;
    x1 = x[1];
    x2 = x[2];
    x3 = x[3];
    gam = mu[1];
    Minf = mu[4];
    p0 = 1/(gam*Minf^2);

    u01 = 1.0;
    u02 = sin(2*pi*(x1-L/2)/L) * cos(2*pi*(x2-L/2)/L) * cos(2*pi*(x3-L/2)/L);
    u03 = - cos(2*pi*(x1-L/2)/L) * sin(2*pi*(x2-L/2)/L) * cos(2*pi*(x3-L/2)/L);
    u04 = 0.0;
    p = p0 + (1/16) * (cos(2*2*pi*(x1-L/2)/L) + cos(2*2*pi*(x2-L/2)/L)) * (cos(2*2*pi*(x3-L/2)/L) + 2.0);
    u05 = 0.5*(u02^2+u03^2+u04^2) / u01 + p / (gam-1);

    u0 = [u01, u02, u03, u04, u05];
    return u0;
end

ext: jl

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
# This file is a part of DataDeps.jl. License is MIT.

# TODO Remove this whole thing once https://github.com/JuliaWeb/HTTP.jl/pull/273

"""
    fetch_http(remotepath, localdir)

Pass in a HTTP[/S] URL  and a directory to save it to,
and it downloads that file, returing the local path.
This is using the HTTP protocol's method of defining filenames in headers,
if that information is present.
"""
function fetch_http(remotepath, localdir)
    @assert(localdir |> isdir)
    filename = get_filename(remotepath)
    localpath = safer_joinpath(localdir, filename)
    Base.download(remotepath, localpath)
end

"""
    safer_joinpath(basepart, parts...)

A variation on `joinpath`, that is more resistant to directory traveral attack
The parts to be joined (excluding the `basepart`),
are not allowed to contain `..`, or begin with a `/`.
If they do then this throws an `DomainError`.
"""
function safer_joinpath(basepart, parts...)
    explain =  "Possible Directory Traversal Attack detected."
    for part in parts
        occursin(part, "..") && throw(DomainError(part, "contains illegal string \"..\". $explain"))
        startswith(part, '/') && throw(DomainError(part, "begins with \"/\". $explain"))
    end
    joinpath(basepart, parts...)
end

"""
    get_filename(remotepath)

Given a remotepath (URL) returns the filename that it should be saved to locally.
"""
function get_filename(remotepath)
    filename = try
        try_get_filename(remotepath)
    catch err
        # Catch *everything* here, as we can always recover and there are many things that can go wrong
        @warn("Could not resolve filename due to error (Fallback engaged)",remotepath, exception=err)
        filename = nothing
    end

    if filename == nothing
        # couldn't get it from the headers
        filename = basename(remotepath)
    end
    filename
end

"""
    try_get_filename(url)

Uses as HEAD request, to attempt to retrieve the filename from the HTTP headers.
Returns a string or nothing if it failes
"""
function try_get_filename(url)
    resp = HTTP.request("HEAD", url,  ["User-Agent"=>"DataDeps.jl (http-lib: HTTP.jl; lang: Julia)"]);
    content_disp = get(Dict(resp.headers),"Content-Disposition","");
    raw = match(r"filename\s*=\s*(.*)", content_disp); # TECHDEBT: Consider if more of this should be moved to process header filename

    process_header_filename(raw)
end

"""
    process_header_filename(raw)

Deal with some of the weird and varied ways filenames can be given.
Not full coverage, but getting the common cases.

Return nothing if input is nothing
"""
function process_header_filename(raw::RegexMatch)::String
    ret = raw[1]
    quoted_match = match(r"\"(.*)\"", ret)

    if quoted_match != nothing
        ret = unescape_string(quoted_match[1]) #It was in quotes, so it will be double escaped
    end

    strip(ret)
end

process_header_filename(::Nothing) = nothing


text:
mutable struct DocumentTermMatrix{T}
    dtm::SparseMatrixCSC{Int, Int}
    terms::Vector{T}
    column_indices::Dict{T, Int}
end

function serialize(io::AbstractSerializer, dtm::DocumentTermMatrix{T}) where {T}
    Serialization.writetag(io.io, Serialization.OBJECT_TAG)
    serialize(io, DocumentTermMatrix{T})
    serialize(io, dtm.dtm)
    serialize(io, dtm.terms)
    nothing
end

function deserialize(io::AbstractSerializer, ::Type{DocumentTermMatrix{T}}) where {T}
    dtm = deserialize(io)
    terms = deserialize(io)
    column_indices = Dict{T,Int}(term => idx for (idx,term) in enumerate(terms))
    DocumentTermMatrix{T}(dtm, terms, column_indices)
end

"""
    columnindices(terms::Vector{String})

Creates a column index lookup dictionary from a vector of terms.
"""
function columnindices(terms::Vector{T}) where T
    column_indices = Dict{T, Int}()
    for i in 1:length(terms)
        term = terms[i]
        column_indices[term] = i
    end
    column_indices
end

"""
    DocumentTermMatrix(crps::Corpus)
    DocumentTermMatrix(crps::Corpus, terms::Vector{String})
    DocumentTermMatrix(crps::Corpus, lex::AbstractDict)
    DocumentTermMatrix(dtm::SparseMatrixCSC{Int, Int},terms::Vector{String})

Represent documents as a matrix of word counts.

Allow us to apply linear algebra operations and statistical techniques.
Need to update lexicon before use.

# Examples
```julia-repl
julia> crps = Corpus([StringDocument("To be or not to be"),
                      StringDocument("To become or not to become")])

julia> update_lexicon!(crps)

julia> m = DocumentTermMatrix(crps)
A 2 X 6 DocumentTermMatrix

julia> m.dtm
2×6 SparseArrays.SparseMatrixCSC{Int64,Int64} with 10 stored entries:
  [1, 1]  =  1
  [2, 1]  =  1
  [1, 2]  =  2
  [2, 3]  =  2
  [1, 4]  =  1
  [2, 4]  =  1
  [1, 5]  =  1
  [2, 5]  =  1
  [1, 6]  =  1
  [2, 6]  =  1
```
"""
function DocumentTermMatrix(crps::Corpus, terms::Vector{T}) where T
    column_indices = columnindices(terms)

    m = length(crps)
    n = length(terms)

    rows = Array{Int}(undef, 0)
    columns = Array{Int}(undef, 0)
    values = Array{Int}(undef, 0)
    for i in 1:m
        doc = crps.documents[i]
        ngs = ngrams(doc)
        for ngram in keys(ngs)
            j = get(column_indices, ngram, 0)
            v = ngs[ngram]
            if j != 0
                push!(rows, i)
                push!(columns, j)
                push!(values, v)
            end
        end
    end
    if length(rows) > 0
        dtm = sparse(rows, columns, values, m, n)
    else
        dtm = spzeros(Int, m, n)
    end
    DocumentTermMatrix(dtm, terms, column_indices)
end
DocumentTermMatrix(crps::Corpus) = DocumentTermMatrix(crps, lexicon(crps))

DocumentTermMatrix(crps::Corpus, lex::AbstractDict) = DocumentTermMatrix(crps, sort(collect(keys(lex))))

DocumentTermMatrix(dtm::SparseMatrixCSC{Int, Int},terms::Vector{T}) where T = DocumentTermMatrix{T}(dtm, terms, columnindices(terms))

"""
    dtm(crps::Corpus)
    dtm(d::DocumentTermMatrix)
    dtm(d::DocumentTermMatrix, density::Symbol)

Creates a simple sparse matrix of DocumentTermMatrix object.

# Examples
```julia-repl
julia> crps = Corpus([StringDocument("To be or not to be"),
                      StringDocument("To become or not to become")])

julia> update_lexicon!(crps)

julia> dtm(DocumentTermMatrix(crps))
2×6 SparseArrays.SparseMatrixCSC{Int64,Int64} with 10 stored entries:
  [1, 1]  =  1
  [2, 1]  =  1
  [1, 2]  =  2
  [2, 3]  =  2
  [1, 4]  =  1
  [2, 4]  =  1
  [1, 5]  =  1
  [2, 5]  =  1
  [1, 6]  =  1
  [2, 6]  =  1

julia> dtm(DocumentTermMatrix(crps), :dense)
2×6 Array{Int64,2}:
 1  2  0  1  1  1
 1  0  2  1  1  1
```
"""
function dtm(d::DocumentTermMatrix, density::Symbol)
    if density == :sparse
        return d.dtm
    else
        return Matrix(d.dtm)
    end
end

function dtm(d::DocumentTermMatrix)
    return d.dtm
end

function dtm(crps::Corpus)
    dtm(DocumentTermMatrix(crps))
end

tdm(crps::DocumentTermMatrix, density::Symbol) = dtm(crps, density)' #'

tdm(crps::DocumentTermMatrix) = dtm(crps)' #'

tdm(crps::Corpus) = dtm(crps)' #'

##############################################################################
#
# Produce the signature of a DTM entry for a document
#
# TODO: Make this more efficieny by reusing column_indices
#
##############################################################################

function dtm_entries(d::AbstractDocument, lex::Dict{T, Int}) where T
    ngs = ngrams(d)
    indices = Array{Int}(undef, 0)
    values = Array{Int}(undef, 0)
    terms = sort(collect(keys(lex)))
    column_indices = columnindices(terms)

    for ngram in keys(ngs)
        if haskey(column_indices, ngram)
            push!(indices, column_indices[ngram])
            push!(values, ngs[ngram])
        end
    end
    return (indices, values)
end

"""
    dtv(d::AbstractDocument, lex::Dict{String, Int})

Produce a single row of a DocumentTermMatrix.

Individual documents do not have a lexicon associated with them,
we have to pass in a lexicon as an additional argument.

# Examples
```julia-repl
julia> dtv(crps[1], lexicon(crps))
1×6 Array{Int64,2}:
 1  2  0  1  1  1
```
"""
function dtv(d::AbstractDocument, lex::Dict{T, Int}) where T
    p = length(keys(lex))
    row = zeros(Int, 1, p)
    indices, values = dtm_entries(d, lex)
    for i in 1:length(indices)
        row[1, indices[i]] = values[i]
    end
    return row
end

function dtv(d::AbstractDocument)
    error("Cannot construct a DTV without a pre-existing lexicon")
end

##############################################################################
#
# The hash trick: use a hash function instead of a lexicon to determine the
# columns of a DocumentTermMatrix-like encoding of the data
#
##############################################################################
"""
    hash_dtv(d::AbstractDocument)
    hash_dtv(d::AbstractDocument, h::TextHashFunction)

Represents a document as a vector with N entries.

# Examples
```julia-repl
julia> crps = Corpus([StringDocument("To be or not to be"),
                      StringDocument("To become or not to become")])

julia> h = TextHashFunction(10)
TextHashFunction(hash, 10)

julia> hash_dtv(crps[1], h)
1×10 Array{Int64,2}:
 0  2  0  0  1  3  0  0  0  0

julia> hash_dtv(crps[1])
1×100 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  …  0  0  0  0  0  0  0  0  0  0  0  0
```
"""
function hash_dtv(d::AbstractDocument, h::TextHashFunction)
    p = cardinality(h)
    res = zeros(Int, 1, p)
    ngs = ngrams(d)
    for ng in keys(ngs)
        res[1, index_hash(ng, h)] += ngs[ng]
    end
    return res
end

hash_dtv(d::AbstractDocument) = hash_dtv(d, TextHashFunction())

"""
    hash_dtm(crps::Corpus)
    hash_dtm(crps::Corpus, h::TextHashFunction)

Represents a Corpus as a Matrix with N entries.
"""
function hash_dtm(crps::Corpus, h::TextHashFunction)
    n, p = length(crps), cardinality(h)
    res = zeros(Int, n, p)
    for i in 1:length(crps)
        doc = crps.documents[i]
        res[i, :] = hash_dtv(doc, h)
    end
    return res
end

hash_dtm(crps::Corpus) = hash_dtm(crps, hash_function(crps))

hash_tdm(crps::Corpus) = hash_dtm(crps)' #'

##############################################################################
#
# Produce entries for on-line analysis when DTM would not fit in memory
#
##############################################################################

mutable struct EachDTV
    crps::Corpus
end

start(edt::EachDTV) = 1

function next(edt::EachDTV, state::Int)
    return (dtv(edt.crps.documents[state], lexicon(edt.crps)), state + 1)
end

done(edt::EachDTV, state::Int) = state > length(edt.crps.documents)

mutable struct EachHashDTV
    crps::Corpus
end

start(edt::EachHashDTV) = 1

function next(edt::EachHashDTV, state::Int)
    (hash_dtv(edt.crps.documents[state]), state + 1)
end

done(edt::EachHashDTV, state::Int) = state > length(edt.crps.documents)

each_dtv(crps::Corpus) = EachDTV(crps)

each_hash_dtv(crps::Corpus) = EachHashDTV(crps)

##
## getindex() methods
##

Base.getindex(dtm::DocumentTermMatrix, k::AbstractString) = dtm.dtm[:, dtm.column_indices[k]]
Base.getindex(dtm::DocumentTermMatrix, i::Any) = dtm.dtm[i]
Base.getindex(dtm::DocumentTermMatrix, i::Any, j::Any) = dtm.dtm[i, j]

"""
    prune!(dtm::DocumentTermMatrix{T}, document_positions; compact::Bool=true, retain_terms::Union{Nothing,Vector{T}}=nothing) where {T}

Delete documents specified by `document_positions` from a document term matrix. Optionally compact the matrix by removing unreferenced terms.
"""
function prune!(dtm::DocumentTermMatrix{T}, document_positions; compact::Bool=true, retain_terms::Union{Nothing,Vector{T}}=nothing) where {T}
    if ((document_positions === nothing) || isempty(document_positions))
        dtm_matrix = dtm.dtm
    else
        docrows_to_retain = [!(idx in document_positions) for idx in 1:size(dtm.dtm, 1)]
        dtm_matrix = dtm.dtm[docrows_to_retain, :]
    end

    if compact
        termcols_to_delete = map(x->x==0, sum(dtm_matrix, dims=(1,)))
        if retain_terms !== nothing
            for idx in 1:length(termcols_to_delete)
                (!termcols_to_delete[idx] || !(dtm.terms[idx] in retain_terms)) && continue
                termcols_to_delete[idx] = false
            end
        end
    else
        termcols_to_delete = Bool[]
    end

    if any(termcols_to_delete)
        dtm.dtm = dtm_matrix[:,[!termcols_to_delete[idx] for idx in 1:length(termcols_to_delete)]]
        dtm.terms = [dtm.terms[idx] for idx in 1:length(dtm.terms) if !termcols_to_delete[idx]]
        dtm.column_indices = Dict{T,Int}(term => idx for (idx,term) in enumerate(dtm.terms))
    else
        dtm.dtm = dtm_matrix
    end

    dtm
end

"""
    merge!(dtm1::DocumentTermMatrix{T}, dtm2::DocumentTermMatrix{T}) where {T}

Merge one DocumentTermMatrix instance into another. Documents are appended to the end. Terms are re-sorted.
For efficiency, this may result in modifications to dtm2 as well.
"""
function merge!(dtm1::DocumentTermMatrix{T}, dtm2::DocumentTermMatrix{T}) where {T}
    (length(dtm2.dtm) == 0) && (return dtm1)

    ncombined_docs = size(dtm1.dtm,1) + size(dtm2.dtm,1)
    termset1 = Set(dtm1.terms)
    termset2 = Set(dtm2.terms)
    termset = union(termset1, termset2)

    if termset1 == termset
        # no new terms added
        combined_terms = dtm1.terms
        ncombined_terms = length(dtm1.terms)
    else
        combined_terms = sort!(collect(termset))
        ncombined_terms = length(combined_terms)
    end

    function permute_terms!(dtm_to_permute, terms)
        (length(dtm_to_permute) == 0) && (return dtm_to_permute)
        terms_perm = map(x->(x===nothing) ? 0 : x, indexin(combined_terms, terms))
        remaining_cols = setdiff(1:ncombined_terms, terms_perm)
        for idx in 1:length(terms_perm)
            if terms_perm[idx] == 0
                terms_perm[idx] = popfirst!(remaining_cols)
            end
        end
        permute!(dtm_to_permute, 1:size(dtm_to_permute,1), terms_perm)
    end
    function expand_columns(S, n)
        (S.n == n) && (return S)
        @assert (n > S.n)
        colptr = S.colptr
        resize!(colptr, n+1)
        colptr[(S.n+2):(n+1)] .= colptr[S.n+1]
        SparseMatrixCSC(S.m, n, colptr, S.rowval, S.nzval)
    end
    function row_append(A, B)
        @assert size(A,2) == size(B,2)
        (length(A) == 0) && (return B)
        (length(B) == 0) && (return A)

        C_colptr = similar(A.colptr)
        C_rowvals = similar(A.rowval, length(A.rowval) + length(B.rowval))
        C_nzval = similar(A.nzval, length(A.nzval) + length(B.nzval))

        offset = 0
        rowval_pos = 0
        nzval_pos = 0
        for col in 1:(length(C_colptr)-1)
            colptr_pos = C_colptr[col] = A.colptr[col] + offset
            # first copy from A
            nvalsA = A.colptr[col+1] - A.colptr[col]
            if nvalsA > 0
                C_rowvals[colptr_pos:(colptr_pos+nvalsA-1)] .= A.rowval[A.colptr[col]:(A.colptr[col+1]-1)]
                C_nzval[colptr_pos:(colptr_pos+nvalsA-1)] .= A.nzval[A.colptr[col]:(A.colptr[col+1]-1)]
                colptr_pos += nvalsA
            end
            # then copy from B
            nvalsB = B.colptr[col+1] - B.colptr[col]
            if nvalsB > 0
                C_rowvals[colptr_pos:(colptr_pos+nvalsB-1)] .= (B.rowval[B.colptr[col]:(B.colptr[col+1]-1)] .+ size(A,1))
                C_nzval[colptr_pos:(colptr_pos+nvalsB-1)] .= B.nzval[B.colptr[col]:(B.colptr[col+1]-1)]
                offset += nvalsB
                colptr_pos += nvalsB
            end
        end
        C_colptr[end] = length(C_rowvals)+1
        SparseMatrixCSC(size(A,1) + size(B,1), size(A,2), C_colptr, C_rowvals, C_nzval)
    end

    dtm1_matrix = (combined_terms === dtm1.terms) ? dtm1.dtm : permute_terms!(expand_columns(dtm1.dtm, ncombined_terms), dtm1.terms)
    dtm2_matrix = permute_terms!(expand_columns(dtm2.dtm, ncombined_terms), dtm2.terms)
    combined_matrix = row_append(dtm1_matrix, dtm2_matrix)

    # set new terms and recompute column_indices
    dtm1.dtm = combined_matrix
    dtm1.terms = combined_terms
    dtm1.column_indices = Dict{T,Int}(term => idx for (idx,term) in enumerate(combined_terms))

    dtm1
end


text:
using Docile, Lexicon, Gnuplot

# config = Config()
config = Config(md_subheader=:skip, mathjax = true)

index = save("Gnuplot.md", Gnuplot, config)
save("index.md", Index([index]), config)

# run(`../mkdocs build`)
# mkdocs gh-deploy --clean

max_stars_repo_name: johncsnyder/Gnuplot


text:
#= test/spfft1.jl
=#

println("spfft2.jl")
@time begin

@testset "spfft2" begin

n1 = 50
n2 = 50
k1 = 11
k2 = 10

# complex transforms

r1 = randperm(n1)[1:k1]
r2 = randperm(n2)[1:k2]

for fcn in (:fft, :bfft)
  sf = Symbol("sp", fcn)
  for T in (ComplexF32, ComplexF64)
    println("  $sf/$T")
    Tr = real(T)

    psf = Symbol("plan_", sf)
    @eval P = $psf($T, $n1, $n2, $r1, $r2)

    x  = Array{T }(undef, n1, n2)
    y  = Array{T }(undef, k1, k2)
    xr = Array{Tr}(undef, size(x))
    yr = Array{Tr}(undef, size(y))

    # f2s
    f2s = Symbol(sf, "_f2s!")
    rand!(x)

    ## (c2c, f2s)
    @eval f = $fcn($x)[$r1,$r2]
    @eval $f2s($y, P, $x)
    @test f ≈ y

    ## (c2r, f2s)
    @eval $f2s($yr, P, $x)
    @test real(f) ≈ yr

    ## (r2c, f2s)
    xr[:] = real(x)
    @eval f = $fcn($xr)[$r1,$r2]
    @eval $f2s($y, P, $xr)
    @test f ≈ y

    # s2f
    s2f = Symbol(sf, "_s2f!")
    rand!(y)

    ## (c2c, s2f)
    fill!(x, 0)
    x[r1,r2] = y
    @eval f = $fcn($x)
    @eval $s2f($x, P, $y)
    @test f ≈ x

    ## (c2r, s2f)
    @eval $s2f($xr, P, $y)
    @test real(f) ≈ xr

    ## (r2c, s2f)
    yr[:] = real(y)
    fill!(x, 0)
    x[r1,r2] = yr
    @eval f = $fcn($x)
    @eval $s2f($x, P, $yr)
    @test f ≈ x
  end
end

# real transforms

m1 = div(n1,2) + 1
m2 = div(n2,2) + 1
r1 = vcat([1,m1], 1 .+ randperm(m1-2)[1:k1-2])  # include edge cases
r2 = vcat([1,m2], 1 .+ randperm(m2-2)[1:k2-2])

for T in (Float32, Float64)
  Tc = Complex{T}

  x  = Array{T }(undef, n1, n2)
  y  = Array{Tc}(undef, k1, k2)
  xc = Array{Tc}(undef, m1, n2)

  # (r2c, f2s)
  println("  sprfft/$T")
  P = plan_sprfft(T, n1, n2, r1, r2)
  rand!(x)
  f = rfft(x)[r1,r2]
  sprfft_f2s!(y, P, x)
  @test f ≈ y

  # (c2r, s2f)
  println("  spbrfft/$T")
  P = plan_spbrfft(T, n1, n2, r1, r2)
  rand!(y)
  y[r1.==1,r2.==1] = real(y[r1.==1,r2.==1])
  even1 = n1 % 2 == 0
  even2 = n2 % 2 == 0
  even1 && (y[r1.==m1,1] = real(y[r1.==m1,1]))
  even2 && (y[1,r2.==m2] = real(y[1,r2.==m2]))
  (even1 && even2) && (y[r1.==m1,r2.==m2] = real(y[r1.==m1,r2.==m2]))
  fill!(xc, 0)
  xc[r1,r2] = y
  f = brfft(xc, n1)
  spbrfft_s2f!(x, P, y)
  @test f ≈ x
end

end
end


text:
using Downloads

@testset "WASI - Cowsay" begin
    mktempdir() do tmp_dir
        # TODO: Find a better way to provide the files
        location = joinpath(tmp_dir, "cowsay.wasm")
        Downloads.download(
            "https://registry-cdn.wapm.io/contents/_/cowsay/0.2.0/target/wasm32-wasi/release/cowsay.wasm",
            location,
        )

        wasm_code = open(read, location, "r") |> Wasmtime.WasmVec

        engine = WasmEngine()
        store = Wasmtime.WasmtimeStore(engine)
        wasm_module = Wasmtime.WasmtimeModule(engine, wasm_code)

        msg = "Meuuuh!"
        stdin_path = joinpath(tmp_dir, "./stdin.txt")
        stdout_path = joinpath(tmp_dir, "./stdout.txt")
        open(stdin_path, "w") do file
            write(file, msg)
        end

        # TODO: Fix program name and argv's
        wasi_config = Wasmtime.WasiConfig("wasi_unstable")
        Wasmtime.wasi_config_inherit_stderr(wasi_config)

        Wasmtime.wasi_config_set_stdout_file(wasi_config, stdout_path)
        Wasmtime.wasi_config_set_stdin_file(wasi_config, stdin_path)

        Wasmtime.wasi_config_inherit_env(wasi_config)
        Wasmtime.wasi_config_set_argv(wasi_config, 1, ["cowsay"])

        Wasmtime.set_wasi!(store, wasi_config)

        linker = Wasmtime.WasmtimeLinker(engine)
        Wasmtime.define_wasi!(linker)

        instance = Wasmtime.WasmtimeInstance(linker, store, wasm_module)
        instance_exports = exports(instance)
        main_func = instance_exports._start

        main_func()

        result = read(stdout_path, String)

        expected = raw"""
         _________
        < Meuuuh! >
         ---------
                \   ^__^
                 \  (oo)\_______
                    (__)\       )\/\
                       ||----w |
                        ||     ||
        """

        @test result == expected
    end
end


text:
using DataFrames
using CSV
#using Plots
using Random
using StatsBase

include("mironov.jl");

function get_TB(sim_res, leafpot, pvs, alpha, beta)

smc_mod = sim_res[2][:,1];

eps_real = 0*smc_mod; eps_im = 0*smc_mod;
for i in 1:length(eps_real)
        eps_i = mironov(1.4e9, smc_mod[i], 0.15);
        eps_real[i] = eps_i[1];
        eps_im[i] = eps_i[2];
end

eps = eps_real + eps_im*1im;

theta = 40*pi/180;
mycostheta = cos(theta);
mysin2theta = sin(theta)^2;

fH = abs2.((mycostheta .- sqrt.(eps .- mysin2theta)) ./ (mycostheta .+ sqrt.(eps .- mysin2theta)));
fV = abs2.((eps*mycostheta .- sqrt.(eps .- mysin2theta)) ./ (eps*mycostheta .+ sqrt.(eps .- mysin2theta)));

#leafpot = mean(sim_res1[3],dims=2);
leafvol_rel = 1 .+ pvs*leafpot;
leafvol_rel[leafvol_rel .< 0] .= 0;
laiM = sim_res[1].LAI_modis;

vod = leafvol_rel .* (alpha .+ laiM*beta);

gamma = exp.(-1*vod / mycostheta);
omega = 0.05
rhfac =  exp(-0.16*mycostheta^2)

Tsoil = sim_res[1].T_SOIL .+ 273.15;
Tcan = (sim_res[1].LW_OUT / 0.97 / K_STEFAN ) .^ 0.25;
predH = gamma .* (1 .- rhfac*fH) .* Tsoil + (1-omega) * (1 .- gamma) .* (1 .+ gamma .* (rhfac*fH)) .* Tcan;
predV = gamma .* (1 .- rhfac*fV) .* Tsoil + (1-omega) * (1 .- gamma) .* (1 .+ gamma .* (rhfac*fV)) .* Tcan;

return predH[:,1], predV[:,1], vod

end

function get_TB_2(smc_mod, Tsoil, Tcan, leafpot, laiM, pvs, alpha, beta)

#smc_mod = sim_res[2][:,1];

eps_real = 0*smc_mod; eps_im = 0*smc_mod;
for i in 1:length(eps_real)
        eps_i = mironov(1.4e9, smc_mod[i], 0.15);
        eps_real[i] = eps_i[1];
        eps_im[i] = eps_i[2];
end

eps = eps_real + eps_im*1im;

theta = 40*pi/180;
mycostheta = cos(theta);
mysin2theta = sin(theta)^2;

fH = abs2.((mycostheta .- sqrt.(eps .- mysin2theta)) ./ (mycostheta .+ sqrt.(eps .- mysin2theta)));
fV = abs2.((eps*mycostheta .- sqrt.(eps .- mysin2theta)) ./ (eps*mycostheta .+ sqrt.(eps .- mysin2theta)));

#leafpot = mean(sim_res1[3],dims=2);
leafvol_rel = 1 .+ pvs*leafpot;
leafvol_rel[leafvol_rel .< 0] .= 0;
#laiM = sim_res[1].LAI_modis;

vod = leafvol_rel .* (alpha .+ laiM*beta);

gamma = exp.(-1*vod / mycostheta);
omega = 0.05
rhfac =  exp(-0.16*mycostheta^2)

#Tsoil = sim_res[1].T_SOIL .+ 273.15;
#Tcan = (sim_res[1].LW_OUT / 0.97 / K_STEFAN ) .^ 0.25;
predH = gamma .* (1 .- rhfac*fH) .* Tsoil + (1-omega) * (1 .- gamma) .* (1 .+ gamma .* (rhfac*fH)) .* Tcan;
predV = gamma .* (1 .- rhfac*fV) .* Tsoil + (1-omega) * (1 .- gamma) .* (1 .+ gamma .* (rhfac*fV)) .* Tcan;

return predH[:,1], predV[:,1], vod

end

using LinearAlgebra

function optim_vod_while(obsHi, obsVi, leafpot, smc_mod, simtab, k, alpha, beta, maxiter, stop_crit, step_size)

eps_real = zeros(length(smc_mod)); eps_im = zeros(length(smc_mod));
for i in 1:length(eps_real)
        eps_i = mironov(1.4e9, smc_mod[i], 0.15);
        eps_real[i] = eps_i[1];
        eps_im[i] = eps_i[2];
end

eps = complex.(eps_real, eps_im);

theta = 40*pi/180;
mycostheta = cos(theta);
mysin2theta = sin(theta)^2;

fH = abs2.((mycostheta .- sqrt.(eps .- mysin2theta)) ./ (mycostheta .+ sqrt.(eps .- mysin2theta)));
fV = abs2.((eps*mycostheta .- sqrt.(eps .- mysin2theta)) ./ (eps*mycostheta .+ sqrt.(eps .- mysin2theta)));

omega = 0.05;
rhfac =  exp(-0.16*mycostheta^2);

Tsoil = simtab.T_SOIL .+ 273.15;
Tcan = (simtab.LW_OUT / 0.97 / K_STEFAN ) .^ 0.25;
laiM = simtab.LAI_modis;

function overall_grad(diff_overall, r_rough, dgamma, gamma_star)

	t1 = (1 .- r_rough) .* Tsoil .* dgamma;
	t2 = (1-omega) * (1 .+ gamma_star .* r_rough) .* Tcan .* (-1* dgamma);
	t3 = (1-omega) * (1 .- gamma_star) .* Tcan .* r_rough .* dgamma;
	ans = mean(diff_overall .* (t1+t2+t3));
	return ans;
end

total_err_old = 1000;
err_diff = 1000;
itercount = 0;

while (itercount < maxiter) & (err_diff > stop_crit)

	itercount = itercount + 1;
	leafvol_rel = 1 .+ k*leafpot;
	leafvol_rel[leafvol_rel .< 0] .= 0;
	vod = leafvol_rel .* (alpha .+ laiM * beta);
	gamma = exp.(-1*vod / mycostheta);

	predH_oi = gamma .* (1 .- rhfac*fH) .* Tsoil + (1-omega) * (1 .- gamma) .* (1 .+ gamma .* (rhfac*fH)) .* Tcan;
	predV_oi = gamma .* (1 .- rhfac*fV) .* Tsoil + (1-omega) * (1 .- gamma) .* (1 .+ gamma .* (rhfac*fV)) .* Tcan;

	dgamma_dk = gamma .* (-1 / mycostheta) .* leafpot;
	dgamma_dalpha = gamma .* (-1 / mycostheta) .* leafvol_rel;
	dgamma_dbeta = gamma .* (-1 / mycostheta) .* leafvol_rel .* laiM;

	diffH = (predH_oi .- obsHi); diffV = (predV_oi .- obsVi);

	grad_k = overall_grad(diffH, rhfac*fH, dgamma_dk, gamma) + overall_grad(diffV, rhfac*fV, dgamma_dk, gamma);
	grad_alpha = overall_grad(diffH, rhfac*fH, dgamma_dalpha, gamma) + overall_grad(diffV, rhfac*fV, dgamma_dalpha, gamma);
	grad_beta = overall_grad(diffH, rhfac*fH, dgamma_dbeta, gamma) + overall_grad(diffV, rhfac*fV, dgamma_dbeta, gamma);

	k = k - step_size*grad_k;
	alpha = alpha - step_size*grad_alpha;
	beta = beta - step_size*grad_beta;

	#k = min(max(0.01,k - step_size*grad_k),0.25);
	#alpha = min(max(0.25,alpha - step_size*grad_alpha),2);
	#beta = min(max(1/100,beta - step_size*grad_beta),1/25);

	total_err_new = (sqrt(mean(diffH .^ 2)) + sqrt(mean(diffV .^ 2)))/2;
	err_diff = total_err_old-total_err_new;
	total_err_old = total_err_new*1;

end

        k = min(max(0.00,k),0.75);
        alpha = min(max(0.0,alpha),2);
        beta = min(max(0,beta ),0.2);

return [k, alpha, beta, total_err_old, itercount];
end

function optim_vod_while2(obsHi, obsVi, leafpot, smc_mod, Tsoil, Tcan, laiM, k, alpha, beta, maxiter, stop_crit, step_size)

eps_real = zeros(length(smc_mod)); eps_im = zeros(length(smc_mod));
for i in 1:length(eps_real)
        eps_i = mironov(1.4e9, smc_mod[i], 0.15);
        eps_real[i] = eps_i[1];
        eps_im[i] = eps_i[2];
end

eps = complex.(eps_real, eps_im);

theta = 40*pi/180;
mycostheta = cos(theta);
mysin2theta = sin(theta)^2;

fH = abs2.((mycostheta .- sqrt.(eps .- mysin2theta)) ./ (mycostheta .+ sqrt.(eps .- mysin2theta)));
fV = abs2.((eps*mycostheta .- sqrt.(eps .- mysin2theta)) ./ (eps*mycostheta .+ sqrt.(eps .- mysin2theta)));

omega = 0.05;
rhfac =  exp(-0.16*mycostheta^2);

#Tsoil = simtab.T_SOIL .+ 273.15;
#Tcan = (simtab.LW_OUT / 0.97 / K_STEFAN ) .^ 0.25;
#laiM = simtab.LAI_modis;

function overall_grad(diff_overall, r_rough, dgamma, gamma_star)

	t1 = (1 .- r_rough) .* Tsoil .* dgamma;
	t2 = (1-omega) * (1 .+ gamma_star .* r_rough) .* Tcan .* (-1* dgamma);
	t3 = (1-omega) * (1 .- gamma_star) .* Tcan .* r_rough .* dgamma;
	ans = mean(diff_overall .* (t1+t2+t3));
	return ans;
end

total_err_old = 1000;
err_diff = 1000;
itercount = 0;

while (itercount < maxiter) & (err_diff > stop_crit)

	itercount = itercount + 1;
	leafvol_rel = 1 .+ k*leafpot;
	leafvol_rel[leafvol_rel .< 0] .= 0;
	vod = leafvol_rel .* (alpha .+ laiM * beta);
	gamma = exp.(-1*vod / mycostheta);

	predH_oi = gamma .* (1 .- rhfac*fH) .* Tsoil + (1-omega) * (1 .- gamma) .* (1 .+ gamma .* (rhfac*fH)) .* Tcan;
	predV_oi = gamma .* (1 .- rhfac*fV) .* Tsoil + (1-omega) * (1 .- gamma) .* (1 .+ gamma .* (rhfac*fV)) .* Tcan;

	dgamma_dk = gamma .* (-1 / mycostheta) .* leafpot;
	dgamma_dalpha = gamma .* (-1 / mycostheta) .* leafvol_rel;
	dgamma_dbeta = gamma .* (-1 / mycostheta) .* leafvol_rel .* laiM;

	diffH = (predH_oi .- obsHi); diffV = (predV_oi .- obsVi);

	grad_k = overall_grad(diffH, rhfac*fH, dgamma_dk, gamma) + overall_grad(diffV, rhfac*fV, dgamma_dk, gamma);
	grad_alpha = overall_grad(diffH, rhfac*fH, dgamma_dalpha, gamma) + overall_grad(diffV, rhfac*fV, dgamma_dalpha, gamma);
	grad_beta = overall_grad(diffH, rhfac*fH, dgamma_dbeta, gamma) + overall_grad(diffV, rhfac*fV, dgamma_dbeta, gamma);

	k = k - step_size*grad_k;
	alpha = alpha - step_size*grad_alpha;
	beta = beta - step_size*grad_beta;

	#k = min(max(0.01,k - step_size*grad_k),0.25);
	#alpha = min(max(0.25,alpha - step_size*grad_alpha),2);
	#beta = min(max(1/100,beta - step_size*grad_beta),1/25);

	total_err_new = (sqrt(mean(diffH .^ 2)) + sqrt(mean(diffV .^ 2)))/2;
	err_diff = total_err_old-total_err_new;
	total_err_old = total_err_new*1;

end

        k = min(max(0.00,k),0.75);
        alpha = min(max(0.0,alpha),2);
        beta = min(max(0,beta ),0.2);

return [k, alpha, beta, total_err_old, itercount];
end

function optim_vod_stages(obsHi, obsVi, leafpot, smc_mod, simtab, k, alpha, beta, maxiter, stop_crit,step_size, nskip)
	stage1 = optim_vod_while(obsHi[1:nskip:end], obsVi[1:nskip:end], leafpot[1:nskip:end], smc_mod[1:nskip:end], simtab[1:nskip:end,:], k, alpha, beta, maxiter, stop_crit, step_size);
	k1, alpha1, beta1, e1, i1 = stage1;
	stage2 = optim_vod_while(obsHi, obsVi, leafpot, smc_mod, simtab, k1, alpha1, beta1, maxiter, stop_crit, step_size);
	return stage2;
end

function optim_vod_stages2(obsHi, obsVi, leafpot, smc_mod, Tsoil, Tcan, laiM, k, alpha, beta, maxiter, stop_crit,step_size, nskip)
	stage1 = optim_vod_while2(obsHi[1:nskip:end], obsVi[1:nskip:end], leafpot[1:nskip:end], smc_mod[1:nskip:end], Tsoil[1:nskip:end], Tcan[1:nskip:end], laiM[1:nskip:end], k, alpha, beta, maxiter, stop_crit, step_size);
	k1, alpha1, beta1, e1, i1 = stage1;
	stage2 = optim_vod_while2(obsHi, obsVi, leafpot, smc_mod, Tsoil, Tcan, laiM, k1, alpha1, beta1, maxiter, stop_crit, step_size);
	return stage2;
end


text:
mutable struct Log <: AbstractLog
  position::Int

  label::AbstractString
  value::Int

  score::Int

  level::Int

  lines::Int
  combo::Int

  stash::Int
  glues::Int

  total::Int

  chain::AbstractString
  piece::Char
end


text:
#=
This is a wrapper for Julia around the Go based C-shared library libstn.go.
=#

import JSON
include("stn.jl")

function test_version()
    """Read in version information from ../codemeta.json and makesure stn.version() returns appropriate string"""
    src = '{"version":""}'
    open("../codemeta.json") do f
        src = readstring(f)
    end
    codemeta = JSON.parse(src)
    expected = "v" * codemeta["version"]
    result = stn.version()
    if expected != result
        println("expected ", expected, " got ", result)
        return 1
    end
    return 0
end

#
# Main processing
#

# Pre-test check
error_count = 0
ok = true

error_count += test_version()

println("Tests completed")

# Summarize our test results
if error_count > 0
    println("Total errors ", error_count)
    exit(1)
end
println("PASS")
println("OK, stn_test.jl")

max_stars_repo_name: rsdoiel/stngo

max_issues_repo_name: rsdoiel/stngo

max_forks_repo_name: rsdoiel/stngo


text:
using DNC
using Flux
using DNC: MemoryAccess, split_ξ, update_state_after_write!, update_state_after_read!, eraseandadd, calcoutput
using DNC: contentaddress, backwardweight, forwardweight, readweight, readweights, writeweight, writeweights
using DNC: usage, memoryretention, allocationweighting

import DNC.readweight, DNC.readweights, DNC.writeweight, DNC.writeweights

function readweight(io::IO, backw, content, forw, readmode)
    rw = DNC.readweight(backw, content, forw, readmode)
    print(io, "readweight:",rw,"\n")
    rw
end

function readweights(io::IO, M, inputs, L, prev_wr)
    k, β, readmode= inputs.kr, inputs.βr, inputs.readmode
    cr = contentaddress(k, M, β)
    b = backwardweight(L, prev_wr)
    f = forwardweight(L, prev_wr)
    wr = readweight(io, b, cr, f, readmode)
end

function writeweight(io::IO, cw, a, gw, ga)
    ww = DNC.writeweight(cw, a, gw, ga)
    print(io, "writeweight:", ww,"\n")
    ww
end

function writeweights(io::IO, M, inputs,
        prev_ww,
        prev_wr,
        prev_usage)
    k, β, ga, gw, e, v, free = inputs.kw, inputs.βw, inputs.ga[1], inputs.gw[1], inputs.e, inputs.v, inputs.f
    cw = contentaddress(k, M, β)
    𝜓 = memoryretention(prev_wr, free)
    u = usage(prev_usage, prev_ww, 𝜓)
    a = allocationweighting(u)
    ww = writeweight(io, cw, a, gw, ga)
end

function (ma::MemoryAccess)(io::IO, inputs)
    R = size(ma.state.wr)[2]
    W = size(ma.M)[2]
    inputs = split_ξ(inputs, ma.inputmaps)
    p, u, ww, wr = ma.state.p, ma.state.u, ma.state.ww, ma.state.wr
    u = usage(u, ww, wr, inputs.f)
    ww= writeweights(io, ma.M, inputs, ww, wr, u)
    ma.M = eraseandadd(ma.M, ww, inputs.e, inputs.v)
    update_state_after_write!(ma.state, ww, u)
    wr = readweights(io, ma.M, inputs, ma.state.L, wr)
    update_state_after_read!(ma.state, wr)
    readvectors = ma.M' * wr
    readvectors
end

function monitoredprediction(io::IO, model, x)
    m = model.cell
    h = m.readvectors
    out = m.controller([x;h])
    v = out[1:m.Y]
    ξ = out[m.Y+1:end]
    r = m.memoryaccess(io, ξ)
    r = reshape(r, size(r)[1]*size(r)[2])
    return calcoutput(v, r, m.Wr)
end

function monitored_loss(model, x, y, mask)
    open("monitor.txt", "a") do io
        ŷ = monitoredprediction(io, model, x)
        loss = sum(mask * Flux.logitbinarycrossentropy.(ŷ, y))
        loss
    end
end


text:
using Thunks
using Test

"Fail test on error without stopping testset. Not totally working..?"
macro safetest(ex)
    return quote
        @test try
            @eval $(esc(ex))
        catch err
            throw(err)
            false
       end
    end
end
@testset "all" begin

@testset "reify" begin
    x = thunk(identity)(2)
    y = thunk(identity)(3)
    z = thunk(+)(x, y)
    @test typeof(x) == typeof(y) == typeof(z) == Thunk
    @test z.evaluated == x.evaluated == y.evaluated == false
    @test z.result == x.result == y.result == nothing
    @test reify(z) == 5
    @test z.evaluated == x.evaluated == y.evaluated == true
    @test z.result == 5
    @test x.result == 2
end

@testset "find symbols in ast" begin
    ex = :([x[1:2] for x in repeat([repeat([a], 3)],3)])
    symbols = Thunks.find_symbols_in_ast(ex)
    @test all([x in symbols for x in [:a, :repeat, :x]])
    ex2 = :(begin
        z = a + 3
        y = collect(1:b)
        abc = x -> rand(3,5)[2:3, 2:3]
    end)
    symbols2 = Thunks.find_symbols_in_ast(ex2)
    @test all([x in symbols2 for x in [:z, :a, :y, :b, :collect, :abc, :rand, :x]])
end

@testset "find assigned symbols in ast" begin
    ex0 = :(collect(1:b; fakekeyword=1, fake=2))
    symbols0 = Thunks.find_assigned_symbols_in_ast(ex0)
    @test all([~in(x, symbols0) for x in [:fakekeyword, :fake]])
    a = 3
    ex = :([x[1:2] for x in repeat([repeat([a], 3)],3)])
    symbols = Thunks.find_assigned_symbols_in_ast(ex)
    @test all([x in symbols for x in [:a, :repeat]])
    @test (~)(:x in symbols)
    b = 5
    ex2 = :(begin
        z = a + 3
        y = collect(1:b, fakekeyword=1)
        abc = x -> rand(3,5)[2:3, 2:3]
    end)
    symbols2 = Thunks.find_assigned_symbols_in_ast(ex2)
    @test all([x in symbols2 for x in [:z, :a, :y, :b, :collect, :abc, :rand]])
    @test all([~in(x, symbols2) for x in [:x, :fakekeyword]])
end

# quick test
# using Thunks; begin x=y=z=1; a = @lazy sum([x,y,z]); reify(a); end

@testset "@lazy" begin
    a = b = c = 2
    t = @lazy sum([a,b,c])
    @test typeof(t) == Thunk
    @test reify(t) == 6
    @lazy begin
        x = 1
        b = 2
        z = 4
        abc = sum([x,b,z])
    end
    @test typeof(abc) == Thunk
    @test reify(abc) == 7
end

function add1(x)
    x + 1
end

@testset "dot broadcast" begin
    a = ones(5)

    @lazy b = a .+ a
    @test all(reify(b) .== (a .*2))

    @lazy c = add1.(a)
    @test all(reify(c) .== (a .+ 1))
end

test() = (1,2,3)

@testset "indexing" begin
    @lazy x = test()[1]
    @test reify(x) == 1
    @lazy y = test()
    y1 = y[1]
    @test reify(y1) == 1
    i = @lazy identity(10)
    x = @lazy collect(1:i)[7:end]
    @test all(reify(x) .== [7,8,9,10])
end

@testset "if" begin
        @lazy x = true ? 1 : 0
        @test reify(x) == 1
        y = @lazy if true
            1
        else
            0
        end
        @test reify(y) == 1
end

@testset "return arity" begin
    @lazy x,y,z = test()
    @test map(reify,(x,y,z)) == (1,2,3)
    @lazy (a,b,c) = test()
    @test map(reify,(x,y,z)) == (1,2,3)
end

function maybe_add1(x; add1=false)
    y = add1 ? x + 1 : x
    y*2
end

@testset "kwargs" begin
    @lazy y = maybe_add1(2; add1=false)
    @test reify(y)==4
    @lazy yy = maybe_add1(2; add1=true)
    @test reify(yy)==6
end

@testset "comprehension Expr" begin
    a = @lazy identity(6)
    z = @lazy [x[1:2] for x in repeat([repeat([a], 3)],3)]
    @test all(sum(reify(z)) .== [18, 18])
end

@testset "Unevaluated" begin
    # will not error
    not_used = Unevaluated(+, ()->(not_defined,1))
    a = Unevaluated(identity, ()->1)
    b = Unevaluated(identity, ()->2)
    c = Unevaluated(+, ()->(a,b))
    @test reify(c) == 3
end

@testset "@noeval" begin
    # will not error
    @noeval begin
        not_used = not_defined + 1
        a = identity(1)
        b = identity(2)
        c = a + b
    end
    # behavior may be unexpected if referential transparency does not hold.
    b = -2
    @test reify(c) == -1
end

@testset "@reversible" begin
    f() = 3
    @reversible begin
        a = f()
        b = identity(2)
        c = a + b
    end
    @test c.evaluated == false
    @test reify(c) == 5
    @test c.evaluated == true
    undo.([a,b,c])
    @test c.evaluated == false
    f() = 1
    @test reify(c) == 3
    @test c.evaluated == true
end

function save_str(path)
    (str) -> open(path, "w") do file
            write(file, str)
        end
end

function read_str(path)
    if isfile(path)
        open(path, "r") do file
            read(file, String)
        end
    else
        nothing
    end
end

function hello_world_thunk()
    @lazy begin
        a = identity("hello")
        b = identity("world")
        c = "$a $b"
    end
    a,b,c
end

@testset "checkpoint" begin
    path, _ = mktemp()
    path *= ".txt"
    a,b,c = hello_world_thunk()
    c = Checkpointable(c, save_str(path), read_str(path))
    @test reify(c) == "hello world"
    @test b.evaluated == true

    a,b,c = hello_world_thunk()
    c = Checkpointable(c, save_str(path), read_str(path))
    @test reify(c) == "hello world"
    @test b.evaluated == false

    p = @lazy mktemp()[1]*".txt"
    a,b,c = hello_world_thunk()
    @lazy begin
        ckpt_f = save_str(p)
        restore_f = read_str(p)
    end
    c = Checkpointable(c, ckpt_f, restore_f)
    @test p.evaluated == false
    @test reify(c) == "hello world"
    @test p.evaluated == true

    a,b,c = hello_world_thunk()
    @lazy begin
        ckpt_f = save_str(p)
        restore_f = read_str(p)
    end
    c = Checkpointable(c, ckpt_f, restore_f)
    @test reify(c) == "hello world"
    @test b.evaluated == false
end

end


text:
struct Embedding
    embedding::AbstractArray{Float64, 2}
end

function random_embedding(npts::Int)
    dist = Normal()
    # Create an uncorrelated source and target
    invariant_embedding_found = false
    while !invariant_embedding_found
        ts = rand(dist, npts, 1)

        embedd = hcat(ts[1:(end - 2)], ts[2:(end - 1)], ts[3:end])

        if invariantset(embedd)
            invariant_embedding_found = true
            return Embedding(embedd)
        end
    end
end

function gaussian_embedding(npts::Int; covariance::Float64 = 0.4)
    dist = Normal()
    # Create an uncorrelated source and target
    invariant_embedding_found = false
    while !invariant_embedding_found
        source = rand(dist, npts, 1)
        dest = covariance .* source[1:end] .+ (1.0 - covariance) .* rand(dist, npts, 1)
        # Embedding
        embedd = hcat(source[1:end-1], source[2:end], dest[2:end])

        if invariantset(embedd)
            invariant_embedding_found = true
            return Embedding(embedd)
        end
    end
end

function gaussian_embedding(npts::Int, covariance::Float64)
    dist = Normal()
    # Create an uncorrelated source and target
    invariant_embedding_found = false
    while !invariant_embedding_found
        source = rand(dist, npts, 1)
        dest = covariance .* source[1:end] .+ (1.0 - covariance) .* rand(dist, npts, 1)
        # Embedding
        embedd = hcat(source[1:end-1], source[2:end], dest[2:end])

        if invariantset(embedd)
            invariant_embedding_found = true
            return Embedding(embedd)
        end
    end
end

function gaussian_embedding_arr(npts::Int; covariance::Float64 = 0.4)
    dist = Normal()
    # Create an uncorrelated source and target
    invariant_embedding_found = false
    while !invariant_embedding_found
        source = rand(dist, npts, 1)
        dest = covariance .* source[1:end] .+ (1.0 - covariance) .* rand(dist, npts, 1)
        # Embedding
        embedd = hcat(source[1:end-1], source[2:end], dest[2:end])

        if invariantset(embedd)
            invariant_embedding_found = true
            return embedd
        end
    end
end

"""
    embed(ts::Vector{Float64}, E::Int, tau::Int)

Embed a time series `ts` in `E` dimensions with embedding lag `tau`.
"""
function embed(ts::Vector{Float64}, E::Int, tau::Int)
  n::Int = length(ts)

  # Initialize
  embedded_ts = zeros(Float64, E, n - ((E - 1) * tau))
  l::Int = n - ((E - 1) * tau)

  for i in 1:E
    start_index = 1 + (i - 1) * tau
    stop_index  = start_index + l - 1
    embedded_ts[i, :] = ts[start_index:stop_index]
  end
  return copy(transpose(embedded_ts))
end

"""
    embed(ts::Vector{Float64}, E::Int, tau::Int)

Embed a time series `ts` in `E` dimensions with embedding lag `tau`.
"""
function embedding(ts::Vector{Float64}, E::Int, tau::Int)
  n::Int = length(ts)

  # Initialize
  embedded_ts = zeros(Float64, E, n - ((E - 1) * tau))
  l::Int = n - ((E - 1) * tau)

  for i in 1:E
    start_index = 1 + (i - 1) * tau
    stop_index  = start_index + l - 1
    embedded_ts[i, :] = ts[start_index:stop_index]
  end
  return Embedding(copy(transpose(embedded_ts)))
end

"""
    embedding_example(n_points::Int, m::Int, tau::Int)

Create an example embedding consisting of `n_points` points in `E` dimension space,
using embedding lag `tau`.
"""
function embedding_example(n_points::Int, E::Int, tau::Int)
  ts = randn(n_points)
  embedding = embed(ts, E, tau)
end

"""
    embedding_example(n_points::Int, m::Int, tau::Int)

Create an example embedding consisting of `n_points` points in `E` dimension space,
using embedding lag `tau`.
"""
function embedding_ex(n_points::Int, E::Int, tau::Int)
  ts = randn(n_points)
  embedding = embedding(ts, E, tau)
end


text:
mutable struct viewer_draw_t <: LCMType
    timestamp::Int64
    num_links::Int32
    link_name::Vector{String}
    robot_num::Vector{Int32}
    position::Matrix{Float32}
    quaternion::Matrix{Float32}
end

@lcmtypesetup(viewer_draw_t,
    link_name => (num_links,),
    robot_num => (num_links,),
    position => (num_links, 3),
    quaternion => (num_links, 4),
)


text:
using Revise
using StrayCopses
using CSV
using StatsBase

using StrayCopses: Node, DecisionTreeContainer, process_node

function prepare_dataset(df, k = 50, n = 3)
    if n == 3
        ids = (1:k, 51:(50+k), 101:(100+k))
    elseif n == 2
        ids = (1:k, 51:(50+k))
    else
        ids = 1:k
    end
    ids = collect(Iterators.flatten(ids))
    X = convert(Matrix, df[ids, 1:size(df, 2) - 1])
    y = convert(Vector, df[ids, size(df, 2)])

    out_classes = unique(y)
    n_classes = length(out_classes)
    in_classes = collect(1:n_classes)
    d = Dict(zip(out_classes, in_classes))
    target = map(z -> d[z], y)
    T = eltype(X)

    return (X = X, y = y, out_classes = out_classes, n_classes = n_classes, target = target, T = T)
end

df = CSV.read(joinpath(@__DIR__, "data", "uci_iris_all.csv"))
X, y, out_classes, n_classes, target, T = prepare_dataset(df);

function build_tree(X, y; max_depth = 1000000, n_features = 4, n_classes = 3)
    T = eltype(X)
    root = Node{T}()
    dtc = DecisionTreeContainer{T}(root, n_features, n_classes, max_depth, 1)
    process_node(dtc, root, X, y)

    return root
end

function lookahead(tree::Node, X)
    res = similar(target)
    for i in axes(target, 1)
        res[i] = predict(tree, X, i)
    end
    res
end

function accuracy(x1, x2)
    sum(x1 .== x2)/length(x1)
end

begin
    tree = build_tree(X, target, max_depth = 1)
    yhat = lookahead(tree, X)
    accuracy(target, yhat)
end

begin
    tree = build_tree(X, target, max_depth = 1_000_000)
    yhat = lookahead(tree, X)
    accuracy(target, yhat)
end

begin
    sc = StrayCopse()
    yhat = predict(md, X)
    md = fit(sc, X, y; max_depth = 3, n_trees = 10)
    accuracy(yhat, y)
end

tree

X0 = [6.5 2.8 4.6 1.5; 6.5 3.0 5.8 2.2]
y0 = [2, 3]

StrayCopses.get_split_indices(X0, 1, 6.5)

StrayCopses.best_split(X0, y0, dtc.n_classes, 1:4)

StrayCopses.feature_best_split(X0, y0, dtc.n_classes, 2)

axes(target)


text:
using YaoToEinsum
using Test, OMEinsum, OMEinsumContractionOrders
using Yao
using Yao.EasyBuild: qft_circuit, variational_circuit, rand_google53
using SymEngine

@testset "YaoToEinsum.jl" begin
    n = 5
    for c in [put(n, 2=>Y), put(n, (5,3)=>SWAP), put(n, (4,2)=>ConstGate.CNOT), put(n, (2,3,1)=>kron(ConstGate.CNOT, X)),
            put(n, 2=>Z), control(n, -3, 2=>X), control(n, 3, 2=>X), control(n, (2, -1), 3=>Y), control(n, (4,1,-2), 5=>Z)]
        @show c
        C = chain([put(n, i=>Rx(rand()*2π)) for i=1:n]..., c)
        code, xs = yao2einsum(C)
        optcode = optimize_code(code, OMEinsumContractionOrders.uniformsize(code, 2), GreedyMethod())
        @test reshape(optcode(xs...; size_info=OMEinsumContractionOrders.uniformsize(code, 2)), 1<<n, 1<<n) ≈ mat(C)
    end
end

@testset "Yao Extensions" begin
    n = 5
    for c in [qft_circuit(n), variational_circuit(n, 2), rand_google53(5; nbits=n)]
        code, xs = yao2einsum(c)
        optcode = optimize_code(code, OMEinsumContractionOrders.uniformsize(code, 2), TreeSA(nslices=3))
        @test reshape(optcode(xs...; size_info=OMEinsumContractionOrders.uniformsize(code, 2)), 1<<n, 1<<n) ≈ mat(c)
    end
end

@testset "boundary conditions" begin
    n = 5
    c = qft_circuit(n)
    initial_state = Dict([i=>rand_state(1) for i=1:n])
    reg = join([initial_state[i] for i=n:-1:1]...)
    reg |> c
    inner = (2,3)
    focus!(reg, inner)
    for final_state in [Dict([i=>rand_state(1) for i in inner]), Dict([i=>1 for i in inner])]
        freg = join(YaoToEinsum.render_single_qubit_state(ComplexF64, final_state[3]), YaoToEinsum.render_single_qubit_state(ComplexF64, final_state[2]))
        code, xs = yao2einsum(c; initial_state=initial_state, final_state=final_state)
        optcode = optimize_code(code, OMEinsumContractionOrders.uniformsize(code, 2), TreeSA(nslices=3))
        @test vec(optcode(xs...; size_info=OMEinsumContractionOrders.uniformsize(code, 2))) ≈ vec(transpose(statevec(freg)) * state(reg))
    end
end

@testset "symbolic" begin
    n = 5
    c = qft_circuit(n)
    initial_state = Dict([i=>zero_state(Basic, 1) for i=1:n])
    code, xs = yao2einsum(c; initial_state=initial_state)
    @test eltype(xs) == AbstractArray{Basic}
end

@testset "fix to basic type" begin
    c = chain(kron(X,X))
    @test yao2einsum(c)[1] isa EinCode
end


text:
using LinearAlgebra
using PyPlot
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCellDG
include("../../../../test/useful_routines.jl")
include("transformation-elasticity-solver.jl")

function angular_position(points)
    cpoints = points[1, :] + im * points[2, :]
    return rad2deg.(angle.(cpoints))
end

function stress_error_vs_angular_position(polyorder, nelmts, penaltyfactor)
    K1, K2 = 247.0, 192.0    # Pa
    mu1, mu2 = 126.0, 87.0   # Pa
    lambda1 = lame_lambda(K1, mu1)
    lambda2 = lame_lambda(K2, mu2)
    stiffness = CutCellDG.HookeStiffness(lambda1, mu1, lambda2, mu2)

    theta0 = -0.067
    transfstress =
        CutCellDG.plane_strain_transformation_stress(lambda1, mu1, theta0)

    meshwidth = [1.0, 1.0]
    numqp = required_quadrature_order(polyorder) + 2

    dx = minimum(meshwidth) / nelmts
    penalty = penaltyfactor / dx * 0.5 * (lambda1 + mu1 + lambda2 + mu2)

    basis = TensorProductBasis(2, polyorder)
    interfacecenter = [0.0, 0.0]
    interfaceradius = 0.5
    outerradius = 2.0
    analyticalsolution = AnalyticalSolution(
        interfaceradius,
        outerradius,
        interfacecenter,
        lambda1,
        mu1,
        lambda2,
        mu2,
        theta0,
    )

    mesh, cellquads, facequads, interfacequads = construct_mesh_and_quadratures(
        meshwidth,
        nelmts,
        basis,
        interfacecenter,
        interfaceradius,
        numqp,
        tinyratio = 0.2,
    )

    nodaldisplacement = nodal_displacement(
        mesh,
        basis,
        cellquads,
        facequads,
        interfacequads,
        stiffness,
        theta0,
        analyticalsolution,
        penalty,
    )

    refseedpoints, spatialseedpoints, seedcellids =
        CutCellDG.seed_zero_levelset_with_interfacequads(interfacequads, mesh)

    normals = CutCellDG.collect_interface_normals(interfacequads, mesh)
    tangents = CutCellDG.rotate_90(normals)

    spatialpoints = spatialseedpoints[1, :, :]
    referencepoints = refseedpoints
    referencecellids = seedcellids

    relspatialpoints = spatialpoints .- interfacecenter
    angularposition = angular_position(relspatialpoints)
    sortidx = sortperm(angularposition)
    angularposition = angularposition[sortidx]

    referencepoints = referencepoints[:, :, sortidx]
    referencecellids = referencecellids[:, sortidx]
    spatialpoints = spatialpoints[:, sortidx]
    normals = normals[:, sortidx]
    tangents = tangents[:, sortidx]

    productstress = product_stress_at_reference_points(
        nodaldisplacement,
        basis,
        stiffness,
        transfstress,
        theta0,
        referencepoints,
        referencecellids,
        mesh,
    )
    parentstress = parent_stress_at_reference_points(
        nodaldisplacement,
        basis,
        stiffness,
        referencepoints,
        referencecellids,
        mesh,
    )
    exactparentstress =
        mapslices(x -> core_stress(analyticalsolution), spatialpoints, dims = 1)
    exactproductstress = mapslices(
        x -> shell_stress(analyticalsolution, x),
        spatialpoints,
        dims = 1,
    )

    productradialtraction =
        CutCellDG.traction_force_at_points(productstress, normals)
    parentradialtraction =
        CutCellDG.traction_force_at_points(parentstress, normals)

    productradialstress = vec(sum(productradialtraction .* normals, dims = 1))
    parentradialstress = vec(sum(parentradialtraction .* normals, dims = 1))

    exactproductradialstress =
        shell_radial_stress(analyticalsolution, interfaceradius)
    exactparentradialstress = core_in_plane_stress(analyticalsolution)

    productradialstresserror =
        abs.(productradialstress .- exactproductradialstress) ./
        abs(exactproductradialstress)
    parentradialstresserror =
        abs.(parentradialstress .- exactparentradialstress) ./
        abs(exactparentradialstress)

    productcircumferentialtraction =
        CutCellDG.traction_force_at_points(productstress, tangents)
    productcircumferentialstress =
        vec(sum(productcircumferentialtraction .* tangents, dims = 1))

    exactproductcircumferentialstress =
        shell_circumferential_stress(analyticalsolution, interfaceradius)
    productcircumferentialstresserror =
        abs.(
            productcircumferentialstress .- exactproductcircumferentialstress,
        ) ./ abs(exactproductcircumferentialstress)

    return angularposition,
    parentradialstresserror,
    productradialstresserror,
    productcircumferentialstresserror

end

polyorder = 3
results9 = stress_error_vs_angular_position(polyorder, 9, 1e2)
results17 = stress_error_vs_angular_position(polyorder, 17, 1e2)
results33 = stress_error_vs_angular_position(polyorder, 33, 1e2)

folderpath = "examples/traction-continuity-with-transformation/cylindrical-particle/cylindrical-bc/"
filelabel = "polyorder-"*string(polyorder)*"-"

fig, ax = PyPlot.subplots()
ax.plot(results9[1], results9[2], label = "9 x 9")
ax.plot(results17[1], results17[2], label = "17 x 17")
ax.plot(results33[1], results33[2], label = "33 x 33")
ax.legend()
ax.set_ylim(0,0.005)
ax.grid()
ax.set_xlabel("Angular position (deg)")
ax.set_ylabel("Normalized error in parent radial stress")
fig.tight_layout()
fig
fig.savefig(folderpath*filelabel*"parent-radial-stress-error.png")

fig, ax = PyPlot.subplots()
ax.plot(results9[1], results9[3], label = "9 x 9")
ax.plot(results17[1], results17[3], label = "17 x 17")
ax.plot(results33[1], results33[3], label = "33 x 33")
ax.legend()
ax.grid()
ax.set_ylim(0,0.005)
ax.set_xlabel("Angular position (deg)")
ax.set_ylabel("Normalized error in product radial stress")
fig.tight_layout()
fig
fig.savefig(folderpath*filelabel*"product-radial-stress-error.png")

fig, ax = PyPlot.subplots()
ax.plot(results9[1], results9[4], label = "9 x 9")
ax.plot(results17[1], results17[4], label = "17 x 17")
ax.plot(results33[1], results33[4], label = "33 x 33")
ax.legend()
ax.grid()
ax.set_ylim(0,0.005)
ax.set_xlabel("Angular position (deg)")
ax.set_ylabel("Normalized error in product circumferential stress")
fig.tight_layout()
fig
fig.savefig(folderpath*filelabel*"product-circumferential-stress-error.png")


text:
using LinQuadOptInterface

const LQOI = LinQuadOptInterface
const MOI  = LQOI.MathOptInterface

const SUPPORTED_OBJECTIVES = [
    LQOI.Linear
    LQOI.SinVar
]

const SUPPORTED_CONSTRAINTS = [
    (LQOI.Linear, LQOI.EQ),
    (LQOI.Linear, LQOI.LE),
    (LQOI.Linear, LQOI.GE),
    (LQOI.SinVar, LQOI.EQ),
    (LQOI.SinVar, LQOI.LE),
    (LQOI.SinVar, LQOI.GE),
    (LQOI.SinVar, LQOI.IV),
    (LQOI.SinVar, MOI.ZeroOne),
    (LQOI.SinVar, MOI.Integer),
    (LQOI.VecVar, MOI.Nonnegatives),
    (LQOI.VecVar, MOI.Nonpositives),
    (LQOI.VecVar, MOI.Zeros),
    (LQOI.VecVar, MOI.SOS1{Float64}),
    (LQOI.VecVar, MOI.SOS2{Float64}),
    (LQOI.VecLin, MOI.Nonnegatives),
    (LQOI.VecLin, MOI.Nonpositives),
    (LQOI.VecLin, MOI.Zeros)
]

mutable struct Optimizer <: LQOI.LinQuadOptimizer
    LQOI.@LinQuadOptimizerBase
    env::Union{Nothing, Env}
    params::Dict{String, Any}

    """
        Optimizer(env = nothing; kwargs...)

    Create a new Optimizer object.

    You can share CPLEX `Env`s between models by passing an instance of `Env` as
    the first argument. By default, a new environment is created for every
    model.
    """
    function Optimizer(env::Union{Nothing, Env} = nothing; kwargs...)
        model = new()
        model.env = env
        # For consistency with MPB, output logs to stdout by default.
        model.params = Dict{String, Any}("CPX_PARAM_SCRIND" => 1)
        for (name, value) in kwargs
            model.params[string(name)] = value
        end
        MOI.empty!(model)
        return model
    end
end

# The existing env is `Nothing`, so create a new one.
LQOI.LinearQuadraticModel(::Type{Optimizer}, ::Nothing) = Model(Env())
# The existing env is `Env`, so pass it through.
LQOI.LinearQuadraticModel(::Type{Optimizer}, env::Env) = Model(env)

function MOI.empty!(model::Optimizer)
    MOI.empty!(model, model.env)
    for (name, value) in model.params
        set_param!(model.inner.env, name, value)
    end
    return
end

MOI.get(::Optimizer, ::MOI.SolverName) = "CPLEX"

LQOI.supported_constraints(::Optimizer) = SUPPORTED_CONSTRAINTS
LQOI.supported_objectives(::Optimizer)  = SUPPORTED_OBJECTIVES

LQOI.backend_type(model::Optimizer, ::MOI.EqualTo{Float64})     = Cchar('E')
LQOI.backend_type(model::Optimizer, ::MOI.LessThan{Float64})    = Cchar('L')
LQOI.backend_type(model::Optimizer, ::MOI.GreaterThan{Float64}) = Cchar('G')

LQOI.backend_type(model::Optimizer, ::MOI.Zeros)                = Cchar('B')
LQOI.backend_type(model::Optimizer, ::MOI.Nonpositives)         = Cchar('L')
LQOI.backend_type(model::Optimizer, ::MOI.Nonnegatives)         = Cchar('E')

function LQOI.change_variable_bounds!(model::Optimizer,
        columns::Vector{Int}, values::Vector{Float64}, senses::Vector{Cchar})
    c_api_chgbds(model.inner, ivec(columns), senses, values)
    return
end

function LQOI.get_variable_lowerbound(model::Optimizer, column::Int)
    return c_api_getlb(model.inner, Cint(column), Cint(column))[1]
end

function LQOI.get_variable_upperbound(model::Optimizer, column::Int)
    return c_api_getub(model.inner, Cint(column), Cint(column))[1]
end

function LQOI.get_number_linear_constraints(model::Optimizer)
    return c_api_getnumrows(model.inner)
end

function LQOI.add_linear_constraints!(model::Optimizer,
        A::LQOI.CSRMatrix{Float64}, sense::Vector{Cchar}, rhs::Vector{Float64})
    c_api_addrows(model.inner, ivec(A.row_pointers), ivec(A.columns),
        A.coefficients, sense, rhs)
    return
end

function LQOI.get_rhs(model::Optimizer, row::Int)
    rhs = Vector{Cdouble}(undef, 1)
    c_api_getrhs(model.inner, rhs, Cint(row), Cint(row))
    return rhs[1]
end

function LQOI.get_linear_constraint(model::Optimizer, row::Int)
    (nzcnt, rmatbeg, rmatind, rmatval) =
            c_api_getrows(model.inner, Cint(row), Cint(row))
    return rmatind[1:nzcnt], rmatval[1:nzcnt]
end

function LQOI.change_matrix_coefficient!(
        model::Optimizer, row::Int, col::Int, coef::Float64)
    c_api_chgcoef(model.inner, Cint(row), Cint(col), coef)
    return
end

function LQOI.change_objective_coefficient!(
        model::Optimizer, col::Int, coef::Float64)
    c_api_chgobj(model.inner, [Cint(col)], [coef])
    return
end

function LQOI.change_rhs_coefficient!(
        model::Optimizer, row::Int, coef::Float64)
    c_api_chgrhs(model.inner, [Cint(row)], [coef])
    return
end

function LQOI.delete_linear_constraints!(
        model::Optimizer, first_row::Int, last_row::Int)
    c_api_delrows(model.inner, Cint(first_row), Cint(last_row))
    return
end

function LQOI.change_variable_types!(
        model::Optimizer, columns::Vector{Int}, vtypes::Vector{Cchar})
    c_api_chgctype(model.inner, ivec(columns), vtypes)
    return
end

function LQOI.change_linear_constraint_sense!(
        model::Optimizer, rows::Vector{Int}, senses::Vector{Cchar})
    c_api_chgsense(model.inner, ivec(rows), senses)
    return
end

function LQOI.set_linear_objective!(
        model::Optimizer, columns::Vector{Int}, coefficients::Vector{Float64})
    n = num_var(model.inner)
    all_coefs = zeros(Float64, n)
    for (col, coef) in zip(columns, coefficients)
        all_coefs[col] += coef
    end
    c_api_chgobj(model.inner, Cint[1:n;], all_coefs)
    return
end

function LQOI.change_objective_sense!(model::Optimizer, symbol)
    if symbol == :min
        c_api_chgobjsen(model.inner, Cint(1))
    else
        c_api_chgobjsen(model.inner, Cint(-1))
    end
    return
end

function LQOI.get_linear_objective!(model::Optimizer, dest)
    c_api_getobj(model.inner, dest, Cint(1), c_api_getnumcols(model.inner))
    return
end

function LQOI.get_objectivesense(model::Optimizer)
    s = c_api_getobjsen(model.inner)
    if s == 1
        return MOI.MIN_SENSE
    else
        return MOI.MAX_SENSE
    end
end

function LQOI.get_number_variables(model::Optimizer)
    return c_api_getnumcols(model.inner)
end

function LQOI.add_variables!(model::Optimizer, N::Int)
    add_vars!(model.inner, zeros(Float64, N), fill(-Inf, N), fill(Inf, N))
    return
end

function LQOI.delete_variables!(model::Optimizer, first_col::Int, last_col::Int)
    c_api_delcols(model.inner, Cint(first_col), Cint(last_col))
    return
end

function LQOI.solve_mip_problem!(model::Optimizer)
    model.inner.has_int = true
    LQOI.solve_linear_problem!(model)
    return
end

function LQOI.solve_linear_problem!(model::Optimizer)
    optimize!(model.inner)
    return
end

function LQOI.get_termination_status(model::Optimizer)
    stat = c_api_getstat(model.inner)

    if stat == 1 # CPX_STAT_OPTIMAL
        return MOI.OPTIMAL
    elseif stat == 3 # CPX_STAT_INFEASIBLE
        return MOI.INFEASIBLE
    elseif stat == 4 # CPX_STAT_INForUNBD
        return MOI.INFEASIBLE_OR_UNBOUNDED
    elseif stat == 2 # CPX_STAT_UNBOUNDED
        return MOI.DUAL_INFEASIBLE
    elseif stat in (12, 21, 22, 36) # CPX_STAT_*ABORT*_OBJ_LIM
        return MOI.OBJECTIVE_LIMIT
    elseif stat in (10, 34) # CPX_STAT_*ABORT_IT_LIM
        return MOI.ITERATION_LIMIT
    elseif stat == 53 # CPX_STAT_CONFLICT_ABORT_NODE_LIM
        return MOI.NODE_LIMIT
    elseif stat in (11, 25, 33, 39) # CPX_STAT_*ABORT*TIME_LIM
        return MOI.TIME_LIMIT
    elseif stat == 5 # CPX_STAT_OPTIMAL_INFEAS
        return MOI.NUMERICAL_ERROR
    # MIP STATUS
    elseif stat in (101, 102) # CPXMIP_OPTIMAL, CPXMIP_OPTIMAL_TOL
        return MOI.OPTIMAL
    elseif stat == 103 # CPXMIP_INFEASIBLE
        return MOI.INFEASIBLE
    elseif stat == 119 # CPXMIP_INForUNBD
        return MOI.INFEASIBLE_OR_UNBOUNDED
    elseif stat == 118 # CPXMIP_UNBOUNDED
        return MOI.DUAL_INFEASIBLE
    elseif stat in (105, 106) # CPXMIP_NODE_LIM*
        return MOI.NODE_LIMIT
    elseif stat in (107, 108, 131, 132) # CPXMIP_*TIME_LIM*
        return MOI.TIME_LIMIT
    else
        return MOI.OTHER_ERROR
    end
end

function LQOI.get_primal_status(model::Optimizer)
    soln_method, soln_type, primal_stat, dual_stat = c_api_solninfo(model.inner)
    if primal_stat == 1
        return MOI.FEASIBLE_POINT
    else
        return MOI.NO_SOLUTION
    end
end

function LQOI.get_dual_status(model::Optimizer)
    if model.inner.has_int
        return MOI.NO_SOLUTION
    end
    soln_method, soln_type, primal_stat, dual_stat = c_api_solninfo(model.inner)
    if dual_stat == 1
        return MOI.FEASIBLE_POINT
    else
        return MOI.NO_SOLUTION
    end
end

function LQOI.get_variable_primal_solution!(model::Optimizer, dest)
    c_api_getx(model.inner, dest)
    return
end

function LQOI.get_linear_primal_solution!(model::Optimizer, dest)
    c_api_getax(model.inner, dest)
    return
end

function LQOI.get_variable_dual_solution!(model::Optimizer, dest)
    c_api_getdj(model.inner, dest)
    return
end

function LQOI.get_linear_dual_solution!(model::Optimizer, dest)
    c_api_getpi(model.inner, dest)
    return
end

function LQOI.get_objective_value(model::Optimizer)
    return c_api_getobjval(model.inner)
end

function LQOI.get_objective_bound(model::Optimizer)
    return get_best_bound(model.inner)
end

function LQOI.get_farkas_dual!(model::Optimizer, dest)
    copy!(dest, get_infeasibility_ray(model.inner))
    return
end

function LQOI.get_unbounded_ray!(model::Optimizer, dest)
    copy!(dest, get_unbounded_ray(model.inner))
    return
end

function LQOI.add_sos_constraint!(model::Optimizer, columns::Vector{Int},
        weights::Vector{Float64}, sos_type::Symbol)
    add_sos!(model.inner, sos_type, columns, weights)
    return
end

function LQOI.get_sos_constraint(model::Optimizer, index::Int)
    (cols, weights, sos_type) = c_api_getsos(model.inner, index - 1)
    if sos_type == Cchar('1')
        return (cols .+ 1, weights, :SOS1)
    else
        @assert sos_type == Cchar('2')
        return (cols .+ 1, weights, :SOS2)
    end
    return
end

function LQOI.delete_sos!(model::Optimizer, begin_index::Int, end_index::Int)
    c_api_delsos(model.inner, begin_index - 1, end_index - 1)
    return
end


text:
using Distributions, StatsBase, Plots, LaTeXStrings; pyplot()

mu0, mu1a, mu1b, mu1c, sd = 15, 16, 18, 20, 2
tauGrid = 5:0.1:25

dist0 = Normal(mu0,sd)
dist1a, dist1b, dist1c  = Normal(mu1a,sd), Normal(mu1b,sd), Normal(mu1c,sd)

falsePositive = ccdf.(dist0,tauGrid)
truePositiveA, truePositiveB, truePositiveC =
    ccdf.(dist1a,tauGrid), ccdf.(dist1b,tauGrid), ccdf.(dist1c,tauGrid)

plot(falsePositive, [truePositiveA truePositiveB truePositiveC],
    c=[:blue :red :green],
    label=[L"H1a: \mu_1 = 16" L"H1b: \mu_1 = 18" L"H1c: \mu_1 = 20"])
plot!([0,1], [0,1], c=:black, ls=:dash, label="H0 = H1 = 15",
    xlims=(0,1), ylims=(0,1), xlabel=L"\alpha", ylabel="Power",
    ratio=:equal, legend=:bottomright)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
struct PressureTrace
    t::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}
    timestamp::DateTime
    data::DataFrame
end

function PressureTrace(filepath, filter=(median_size = 15, lowpass_mul=100, filt=Butterworth(2)))
    open(filepath) do f
        header = Dict(Pair(split(readline(f), '\t', limit=3)[1:2]...) for i = 1:23)

        t = range(parse(Float64, header["X0"]),
        length = parse(Int, header["Samples"]),
          step = parse(Float64, header["Delta_X"]))

        date = Date(header["Date"], "yyyy/mm/dd")
        time = Time(first(split(header["Time"], '.')))
        timestamp = DateTime(date, time)

        seekstart(f)
        data = CSV.File(f, datarow=24, drop=[1], delim='\t',
            header="PT" .* string.(0:12)) |> DataFrame

        if !isnothing(filter)
            freq = 1/step(t)
            lowpass = digitalfilter(Lowpass(freq/filter.lowpass_mul, fs=freq), filter.filt)
            mapcols!(data) do s
                s̄ = mapwindow(median, s, filter.median_size)
                ŝ = filtfilt(lowpass, s)
                [median((s[i], s̄[i], ŝ[i])) for i in eachindex(s)]
            end
        end
        PressureTrace(t, timestamp, data)
    end
end

"""
    xt_data(ptrace, ptlocs, driver, driven, trigger, ref_PT)

Calculate the time of shock arrival at each pressure transducer.
Uses the driven gas properties and the Δx/Δt between the trigger
and reference pressure transducers to calculate an approximate Mach
number, which is used to determine driver conditions and shocked
gas pressure for determination of shock arrival time at each transducer.

    Arguments:
    - `ptrace::PressureTrace`
    - `ptlocs::DataFrame`
    - `driver::Fluid`, e.g. `Species("Helium")`
    - `driven::Fluid`, e.g. `Mixture(["Helium", "Acetone"], zs=[0.95, 0.05])`
    - `trigger::Pair{Symbol, Unitful.Quantity}`, e.g. `:PT3 => 10u"psi"`
    - `ref_PT::Symbol`, e.g. `:PT6`

# Examples
```julia-repl
julia> ptlocs = CSV.File("PT_locations.csv") |> DataFrame
12×3 DataFrame
 Row │ name    x_m       σ_m
     │ String  Float64   Float64
─────┼─────────────────────────────-
   1 │ PT1     0.506412  0.00254
   2 │ PT2     1.36207   0.0035921
   3 │ PT3     1.8415    0.0035921
   4 │ PT4     2.05422   0.0035921
   5 │ PT5     2.33997   0.00299529
   6 │ PT6     3.84492   0.00392726
   7 │ PT7     4.49897   0.00392726
   8 │ PT8     5.10222   0.00392726
   9 │ PT9     5.42607   0.00392726
  10 │ PT10    5.98329   0.00392726
  11 │ PT11    6.50716   0.00392726
  12 │ PT12    6.75481   0.00392726

julia> xt_data(PressureTrace("run3/ptrace.lvm"), ptlocs,
        Species("N2"), Species("Ar"), :PT3 => 10u"psi", :PT6)
((driver = Species(N2, 298.1 K, 4.275e+06 Pa),
  driven = Species(Ar, 298.1 K, 1.013e+05 Pa),
 shocked = Species(Ar, 712.5 K, 6.039e+05 Pa),
      u2 = 429.3754150350808 m s^-1),
      12×2 DataFrame
 Row │ x         t_shock
     │ Float64   Float64
─────┼────────────────────
   1 │ 0.506412  0.008039
   2 │ 1.36207   0.009323
   3 │ 1.8415    0.010003
   4 │ 2.05422   0.010307
   5 │ 2.33997   0.010695
   6 │ 3.84492   0.012797
   7 │ 4.49897   0.013711
   8 │ 5.10222   0.014637
   9 │ 5.42607   0.015019
  10 │ 5.98329   0.015796
  11 │ 6.50716   0.016535
  12 │ 6.75481   0.016881)

```
"""
function xt_data(ptrace::PressureTrace, ptlocs::DataFrame, driver::Fluid, driven::Fluid, trigger::Pair, ref_PT::Symbol)
    trigger_PT, trigger_thresh = trigger

    # Get PT indices (PT1 => 1, PT9 => 9, etc) within ptlocs
    trigger_idx = parse(Int, string(trigger_PT)[3:end])
    ref_idx = parse(Int, string(ref_PT)[3:end])

    # find trigger time index
    trigger_sensor = >(ustrip(trigger_thresh |> u"psi"))
    trigger_ptrace_idx = findfirst(trigger_sensor, ptrace.data[!, trigger_PT])

    ref_ptrace_idx = findfirst(trigger_sensor, ptrace.data[!, ref_PT])
    Δt = ptrace.t[ref_ptrace_idx] - ptrace.t[trigger_ptrace_idx]
    Δx = ptlocs[ref_idx, :x_m] - ptlocs[trigger_idx, :x_m]

    # calculate Mach number between trigger and reference PTs
    W₀ = (Δx / Δt)*u"m/s"
    M₀ = W₀/soundspeed(driven)

    # determine gas states
    states = shockcalc(driver, driven, M₀)
    shocked_psi = pressure(states.shocked) |> u"psi" |> ustrip
    t_shock = [isnothing(i) ? NaN : ptrace.t[i] for i in findfirst.(>(0.5*shocked_psi), eachcol(ptrace.data))]

    return (;ptrace, states, shock=(;W₀, M₀), xt=DataFrame(:x => ptlocs[!, :x_m], :t_shock => t_shock))
end

"""
An alternate `ptrace` syntax allows for Base Julia-only arguments, namely filepaths, strings, and symbols.

# Examples
```julia-repl
julia> xt_data("run1/ptrace.lvm", "PTlocations.csv", "N2", "Ar", :PT3 => 10, :PT6)
((driver = Species(N2, 298.1 K, 1.633e+06 Pa),
  driven = Species(Ar, 298.1 K, 1.013e+05 Pa),
 shocked = Species(Ar, 563.7 K, 4.085e+05 Pa),
      u2 = 316.0718112278267 m s^-1),
  12×2 DataFrame
   Row │ x         t_shock
       │ Float64   Float64
  ─────┼────────────────────
     1 │ 0.506412  0.007734
     2 │ 1.36207   0.009218
     3 │ 1.8415    0.010002
     4 │ 2.05422   0.010372
     5 │ 2.33997   0.010854
     6 │ 3.84492   0.013367
     7 │ 4.49897   0.014464
     8 │ 5.10222   0.023391
     9 │ 5.42607   0.016024
    10 │ 5.98329   0.016962
    11 │ 6.50716   0.017844
    12 │ 6.75481   0.018261)
```
"""
function xt_data(ptrace_path, ptloc_path, drivergas, drivengas, trigger, ref_PT;
    ptrace_filter = (median_size = 15, lowpass_mul=100, filt=Butterworth(2)))

    ptrace = PressureTrace(ptrace_path, ptrace_filter)
    ptlocs = CSV.File(ptloc_path) |> DataFrame
    driver = Species(drivergas)
    driven = Species(drivengas)
    trigger_psi = first(trigger) => last(trigger)*u"psi"
    xt_data(ptrace, ptlocs, driver, driven, trigger_psi, ref_PT)
end


text:
"""
    generate_overload(sig)

Takes a signature tuple type, for a primal function that has an `rrule` and generates
appropriate overloads for Nabla's `Node` types to allow performing AD.
This is the hook function for `ChainRulesCore.on_new_rule(hook, rrule)`.

For example, if `generate_overload` is called on `Tuple{typeof{identity}, Any}`
then approximately the following code is `@eval`ed:

```julia
function Base.identity(x1::Node{<:Any}; kwargs...)
    args = (x1,)
    (primal_val, pullback) = rrule(op, unbox.(args)...; kwargs...)
    tape = get_tape(args)
    branch = Branch(primal_val, op, args, kwargs.data, tape, length(tape) + 1, pullback)
    push!(tape, branch)
    return branch
end

@inline function preprocess(
    op::typeof(identity), y::Branch, ȳ, x1::Union{Any, Node{<:Any}}
)
    pullback = getfield(y, :pullback)
        @assert pullback !== nothing "pullback not set, ..."
    return pullback(ȳ)
end

@inline function ∇(
    op::typeof(identity), ::Type{Arg{N}}, p, ::Any, ::Any, x1::Union{Any, Node{<:Any}};
    kwargs...
) where N
    return p[N + 1]  # skip dself (N==1) as we don't support functors
end
```

The real code evaluated is a little more complex with macro-hygiene and handling for
various complicated type-signatures, including multiple arguments.

It does not generate any code for `rrules` for primal functions that Nabla does not support.
See [`should_use_rrule`](@ref) for more details on what rules we do not use.

This function returns true or false as to whether or not code was generated. While this has
no actual effect in itself, it can be useful for checking how many rules Nabla supports.
"""
function generate_overload(sig)
    should_use_rrule(sig) || return false

    original_signature_def = ExprTools.signature(sig; extra_hygiene=true)
    unionized_signature_def = copy(original_signature_def)
    unionized_signature_def[:args] = unionise_sig.(original_signature_def[:args])

    fdef = quote
        @inline $(preprocess_declaration(unionized_signature_def))
        @inline $(∇_declaration(unionized_signature_def))
        $(overload_declarations!(original_signature_def)...)
    end
    eval(fdef)

    return true
end

"""
    should_use_rrule(sig)

Should we make use of the chainrules `rrule` for the primal function with the given
signature tuple type (`sig`).

We do not use rules for:
    - builtin functions
    - functors / closures
    - functions without any positional arguments
    - functions from the `NaNMath` module
    - functions for working with complex numbers.
    - Non-differentiable functions that we define directly on `Node`s better (like `size`)
    - Non-differentiable functions that are never used in practice and that cause a lot of
      compiler invalidations and so cause a large increase in loading time.

Finally this excludes function that at time of last update Nabla had its own rules for
because ChainRules didn't support them.
Generally, for this category once they are added to ChainRules, we should change to using
them from there. This requires also deleting the code from Nabla that provides those rules
currently, so that there is no clash.
"""
function should_use_rrule(sig)
    opT, argTs = Iterators.peel(ExprTools.parameters(sig))
    opT <: Core.Builtin && return false  # can't do operator overloading for builtins

    isabstracttype(opT) || fieldcount(opT) == 0 || return false  # not handling functors
    isempty(argTs) && return false  # we are an operator overloading AD, need operands

    # Don't care about NaNMath
    opT isa DataType && nameof(opT.name.module) == :NaNMath  && return false

    # Ignore functions that have complex ranges. This may change when Nabla supports complex
    # numbers.
    opT ∈ typeof.((
        SpecialFunctions.hankelh1, SpecialFunctions.hankelh2,
        log1p, rem2pi, mod, atan, rem,
    ))  && return false
    opT <: Type{<:Complex} && return false  # skip complex constructor

    # Ignore these functions because they have better Nabla specific versions.
    opT ∈ typeof.((
        isapprox, axes, size, length, isassigned, one, zero,
        Base.Broadcast.combine_styles,  #TODO should i keep this?
    )) && return false

    # Ignore these functions because in practice they are never used and defining them cause
    # a ton of compiler invalidations, making loading slow.
    opT ∈ typeof.((
        string, repr, print, println, write, readlines, eachline, Core.print, Core.println,
        isequal, ==, in, haskey,
        isnothing, ismissing, isfile,
        isbitstype, isbits, isabstracttype, isconcretetype,
        startswith, endswith, join, joinpath, normpath, chomp,
        schedule,  # this one is huge, causes over 2500 invalidations
    )) && return false

    # Rules currently implemented directly in Nabla, but that could use ChainRules in future
    sig <: Union{
        Tuple{typeof(+),AbstractArray,LinearAlgebra.UniformScaling},
        Tuple{typeof(+),LinearAlgebra.UniformScaling,AbstractArray},
        Tuple{typeof(/),Number,AbstractArray},
        Tuple{typeof(LinearAlgebra.BLAS.symm),Char,Char,AbstractArray,AbstractArray},
        Tuple{typeof(LinearAlgebra.BLAS.symm),Char,Char,Number,AbstractArray,AbstractArray},
        Tuple{typeof(LinearAlgebra.BLAS.symv),Char,AbstractArray,AbstractArray},
        Tuple{typeof(LinearAlgebra.BLAS.symv),Char,Number,AbstractArray,AbstractArray},
        Tuple{typeof(LinearAlgebra.BLAS.trmm),Char,Char,Char,Char,Number,AbstractArray,AbstractArray},
        Tuple{typeof(LinearAlgebra.BLAS.trmv),Char,Char,Char,AbstractArray,AbstractArray},
        Tuple{typeof(LinearAlgebra.BLAS.trsm),Char,Char,Char,Char,Number,AbstractArray,AbstractArray},
        Tuple{typeof(LinearAlgebra.BLAS.trsv),Char,Char,Char,AbstractArray,AbstractArray},
        Tuple{typeof(Statistics.mean),Function,AbstractArray},
        Tuple{typeof(\),AbstractArray,Number},
        Tuple{typeof(broadcast),Any,Vararg},
        Tuple{typeof(copy),Any},
        Tuple{typeof(float),Any},
        Tuple{typeof(getindex),Ref},
        Tuple{typeof(kron),AbstractArray,AbstractArray},
        Tuple{typeof(map),Function,Vararg},
        Tuple{typeof(mapfoldl),Any,Union{typeof(+), typeof(Base.add_sum)},Union{Number,AbstractArray}},
        Tuple{typeof(mapfoldr),Any,Union{typeof(+), typeof(Base.add_sum)},Union{Number,AbstractArray}},
        Tuple{typeof(mapreduce),Any,Union{typeof(+), typeof(Base.add_sum)},AbstractArray},
        Tuple{typeof(sum),Function,AbstractArray},
        Tuple{typeof(sum),typeof(abs2),AbstractArray},
    } && return false

    return true  # no exclusion applies
end
"""
    overload_declarations!(original_signature_def)

Given a `signature_def` dictionary as returned by `ExprTools.signature` this returns
the ASTs for the overloads of the primal functions to accept `Nabla.Node`s.
The `signature_def` should *not* have been unionized, as this function will instead generate
1 method for each position a node could be in.
Note: this mutate `signature_def` and so should not be called if others functions also need
to use it after.
"""
function overload_declarations!(signature_def)
    # Our manual macro-hygiene is not complete here.
    # the argument names and `op`, `tape` `args`, `kwargs` etc could conflict with
    # where-params. but for sake of outputting readable code we are not gensyming everything
    # chance of conflict seems low as where-params are normally upper-case.
    @assert(signature_def[:name].head == :(::))
    @assert(signature_def[:name].args[1] == :op)

    original_signature_args = signature_def[:args]
    signature_def[:kwargs] = [:(kwargs...)]
    signature_def[:body] = quote
        args = $(ExprTools.args_tuple_expr(original_signature_args))
        # uncommenting the below to is useful for debugging what rrule is being hit.
        # @show InteractiveUtils.@which rrule(op, unbox.(args)...)
        primal_val, pullback = rrule(op, unbox.(args)...; kwargs...)
        tape = get_tape(args)

        branch = Branch(primal_val, op, args, kwargs.data, tape, length(tape) + 1, pullback)
        push!(tape, branch)
        return branch
    end

    # we need to generate a version of this for each place that an arg could be a Node
    n_args = length(original_signature_args)
    definitions = Expr[]
    for swap_mask in Iterators.product(ntuple(_->(true, false), n_args)...)
        any(swap_mask) || continue  # don't generate if not swapping anything.
        signature_def[:args] = map(swap_mask, original_signature_args) do swap, orig_arg
            if swap
                @assert Meta.isexpr(orig_arg, :(::), 2)
                Expr(:(::), orig_arg.args[1], node_type(orig_arg.args[2]))
            else
                orig_arg
            end
        end
        push!(definitions, ExprTools.combinedef(signature_def))
    end

    return definitions
end

"""
    preprocess_declaration(unionized_signature_def)

Generates AST for overloads for [`Nabla.preprocess`](@ref) that will call the pullback
stored on the `Branch`.
Roughly speaking generated code like:
`preprocess(f::opT, y::Branch, ȳ, xs...)) = y.pullback(ȳ)`
We need the pullback value to use to compute the sensitivies of the primal inputs, that will
be queries by `∇(::opT, ::Type{Arg{N}}, p, y, ȳ, xs...)` where `p` is that pullback value
return by the `preprocess` function.

Note that the `unionised_signature_def` must already have been unionised to accept `Node`s.
"""
function preprocess_declaration(signature_def)
    op = signature_def[:name]
    args = signature_def[:args]
    y = gensym(:y)
    ȳ = gensym(:ȳ)

    # preprocess has a similar definition, signature-wise, to what is in signature_def
    preprocess_def = Dict{Symbol, Any}(
        :name => :preprocess,
        :args => [op, :($y::Branch), ȳ, args...],
        :body => quote
            pullback = getfield($y, :pullback)  # avoid issues with getproperty overloading
            @assert(pullback !== nothing, "pullback not set, probably because different code path used for preprocess vs for ∇. Probably need to delete a defination for ∇")
            return pullback($ȳ)
        end,
    )

    where_params = get(signature_def, :whereparams, nothing)
    if where_params !== nothing
        preprocess_def[:whereparams] = where_params
    end
    return ExprTools.combinedef(preprocess_def)
end

"""
    ∇_declaration(unionised_signature_def)

Generates that AST for the overload of the `∇` function which returns the gradient for
specified arguments.
Basically this generates things like:
`∇(::opT, ::Type{Arg{N}}, p, y, ȳ, xs...) where N = p[N+1]  # Skip dself`
where `p` is the pullback computed by [`preprocess`](@ref)

Note that the `unionised_signature_def` must already have been unionised to accept `Node`s.
"""
function ∇_declaration(signature_def)
    # For readability lets name all the parts, NB: this is being a bit too cute.
    op = signature_def[:name]
    args = signature_def[:args]

    N = gensym(:N)
    p = gensym(:p)
    y = :(::Any)
    ȳ = :(::Any)

    ∇_def = Dict{Symbol, Any}(
        :name => :∇,
        :args => [op, :(::Type{Arg{$N}}), p, y, ȳ, args...],
        :whereparams => [N; get(signature_def, :whereparams, [])],
        :body => quote $p[$N+1] end,  # skip dself
        :kwargs => [:(kwargs...)],
    )
    return ExprTools.combinedef(∇_def)
end

# Find a tape, ds might be Nodes or might be something else.
# All nodes should have the same tape, so the first one will do
get_tape(ds) = first(tape(d) for d in ds if d isa Node)


text:
# include("rules/patterns.jl")
export Pattern
export PatLiteral
export PatVar
export PatTerm
export PatAllTerm
export PatTypeAssertion
export PatSplatVar
export PatEquiv
export patvars
export @pat

# include("rules/rule_types.jl")
export Rule
export SymbolicRule
export RewriteRule
export BidirRule
export EqualityRule
export UnequalRule
export DynamicRule
export MultiPatRewriteRule
export canprune

# include("rules/rule_dsl.jl")
export Rule
export Theory
export @rule
export @theory
export @pruningrule


text:
using Lighthouse
using Documenter

makedocs(; modules=[Lighthouse], sitename="Lighthouse",
         authors="Beacon Biosignals and other contributors",
         pages=["API Documentation" => "index.md",
                "Plotting" => "plotting.md"],
         # makes docs fail hard if there is any error building the examples,
         # so we don't just miss a build failure!
         strict = true)

deploydocs(repo="github.com/beacon-biosignals/Lighthouse.jl.git",
           devbranch="main", push_preview=true)


text:
"""
	$(SIGNATURES)

Ask user for a choice from a set of strings.
Returns chosen string or "" if no match.
"""
function ask_for_choice(questionStr :: String, choiceV :: Vector{String} = ["Yes", "No"])
	print(questionStr * "  ")
	print(choiceV);
	print("  ");
	ans1 = readline(stdin);
	idx = findfirst(isequal.(lowercase.(choiceV),  lowercase(ans1)));
	if isnothing(idx)
		outStr = "";
	else
		outStr = choiceV[idx];
	end
	return outStr
end

"""
	$(SIGNATURES)

Ask a yes/no question. Repeats until user makes a valid choice.
Returns `true` for "yes".
"""
function ask_yes_no(questionStr :: String)
	done = false;
	res = false;
	while !done
		outStr = ask_for_choice(questionStr, ["Yes", "No"]);
		if outStr == "Yes"
			done = true;
			res = true;
		elseif outStr == "No"
			done = true;
			res = false;
		end
	end
	return res
end

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