text:
module GridInterpolations

using StaticArrays
using Printf
using LinearAlgebra

export AbstractGrid, RectangleGrid, SimplexGrid, dimensions, length, label, ind2x, ind2x!, interpolate, maskedInterpolate, interpolants, vertices

abstract type AbstractGrid{D} end # D is the dimension

mutable struct RectangleGrid{D} <: AbstractGrid{D}
    cutPoints::Vector{Vector{Float64}}
    cut_counts::Vector{Int}
    cuts::Vector{Float64}
    index::Vector{Int}
    weight::Vector{Float64}
    index2::Vector{Int}
    weight2::Vector{Float64}

    function RectangleGrid{D}(cutPoints...) where D
        cut_counts = Int[length(cutPoints[i]) for i = 1:length(cutPoints)]
        cuts = vcat(cutPoints...)
        myCutPoints = Array{Vector{Float64}}(undef, length(cutPoints))
        numDims = length(cutPoints)
        @assert numDims == D
        for i = 1:numDims
            if length(Set(cutPoints[i])) != length(cutPoints[i])
                error(@sprintf("Duplicates cutpoints are not allowed (duplicates observed in dimension %d)",i))
            end
            if !issorted(cutPoints[i])
                error("Cut points must be sorted")
            end
            myCutPoints[i] = cutPoints[i]
        end
        index = zeros(Int, 2^numDims)
        weight = zeros(Float64, 2^numDims)
        index[1] = 1
        weight[1] = 1.0
        index2 = zeros(Int, 2^numDims)
        weight2 = zeros(Float64, 2^numDims)
        index2[1] = 1
        weight2[1] = 1.0
        return new(myCutPoints, cut_counts, cuts, index, weight, index2, weight2)
    end
end

RectangleGrid(cutPoints...) = RectangleGrid{length(cutPoints)}(cutPoints...)

mutable struct SimplexGrid{D} <: AbstractGrid{D}
    cutPoints::Vector{Vector{Float64}}
    cut_counts::Vector{Int}
    cuts::Vector{Float64}
    index::Vector{Int}
    weight::Vector{Float64}
    x_p::Vector{Float64} # residuals
    ihi::Vector{Int} # indices of cuts above point
    ilo::Vector{Int} # indices of cuts below point
    n_ind::Vector{Int}

    function SimplexGrid{D}(cutPoints...) where D
        cut_counts = Int[length(cutPoints[i]) for i = 1:length(cutPoints)]
        cuts = vcat(cutPoints...)
        myCutPoints = Array{Vector{Float64}}(undef, length(cutPoints))
        numDims = length(cutPoints)
        @assert numDims == D
        for i = 1:numDims
            if length(Set(cutPoints[i])) != length(cutPoints[i])
                error(@sprintf("Duplicates cutpoints are not allowed (duplicates observed in dimension %d)",i))
            end
            if !issorted(cutPoints[i])
                error("Cut points must be sorted")
            end
            myCutPoints[i] = cutPoints[i]
        end
        index = zeros(Int, numDims+1) # d+1 points for simplex
        weight = zeros(Float64, numDims+1)
        x_p = zeros(numDims) # residuals
        ihi = zeros(Int, numDims) # indicies of cuts above point
        ilo = zeros(Int, numDims) # indicies of cuts below point
        n_ind = zeros(Int, numDims)
        return new(myCutPoints, cut_counts, cuts, index, weight, x_p, ihi, ilo, n_ind)
    end
end

SimplexGrid(cutPoints...) = SimplexGrid{length(cutPoints)}(cutPoints...)

Base.length(grid::RectangleGrid) = prod(grid.cut_counts)
Base.length(grid::SimplexGrid) = prod(grid.cut_counts)

dimensions(grid::AbstractGrid{D}) where D = D
Base.ndims(grid::AbstractGrid{D}) where D = D

label(grid::RectangleGrid) = "multilinear interpolation grid"
label(grid::SimplexGrid) = "simplex interpolation grid"

# showall returns a valid constructor incantation - it will be called when repr is called on a grid
function Base.show(io::IO, grid::AbstractGrid)
    if get(io, :compact, false)
        print(io, "$(typeof(grid)) with $(length(grid)) points")
    else
        print(io, "$(typeof(grid))(")
        for v in grid.cutPoints
            show(io, v)
            print(io, ',')
        end
        print(io, ')')
    end
end

function ind2x(grid::AbstractGrid, ind::Int)
    ndims = dimensions(grid)
    x = Array{Float64}(undef, ndims)
    ind2x!(grid, ind, x)
    x::Array{Float64}
end

function ind2x!(grid::AbstractGrid, ind::Int, x::AbstractArray)
    # Populates x with the value at ind.
	# In-place version of ind2x.
	# Example:
	#   rgrid = RectangleGrid([2,5],[20,50])
	#   x = [0,0]
	#   ind2x!(rgrid,4,x)  # x now contains [5,50]
    #   @show x            # displays [5,50]
    ndims = dimensions(grid)
    stride = grid.cut_counts[1]
    for i=2:ndims-1
        stride *= grid.cut_counts[i]
    end

    for i=(ndims-1):-1:1
        rest = rem(ind-1, stride) + 1
        x[i + 1] = grid.cutPoints[i + 1][div(ind - rest, stride) + 1]
        ind = rest
        stride = div(stride, grid.cut_counts[i])
    end
    x[1] = grid.cutPoints[1][ind]
    nothing
end

# masked interpolation ignores points that are masked
function maskedInterpolate(grid::AbstractGrid, data::DenseArray, x::AbstractVector, mask::BitArray{1})
    index, weight = interpolants(grid, x)
    val = 0
    totalWeight = 0
    for i = 1:length(index)
        if mask[index[i]]
            continue
        end
        val += data[index[i]] * weight[i]
        totalWeight += weight[i]
    end
    return val / totalWeight
end

interpolate(grid::AbstractGrid, data::Matrix, x::AbstractVector) = interpolate(grid, map(Float64, data[:]), x)

function interpolate(grid::AbstractGrid, data::DenseArray, x::AbstractVector)
    index, weight = interpolants(grid, x)
    dot(data[index], weight)
end

function interpolants(grid::RectangleGrid, x::AbstractVector)
    if any(isnan, x)
        throw(DomainError("Input contains NaN!"))
    end
    cut_counts = grid.cut_counts
    cuts = grid.cuts

    # Reset the values in index and weight:
    fill!(grid.index,0)
    fill!(grid.index2,0)
    fill!(grid.weight,0)
    fill!(grid.weight2,0)
    grid.index[1] = 1
    grid.index2[1] = 1
    grid.weight[1] = 1.
    grid.weight2[1] = 1.

    l = 1
    subblock_size = 1
    cut_i = 1
    n = 1
    for d = 1:length(x)
        coord = x[d]
        lasti = cut_counts[d]+cut_i-1
        ii = cut_i

        if coord <= cuts[ii]
            i_lo, i_hi = ii, ii
        elseif coord >= cuts[lasti]
            i_lo, i_hi = lasti, lasti
        else
            while cuts[ii] < coord
                ii = ii + 1
            end
            if cuts[ii] == coord
                i_lo, i_hi = ii, ii
            else
                i_lo, i_hi = (ii-1), ii
            end
        end

        if i_lo == i_hi
            for i = 1:l
                grid.index[i] += (i_lo - cut_i)*subblock_size
            end
        else
            low = (1 - (coord - cuts[i_lo])/(cuts[i_hi]-cuts[i_lo]))
            for i = 1:l
                grid.index2[i  ] = grid.index[i] + (i_lo-cut_i)*subblock_size
                grid.index2[i+l] = grid.index[i] + (i_hi-cut_i)*subblock_size
            end
            copyto!(grid.index,grid.index2)
            for i = 1:l
                grid.weight2[i  ] = grid.weight[i]*low
                grid.weight2[i+l] = grid.weight[i]*(1-low)
            end
            copyto!(grid.weight,grid.weight2)
            l = l*2
            n = n*2
        end
        cut_i = cut_i + cut_counts[d]
        subblock_size = subblock_size*(cut_counts[d])
    end

    if l<length(grid.index)
        # This is true if we don't need to interpolate all dimensions because we're on a boundary:
        return grid.index[1:l]::Vector{Int}, grid.weight[1:l]::Vector{Float64}
    end
    grid.index::Vector{Int}, grid.weight::Vector{Float64}
end

function interpolants(grid::SimplexGrid, x::AbstractVector)

    weight = grid.weight
    index  = grid.index

    x_p = grid.x_p # residuals
    ihi = grid.ihi # indicies of cuts above point
    ilo = grid.ilo # indicies of cuts below point
    n_ind = grid.n_ind

    cut_counts = grid.cut_counts
    cuts = grid.cuts

    cut_i = 1

    for i = 1:dimensions(grid)
        # find indicies of coords if match
        coord = x[i]
        lasti = cut_counts[i]+cut_i-1
        ii = cut_i
        # check bounds, snap to closest if out
        if coord <= cuts[ii]
            ihi[i] = ii
            ilo[i] = ii
            x_p[i] = 0.0
        elseif coord >= cuts[lasti]
            ihi[i] = lasti
            ilo[i] = lasti
            x_p[i] = 0.0
        else
            # increment through cut points if in bounds
            while cuts[ii] < coord
                ii += 1
            end
            # if on cut assign cut indecies
            if cuts[ii] == coord
                ilo[i] = ii
                ihi[i] = ii
                x_p[i] = 0.0
            else
                # if between cuts assign lo and high indecies and translate
                ilo[i] = ii-1
                ihi[i] = ii
                lo = cuts[ilo[i]]
                hi = cuts[ihi[i]]
                x_p[i] = (x[i] - lo) / (hi - lo)
            end
        end
        cut_i = cut_i + cut_counts[i]
    end

    # initialize sort indecies
    for i = 1:length(n_ind); n_ind[i] = i; end
    # sort translated and scaled x values
    sortperm!(n_ind, x_p, rev=true) ############################################# killer of speed
    x_p = x_p[n_ind]
    n_ind = n_ind .- 1

    # get weight
    for i = 1:(length(x_p)+1)
        if i == 1
            weight[i] = 1 - x_p[i]
        elseif i == length(x_p)+1
            weight[i] = x_p[i-1]
        else
            weight[i] = x_p[i-1] - x_p[i]
        end
    end

    # get indices
    fill!(index, 0)
    i_index = 0
    for i = 1:(length(x_p)+1)
        siz = 1
        ct = 0
        good_count = 1
        if i > 1
            i_index = i_index + 2^(n_ind[i-1])
        end
        for k = 1:length(x)
            onHi = ((i_index & good_count) > 0)
            good_count <<= 1
            if onHi
                index[i] += (ihi[k] - 1 - ct) * siz
            else
                index[i] += (ilo[k] - 1 - ct) * siz
            end
            siz = siz*cut_counts[k]
            ct += cut_counts[k]
        end
        index[i] += 1
    end

    weight = weight ./ sum(weight)

    return index::Vector{Int}, weight::Vector{Float64}
end

"Return a vector of SVectors where the ith vector represents the vertex corresponding to the ith index of grid data."
function vertices(grid::AbstractGrid)
    n_dims = dimensions(grid)
    mem = Array{Float64,2}(undef, n_dims, length(grid))

    for idx = 1 : length(grid)
        this_idx::Int = idx-1

        # Get the correct index into each dimension
        # and populate vertex index with corresponding cut point
        for j = 1 : n_dims
            cut_idx::Int = this_idx % grid.cut_counts[j]
            this_idx = div(this_idx,grid.cut_counts[j])
            mem[j, idx] = grid.cutPoints[j][cut_idx+1]
        end
    end

    #=
    This relies on the memory layout of Matrix to stay the same, so is a
    possible source of future errors. However, it is documented
    (http://juliaarrays.github.io/StaticArrays.jl/stable/pages/
    api.html#Arrays-of-static-arrays-1), and tests should catch these errors.
    =#
    return reshape(reinterpret(SVector{n_dims, Float64}, mem), (length(grid),))
end

#################### sortperm! is included in Julia v0.4 ###################

using Base.Order # for sortperm!, should be availiable in v 0.4
using Base.Sort # for sortperm!
const DEFAULT_UNSTABLE = QuickSort

function sortperm!(x::Vector{I}, v::AbstractVector; alg::Algorithm=DEFAULT_UNSTABLE,
                   lt::Function=isless, by::Function=identity, rev::Bool=false, order::Ordering=Forward) where {I<:Integer}
    sort!(x, alg, Perm(ord(lt,by,rev,order),v))
end

function Base.iterate(iter::RectangleGrid, state::Int64=1)
    return state<=length(iter) ? (ind2x(iter, state), state+1) : nothing
end

Base.getindex(grid::RectangleGrid, key::CartesianIndex) = ind2x(grid, LinearIndices(Dims((grid.cut_counts...,)))[key])
Base.getindex(grid::RectangleGrid, indices...) = ind2x(grid, LinearIndices(Dims((grid.cut_counts...,)))[indices...])

end # module


text:
## center_of_mass_velocity
using RigidBodyDynamics
function center_of_mass_velocity(state::MechanismState)
    h = momentum(state)
    FreeVector3D(h.frame, linear(h)) / mass(state.mechanism)
end

## hrep conversion
function RigidBodyDynamics.Contact.HRep(hr::Polyhedra.HRepresentation)
    M = Polyhedra.nhalfspaces(hr)
    N = Polyhedra.fulldim(hr)
    return map(Polyhedra.halfspaces(hr)) do halfspace
        RigidBodyDynamics.Contact.HalfSpace(SVector{N}(halfspace.a), halfspace.β)
    end |> SVector{M} |> RigidBodyDynamics.Contact.HRep
end

function Polyhedra.hrep(hr::RigidBodyDynamics.Contact.HRep)
    # hrep(Vector([Polyhedra.HalfSpace(halfspace.outward_normal, halfspace.offset) for halfspace in hr.halfspaces]))
    hrep([Polyhedra.HalfSpace(halfspace.outward_normal, halfspace.offset) for halfspace in hr.halfspaces])
end

## setelement! for ContactModel
using MeshCatMechanisms, RigidBodyDynamics
function MeshCatMechanisms.setelement!(mvis::MechanismVisualizer, contact_model::ContactModel, args...; base_name="collision element ")
    mechanism = mvis.state.mechanism
    for group in contact_model.collision_groups
        for element in group
            geometry = element.geometry
            frame = element.transform.from
            if geometry isa MeshCatMechanisms.GeometryLike
                setelement!(mvis, frame, element.geometry, args...)
            elseif geometry isa Contact.HRep
                setelement!(mvis, frame, GLNormalMesh(polyhedron(hrep(geometry))))
                setobject!(mvis[frame][:triad], Triad(0.1))
            end
        end
    end
end

## PosePlan
using QPControl.Trajectories: Piecewise, Constant
using CentroidalTrajOpt: closest_pose

function QPWalkingControl.PosePlan(pose0::Transform3D,
        point_trajectory::Piecewise{<:Constant{<:SVector{3}}},
        contact_indicators::Piecewise{<:Constant{<:AbstractVector{Bool}}},
        region_data::AbstractVector{<:ContactRegion})
    @assert point_trajectory.breaks == contact_indicators.breaks
    T = eltype(pose0)
    plan = PosePlan{T}()
    previous_region_index = findfirst(contact_indicators.subfunctions[1].value)
    previous_pose = pose0
    start_time = nothing
    n = length(point_trajectory.subfunctions)
    for i = 2 : n
        region_index = findfirst(contact_indicators.subfunctions[i].value)
        if region_index == previous_region_index
            continue
        elseif region_index === nothing
            start_time = point_trajectory.breaks[i]
            previous_region_index = region_index
            continue
        else
            point = point_trajectory.subfunctions[i].value
            region = region_data[region_index]
            normal = region.transform.linear * SVector(0, 0, 1)
            final_pose = Transform3D(previous_pose.from, previous_pose.to, closest_pose(AffineMap(previous_pose), point, normal))
            duration = point_trajectory.breaks[i] - start_time
            push!(plan, start_time, duration, final_pose)
            previous_pose = final_pose
            previous_region_index = region_index
        end
    end
    return plan
end

function check_kinematic_constraints_satisfied(contacts, c0, max_com_to_contact_distance, min_inter_contact_distance)
    ## Basic initial state feasibility checks
    for i in eachindex(contacts)
        p_i_0 = last(contacts[i])
        norm(c0 - p_i_0) <= max_com_to_contact_distance || error()
        for j in eachindex(contacts)
            p_j_0 = last(contacts[j])
            if i != j
                norm(p_i_0 - p_j_0) >= min_inter_contact_distance || error()
            end
        end
    end
end

ext: jl


text:
"""`semimajoraxis(x, y, z, u, v, w, m1, m2)`

Calculates semimajor axis for the two
body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies
with masses `m1` and `m2`."""
function semimajoraxis(x, y, z, u, v, w, m1, m2)
    r = sqrt(x^2+y^2+z^2)
    vsq = u^2+v^2+w^2
    M = m1+m2
    return 1/( (2/r)-(vsq/M) )
end

"""`eccentricity(x, y, z, u, v, w, m1, m2)`

Calculates eccentricity for the two
body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies
with masses `m1` and `m2`."""
function eccentricity(x, y, z, u, v, w, m1, m2)
    hsq = (y*w-z*v)^2 + (z*u-x*w)^2 + (x*v-y*u)^2
    a = semimajoraxis(x, y, z, u, v, w, m1, m2)
    M = m1+m2
    quotient = hsq/( M*a )
    return sqrt(1 - quotient)
end

"""`inclination(x, y, z, u, v, w)`

Calculates inclination for the two
body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies."""
function inclination(x, y, z, u, v, w)
    hz = x*v-y*u
    h = sqrt( (y*w-z*v)^2 + (z*u-x*w)^2 + (x*v-y*u)^2 )
    return acos(hz/h)
end

"""`aei(x, y, z, u, v, w, m1, m2)`

Returns semimajor axis `a`, eccentricity `e` and inclination `I` for the two-body
problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)`vectors
between two bodies with masses `m1` and `m2`."""
function aei(x, y, z, u, v, w, m1, m2)
    r = sqrt(x^2+y^2+z^2)
    vsq = u^2+v^2+w^2
    M = m1+m2
    a = 1/( (2/r)-(vsq/M) )
    hsq = (y*w-z*v)^2 + (z*u-x*w)^2 + (x*v-y*u)^2
    h = sqrt(hsq)
    quotient = hsq/( M*a )
    e = sqrt(1 - quotient)
    hz = x*v-y*u
    inc = acos(hz/h)
    return a, e, inc
end

"""`ae(x, y, z, u, v, w)`

Returns semimajor axis `a` and eccentricity `e` for the two-body problem defined
by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two
bodies with masses `m1` and `m2`."""
function ae(x, y, z, u, v, w, m1, m2)
    r = sqrt(x^2+y^2+z^2)
    vsq = u^2+v^2+w^2
    M = m1+m2
    a = 1/( (2/r)-(vsq/M) )
    hsq = (y*w-z*v)^2 + (z*u-x*w)^2 + (x*v-y*u)^2
    quotient = hsq/( M*a )
    e = sqrt(1 - quotient)
    return a, e
end

"""`longascnode(x, y, z, u, v, w)`

Calculates longitude of ascending node for the two
body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies."""
function longascnode(
    x, y, z,
    u, v, w
    )

    #the longitude of ascending node is computed as
    #the angle between x-axis and the vector n = (-hy,hx,0)
    #where hx, hy, are resp., the x and y comps. of ang mom vector per unit mass, h

    res = atan( y*w-z*v, x*w-z*u)

    if constant_term(res) ≥ zero(constant_term(res))
        return res
    else
        return res+2pi
    end

end

"""`argperi(x, y, z, u, v, w, m1, m2)`

Calculates the instantaneous (osculating) argument of pericentre for the two
body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies
with masses `m1` and `m2`."""
function argperi(
    x, y, z,
    u, v, w,
    m1, m2
    )

    #n = (z-axis unit vector)×h = (-hy, hx, 0)
    n = [x*w-z*u, y*w-z*v, zero(x)]
    e = lrlvec(x,y,z,u,v,w,m1,m2) #cartesian comps. of Laplace-Runge-Lenz vector
    n = n/sqrt(n[1]^2+n[2]^2+n[3]^2)
    e = e/sqrt(e[1]^2+e[2]^2+e[3]^2)
    cosω = dot(n, e)

    if constant_term(e[3]) >= zero(constant_term(x))
        return acos(cosω)
    else
        return 2pi-acos(cosω)
    end

end

"""`longperi(x, y, z, u, v, w, m1, m2)`

Calculates the instantaneous (osculating) longitude of pericentre for the two
body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies
with masses `m1` and `m2`."""
function longperi(
    x, y, z,
    u, v, w,
    m1, m2
    )

    return argperi(x, y, z, u, v, w, m1, m2)+longascnode(x, y, z, u, v, w)
end

"""`trueanomaly(x, y, z, u, v, w, m1, m2)`

Calculates the instantaneous (osculating) true anomaly for the two
body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies
with masses `m1` and `m2`."""
function trueanomaly(
    x, y, z,
    u, v, w,
    m1, m2
    )

    R = sqrt(x^2+y^2+z^2)
    V2 = u^2+v^2+w^2
    V = sqrt(V2)

    a = 1.0/( (2.0/R)-(V2/(m1+m2)) )

    h2 = (y*w-z*v)^2 + (z*u-x*w)^2 + (x*v-y*u)^2
    h = sqrt(h2)

    quotient=h2/( (m1+m2)*a )

    e = sqrt(1.0 - quotient)

    Rdot = (x*u+y*v+z*w)/R

    sin_f = Rdot*a*(1.0-e^2)/( e*h )
    cos_f = (1.0/e)*( -1.0+a*(1.0-e^2)/R )

    res = atan(sin_f, cos_f)

    if res ≥ zero(res)
        return res
    else
        return res+2pi
    end
end

"""`ecanomaly(x, y, z, u, v, w, m1, m2)`

Calculates the instantaneous (osculating) eccentric anomaly for the two
body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies
with masses `m1` and `m2`."""
function ecanomaly(
    x, y, z,
    u, v, w,
    m1, m2
    )

    a = semimajoraxis(x, y, z, u, v, w, m1, m2)
    e = eccentricity(x, y, z, u, v, w, m1, m2)
    R = sqrt(x^2+y^2+z^2)
    Rdot = (x*u+y*v+z*w)/R
    mu = m1+m2

    cosE = (a-R)/(a*e)
    sinE = Rdot*R/( e*sqrt(mu*a) )

    res = atan(sinE, cosE)

    return res

end

"""`meananomaly(x, y, z, u, v, w, m1, m2)`

Calculates the instantaneous (osculating) mean anomaly for the two
body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies
with masses `m1` and `m2`."""
function meananomaly(
    x, y, z,
    u, v, w,
    m1, m2
    )

    E = ecanomaly(x, y, z, u, v, w, m1, m2)
    e = eccentricity(x, y, z, u, v, w, m1, m2)

    return E-e*sin(E)
end

"""`timeperipass(t, x, y, z, u, v, w, m1, m2)`

Calculates the instantaneous (osculating) time of pericentre passage, at time `t`, for the two
body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies
with masses `m1` and `m2`."""
function timeperipass(
    t, x, y, z,
    u, v, w,
    m1, m2
    )

    me_an = meananomaly(x, y, z, u, v, w, m1, m2) #mean anomaly
    a = semimajoraxis(x, y, z, u, v, w, m1, m2)
    me_mo = meanmotion(m1+m2, a) #mean motion

    return t-me_an/me_mo
end

"""`rungelenzx(x, y, z, u, v, w, m1, m2)`

Calculates the instantaneous (osculating) x-component of the Laplace-Runge-Lenz vector for the two
body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies
with masses `m1` and `m2`."""
function rungelenzx(
    x, y, z,
    u, v, w,
    m1, m2
    )

    # v = [u, v, w]
    # h = [y.*w-z.*v, z.*u-x.*w, x.*v-y.*u]
    mu = m1+m2
    # r = [x, y, z]
    rmag = sqrt(x^2+y^2+z^2)

    lrl_x = ( -(z*u-x*w)*w+(x*v-y*u)*v )/(mu)-x/rmag

    return lrl_x
end

"""`rungelenzy(x, y, z, u, v, w, m1, m2)`

Calculates the instantaneous (osculating) y-component of the Laplace-Runge-Lenz vector for the two
body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies
with masses `m1` and `m2`."""
function rungelenzy(
    x, y, z,
    u, v, w,
    m1, m2
    )

    # v = [u, v, w]
    # h = [y.*w-z.*v, z.*u-x.*w, x.*v-y.*u]
    mu = m1+m2
    # r = [x, y, z]
    rmag = sqrt(x^2+y^2+z^2)

    lrl_y = ( -(x*v-y*u)*u+(y*w-z*v)*w )/(mu)-y/rmag

    return lrl_y
end

"""`rungelenzz(x, y, z, u, v, w, m1, m2)`

Calculates the instantaneous (osculating) z-component of the Laplace-Runge-Lenz vector for the two
body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies
with masses `m1` and `m2`."""
function rungelenzz(
    x, y, z,
    u, v, w,
    m1, m2
    )

    # v = [u, v, w]
    # h = [y.*w-z.*v, z.*u-x.*w, x.*v-y.*u]
    mu = m1+m2
    # r = [x, y, z]
    rmag = sqrt(x^2+y^2+z^2)

    lrl_z = ( -(y*w-z*v)*v+(z*u-x*w)*u )/(mu)-z/rmag

    return lrl_z
end

"""`rungelenzmag(x, y, z, u, v, w, m1, m2)`

Calculates the instantaneous (osculating) magnitude of the Laplace-Runge-Lenz vector for the two
body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies
with masses `m1` and `m2`."""
function rungelenzmag(
    x, y, z,
    u, v, w,
    m1, m2
    )

    # v = [u, v, w]
    # h = [y.*w-z.*v, z.*u-x.*w, x.*v-y.*u]
    mu = m1+m2
    # r = [x, y, z]
    rmag = sqrt(x^2+y^2+z^2)

    lrl_x = ( -(z*u-x*w)*w+(x*v-y*u)*v )/(mu)-x/rmag
    lrl_y = ( -(x*v-y*u)*u+(y*w-z*v)*w )/(mu)-y/rmag
    lrl_z = ( -(y*w-z*v)*v+(z*u-x*w)*u )/(mu)-z/rmag

    return sqrt(lrl_x^2+lrl_y^2+lrl_z^2)
end

"""`lrlvec(x, y, z, u, v, w, m1, m2)`

Calculates the instantaneous (osculating) cartesian components of the Laplace-Runge-Lenz vector for the two
body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies
with masses `m1` and `m2`."""
function lrlvec(
    x, y, z,
    u, v, w,
    m1, m2
    )

    mu = m1+m2
    rmag = sqrt(x^2+y^2+z^2)

    lrl_x = ( -(z*u-x*w)*w+(x*v-y*u)*v )/(mu)-x/rmag
    lrl_y = ( -(x*v-y*u)*u+(y*w-z*v)*w )/(mu)-y/rmag
    lrl_z = ( -(y*w-z*v)*v+(z*u-x*w)*u )/(mu)-z/rmag

    return [lrl_x, lrl_y, lrl_z]
end

# compute eccentric anomaly (E) from eccentricity (e) and mean anomaly (M)
function eccentricanomaly(e::T, M::T) where {T <: Number}
    M0 = mod2pi(M)
    E0 = M0 + sign(sin(M0))*0.85*e #Murray-Dermotts' initial estimate
    # successive approximations via Newtons' method
    for i in 0:4
        #TODO: implement modified Newton's method for Kepler's equation (Murray-Dermott)
        Eans = E0 - (E0-e*sin(E0)-M0)/(1.0-e*cos(E0))
        E0 = Eans
    end
    return E0
end

# compute true anomaly (f) from eccentricity (e) and eccentric anomaly (E)
function trueanomaly(e,E)
    Enew = mod2pi(E)
    return 2.0*atan(  sqrt((1.0+e)/(1.0-e))*tan(Enew/2)  )
end

# compute true anomaly from eccentricity and mean anomaly
function meanan2truean(e,M)
    return trueanomaly(e, eccentricanomaly(e, M))
end

# get mean motion from mass parameter (mu) and semimajor axis (a)
function meanmotion(mu,a)
    return sqrt(mu/(a^3))
end

# get mean anomaly from mean motion (n), time (t) and time of pericenter passage (taup)
function meananomaly(n, t, taup)
    return n*(t-taup)
end

# compute true anomaly from time, a, e, mu and taup
function time2truean(a, e, mu, t, taup)
    return meanan2truean(e, meananomaly(meanmotion(mu, a), t, taup))
end


text:
module RTLSDR

import Base.open, Base.close

export
    RtlSdr,
    open,
    close,
    read_bytes,
    read_samples,
    packed_bytes_to_iq

export
    get_rate,
    set_rate,
    get_freq,
    set_freq,
    set_agc_mode,
    set_tuner_gain_mode

include("c_interface.jl")

mutable struct RtlSdr
    valid_ptr::Bool
    dongle_ptr::Ptr{rtlsdr_dev}

    function RtlSdr()
        dp = rtlsdr_open()

        r = new(true, dp)

        # default sample rate and center freq
        set_rate(r, 2.0e6)
        set_freq(r, 88.5e6)

        return r
    end
end

function open(r::RtlSdr, device_index::Int)
end

function close(r::RtlSdr)
    rtlsdr_close(r.dongle_ptr)
    r.valid_ptr = false
end

"""
`set_rate(r::RtlSdr, sample_rate_Hz)
"""
function set_rate(r::RtlSdr, sample_rate_Hz)
    @assert r.valid_ptr
    rtlsdr_set_sample_rate(r.dongle_ptr, sample_rate_Hz)
end
function get_rate(r::RtlSdr)
    @assert r.valid_ptr
    rate = rtlsdr_get_sample_rate(r.dongle_ptr)
    return Int(rate)
end

"""
`set_freq(r::RtlSdr, freq_Hz)`

Interface for `rtlsdr_set_center_freq`.
"""
function set_freq(r::RtlSdr, freq_Hz)
    @assert r.valid_ptr
    rtlsdr_set_center_freq(r.dongle_ptr, freq_Hz)
end
function get_freq(r::RtlSdr)
    @assert r.valid_ptr
    freq = rtlsdr_get_center_freq(r.dongle_ptr)
    return Int(freq)
end

"""
Interface for `rtlsdr_set_agc_mode`
"""
function set_agc_mode(r::RtlSdr, mode)
    @assert r.valid_ptr
    rtlsdr_set_agc_mode(r.dongle_ptr, mode)
end

"""
Interface for `rtlsdr_set_tuner_gain_mode`
"""
function set_tuner_gain_mode(r::RtlSdr, manual)
    @assert r.valid_ptr
    rtlsdr_set_tuner_gain_mode(r.dongle_ptr, manual)
end

"""
`read_bytes(r::RtlSdr, num_bytes)`

High-level interface for `rtlsdr_read_sync`.

`num_bytes` must be a multiple of 512.

Returns a vector of length `num_bytes` of Uint8 (bytes).
"""
function read_bytes(r::RtlSdr, num_bytes)
    @assert r.valid_ptr
    return read_bytes(r.dongle_ptr, num_bytes)
end

"""
`read_samples(r::RtlSdr, num_samples)`

Returns a vector of length `num_samples` with complex numbers.
"""
function read_samples(r::RtlSdr, num_samples)
    @assert r.valid_ptr
    num_bytes = 2num_samples
    raw_data = read_bytes(r.dongle_ptr, num_bytes)
    return packed_bytes_to_iq(raw_data)
end

"""
`packed_bytes_to_iq(bytes)`
"""
function packed_bytes_to_iq(bytes)
    num_bytes = length(bytes)
    num_iq = round(Int, num_bytes/2.0, RoundDown)
    iq_vals = zeros(Complex{Float64}, num_iq)
    den = 255.0/2.0

    for i = 1:num_iq
        iq_vals[i] = bytes[2i-1]/den - 1.0 + im*(bytes[2i]/den - 1.0)
    end

    return iq_vals
end

end # module


text:
using BandedMatrices, LinearAlgebra
export ProductInversion

"""
    ProductInversion(r, μ; verbose=false) <: Decomposition
Decompose a timeseries `s` into a **product** `x * r`, given that you have a good
estimate of the second factor `r` (the "input") and you need `x` (the "multiplier")
but you can't do simply `x =  s ./ r` because `r` contains zeros.

This method works well when the characteristic timescales of `x` are comparable,
or larger than those of `r` but not much smaller.

The second argument `μ` is a regularization parameter. In short, we estimate `r`
by minimizing a cost with two components: that `x` is close to `s/r` and that `x`
is smooth. `μ` is the multiplier of the smoothness cost.

You can give a vector as `μ`. The process will be repeated for all `μ` and the
[`rmse`](@ref) between `s` and the estimated `x * r` will be computed. The
`x` that gives the least error will be returned finally.
If `verbose = true`, the method also prints the pairs `(μ, err)` for each `μ`.

Use the low level `SignalDecomposition.matrix_invert(s, r, μ::Real) → x, err`
to get the error values.
"""
struct ProductInversion{R, M} <: Decomposition
    r::R
    μ::M
    x::R
    xdummy::R
    verbose::Bool
end
ProductInversion(r, μ; verbose=false) = ProductInversion(r, μ, copy(r), copy(r), verbose)

# TODO: performant method `update!(method, s, μ)`

function decompose(t, s, method::ProductInversion{R, <: Real}) where R
    x, err = matrix_invert(s, method.r, method.μ)
    return x, method.r
end

function decompose(t, s, method::ProductInversion{R, <: AbstractArray}) where R
    isequispaced(t) || error("Input time axis must be equispaced for method ProductInversion.")
    leasterr = Inf
    bestx = copy(s)
    bestμ = method.μ[1]
    for μ in method.μ
        x, err = matrix_invert(s, method.r, μ, method.x, method.xdummy)
        method.verbose && @show (μ, err)
        if err < leasterr
            bestx .= x
            leasterr = err
            bestμ = μ
        end
    end
    method.verbose && @show (bestμ, leasterr)
    return bestx, method.r
end

function matrix_invert(s::AbstractVector{T}, r, μ::Real, x = copy(r), xdummy = copy(r)) where {T}
    # TODO: This can be optimized by pre-allocating the diags dictinary
    # and simply updating the values for each μ
    # TODO : add @inbounds
    L = length(s)
    diags = Dict(
        0 => T[6μ + r[i]^2 for i in 1:L],
        1 => fill(T(-4μ), L-1),
        -1 => fill(T(-4μ), L-1),
        2 => fill(T(μ), L-2),
        -2 => fill(T(μ), L-2),
    )
    # take care of border elements
    diags[1][1] = -8μ
    diags[2][1:2] .= 2μ
    diags[-1][end] = -8μ
    diags[-2][end-1:end] .= 2μ
    B = BandedMatrix(diags...)
    @. xdummy = s * r
    ldiv!(x, B, xdummy)
    @. xdummy = x * r
    err = rmse(s, xdummy)
    return x, err
end


text:
#----------------------------------------------------------------------------
# Destructuring.Function;
#----------------------------------------------------------------------------
import ..Destructuring.Function: destructure, code

"""
`destructure(pattern, value, body)`

Creates a let expression where every binding name in `pattern` is attributed
a correspond value found in `value`, with `body` begin the body of the let expression.

examples:
```julia

eval(destructure(:(x+y), :(:(1+2)), quote x,y end))  #> (1,2)

eval(destructure(:(*{x}), :(:(1,2,3)), quote x end)) #> [1,2,3]

pattern = :(for x in coll; b end)
value   = :(for i in 1:10; print(i) end)
eval(destructure(pattern, :value, quote x,coll,b end)) #> (:x, :(1:10), :(print(x))
```

"""
destructure;

"""
`code(vars, ptree, body) -> Expr`

"""
code;

#----------------------------------------------------------------------------
# Destructuring.Applications;
#----------------------------------------------------------------------------
import ..Destructuring.Applications: @letds, @macrods, @anonds, @funds

"""
```julia
@letds p1=v1 p2=v2 ... pn=vn begin
  ...
end
```

Let expression with destructuring properties.

examples:
```julia
@letds f(x::T1,y::T2) = :(add(a::Int,b::Int)) begin
    f == :add && x == :a && y == :b && T1 == T2 == :Int
end

@letds [*{[*{x,y}]}] = :[[1,2,3,4],[a,b,c,d]] begin
  x == Any[Any[1,3], Any[:a, :c]] &&
  y == Any[Any[2,4], Any[:b, :d]]
end

@letds [*{x}, 3, *{y}] = :[1,2,3,4,5] begin
  x == Any[1,2] && y == Any[3,4]
end
```

"""
:@letds;

"""
```julia
@macrods name(p1,p2,...,pn) begin
  ...
end
```

Defines a macro with destructuring properties.

examples:
```julia
@macrods swap((x,y)) :(\$y,\$x)
@macrods split(:C{at} => [*{a}, at, *{b}]) (a,b)

@swap((1,2)) == (2,1)
@split(3 => [1,2,3,4,5]) == (Any[1,2], Any[4,5])

```

"""
:@macrods;

"""
```julia
@anonds (p1,p2,...,pn) -> body
```

Defines an anonymous function with destructuring properties.

examples:
```julia
extremes = @anonds [first, *{middle}, last] -> (first, last)
extremes(:[1,2,3,4]) == (1,4)

first_index = @anonds (:C{i} => [:?{x}, i, *{a}]) -> length(x)+1
last_index  = @anonds (:C{i} => [*{x},  i, *{a}]) -> length(x)+1

ex = :(0 => [1,1,1,0,1,1,0])
first_index(ex) == 4
last_index(ex)  == 7

```

"""
:@anonds;

"""
```julia
@funds name(p1,p2,...,pn) = body
```

Defines an function with destructuring properties.

examples:
```julia
@funds odd_even([*{odd,even}]) = (odd, even)
odd_even(:[1,2,3,4,5,6]) == (Any[1,3,5],Any[2,4,6])

@funds extract_int([*{a}, :T{x, Integer}, *{b})) = x
extract_int(:[a,b,4,d,e,f,g]) == 4

```

"""
:@funds;


text:
using LinSepVerify
using Test

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


text:
using Test

@testset "Directed" begin
  include("Directed.jl")
end

@testset "Undirected" begin
  include("Undirected.jl")
end

@testset "Algebraic" begin
  include("Algebraic.jl")
end

@testset "Algorithms" begin
  include("Algorithms.jl")
end

@testset "Expressions" begin
  include("Expressions.jl")
end

@testset "GraphML" begin
  include("GraphML.jl")
end

@testset "JSON" begin
  include("JSON.jl")
end


text:
# Autogenerated wrapper script for GDAL_jll for i686-w64-mingw32-cxx11
export gdal_contour_path, gdal_grid_path, gdal_rasterize_path, gdal_translate_path, gdaladdo_path, gdalbuildvrt_path, gdaldem_path, gdalinfo_path, gdallocationinfo_path, gdalmanage_path, gdalsrsinfo_path, gdaltindex_path, gdaltransform_path, gdalwarp_path, libgdal, nearblack_path, ogr2ogr_path, ogrinfo_path, ogrlineref_path, ogrtindex_path

using GEOS_jll
using PROJ_jll
using Zlib_jll
using SQLite_jll
using OpenJpeg_jll
using Expat_jll
using Zstd_jll
using Libtiff_jll
using libgeotiff_jll
using LibCURL_jll
using LibSSH2_jll
using MbedTLS_jll
using nghttp2_jll
JLLWrappers.@generate_wrapper_header("GDAL")
JLLWrappers.@declare_executable_product(gdal_contour_path)
JLLWrappers.@declare_executable_product(gdal_grid_path)
JLLWrappers.@declare_executable_product(gdal_rasterize_path)
JLLWrappers.@declare_executable_product(gdal_translate_path)
JLLWrappers.@declare_executable_product(gdaladdo_path)
JLLWrappers.@declare_executable_product(gdalbuildvrt_path)
JLLWrappers.@declare_executable_product(gdaldem_path)
JLLWrappers.@declare_executable_product(gdalinfo_path)
JLLWrappers.@declare_executable_product(gdallocationinfo_path)
JLLWrappers.@declare_executable_product(gdalmanage_path)
JLLWrappers.@declare_executable_product(gdalsrsinfo_path)
JLLWrappers.@declare_executable_product(gdaltindex_path)
JLLWrappers.@declare_executable_product(gdaltransform_path)
JLLWrappers.@declare_executable_product(gdalwarp_path)
JLLWrappers.@declare_library_product(libgdal, "libgdal-28.dll")
JLLWrappers.@declare_executable_product(nearblack_path)
JLLWrappers.@declare_executable_product(ogr2ogr_path)
JLLWrappers.@declare_executable_product(ogrinfo_path)
JLLWrappers.@declare_executable_product(ogrlineref_path)
JLLWrappers.@declare_executable_product(ogrtindex_path)
function __init__()
    JLLWrappers.@generate_init_header(GEOS_jll, PROJ_jll, Zlib_jll, SQLite_jll, OpenJpeg_jll, Expat_jll, Zstd_jll, Libtiff_jll, libgeotiff_jll, LibCURL_jll, LibSSH2_jll, MbedTLS_jll, nghttp2_jll)
    JLLWrappers.@init_executable_product(
        gdal_contour_path,
        "bin\\gdal_contour.exe",
    )

    JLLWrappers.@init_executable_product(
        gdal_grid_path,
        "bin\\gdal_grid.exe",
    )

    JLLWrappers.@init_executable_product(
        gdal_rasterize_path,
        "bin\\gdal_rasterize.exe",
    )

    JLLWrappers.@init_executable_product(
        gdal_translate_path,
        "bin\\gdal_translate.exe",
    )

    JLLWrappers.@init_executable_product(
        gdaladdo_path,
        "bin\\gdaladdo.exe",
    )

    JLLWrappers.@init_executable_product(
        gdalbuildvrt_path,
        "bin\\gdalbuildvrt.exe",
    )

    JLLWrappers.@init_executable_product(
        gdaldem_path,
        "bin\\gdaldem.exe",
    )

    JLLWrappers.@init_executable_product(
        gdalinfo_path,
        "bin\\gdalinfo.exe",
    )

    JLLWrappers.@init_executable_product(
        gdallocationinfo_path,
        "bin\\gdallocationinfo.exe",
    )

    JLLWrappers.@init_executable_product(
        gdalmanage_path,
        "bin\\gdalmanage.exe",
    )

    JLLWrappers.@init_executable_product(
        gdalsrsinfo_path,
        "bin\\gdalsrsinfo.exe",
    )

    JLLWrappers.@init_executable_product(
        gdaltindex_path,
        "bin\\gdaltindex.exe",
    )

    JLLWrappers.@init_executable_product(
        gdaltransform_path,
        "bin\\gdaltransform.exe",
    )

    JLLWrappers.@init_executable_product(
        gdalwarp_path,
        "bin\\gdalwarp.exe",
    )

    JLLWrappers.@init_library_product(
        libgdal,
        "bin\\libgdal-28.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        nearblack_path,
        "bin\\nearblack.exe",
    )

    JLLWrappers.@init_executable_product(
        ogr2ogr_path,
        "bin\\ogr2ogr.exe",
    )

    JLLWrappers.@init_executable_product(
        ogrinfo_path,
        "bin\\ogrinfo.exe",
    )

    JLLWrappers.@init_executable_product(
        ogrlineref_path,
        "bin\\ogrlineref.exe",
    )

    JLLWrappers.@init_executable_product(
        ogrtindex_path,
        "bin\\ogrtindex.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    points, lines, triangles, tetrahedra, domains=read_nastran(filename)

Read simplicial nastran mesh.

# Note
There is currently no sanity check for non-simplicial elements. These elements are just skipped while reading.
"""
function read_nastran(filename)
  number_of_points=0
  number_of_lines=0
  number_of_triangles=0
  number_of_tetrahedra=0
  name_tags=Dict()

  #first read to get number of points
  open(filename) do fid
    while !eof(fid)
      line=readline(fid)
      if line[1]=='$' #line is a comment
        #comments are freqeuntly used to leave more information
        #especially the domain tags
        #this is not standardized so it depends on the software used to write the file

        #ANSA-style tag comment
        #example:
        #$ANSA_NAME_COMMENT;2;PSHELL;TA_IN_Inlet;;NO;NO;NO;NO;
        if length(line)>=18 && line[2:18]=="ANSA_NAME_COMMENT"
          data=split(line[2:end],";")
          if data[3]=="PSOLID" || data[3]=="PSHELL"
            name_tags[data[2]]=data[4]
          end
        end

        #Centaur-style tag comment
        #example:
        #$HMNAME COMP                   1"Walls"
        if length(line)>=12 && line[2:12]=="HMNAME COMP"
          data=split(strip(line[13:end]),'"')
          name_tags[data[1]]=data[2]
        end

        continue
      end
      if length(line)<8 #line is too short probably it is the ENDDATA tag.
        continue
      end
      head=line[1:8]
      if head =="GRID    " || head == "GRID*   " || head[1:5] == "GRID," #TODO allow for whitespace between commas in free format
        number_of_points+=1
      elseif (head[1:6] == "CTRIA3") || (head[1:6] == "CTRIA6")
        number_of_triangles+=1
      elseif head[1:6] == "CTETRA"
        number_of_tetrahedra+=1
      end
    end
  end
  tri_idx=0
  tet_idx=0
  #TODO: sanity check for non-consecutive indexing
  points=Array{Float64}(undef,3,number_of_points)
  lines=Array{Array{UInt32,1}}(undef,number_of_lines)
  triangles=Array{Array{UInt32,1}}(undef, number_of_triangles)
  tetrahedra=Array{Array{UInt32,1}}(undef, number_of_tetrahedra)
  domains=Dict()
  surface_tags=Dict()
  volume_tags=Dict()
  #second read to get data
  open(filename) do fid
    while !eof(fid)
      line=readline(fid)
      if line[1]=='$' #line is a comment
        continue
      end
      # data = parse_nas_line(line)
      # if length(data)==0 #line is to short probaboly it is the ENDDATA tag.
      #   continue
      # end

      if length(line)<8 #line is too short probably it is the ENDDATA tag.
        continue
      end

      free=false #sentinel value for checking for free format
      if occursin(",",line)
        free=true
      end

      head=line[1:8]

      if head =="GRID    " #grid point short format
        data = parse_nas_line(line)
        idx=parse(UInt,data[2])
        points[1,idx]=parse_nas_number(data[4])
        points[2,idx]=parse_nas_number(data[5])
        points[3,idx]=parse_nas_number(data[6])

      elseif head =="GRID*   " #grid point long format
        data = parse_nas_line(line,format=:long)
        idx=parse(UInt,data[2])
        points[1,idx]=parse_nas_number(data[4])
        points[2,idx]=parse_nas_number(data[5])
        line=readline(fid) #long format has two lines
        data = parse_nas_line(line,format=:long)
        points[3,idx]=parse_nas_number(data[2])

      elseif head[1:5] =="GRID,"#grid point free format
        data = parse_nas_line(line,format=:free)
        idx=parse(UInt,data[2])
        points[1,idx]=parse_nas_number(data[4])
        points[2,idx]=parse_nas_number(data[5])
        points[3,idx]=parse_nas_number(data[6])

      elseif (head[1:6] == "CTRIA3") || (head[1:6] == "CTRIA6") #Triangle
        tri_idx+=1
        if free
          data = parse_nas_line(line,format=:free)
        else
          data = parse_nas_line(line)
        end
          #TODO: find nastran file specification and see whether there is
          #something like long format for stuff other than GRID
        idx=tri_idx
        dom=strip(data[3])
        if dom in keys(name_tags)
          dom=name_tags[dom]
        else
          dom="surf"*lpad(dom,4,"0")
        end
        if dom in keys(domains)
          append!(domains[dom]["simplices"],idx)
        else
          domains[dom]=Dict("dimension"=>2,"simplices"=>UInt32[idx])
        end
        tri=[parse(UInt32,data[4]),parse(UInt32,data[5]),parse(UInt32,data[6])]
        triangles[idx]=tri

      elseif head[1:6] == "CTETRA" #Tetrahedron
        tet_idx+=1
        if free
          data = parse_nas_line(line,format=:free)
        else
          data = parse_nas_line(line)
        end
        idx=tet_idx#parse(UInt,data[2])-number_of_triangles
        dom=strip(data[3])
        if dom in keys(name_tags)
          dom=name_tags[dom]
        else
          dom="vol"*lpad(dom,4,"0")
        end
        if dom in keys(domains)
          append!(domains[dom]["simplices"],idx)
        else
          domains[dom]=Dict("dimension"=>3,"simplices"=>UInt32[idx])
        end
        tet=[parse(UInt32,data[4]),parse(UInt32,data[5]),parse(UInt32,data[6]),parse(UInt32,data[7])]
        tetrahedra[idx]=tet
      end
    end
  end

  #sanity check when using higher order simplices (CTRIA6 or CTETRA with ten nodes)
  #Comsol uses this stuff -.-
  used_idx=sort(unique(vcat(tetrahedra...)))
  if used_idx[end]==length(used_idx)
    if length(triangles)!=0
      tri_max_idx=maximum(vcat(triangles...))
    else
      tri_max_idx=0
    end
    if tri_max_idx > used_idx[end]
      println("Warning: Nastran mesh needs manual reordering of point indeces. Not all points are actually used to define simplices.")
    else
      points=points[:,1:used_idx[end]]
    end
  else
    println("Warning: Nastran mesh needs manual reordering of point indeces. Not all points are actually used to define simplices.")
  end
  return points, lines, triangles, tetrahedra, domains
end

function parse_nas_line(txt;format=:short)
  #split line
  if format==:short
    number_of_entries=length(txt)÷8
    data=Array{String}(undef,number_of_entries)
    for idx=1:number_of_entries
      data[idx]=txt[(1+8*(idx-1)):(idx*8)]
    end
  elseif format==:long
    number_of_entries=(length(txt)-8)÷16+1
    data=Array{String}(undef,number_of_entries)
    data[1]=txt[1:8]
    for idx=2:number_of_entries
      data[idx]=txt[(1+16*(idx-1)-8):((idx*16)-8)]
    end
  elseif format==:free
    data=split(replace(txt," "=>""),",")
  end
  return data
end

function parse_nas_number(txt)
  # nastran numbers may or may not have sign symbols
  # and leading 0
  # also there might be lower or upper Case indicating
  # scientific notation or just a directly a sign symbol
  # This means, everything that is correctly parsed as
  # as a julia Float is a nastran float
  # but also something exotic like "-.314+7" is valid
  # nastran and would be the same as "-0.314E+7"
  # this routine correctly parses nastran numbers

  txt=strip(txt)
  E=false
  ins=0
  if !occursin("E",txt) && !occursin("e",txt)
    for (idx,c) in enumerate(txt)
      if c=="E"
        E=true
      end
      if idx!=1 && ( (!E && c=='-') || (!E && c=='+'))
        ins=idx
        break
      end
    end
    if ins >1
      txt=txt[1:(ins-1)]*"E"*txt[ins:end]
    end
  end
  return parse(Float64,txt)
end

# function relabel_nas_domains!(domains,filename)
#   open(filename) do fid
#     while !eof(fid)
#       line=readline(fid)
#       if line[1]=='#' #line is a comment
#         continue
#       end
#       data=split(line,",")
#       old_idx=strip(data[1])
#       new_idx=strip(data[2])
#       if old_idx in keys(domains)
#         domains[new_idx]=domains[old_idx]
#         delete!(domains,old_idx)
#       else
#         println("Warning: key `$old_idx` not in domains!")
#       end
#     end
#   end
#   return domains
# end


text:
"""
    samin(obj_fn, x, lb, ub)

Minimization of a function, over box constraints, via simulated annealing.
Call samin() for an example, and/or edit(samin,()) to see the code.

History: Bassed on Octave code samin.cc, also by Creel,
which was originally based on Gauss code by E.G. Tsionas. A source
for the Gauss code is http://web.stanford.edu/~doubleh/otherpapers/sa.txt
The original Fortran code by W. Goffe is at
http://www.degruyter.com/view/j/snde.1996.1.3/snde.1996.1.3.1020/snde.1996.1.3.1020.xml?format=INT
Tsionas and Goffe agreed to MIT licensing of samin.jl in email
messages to Creel.

This Julia code uses the same names for control variables,
for the most part. A notable difference is that the initial
temperature can be found automatically to ensure that the active
bounds when the temperature begins to reduce cover the entire
parameter space (defined as a n-dimensional rectangle that is the
Cartesian product of the(lb_i, ub_i), i = 1,2,..n. The code also
allows for parameters to be restricted, by setting lb_i = ub_i,
for the appropriate i.

References:
Goffe, William L. (1996) "SIMANN: A Global Optimization Algorithm
  using Simulated Annealing " Studies in Nonlinear Dynamics & Econometrics
  Oct96, Vol. 1 Issue 3.

Goffe, et. al. (1994) "Global Optimization of Statistical Functions
  with Simulated Annealing", Journal of Econometrics,
  V. 60, N. 1/2.

Usage details: x, obj, convergence, details = samin(f,
                                          x_init,
                                          lb,
                                          ub,
                                          nt,
                                          ns,
                                          rt,
                                          maxevals,
                                          neps,
                                          functol,
                                          paramtol,
                                          verbosity)
Arguments:
REQUIRED
* f: objective function
* x_init: starting value
* lb:  vector of lower bounds
* ub: vector of upper bounds
KEYWORDS
* nt:  integer: (default = 5) reduce temperature every nt*ns*dim(x_init) evaluations
* ns:  integer: (default = 5) adjust bounds every ns*dim(x_init) evaluations
* rt:  (0 < rt <1): (default = 0.5) geometric temperature reduction factor: when temp changes, new temp is t=rt*t
* maxevals:  integer: limit on function evaluations
* neps:  integer: (default = 5) number of previous best values the final result is compared to
* functol: (> 0): (default = 1e-8) the required tolerance level for function value
                   comparisons
* paramtol: (> 0): (default = 1e-5) the required tolerance level for parameters
* verbosity: scalar: 0, 1, 2 or 3 (default = 1).
    * 0 = no screen output
    * 1 = only final results to screen
    * 2 = summary every temperature change, without param values
    * 3 = summary every temperature change, with param values
* coverage_ok: (0 or 1) (default = 0) 0: increase temperature until parameter space is
        covered by the trial values. 1: start decreasing temperature immediately
Returns:
* x: the minimizer
* obj: the value of f() at x
* convergence:
    0 if no convergence within maxevals function evaluations
    1 if normal convergence to a point interior to the parameter space
    2 if convergence to point very near bounds of parameter space
      (suggest re-running with looser bounds)
* details: a px3 matrix. p is the number of times improvements were found.
           The columns record information at the time an improvement was found
           * first: cumulative number of function evaluations
           * second: temperature
           * third: function value

"""

function samin()
    println("samin(), with no arguments, runs an example of usage")
    println("execute edit(samin,()) to see the example")

    junk=2. # shows use of obj. fun. as a closure
    function sse(x)
        objvalue = junk + sum(x.*x)
    end
    k = 5
    x = rand(k,1)
    lb = -ones(k,1)
    ub = -lb
    # converge to global opt, see final parameters
    println("normal convergence, see only final results")
    @time xopt = samin(sse, x, lb, ub, verbosity=1)
    # no convergence within iteration limit
    println("no convergence within iter limit")
    @time xopt = samin(sse, x, lb, ub, maxevals=10, verbosity=1)
    # initial point out of bounds
    println("initial point out of bounds")
    lb = 0.5*ub
    x[1,1] = 0.2
    xopt = samin(sse, x, lb, ub, verbosity=1)
    # optimum on bound of parameter space
    println("optimum on bounds of parameter space")
    x = 0.5 .+ 0.5*rand(k,1)
    xopt = samin(sse, x, lb, ub, verbosity=1)
end

function samin(obj_fn, x, lb, ub; nt=5, ns=5, rt=0.5, maxevals=1e6, neps=5, functol=1e-8, paramtol=1e-5, verbosity=1, coverage_ok=0)
    n = size(x,1) # dimension of parameter
    #  Set initial values
    nacc = 0 # total accepted trials
    t = 2.0 # temperature - will initially rise or fall to cover parameter space. Then it will fall
    converge = 0 # convergence indicator 0 (failure), 1 (normal success), or 2 (convergence but near bounds)
    # most recent values, to compare to when checking convergend
    fstar = typemax(Float64)*ones(neps)
    # Initial obj_value
    xopt = copy(x)
    f = obj_fn(x)
    fopt = copy(f) # give it something to compare to
    func_evals = 0 # total function evaluations (limited by maxeval)
    details = [func_evals t fopt]
    bounds = ub - lb
    # check for out-of-bounds starting values
    for i = 1:n
        if(( x[i] > ub[i]) || (x[i] < lb[i]))
            @printf("samin: initial parameter %d out of bounds\n", i)
            converge = 0
            return xopt, fopt, converge, details
        end
    end
    # main loop, first increase temperature until parameter space covered, then reduce until convergence
    while (converge==0)
        # statistics to report at each temp change, set back to zero
        nup = 0
        nrej = 0
        nnew = 0
        ndown = 0
        lnobds = 0
        # repeat nt times then adjust temperature
        for m = 1:nt
            # repeat ns times, then adjust bounds
            nacp = zeros(n)
            for j = 1:ns
                # generate new point by taking last and adding a random value
                # to each of elements, in turn
                for h = 1:n
                    # new Sept 2011, if bounds are same, skip the search for that vbl.
                    # Allows restrictions without complicated programming
                    if (lb[h] != ub[h])
                        xp = copy(x)
                        xp[h] += (2.0 * rand() - 1.0) * bounds[h]
                        if((xp[h] < lb[h]) || (xp[h] > ub[h]))
                            xp[h] = lb[h] + (ub[h] - lb[h]) * rand()
                            lnobds += 1
                        end
                        # Evaluate function at new point
                        fp = obj_fn(xp)
                        func_evals += 1
                        #  Accept the new point if the function value decreases
                        if (fp <= f)
                            x = copy(xp)
                            f = copy(fp)
                            nacc += 1 # total number of acceptances
                            nacp[h] += 1 # acceptances for this parameter
                            nup += 1
                            #  If lower than any other point, record as new optimum
                            if(fp < fopt)
                                xopt = copy(xp)
                                fopt = copy(fp)
                                nnew +=1
                                details = [details; [func_evals t fp]]
                            end
                        # If the point is higher, use the Metropolis criteria to decide on
                        # acceptance or rejection.
                        else
                            p = exp(-(fp - f) / t)
                            if(rand() < p)
                                x = copy(xp)
                                f = copy(fp)
                                nacc += 1
                                nacp[h] += 1
                                ndown += 1
                            else
                                nrej += 1
                            end
                        end
                        # If maxevals exceeded, terminate the algorithm
                        if(func_evals >= maxevals)
                            if(verbosity >= 1)
                                println("================================================")
                                println("SAMIN results")
                                println("NO CONVERGENCE: MAXEVALS exceeded")
                                @printf("\n     Obj. value:  %16.5f\n\n", fopt)
                                if(verbosity >=2)
                                    println("       parameter      search width")
                                    for i=1:n
                                        @printf("%16.5f  %16.5f \n", xopt[i], bounds[i])
                                    end
                                end
                                println("================================================")
                            end
                            converge = 0
                            return xopt, fopt, converge, details
                        end
                    end
                end
            end
            #  Adjust bounds so that approximately half of all evaluations are accepted
            test = 0
            for i = 1:n
                if (lb[i] != ub[i])
                    ratio = nacp[i] / ns
                    if(ratio > 0.6) bounds[i] = bounds[i] * (1.0 + 2.0 * (ratio - 0.6) / 0.4) end
                    if(ratio < .4) bounds[i] = bounds[i] / (1.0 + 2.0 * ((0.4 - ratio) / 0.4)) end
                    # keep within initial bounds
                    if(bounds[i] > (ub[i] - lb[i]))
                        bounds[i] = ub[i] - lb[i]
                        test += 1
                    end
                else
                    test += 1 # make sure coverage check passes for the fixed parameters
                end
            end
            nacp = 0 # set back to zero
            # check if we cover parameter space, if we have yet to do so
            if (coverage_ok != 1) coverage_ok = (test == n) end
        end

        # intermediate output, if desired
        if(verbosity >= 2)
            println("samin: intermediate results before next temperature change")
            println("temperature: ", round(t,5))
            println("current best function value: ", round(fopt,5))
            println("total evaluations so far: ", func_evals)
            println("total moves since last temperature reduction: ", nup + ndown + nrej)
            println("downhill: ", nup)
            println("accepted uphill: ", ndown)
            println("rejected uphill: ", nrej)
            println("out of bounds trials: ", lnobds)
            println("new minima this temperature: ", nnew)
            println()
            if (verbosity > 2)
                println("       parameter      search width")
                for i=1:n
                    @printf("%16.5f  %16.5f \n", xopt[i], bounds[i])
                end
            end
            println()
        end
        # Check for convergence, if we have covered the parameter space
        if (coverage_ok)
            # last value close enough to last neps values?
            fstar[1] = f
            test = 0
            for i=1:neps
                test += (abs(f - fstar[i]) > functol)
            end
            test = (test > 0) # if different from zero, function conv. has failed
            # last value close enough to overall best?
            if (((fopt - f) <= functol) && (!test))
                # check for bound narrow enough for parameter convergence
                for i = 1:n
                    if (bounds[i] > paramtol)
                        converge = 0 # no conv. if bounds too wide
                        break
                    else
                        converge = 1
                    end
                end
            end
            # check if optimal point is near boundary of parameter space, and change message if so
            if (converge == 1) && (lnobds > 0)
                converge = 2
            end
            # Like to see the final results?
            if (converge > 0)
                if (verbosity >= 1)
                    println("================================================")
                    println("SAMIN results")
                    if (converge == 1)
                        println("==> Normal convergence <==")
                    end
                    if (converge == 2)
                        println("==> WARNING <==: Last point satisfies convergence criteria,")
                        println("but is near boundary of parameter space.")
                        println(lnobds, " out of  ", (nup+ndown+nrej), " evaluations were out of bounds in the last round.")
                        println("Expand bounds and re-run, unless this is a constrained minimization.")
                    end
                    println("total number of objective function evaluations: ", func_evals)
                    @printf("\n     Obj. value:  %16.10f\n\n", fopt)
                    println("       parameter      search width")
                    for i=1:n
                        @printf("%16.5f  %16.5f \n", xopt[i], bounds[i])
                    end
                    println("================================================")
                end
            end
            # Reduce temperature, record current function value in the
            # list of last "neps" values, and loop again
            t *= rt
            for i = neps:-1:2
                fstar[i] = fstar[i-1]
            end
            f = copy(fopt)
            x = copy(xopt)
        else  # coverage not ok - increase temperature quickly to expand search area
            t *= 10.0
            for i = neps:-1:2
                fstar[i] = fstar[i-1]
            end
            f = fopt
            x = xopt
        end
    end
    return xopt, fopt, converge, details
end


text:
# Performance tests for InterestRates module

println("Running perftests...")

using Test
using InterestRates

function run_perftests()

    dt_curve = Date(2015,08,08)

    vert_x = [11, 15, 19, 23, 25, 200, 500]
    vert_y = [0.10, 0.12, 0.20, 0.19, 0.21, 0.21, 0.22]
    NS_params = [0.1, 0.2, 0.3, 0.5]
    SVEN_params = [0.1, 0.2, 0.3, 0.4, 0.5, 0.8]

    days_vec = [Dates.Day(i) for i=1:1500]

    mat_vec = dt_curve + days_vec

    curve_linear = InterestRates.IRCurve("dummy-linear", InterestRates.Actual360(),
        InterestRates.ExponentialCompounding(), InterestRates.Linear(), dt_curve,
        vert_x, vert_y)

    curve_flatforward = InterestRates.IRCurve("dummy-flatforward", InterestRates.Actual360(),
        InterestRates.ExponentialCompounding(), InterestRates.FlatForward(), dt_curve,
        vert_x, vert_y)

    curve_ns = InterestRates.IRCurve("dummy-ns", InterestRates.Actual360(),
        InterestRates.ExponentialCompounding(), InterestRates.NelsonSiegel(), dt_curve,
        NS_params)

    curve_sven = InterestRates.IRCurve("dummy-sven", InterestRates.Actual360(),
        InterestRates.ExponentialCompounding(), InterestRates.Svensson(), dt_curve,
        SVEN_params)

    curve_spline_rates = InterestRates.IRCurve("dummy-spline_rates", InterestRates.Actual360(),
        InterestRates.ExponentialCompounding(), InterestRates.CubicSplineOnRates(), dt_curve,
        vert_x, vert_y)

    curve_spline_df = InterestRates.IRCurve("dummy-spline_df", InterestRates.Actual360(),
        InterestRates.ExponentialCompounding(), InterestRates.CubicSplineOnDiscountFactors(), dt_curve,
        vert_x, vert_y)

    curve_step = InterestRates.IRCurve("dummy-step", InterestRates.Actual360(),
        InterestRates.ExponentialCompounding(), InterestRates.StepFunction(), dt_curve,
        vert_x, vert_y)

    c_array = [curve_linear, curve_flatforward, curve_ns, curve_sven, curve_spline_rates, curve_spline_df, curve_step]

    # Warm up
    for c in c_array
        zero_rate(c, mat_vec)
        ERF(c, mat_vec)
        discountfactor(c, mat_vec)
    end

    # Check results
    for c in c_array
        zr = zero_rate(c, mat_vec)
        for i in 1:length(zr)
            @test zr[i] == zero_rate(c, mat_vec[i])
        end

        erf = ERF(c, mat_vec)
        for i in 1:length(erf)
            @test erf[i] == ERF(c, mat_vec[i])
        end

        df = discountfactor(c, mat_vec)
        for i in 1:length(df)
            @test df[i] == discountfactor(c, mat_vec[i])
        end
    end

    sp = InterestRates.splinefit(vert_x, vert_y)
    InterestRates.splineint(sp, convert(Vector{Int}, 1:30))

    # Perftests
    for c in c_array
        println("$(InterestRates.curve_get_method(c))")
        @time for i=1:100 zero_rate(c, mat_vec) end
        @time for i=1:100 ERF(c, mat_vec) end
        @time for i=1:100 discountfactor(c, mat_vec) end
    end
end

run_perftests()

println("Perftests end")


text:
function test()
    if (i = 2) < 3
        @assert(i == 2)
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    test()
end


text:
function solve{islinear,isstochastic,MeshType<:FEMMesh,F1,F2,F3,F4,F5,F6,F7,DiffType}(
  prob::PoissonProblem{islinear,isstochastic,MeshType,F1,F2,F3,F4,F5,F6,F7,DiffType},
  alg::FEMDiffEqPoisson;solver::Symbol=:Direct,autodiff::Bool=false,
  method=:trust_region,show_trace=false,iterations=1000,kwargs...)
  #Assemble Matrices
  A,M,area = assemblematrix(prob.mesh,lumpflag=true)
  #Unroll some important constants
  @unpack dt,bdnode,node,elem,N,NT,freenode,dirichlet,neumann = prob.mesh
  @unpack f,Du,f,gD,gN,analytic,u0,numvars,σ,noisetype,D = prob

  #Setup f quadrature
  mid = Array{eltype(node)}(size(@view node[vec(@view elem[:,2]),:])...,3)
  mid[:,:,1] .= (view(node,vec(@view elem[:,2]),:) .+ @view node[vec(@view elem[:,3]),:])./2
  mid[:,:,2] .= (view(node,vec(@view elem[:,3]),:) .+ @view node[vec(@view elem[:,1]),:])./2
  mid[:,:,3] .= (view(node,vec(@view elem[:,1]),:) .+ @view node[vec(@view elem[:,2]),:])./2

  #Setup u
  u = copy(u0)
  u[bdnode] = gD(@view node[bdnode,:])

  #isstochastic Part
  if isstochastic
    rands = getNoise(u,node,elem,noisetype=noisetype)
    dW = next(rands)
    rhs = (u) -> quadfbasis(f,gD,gN,A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars) + quadfbasis(σ,gD,gN,A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars).*dW
  else #Not isstochastic
    rhs = (u) -> quadfbasis(f,gD,gN,A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars)
  end
  Dinv = map((x)->inv(x),D) # Special form so that units work
  #Solve
  if islinear
    if solver==:Direct
      u[freenode,:]=A[freenode,freenode]\(Dinv.*rhs(u))[freenode]
    elseif solver==:CG
      for i = 1:size(u,2)
        u[freenode,i],ch=cg!(u[freenode,i],A[freenode,freenode],Dinv.*rhs(u)[freenode,i]) # Needs diffusion constant
      end
    elseif solver==:GMRES
      for i = 1:size(u,2)
        u[freenode,i],ch=gmres!(u[freenode,i],A[freenode,freenode],Dinv.*rhs(u)[freenode,i]) # Needs diffusion constants
      end
    end
    #Adjust result
    if isempty(dirichlet) #isPureneumann
      patchArea = accumarray(vec(elem),[area;area;area]/3, [N 1])
      uc = sum(u.*patchArea)/sum(area)
      u = u - uc   # Impose integral of u = 0
    end
  else #Nonlinear
    rhs! = (resid,u) -> begin
      u = reshape(u,N,numvars)
      resid = reshape(resid,N,numvars)
      resid[freenode,:]=D.*(A[freenode,freenode]*u[freenode,:])-rhs(u)[freenode,:]
      u = vec(u)
      resid = vec(resid)
    end
    u = vec(u)
    nlres = NLsolve.nlsolve(rhs!,u,autodiff=autodiff,method=method,show_trace=show_trace,iterations=iterations)
    u = nlres.zero
    if numvars > 1
      u = reshape(u,N,numvars)
    end
  end

  #Return
  if typeof(analytic)!=Void # True solution exists
    return(FEMSolution([u],analytic(node),Du,prob))
  else #No true solution
    return(FEMSolution([u],prob))
  end
end

## Evolution Equation Solvers
#Note
#rhs(u,i) = Dm[freenode,freenode]*u[freenode,:] + dt*f(node,(i-.5)*dt)[freenode] #Nodel interpolation 1st 𝒪
function solve{islinear,isstochastic,MeshType<:FEMMesh,F,F2,F3,F4,F5,F6,F7,DiffType}(
  prob::HeatProblem{islinear,isstochastic,MeshType,F,F2,F3,F4,F5,F6,F7,DiffType},
  alg::AbstractHeatFEMAlgorithm;
  solver::Symbol=:LU,save_everystep::Bool = false,timeseries_steps::Int = 100,
  autodiff::Bool=false,method=:trust_region,show_trace=false,iterations=1000,
  progress_steps::Int=1000,progressbar::Bool=false,progressbar_name="FEM",kwargs...)
  #Assemble Matrices
  A,M,area = assemblematrix(prob.mesh,lumpflag=true)

  #Unroll some important constants
  @unpack dt,tspan,bdnode,node,elem,N,NT,freenode,dirichlet,neumann = prob.mesh
  @unpack f,u0,Du,gD,gN,analytic,numvars,σ,noisetype,D = prob

  if dt != 0
    numiters = round(Int64,tspan[2]/dt)
  else
    numiters = 0
  end

  #Set Initial
  u = copy(u0)
  t = tspan[1]

  #Setup timeseries

  timeseries = Vector{typeof(u)}(0)
  push!(timeseries,copy(u))
  ts = Float64[t]

  sqrtdt= sqrt(dt)
  #Setup f quadraturef
  mid = Array{eltype(node)}(size(node[vec(elem[:,2]),:])...,3)
  mid[:,:,1] = (node[vec(elem[:,2]),:]+node[vec(elem[:,3]),:])/2
  mid[:,:,2] = (node[vec(elem[:,3]),:]+node[vec(elem[:,1]),:])/2
  mid[:,:,3] = (node[vec(elem[:,1]),:]+node[vec(elem[:,2]),:])/2

  #Setup for Calculations
  Minv = sparse(inv(M)) #sparse(Minv) needed until update

  #Heat Equation Loop
  u,timeseries,ts=femheat_solve(FEMHeatIntegrator{islinear,typeof(alg),isstochastic,typeof(f),typeof(gD),typeof(gN),typeof(σ),eltype(u),typeof(u),typeof(node),typeof(area),typeof(t),typeof(D),typeof(Minv),typeof(A)}(N,NT,dt,t,Minv,D,A,freenode,f,gD,gN,u,node,elem,area,bdnode,mid,dirichlet,neumann,islinear,numvars,sqrtdt,σ,noisetype,numiters,save_everystep,timeseries,ts,solver,autodiff,method,show_trace,iterations,timeseries_steps,progressbar,progress_steps,progressbar_name))

  if typeof(analytic)!=Void #True Solution exists
    return(FEMSolution(timeseries,analytic(tspan[2],node),Du,ts,prob))
  else #No true solution
    return(FEMSolution(timeseries,ts,prob))
  end
end

"""
`quadfbasis(f,gD,gN,A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars;gNquad𝒪=2)`

Performs the order 2 quadrature to calculate the vector from the term ``<f,v>`` for linear elements.
"""
function quadfbasis(f,gD,gN,A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars;gNquad𝒪=2)
  b = zeros(area[1].*u) #size(bt1,2) == numvars #area is for units
  if islinear
    bt1 = area.*(f(mid[:,:,2])+f(mid[:,:,3]))/6
    bt2 = area.*(f(mid[:,:,3])+f(mid[:,:,1]))/6
    bt3 = area.*(f(mid[:,:,1])+f(mid[:,:,2]))/6
  else
    u1 = (u[vec(elem[:,2]),:]+u[vec(elem[:,3]),:])/2
    u2 = (u[vec(elem[:,3]),:]+u[vec(elem[:,1]),:])/2
    u3 = (u[vec(elem[:,1]),:]+u[vec(elem[:,2]),:])/2
    bt1 = area.*(f(mid[:,:,2],u2)+f(mid[:,:,3],u3))/6
    bt2 = area.*(f(mid[:,:,3],u3)+f(mid[:,:,1],u1))/6
    bt3 = area.*(f(mid[:,:,1],u1)+f(mid[:,:,2],u2))/6
  end

  for i = 1:numvars # accumarray the bt's
    for j = 1:NT
      b[elem[j,1],i] += bt1[j,i]
    end
    for j = 1:NT
      b[elem[j,2],i] += bt2[j,i]
    end
    for j = 1:NT
      b[elem[j,3],i] += bt3[j,i]
    end
  end

  if(!isempty(dirichlet))
    uz = zeros(b)
    uz[bdnode,:] = gD(node[bdnode,:])
    b = b-A*uz
  end

  if(!isempty(neumann))
    el = sqrt.(float(sum((node[neumann[:,1],:] - node[neumann[:,2],:]).^2,2)))
    λgN,ωgN = quadpts1(gNquad𝒪)
    ϕgN = λgN                # linear bases
    nQuadgN = size(λgN,1)
    ge = zeros(size(neumann,1),2,numvars) # Does this need units?

    for pp = 1:nQuadgN
        # quadrature points in the x-y coordinate
        ppxy = λgN[pp,1]*node[neumann[:,1],:] +
               λgN[pp,2]*node[neumann[:,2],:]
        gNp = gN(ppxy)
        for igN = 1:2, var = 1:numvars
            ge[:,igN,var] = ge[:,igN,var] + vec(ωgN[pp]*ϕgN[pp,igN]*gNp[:,var])
        end
    end
    ge = ge.*repeat(el,outer=[1,2,numvars]) # tuple in v0.5?

    for i=1:numvars
      b[:,i] = b[:,i] + accumarray(vec(neumann), vec(ge[:,i]),[N,1])
    end
  end
  if numvars == 1
    b = vec(b)
  end
  return(b)
end


text:
export
    FourierNeuralOperator,
    MarkovNeuralOperator

"""
    FourierNeuralOperator(;
        ch=(2, 64, 64, 64, 64, 64, 128, 1),
        modes=(16, ),
        σ=gelu
    )

Fourier neural operator learns a neural operator with Dirichlet kernel to form a Fourier transformation.
It performs Fourier transformation across infinite-dimensional function spaces and learns better than neural operator.
"""
function FourierNeuralOperator(;
    ch=(2, 64, 64, 64, 64, 64, 128, 1),
    modes=(16, ),
    σ=gelu
)
    Transform = FourierTransform

    return Chain(
        Dense(ch[1], ch[2]),
        OperatorKernel(ch[2]=>ch[3], modes, Transform, σ),
        OperatorKernel(ch[3]=>ch[4], modes, Transform, σ),
        OperatorKernel(ch[4]=>ch[5], modes, Transform, σ),
        OperatorKernel(ch[5]=>ch[6], modes, Transform),
        Dense(ch[6], ch[7], σ),
        Dense(ch[7], ch[8]),
    )
end

"""
    MarkovNeuralOperator(;
        ch=(1, 64, 64, 64, 64, 64, 1),
        modes=(24, 24),
        σ=gelu
    )

Markov neural operator learns a neural operator with Fourier operators.
With only one time step information of learning, it can predict the following few steps with low loss
by linking the operators into a Markov chain.
"""
function MarkovNeuralOperator(;
    ch=(1, 64, 64, 64, 64, 64, 1),
    modes=(24, 24),
    σ=gelu
)
    Transform = FourierTransform

    return Chain(
        Dense(ch[1], ch[2]),
        OperatorKernel(ch[2]=>ch[3], modes, Transform, σ),
        OperatorKernel(ch[3]=>ch[4], modes, Transform, σ),
        OperatorKernel(ch[4]=>ch[5], modes, Transform, σ),
        OperatorKernel(ch[5]=>ch[6], modes, Transform, σ),
        Dense(ch[6], ch[7]),
    )
end


text:
#= Electronic Structure
Implements types for dealing with electronic configurations of atoms.

TODO: See if tuple implementation of orbitals needs to be improved on.

=#

struct Electrons
    root::Union{String, Nothing}
    orbitals::Tuple
end

function parse_config(config::String)::Electrons
    config_pattern = r"(?:\[(\w+)\])?((?:\w|\s?)+)"
    captures = match(config_pattern, config).captures
    orbitals = sort(split(captures[2]), by = o -> o[1])

    Electrons(
        captures[1],
        tuple(
            map(o -> tuple(o[1:2], parse(Int, o[3:end])), orbitals)...
        )
    )
end

Electrons(number::Int) = begin
    config = periodic_table_electrons[number]
    parse_config(config)
end

Electrons(symbol::String) =  begin
    number = periodic_table[symbol]
    Electrons(number)
end

Base.show(io::IO, e::Electrons) = print(io, "[", e.root, "] ", e.orbitals)

# Base.show(io::IO, ::MIME"text/plain", e::Electrons) =
#            print(io, "$(typeof(a)):\n   ", e)


text:
Base.@kwdef mutable struct Config <: AMGXObject
    handle::API.AMGX_config_handle = API.AMGX_config_handle(C_NULL)
    ref_count::Threads.Atomic{Int} = Threads.Atomic{Int}(0)
    function Config(handle::API.AMGX_config_handle, ref_count::Threads.Atomic{Int})
        config = new(handle, ref_count)
        finalizer(warn_not_destroyed_on_finalize, config)
        return config
    end
end
get_api_destroy_call(::Type{Config}) = API.AMGX_config_destroy
dec_refcount_parents(config::Config) = nothing

function create!(config::Config, content::String)
    cfg_handle_ptr = Ref{API.AMGX_config_handle}()
    @checked API.AMGX_config_create(cfg_handle_ptr, content)
    config.handle = cfg_handle_ptr[]
    return config
end
Config(content::String) = create!(Config(), content)

function create!(config::Config, d::Dict)
    str = sprint(JSON.print, d)
    create!(config, str)
end
Config(d::Dict) = create!(Config(), d)


text:
import Pkg

function runapp(name::AbstractString)
    folder = joinpath(@__DIR__, "src", name)
    if isdir(folder)
        @info "Running $name"
        Pkg.activate(folder)
        push!(LOAD_PATH, joinpath(@__DIR__, "src"))
        include( joinpath(folder, "src", "$name.jl") )
        return true
    end
    return false
end

if !isempty(ARGS)
    name = ARGS[1]
    if !runapp(name)
        @error "App '$name' does not exist"
    end
else
    @error """No name argument provided

    Usage: julia run.jl <appname>
    Example: julia run.jl Asteroids

    """
end

max_stars_repo_name: jorge-brito/LynxGallery

max_issues_repo_name: jorge-brito/LynxGallery

max_forks_repo_name: jorge-brito/LynxGallery


text:
module GSDicts

using HTTP
using GoogleCloud
using JSON
import BigArrays.BackendBase: AbstractBigArrayBackend, get_info, get_scale_name

export GSDict, get_credential_filename

const GZIP_MAGIC_NUMBER = UInt8[0x1f, 0x8b, 0x08]

struct GSDict <: AbstractBigArrayBackend
    bucketName  	::String
    keyPrefix   	::String
    session	::GoogleCloud.session.GoogleSession
end

"""
    GSDict( path::String; gzip::Bool = GZIP )
construct an associative datastructure based on Google Cloud Storage
format.
"""
function GSDict(path::String;
                credentialFileName = get_credential_filename())
    @assert startswith(path, "gs://")
    bucketName, keyPrefix = splitgs(path)
    bucketName = replace(bucketName, "gs://"=>"")

    session = GoogleSession( credentialFileName, ["devstorage.full_control"])
    set_session!(storage, session)    # storage is the API root, exported from GoogleCloud.jl
    GSDict( bucketName, keyPrefix, session )
end

##################### properties ##############################
function get_info( self::GSDict )
    storage(:Object, :get, self.bucketName, joinpath(self.keyPrefix, "info"))
end

function get_scale_name(self::GSDict)
    basename( self.keyPrefix )
end
"""
get the voxel offset in info
"""
function get_voxel_offset(h::GSDict)
    h.configDict[:offset]
end

###################### Base functions #########################
@inline function Base.delete!( d::GSDict, key::String )
    storage(:Object, :delete, d.bucketName, joinpath(d.keyPrefix, key))
end

function Base.setindex!( d::GSDict, value::Vector{UInt8}, key::AbstractString )
	#authorize( d.googleSession )
    if all(value[1:3] .== GZIP_MAGIC_NUMBER)
        return storage(:Object, :insert, d.bucketName;
                        name=joinpath(d.keyPrefix, key),
                        data=value, gzip=true,
                        content_type="application/octet-stream",
                        fields="")
    else
        response = storage(:Object, :insert, d.bucketName;
                        name=joinpath(d.keyPrefix, key),
                        data=value, gzip=false,
                        content_type="application/octet-stream",
                        fields="")
    end
end

function Base.setindex!( d::GSDict, value::Dict, key::AbstractString )
	#authorize( d.googleSession )
    response = storage(:Object, :insert, d.bucketName;
                       name=joinpath(d.keyPrefix, key),
                       data=JSON.json(value), gzip=false,
                       content_type="application/json")
end

function Base.setindex!( d::GSDict, value::AbstractString, key::AbstractString )
	#authorize( d.googleSession )
    response = storage(:Object, :insert, d.bucketName;
                       name=joinpath(d.keyPrefix, key),
                       data=value, gzip=false,
                       content_type="text/plain")
end

function Base.getindex( d::GSDict, key::AbstractString)
    try

        return storage(:Object, :get, d.bucketName, joinpath(d.keyPrefix, key))
    catch err
        if isa(err, HTTP.ExceptionRequest.StatusError) && err.status==404
            throw(KeyError("NoSuchKey in Google Cloud Storage: $(key)"))
        else
            println("get an unknown error: ", err)
            println("error type is: ", typeof(err))
            rethrow
        end
    end
end

function Base.keys( d::GSDict )
    ds = storage(:Object, :list, d.bucketName; prefix=d.keyPrefix, fields="items(name)")
    ret = Vector{String}()
    for i in eachindex(ds)
        chunkFileName = replace(ds[i][:name], "$(rstrip(d.keyPrefix, '/'))/" => "" )
        push!(ret, chunkFileName)
    end
    return ret
end

function Base.haskey( d::GSDict, key::String )
    @warn("this haskey function will download the object rather than just check whether it exist or not")
    response = storage(:Object, :get, d.bucketName, joinpath(d.keyPrefix, key))
    !GoogleCloud.api.iserror(response)
end

################### utility functions #################

function get_credential_filename()
    if isfile(expanduser("~/.google_credentials.json"))
        return expanduser("~/.google_credentials.json")
    elseif isfile(joinpath(dirname(@__FILE__), "../.google_credentials.json"))
        return joinpath(dirname(@__FILE__), "../.google_credentials.json")
    elseif isfile(expanduser("~/.cloudvolume/secrets/google-secret.json"))
        return expanduser("~/.cloudvolume/secrets/google-secret.json")
    elseif isfile("/secrets/google-secret.json")
        return "/secrets/google-secret.json"
    else
        # to enable building of this package
        @warn("google credential file is not in default place!")
        return nothing
    end
end

"""
    split gs path to bucket name and key
"""
function splitgs( path::String )
    path = replace(path, "gs://"=>"")
    bucketName, key = split(path, "/", limit=2)
    key = rstrip(key, '/')
    return String(bucketName), String(key)
end

end # end of module


text:
struct fitresult
        beta
        residuals
        mse
        converged::Bool

        function fitresult(modelfun, X, Y, optimresult::AbsRes)
            beta = Optim.minimizer(optimresult)
            residuals = (Y - modelfun(X, beta)).^2
            mse = Optim.minimum(optimresult)
            converged = Optim.converged(optimresult)
            new(beta, residuals, mse, converged)
        end
end

function nlinfit(modelfun, X, Y, beta0, method::AbsOptim=NewtonTrustRegion())
        if true # test__modelfun(modelfun, X[1:2,:], beta0)

            # Residual function given the model:
            residual(beta_) = sum((Y - modelfun(X, beta_)).^2)

            optimresult = optimize(residual, beta0, method)
            result = fitresult(modelfun, X, Y, optimresult)

            return result

        else

            error("Something is wrong!")

        end

end

function test__modelfun(modelfun, x, b)
        try
            modelfun(x, b)
            return true
        catch
            error("modelfun is not correctly defined")
        end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct PeerExpressRouteCircuitConnectionsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Gets the specified Peer Express Route Circuit Connection from the specified express route circuit.
Param: resourceGroupName::String (required)
Param: circuitName::String (required)
Param: peeringName::String (required)
Param: connectionName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PeerExpressRouteCircuitConnection
"""
function _swaggerinternal_peerExpressRouteCircuitConnectionsGet(_api::PeerExpressRouteCircuitConnectionsApi, resourceGroupName::String, circuitName::String, peeringName::String, connectionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PeerExpressRouteCircuitConnection, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/expressRouteCircuits/{circuitName}/peerings/{peeringName}/peerConnections/{connectionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "circuitName", circuitName)  # type String
    Swagger.set_param(_ctx.path, "peeringName", peeringName)  # type String
    Swagger.set_param(_ctx.path, "connectionName", connectionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function peerExpressRouteCircuitConnectionsGet(_api::PeerExpressRouteCircuitConnectionsApi, resourceGroupName::String, circuitName::String, peeringName::String, connectionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_peerExpressRouteCircuitConnectionsGet(_api, resourceGroupName, circuitName, peeringName, connectionName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function peerExpressRouteCircuitConnectionsGet(_api::PeerExpressRouteCircuitConnectionsApi, response_stream::Channel, resourceGroupName::String, circuitName::String, peeringName::String, connectionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_peerExpressRouteCircuitConnectionsGet(_api, resourceGroupName, circuitName, peeringName, connectionName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets all global reach peer connections associated with a private peering in an express route circuit.
Param: resourceGroupName::String (required)
Param: circuitName::String (required)
Param: peeringName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PeerExpressRouteCircuitConnectionListResult
"""
function _swaggerinternal_peerExpressRouteCircuitConnectionsList(_api::PeerExpressRouteCircuitConnectionsApi, resourceGroupName::String, circuitName::String, peeringName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PeerExpressRouteCircuitConnectionListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/expressRouteCircuits/{circuitName}/peerings/{peeringName}/peerConnections", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "circuitName", circuitName)  # type String
    Swagger.set_param(_ctx.path, "peeringName", peeringName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function peerExpressRouteCircuitConnectionsList(_api::PeerExpressRouteCircuitConnectionsApi, resourceGroupName::String, circuitName::String, peeringName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_peerExpressRouteCircuitConnectionsList(_api, resourceGroupName, circuitName, peeringName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function peerExpressRouteCircuitConnectionsList(_api::PeerExpressRouteCircuitConnectionsApi, response_stream::Channel, resourceGroupName::String, circuitName::String, peeringName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_peerExpressRouteCircuitConnectionsList(_api, resourceGroupName, circuitName, peeringName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export peerExpressRouteCircuitConnectionsGet, peerExpressRouteCircuitConnectionsList


text:
import PredictMD
import Test

Test.@test_throws ErrorException PredictMD.simple_linear_regression([1, 2], [])
Test.@test_throws ErrorException PredictMD.simple_linear_regression([], [])


text:
# code from https://github.com/FluxML/Zygote.jl/pull/294

import Flux.Zygote

reversesumscan_(x::AbstractArray; dims::Integer) =
    reverse(cumsum(reverse(x, dims=dims), dims=dims), dims=dims)

Zygote.@adjoint function cumprod(A::AbstractArray; dims::Integer)
    return cumprod(A; dims=dims), function (∆)
        dim_size = size(A, dims)
        if dim_size == 1
            return (∆, nothing)
        end
        ndims = length(size(A))

        # Simple case with nonzero elements in the input
        if all(A .!= 0)
            output = cumprod(A, dims=dims)
            return (reversesumscan_(∆ .* output, dims=dims) ./ A, nothing)
        end

        grad_input = similar(A)

        pre_idx = Any[Colon() for _=1:dims-1]
        post_idx = [Colon() for _=dims+1:ndims]
        for k=1:dim_size

            if k == 1
                idx_kplus1 = vcat(pre_idx, [k+1:dim_size], post_idx)
                prods_from_k_plus_1 = cumprod(A[idx_kplus1...], dims=dims)

                # We add 1s to the omitted_products of the first element
                ones_size = [size(A)...]
                ones_size[dims] = 1
                ones_ = ones(eltype(A), ones_size...)

                omitted_products = cat(ones_, prods_from_k_plus_1, dims=dims)

            elseif k == dim_size
                idx_kminus1 = vcat(pre_idx, [1:k-1], post_idx)
                prods_until_k = prod(A[idx_kminus1...], dims=dims)
                omitted_products = prods_until_k

            else
                idx_kplus1 = vcat(pre_idx, [k+1:dim_size], post_idx)
                prods_from_k_plus_1 = cumprod(A[idx_kplus1...], dims=dims)

                idx_kminus1 = vcat(pre_idx, [1:k-1], post_idx)
                prods_until_k = prod(A[idx_kminus1...], dims=dims)

                omitted_products = prods_until_k .* prods_from_k_plus_1
                omitted_products = cat(prods_until_k, omitted_products, dims=dims)
            end

            @assert size(omitted_products, dims) == dim_size - k + 1 "$(size(omitted_products, dims)) != $(dim_size - k - 1)"

            idx_ktodimsize = vcat(pre_idx, [k:dim_size], post_idx)
            idx_k = vcat(pre_idx, [k:k], post_idx)

            grad_input[idx_k...] .= sum(∆[idx_ktodimsize...] .* omitted_products, dims=dims)
        end
        return (grad_input, nothing)
    end
end


text:
abstract type AbstractWall end
abstract type AbstractParticle end
abstract type AbstractPoint end

mutable struct Point2D{T} <: AbstractPoint
    x::T
    y::T
end

mutable struct Point3D{T} <: AbstractPoint
    x::T
    y::T
    z::T
end

struct Particle2D{T} <: AbstractParticle
    pos::Point2D{T}
    radius::T
    materialID::Int
end

struct Particle3D{T} <: AbstractParticle
    pos::Point3D{T}
    radius::T
    materialID::Int
end

struct LineWall{T} <: AbstractWall
    nodes::Vector{Point2D{T}}     # 2 points, defining start and end
    thickness::T                  # gives line "area"
    n::Vector{T}                  # normal unit vector of wall
    t::Vector{T}                  # tangent unit vector of wall
end

struct CircleWall{T} <: AbstractWall
    center::Point2D{T}        # point that defines center
    radius::T                 # radius of circle
    thickness::T              # gives circle curve "area"
end

struct ArcWall{T} <: AbstractWall
    nodes::Vector{Point2D{T}}   # 3 nodes, 1) start, 2) center, and 3) end going counterclockwise
    thickness::T                # gives curve "area"
end

# initializes new wall defined by 2 points (Start,End)
function LineWall(nodes::Vector{Point2D{T}},thickness::T) where T<:Real
    if length(nodes) != 2
        throw(DimensionMismatch("line is defined by 2 points"))
    end

    x0 = [nodes[1].x,nodes[2].x]
    x1 = [nodes[1].y,nodes[2].y]

    n = zeros(Float64,2)
    t = zeros(Float64,2)

    # compute wall length
    WL = sqrt((nodes[2].x - nodes[1].x)^2 + (nodes[2].y - nodes[1].y)^2)

    # comute tangent
    t[1] = (nodes[2].x - nodes[1].x)/WL
    t[2] = (nodes[2].y - nodes[1].y)/WL

    n[1] = -t[2]
    n[2] =  t[1]

    return LineWall(nodes,thickness,n,t)
end
LineWall(n1,n2,thickness) = LineWall([n1,n2],thickness)

function ArcWall(p1,p2,p3,thickness::Float64)
    r1 = sqrt((p2.x-p1.x)^2 + (p2.y-p1.y)^2)
    r2 = sqrt((p2.x-p3.x)^2 + (p2.y-p3.y)^2)

    TOL = 1e-12
    if abs(r1-r2) > TOL
        error("points do not make a circle with constant radius, r1=$(round(r1,2)), r2=$(round(r2,2))")
    end

    return ArcWall([p1,p2,p3],thickness)
end

mutable struct Cell
    particleIDList::Vector{Int}     # list of indices of particles that are within cell
    square::Vector{Float64}         # bounding square that defines position of cell
    neighborList::Vector{Int}       # list of cell neighbors, used for fast access
end

struct CellList
    cells::Vector{Cell}             # array of cells -- this is main component
    bounds::Vector{Float64}         # total bounds of cell list in [x0, x1, y0, y1] format
    sideLength::Float64             # length of one cell inside cell list
end

mutable struct scratch_data
    n_particles::Int
    polygon::Vector{Point2D{Float64}}
    extremePoint::Point2D{Float64}
    pointOnWall::Point2D{Float64}
    cfX::Vector{Float64}
    cfY::Vector{Float64}
    afX::Vector{Float64}
    afY::Vector{Float64}
    gfX::Vector{Float64}
    gfY::Vector{Float64}

    function scratch_data(n_particles::Int)
        polygon      = [Point2D(0.0,0.0),Point2D(0.0,0.0),Point2D(0.0,0.0),Point2D(0.0,0.0)]
        extremePoint = Point2D(100_000.0,0.0)
        pointOnWall  = Point2D(0.0,0.0)
        cfX          = zeros(n_particles)
        cfY          = zeros(n_particles)
        afX          = zeros(n_particles)
        afY          = zeros(n_particles)
        gfX          = zeros(n_particles)
        gfY          = zeros(n_particles)
        return new(n_particles,polygon,extremePoint,pointOnWall,cfX,cfY,afX,afY,gfX,gfY)
    end
end


text:
# This file is part of the TaylorIntegration.jl package; MIT licensed

module TaylorIntegration

using TaylorSeries

import Compat.view

export  taylorinteg

include("integration_methods.jl")

end #module


text:
using LAP_julia: load_H, make_flow_from_H, gen_rand_points

function get_valid_landmarks(flow, orig_inds)
    rng = extrema.(indices_spatial(flow))
    orig_points = transpose(LAP_julia.inds_to_points(orig_inds))
    moved_points = LAP_julia.move_landmarks(orig_points, flow.*(-1))
    stayed_in_mask = is_in_bounds.(moved_points[:,1], rng[1]...) .& is_in_bounds.(moved_points[:,2], rng[2]...)
    filtered_orig_points = orig_points[stayed_in_mask, :]
    filtered_orig_inds = points_to_inds(transpose(filtered_orig_points))
end

## create a birl like structure for with the oxford affine dataset

base_path = "/Users/MrTrololord/Downloads/ox_affine_separate/"
folders = filter(x -> !occursin(".", x), readdir(joinpath(base_path)))

folder = folders[2]

for folder in folders

#converted to png

    # my_path = joinpath(base_path, folder)
    # for file in filter(x -> occursin("ppm", x), readdir(my_path))
    #     @info file
    #     img_path = joinpath(my_path, file)
    #     file_name = split(file, ".")[1]
    #     out_path = joinpath(my_path, file_name*".png")
    #     run(`convert $img_path $out_path`)
    # end

    # # all pngs
    # image_paths = []
    my_path = joinpath(base_path, folder)
    # @info my_path
    # for file in filter(x -> occursin(".png", x), readdir(my_path))
    #     push!(image_paths, joinpath(my_path, file))
    # end

    # # delete png base path
    # deleting = "/Users/MrTrololord/Downloads/"
    # image_paths_relative = map(x -> replace(x, deleting => ""), target_image_paths)

    # get target images --> img1
    target_img_paths = []
    for file in filter(x -> occursin("img1.png", x), readdir(my_path))
        push!(target_img_paths, joinpath(my_path, file))
    end

    target_img_paths

    # paths of target_source_transform triads
    target_source_transform_paths = []
    for target in filter(x -> occursin("img1.png", x), readdir(my_path))
        target_path = joinpath(my_path, target)
        for source in filter(x -> !occursin("img1.png", x) && occursin("png", x), readdir(my_path))
            source_path = joinpath(my_path, source)
            img_number = parse(Int, source[end-4])
            tranform_file = filter(x -> !occursin("img", x) && occursin(string(img_number), x), readdir(my_path))[1]
            tranform_path = joinpath(my_path, tranform_file)
            push!(target_source_transform_paths, [target_path, source_path, tranform_path])
        end
    end
    target_source_transform_paths = permutedims(hcat(target_source_transform_paths...))

    using DataFrames, CSV
    COUNT = 500

    # filter target landmarks not it source img
    # create target_landmarks
    target_img = load(joinpath(base_path, folder, "img1.png"))
    target_lnd_path = joinpath(base_path, folder, "target.csv")
    target_inds = LAP_julia.gen_rand_points(target_img, COUNT, "Gridded")
    for k in 2:6
        # get flow
        transform_path = joinpath(base_path, folder, "H1to$(k)p")
        cur_H = LAP_julia.load_H(transform_path)
        cur_flow = LAP_julia.make_flow_from_H(cur_H, size(target_img))
        # shift and filter landmarks
        target_inds = LAP_julia.get_valid_landmarks(cur_flow.*(-1), target_inds)
    end
    @info size(target_inds)
    @info size(target_img)
    LAP_julia.save_landmarks(target_inds, target_lnd_path)

    # warp target landmarks to source. and create
    for (k, (target_path, source_path, transform_path)) in enumerate(eachrow(target_source_transform_paths))
        @info (k, target_path, source_path, transform_path)
        target_lnd_path = joinpath(base_path, folder, "target.csv")
        cur_H = LAP_julia.load_H(transform_path)
        target_img = load(target_path)
        cur_flow = LAP_julia.make_flow_from_H(cur_H, size(target_img))
        img_number = transform_path[end-1]
        save_path = joinpath(base_path, folder, "source$img_number.csv")
        @info save_path size(cur_flow)
        LAP_julia.save_shift_landmarks(target_lnd_path, cur_flow.*(-1), save_path)
    end

    # create loc table
    tmp_df = CSV.read("/Volumes/davidkunz/Documents/bakalarka/birl_small/dataset_small_repaired.csv") |> DataFrame

    df = deepcopy(tmp_df)
    delete!(df, 1:size(tmp_df,1))

    # rename!(df, Symbol.(replace.(string.(names(df)), Ref(" "=> "_"))))
    # rename!(df, Symbol.(names(df)))

    eltypes(df)
    abc = [:a , :b, :c, :d, :e, :f, :g, :h, :i, :j, :k]
    new_df = DataFrame(a = eltypes(df)[1][],
                       b = eltypes(df)[2][],
                       c = eltypes(df)[3][],
                       d = eltypes(df)[4][],
                       e = eltypes(df)[5][],
                       f = eltypes(df)[6][],
                       g = eltypes(df)[7][],
                       h = eltypes(df)[8][],
                       i = eltypes(df)[9][],
                       j = eltypes(df)[10][],
                       k = eltypes(df)[11][])

    i = 0
    for k in 2:6
        i = i + 1
        row = Dict{String}{Any}()
        row["Column1"] = i
        # images
        row["Target image"] = joinpath("..", folder, "img1.png")
        row["Source image"] = joinpath("..", folder, "img$k.png")
        # lnds
        row["Target landmarks"] = joinpath("..", folder, "target.csv")
        row["Source landmarks"] = joinpath("..", folder, "source$k.csv")
        # sizes
        source_img = load(joinpath(base_path, folder, "img$k.png"))
        row["Image size [pixels]"] = string(size(source_img))
        row["Image diagonal [pixels]"] = round(sqrt(sum(size(source_img).^2)), digits=1)
        row["status"] = "training"
        # empty values
        row["Warped target landmarks"] = missing
        row["Warped source landmarks"] = missing
        row["Execution time [minutes]"] = missing
        global row

        @info rename_dict_keys(names(tmp_df), abc, row)
        push!(new_df, rename_dict_keys(names(tmp_df), abc, row))
    end

    # rename back
    renamed = rename!(new_df, names(tmp_df))

    CSV.write(joinpath(base_path, folder, "loc_table.csv"), renamed)
end

function rename_dict_keys(from_name, to_name, dict)
    new_dict = Dict()
    rename_dict = Dict(from_name .=> to_name)
    for (key, value) in rename_dict
        new_dict[value] = dict[key]
    end
    return new_dict
end