text:
struct DICOMData
    meta::Dict{Tuple{UInt16,UInt16},Any}
    endian::Symbol
    isexplicit::Bool
    vr::Dict{Tuple{UInt16,UInt16},String}
end

function Base.getproperty(dcm::DICOMData, sym::Symbol)
    if sym ∈ fieldnames(DICOMData)
        return getfield(dcm, sym)
    else
        return lookup(dcm, sym)
    end
end

Base.setproperty!(dcm::DICOMData, sym::Symbol, val) = setindex!(dcm, val, sym)

Base.setindex!(dcm::DICOMData, val, sym::Symbol) =
    setindex!(dcm.meta, val, fieldname_dict[sym])
Base.setindex!(dcm::DICOMData, val, str::String) =
    setindex!(dcm.meta, val, fieldname_dict[Symbol(str)])
Base.setindex!(dcm::DICOMData, val, tag::Tuple{UInt16,UInt16}) =
    setindex!(dcm.meta, val, tag)

Base.getindex(dcm::DICOMData, sym::Symbol) = lookup(dcm, sym)
Base.getindex(dcm::DICOMData, str::String) = lookup(dcm, str)
Base.getindex(dcm::DICOMData, tag::Tuple{UInt16,UInt16}) = dcm.meta[tag]

Base.get(dcm::DICOMData, key::Tuple{UInt16,UInt16}, default) = get(dcm.meta, key, default)

Base.keys(dcm::DICOMData) = keys(dcm.meta)
Base.haskey(dcm::DICOMData, sym::Symbol) = haskey(dcm.meta, fieldname_dict[sym])
Base.haskey(dcm::DICOMData, str::String) = haskey(dcm.meta, fieldname_dict[Symbol(str)])
Base.haskey(dcm::DICOMData, tag::Tuple{UInt16,UInt16}) = haskey(dcm.meta, tag)

function lookup(dcm::DICOMData, str::String)
    sym = Symbol(filter(x -> !isspace(x), str))
    return lookup(dcm, sym)
end

function lookup(dcm::DICOMData, sym::Symbol)
    return get(dcm.meta, fieldname_dict[sym], nothing)
end


text:
#
# Convex optimization solver
#
# author: Atsushi Sakai
#

using Test
const __MAIN__ = length(PROGRAM_FILE)!=0 && occursin(PROGRAM_FILE, @__FILE__)

export solve_quadratic_programming
export solve_linear_programming_with_simplex_method

function solve_linear_programming_with_simplex_method(c;
	   	G=nothing, h=nothing, A=nothing, b=nothing)
	"""

	solve linear programming with simplex method
          x = argmin(c.T*x) s.t Gx <= h, Ax = b x>=0

	"""
	nx = length(c)
	n, m = size(G)

    if A == nothing || b == nothing
	    a = vcat(hcat(G, Diagonal{Float64}(I,n), h),
			  hcat(-c', zeros(1, n+1)))
    else
        GA = vcat(G,A,-A)
        hb = vcat(h,b,-b)
        n += 2*size(A)[1]
	    a = vcat(hcat(GA, Diagonal{Float64}(I,n), hb),
			     hcat(-c', zeros(1, n+1)))
        display(a)
    end

	nr, nc = size(a)
	x=y=1

	while true
        # select row
		min = typemax(Float64)
		for k = 1:nc
			if (a[nr, k] < min)
				min = a[nr, k]
				y = k
			end
		end
		if (min >= 0) break; end

		# select col
		min = typemax(Float64)
		for k = 1:nr
			p = a[k, nc] / a[k, y]
			if (a[k, y] > 0 && p < min)
				min = p
				x = k
			end
		end

		p = a[x, y]

		for k = 1:nc
			a[x, k] /= p
		end

		for k = 1:nr
			if k != x
				d = a[k, y]
				for j = 1:nc
					a[k, j] -= d * a[x, j]
				end
			end
		end
	end

	x_hat = zeros(nx)
	for i in 1:nx
		for j in 1:nr
			if a[j, i] == 1.0
				x_hat[i] = a[j, end]
			end
		end
	end
	obj = a[nr, nc]
	# println(obj)
    println(x_hat)

	return x_hat, Dict([ ("objectives", obj)])

end

function solve_quadratic_programming(P, q)
	"""
	solve quadratic programming with only equality constraints
          x = argmin(0.5*x*P*x + q.T*x)
	"""

	x_hat = P \ -q'

	return x_hat
end

function solve_quadratic_programming(P, q, A, b)
	"""
	solve quadratic programming with only equality constraints
          x = argmin(0.5*x*P*x + q.T*x)
          s.t Ax = b
	"""

	n, m = size(A)

	K = [P A';
		A zeros(n, n)]
	# display(K)
	d = [-q';
		 b]
	# display(d)
	x_hat = K \ d

	return x_hat[1:2,:]
end

ext: jl


text:
using Plots
using LinearAlgebra
using Polynomials

using SparseArrays
using Arpack

using ForwardDiff

using BenchmarkTools

function tria6(l1,l2,l3,xs,ys)
    L1 = l1; L2=l2; L3 = l3;
    N1 = (2*L1-1)*L1; N2 = (2*L2-1)*L2; N3 = (2*L3-1)*L3;
    N4 = 4*L1*L2; N5 = 4*L2*L3; N6 = 4*L3*L1;
    x = dot(xs',[N1 N2 N3 N4 N5 N6])
    y = dot(ys',[N1 N2 N3 N4 N5 N6])
    return x,y
end
xs = 10*[-1 0.3 0 -.5 0.2 -.5]
ys = 10*[0 0 1 0.2 .5 .8]
scatter()
for l1 = 0:0.02:1
    for l2 = 0:0.02:1-l1
        scatter!(tria6(l1,l2,1-l1-l2,xs,ys))
    end
end
scatter!()

function tria6jacobian(l1,l2,l3,xs,ys)
    #calculates jacobian matrix dN1/dl1 ....
    L1 = l1; L2=l2; L3=l3;
    N1 = (2*L1-1)*L1; N2 = (2*L2-1)*L2; N3 = (2*L3-1)*L3;
    N4 = 4*L1*L2
    N5 = 4*L2*L3
    N6 = 4*L3*L1
    xy =[xs;ys];
    Ndot = [4*l1-1 0 0;
            0 4*l2-1 0;
            0 0 4*l3-1;
            4*l2 4*l1 0;
            0 4*l2 4*l3;
            4*l3 0 4*l1]
    jacobian = xy*Ndot;
    return jacobian
end
YM = 230e9; #YOUNGS MODULUS
ν = 0.33 # POISON RATIO
d11 = YM*(1-ν)/(1+ν)/(1-2*ν); d12 = ν*d11/(1-ν); d33 = YM/(2*(1+ν));
E = [d11 d12 0;
     d12 d11 0;
     0 0 d33]
function calcBEB(l1,l2,l3,xs,ys)
    L1 = l1; L2= l2; L3 = l3;
    jacobian = tria6jacobian(l1,l2,l3,xs,ys)
    JP = [1 1 1; jacobian] #according to p368 of Fellipa C. Colorado isop
    P = JP\[0 0; 1 0; 0 1]
    Ndot = [4*l1-1 0 0;
            0 4*l2-1 0;
            0 0 4*l3-1;
            4*l2 4*l1 0;
            0 4*l2 4*l3;
            4*l3 0 4*l1]
    Nderivxy = P'*Ndot'; # is [dN1/dx ... dN6/dx;
                        #     dN1/dy ... dN6/dy]
    N1 = (2*L1-1)*L1; N2 = (2*L2-1)*L2; N3 = (2*L3-1)*L3;
    N4 = 4*L1*L2
    N5 = 4*L2*L3
    N6 = 4*L3*L1
    N = [N1 0 N2 0 N3 0 N4 0 N5 0 N6 0;
         0 N1 0 N2 0 N3 0 N4 0 N5 0 N6]
    N1x = Nderivxy[1,1]; N2x = Nderivxy[1,2]; N3x = Nderivxy[1,3]; N4x = Nderivxy[1,4]; N5x = Nderivxy[1,5]; N6x = Nderivxy[1,6];
    N1y = Nderivxy[2,1]; N2y = Nderivxy[2,2]; N3y = Nderivxy[2,3]; N4y = Nderivxy[2,4]; N5y = Nderivxy[2,5]; N6y = Nderivxy[2,6];

    B = [N1x 0 N2x 0 N3x 0 N4x 0 N5x 0 N6x 0;
         0  N1y 0 N2y 0 N3y 0 N4y 0 N5y 0 N6y;
         N1y N1x N2y N2x N3y N3x N4y N4x N5y N5x N6y N6x]
    BEB = B'*E*B
    return BEB
end

function tria6Ke(xs,ys)
    Ke  = zeros(12,12)
    g1 = [0.108103, 0.445948, 0.445948]; w1 =0.223382;
    g2 = [0.445948, 0.108103, 0.445948]; w2 =0.223382;
    g3 = [0.445948, 0.445948, 0.108103]; w3 =0.223382;
    g4 = [0.816848, 0.0915762, 0.0915762]; w4 = 0.109952;
    g5 = [0.0915762, 0.816848, 0.0915762]; w5 = 0.109952;
    g6 = [0.0915762, 0.0915762, 0.816848]; w6 = 0.109952;
    g = [[g1] [g2] [g3] [g4] [g5] [g6]]; w = [w1 w2 w3 w4 w5 w6];

    for i=1:6
        Ke = Ke + w[i]*calcBEB(g[i][1],g[i][2],g[i][3],xs,ys)
    end
    return Ke
end

#f = open("C:\\User\\billykon\\Desktop\\Mesh_1.unv","r")

#function edgeToLs(xedge,yedge,TOPOLOGY):

max_stars_repo_name: billyon/first

max_issues_repo_name: billyon/first

max_forks_repo_name: billyon/first


text:
using Permutations
using DelimitedFiles
using Plots
using DeferredAcceptance
using StatsBase

students = readdlm("examples/sysopt/students.dat", Int)
schools = readdlm("examples/sysopt/schools.dat", Int)
n, m = size(schools)
capacities = ones(Int64, m)

descr = "Hybrid market with $n students, $m schools"
println(descr)

schools_STB = singletiebreaking(schools)
schools_MTB = multipletiebreaking(schools)
schools_XTB = hybridtiebreaking(schools, 0.5)
schools_WXTB = welfaretiebreaking(students, schools, 0.5)

# Update rank risks
println("Starting STB")
cdf_STB = DA_rank_dist(students, schools_STB, capacities)
println("Starting MTB")
cdf_MTB = DA_rank_dist(students, schools_MTB, capacities)
println("Starting XTB")
cdf_XTB = DA_rank_dist(students, schools_XTB, capacities)
println("Starting WXTB")
cdf_WXTB = DA_rank_dist(students, schools_WXTB, capacities)

TTC_results = TTC_match(students, capacities)[2]
TTC_cmap = countmap(TTC_results)
cdf_TTC = cumsum([get(TTC_cmap, i, 0) for i in 1:m])

LPso_results = vec([students[c, s] for (s, c) in enumerate(readdlm("examples/sysopt/LPSystemOpt.dat", Int))])
LPso_cmap = countmap(LPso_results)
cdf_LPso = cumsum([get(LPso_cmap, i, 0) for i in 1:m])

IPso_results = vec([students[c, s] for (s, c) in enumerate(readdlm("examples/sysopt/IPStableOpt.dat", Int))])
IPso_cmap = countmap(IPso_results)
cdf_IPso = cumsum([get(IPso_cmap, i, 0) for i in 1:m])

p = plot([cdf_STB, cdf_MTB, cdf_XTB, cdf_WXTB, cdf_TTC, cdf_LPso, cdf_IPso],
     label = ["DA-STB" "DA-MTB" "DA-XTB, λ=0.5" "DA-WXTB, λ=0.5" "RSD to top trading cycles" "System welfare opt (TUM LP)" "Stable welfare opt (IP)"],
     lc = [:teal :rebeccapurple :firebrick :olivedrab :orangered :gold :black],
        ls = [:dash :dot :dashdot :dash :dot :solid :dashdot],
        legend = :bottomright,
        title = descr, titlefontsize=11,
        xlabel = "rank", ylabel= "number of students")

display(p)

# savefig(p, string("examples/sysopt/plot.pdf"))
# savefig(p, string("examples/sysopt/plot.png"))

max_stars_repo_name: maxkapur/DeferredAcceptance

max_issues_repo_name: maxkapur/DeferredAcceptance

max_forks_repo_name: maxkapur/DeferredAcceptance


text:
#=
MMFF94Parameters:
- Julia version:
- Author: Dan
- Date: 2021-07-30
=#
"""
Radii derived from [MMFF94 Validation Suite](http://www.ccl.net/cca/data/MMFF94/).
"""
const MMFF94radii = Float64[
     0.33, 0.0,
	 1.34, 0.90, 0.81, 0.77, 0.73, 0.72, 0.74, 0.0,
	 1.54, 1.30, 1.22, 1.15, 1.09, 1.03, 1.01, 0.0,
	 1.96, 1.74,
	 1.44, 1.36, 0.00, 0.00, 0.00,
	 0.00, 0.00, 0.00, 1.38, 1.31,
	 1.19, 1.20, 1.20, 1.16, 1.15, 0.0,
	 2.11, 1.92,
	 1.62, 1.48, 0.00, 0.00, 0.00,
	 0.00, 0.00, 0.00, 1.53, 1.48,
	 1.46, 1.40, 1.41, 1.35, 1.33, 0.0
]

"""
Stretch parameters derived from [MMFF94 Validation Suite](http://www.ccl.net/cca/data/MMFF94/).
"""
const MMFF94StrechParameters = Float64[
     2.20, 0.0,
	 0.97, 1.47, 2.01, 2.5, 3.07, 3.5, 4.10, 0.0,
	 1.01, 1.23, 1.47, 1.74, 2.06, 2.44, 2.83, 0.0,
	 0.91, 1.04,
	 1.3, 1.5, 1.6, 1.6, 1.5,
	 1.8, 1.8, 1.8, 1.9, 1.6,
	 1.82, 2.02, 2.20, 2.48, 2.74, 0.0,
	 0.89, 0.99,
	 1.3, 1.4, 1.6, 1.8, 1.9,
	 2.2, 2.2, 2.2, 1.9, 1.7,
	 1.49, 1.72, 1.82, 2.01, 2.21, 0.0
]


text:
"""
    GridWorldBase{O} <: AbstractArray{Bool, 3}

A basic representation of grid world.
The first dimension uses multi-hot encoding to encode objects in a tile.
The second and third dimensions correspond to the height and width of the grid respectively.
"""
struct GridWorldBase{O} <: AbstractArray{Bool, 3}
    grid::BitArray{3}
    objects::O
end

function GridWorldBase(objects::Tuple{Vararg{AbstractObject}}, height::Int, width::Int)
    grid = falses(length(objects), height, width)
    GridWorldBase(grid, objects)
end

get_grid(world::GridWorldBase) = world.grid
get_objects(world::GridWorldBase) = world.objects

get_num_objects(grid::AbstractArray{Bool, 3}) = size(grid, 1)
get_height(grid::AbstractArray{Bool, 3}) = size(grid, 2)
get_width(grid::AbstractArray{Bool, 3}) = size(grid, 3)

#####
# Indexing of GridWorldBase objects
#####

@forward GridWorldBase.grid Base.size, Base.getindex, Base.setindex!

@generated function Base.to_index(::GridWorldBase{O}, object::X) where {X<:AbstractObject, O}
    i = findfirst(X .=== O.parameters)
    isnothing(i) && error("unknow object $object")
    :($i)
end

Base.getindex(world::GridWorldBase, object::AbstractObject, args...) = getindex(get_grid(world), Base.to_index(world, object), args...)
Base.setindex!(world::GridWorldBase, value::Bool, object::AbstractObject, args...) = setindex!(get_grid(world), value, Base.to_index(world, object), args...)

#####
# methods for agent view
#####

function get_grid_inds((i, j), (m, n))
    temp1 = m ÷ 2
    temp2 = i - temp1
    temp3 = n ÷ 2
    temp4 = j - temp3
    return CartesianIndices((temp2 : temp2 + m - 1, temp4 : temp4 + n - 1))
end

function get_grid_inds((i, j), (m, n), ::Up)
    temp1 = n - 1
    temp2 = temp1 ÷ 2
    temp3 = j - temp2
    return CartesianIndices((i - m + 1 : i, temp3 : temp3 + temp1))
end

function get_grid_inds((i, j), (m, n), ::Down)
    temp1 = n - 1
    temp2 = temp1 ÷ 2
    temp3 = j - temp2
    return CartesianIndices((i : i + m - 1, temp3 : temp3 + temp1))
end

function get_grid_inds((i, j), (m, n), ::Left)
    temp1 = n - 1
    temp2 = temp1 ÷ 2
    temp3 = i - temp2
    return CartesianIndices((temp3 : temp3 + temp1, j - m + 1 : j))
end

function get_grid_inds((i, j), (m, n), ::Right)
    temp1 = n - 1
    temp2 = temp1 ÷ 2
    temp3 = i - temp2
    return CartesianIndices((temp3 : temp3 + temp1, j : j + m - 1))
end

map_ind((i,j), (m, n), ::Up) = (m-i+1, n-j+1)
map_ind((i,j), (m, n), ::Down) = (i,j)
map_ind((i,j), (m, n), ::Left) = (m-j+1, i)
map_ind((i,j), (m, n), ::Right) = (j, n-i+1)

function get_grid(src_grid::AbstractArray{Bool, 3}, pos::CartesianIndex{2}, dir::AbstractDirection, window_size, layers)
    dest_grid = falses(length(layers), window_size...)
    get_grid!(dest_grid, src_grid, pos, dir, window_size, layers)
    return dest_grid
end

function get_grid!(dest_grid::AbstractArray{Bool, 3}, src_grid::AbstractArray{Bool, 3}, pos::CartesianIndex{2}, dir::AbstractDirection, window_size, layers)
    window_inds = get_grid_inds(pos.I, window_size, dir)

    valid_inds = CartesianIndices((get_height(src_grid), get_width(src_grid)))

    @views for ind in CartesianIndices(window_inds)
        pos = window_inds[ind]
        if pos ∈ valid_inds
            dest_grid[:, map_ind(ind.I, window_size, dir)...] .= src_grid[layers, pos]
        end
    end

    return nothing
end

function get_grid(src_grid::AbstractArray{Bool, 3}, pos::CartesianIndex{2}, window_size, layers)
    dest_grid = falses(length(layers), window_size...)
    get_grid!(dest_grid, src_grid, pos, window_size, layers)
    return dest_grid
end

function get_grid!(dest_grid::AbstractArray{Bool, 3}, src_grid::AbstractArray{Bool, 3}, pos::CartesianIndex{2}, window_size, layers)
    window_inds = get_grid_inds(pos.I, window_size)

    valid_inds = CartesianIndices((get_height(src_grid), get_width(src_grid)))

    @views for ind in CartesianIndices(window_inds)
        pos = window_inds[ind]
        if pos ∈ valid_inds
            dest_grid[:, ind.I...] .= src_grid[layers, pos]
        end
    end

    return nothing
end

#####
# utils
#####

function Random.rand(rng::Random.AbstractRNG, f::Function, inds::Union{Vector{CartesianIndex{2}}, CartesianIndices{2}}; max_try::Int = 1000)
    for _ in 1:max_try
        pos = rand(rng, inds)
        if f(pos)
            return pos
        end
    end
    @warn "number of tries exceeded max_try = $max_try"
    return inds[1]
end

function Random.rand(rng::Random.AbstractRNG, f::Function, world::GridWorldBase; max_try = 1000)
    inds = CartesianIndices((get_height(world), get_width(world)))
    rand(rng, f, inds, max_try = max_try)
end


text:
# functions for interpolating between different SBP operators

#TODO: interpolate q or q_vec
"""
  This function interpolates a solution field from one SBP operator to
  another.  This requires the meshes be identical except for the difference
  in SBP operator.  In particular, the element numbering scheme must be the
  same for this to work.

  The field in question can be any field defined on the nodes of an SBP
  operator.  The interpolation is generally accurate to the degree of the
  SBP operator.

  **Inputs**

   * sbp_old: the SBP operator for the mesh with the solution to be interpolated
   * q_old: the 3D array containing the the solution to be interpolated
   * sbp_new: the SBP operator for the new mesh

  **Inputs/Outputs**

   * q_new: the array to put the new solution into (overwritten)
"""
function interpField(sbp_old::AbstractOperator, q_old::Abstract3DArray,
                        sbp_new::AbstractOperator, q_new::Abstract3DArray)

  @assert size(q_old, 3) == size(q_new, 3)
  @assert size(q_old, 1) == size(q_new, 1)

  # construct interpolation operator
  node_coords = calcnodes(sbp_new)
  interp_op = SummationByParts.buildinterpolation(sbp_old, node_coords)

  # apply interpolation operator
  applyInterpolation(interp_op, q_old, q_new)

  return nothing
end

"""
  This function interpolates a vector from one mesh to another.  The two
  meshes must have the elements numbered in the same order.

  Currently, the vector must have `mesh.numDofPerNode` degrees of freedom
  per node.

  **Inputs**

   * mesh_old: the old mesh
   * sbp_old: the old SBP operator
   * q_old: the vector containing the data, length `mesh_old.numDof`
   * mesh_new: the new mesh
   * sbp_new: the new SBP operator

  **Inputs/Outputs**

   * q_new: vector to overwrite with the interpolated values
"""
function interpField(mesh_old::AbstractMesh, sbp_old::AbstractOperator,
                     q_old::AbstractVector, mesh_new::AbstractMesh,
                     sbp_new::AbstractOperator, q_new::AbstractVector)

  @assert length(q_old) == mesh_old.numDof
  @assert length(q_new) == mesh_new.numDof
  # construct interpolation operator
  node_coords = calcnodes(sbp_new)
  interp_op = SummationByParts.buildinterpolation(sbp_old, node_coords)

  applyInterpolation(interp_op, mesh_old, q_old, mesh_new, q_new)

  return nothing
end

function applyInterpolation(interp_op::AbstractMatrix, q_old::Abstract3DArray,
                            q_new::Abstract3DArray)

  @assert size(q_old, 3) == size(q_new, 3)
  @assert size(q_old, 1) == size(q_new, 1)

  numEl = size(q_old, 3)

  for i=1:numEl
    qold_i = sview(q_old, :, :, i)
    qnew_i = sview(q_new, :, :, i)

    # to interpolate all dofs at the same time (in the layout of the q arrays)
    # do q_old[:, :, i]*interp_op.' = q_new[:, :, i]
    smallmatmatT!(qold_i, interp_op, qnew_i)
  end

  return nothing
end

function applyInterpolation(interp_op::AbstractMatrix,
                        mesh_old::AbstractMesh, q_old::AbstractVector{Told},
                        mesh_new::AbstractMesh, q_new::AbstractVector{Tnew}
                        ) where {Told, Tnew}

  q_el_old = zeros(Told, mesh_old.numDofPerNode, mesh_old.numNodesPerElement)
  q_el_new = zeros(Tnew, mesh_new.numDofPerNode, mesh_new.numNodesPerElement)
  for i=1:mesh_old.numEl

    @simd for j=1:mesh_old.numNodesPerElement
      @simd for k=1:mesh_old.numDofPerNode
        q_el_old[k, j] = q_old[mesh_old.dofs[k, j, i]]
      end
    end

    smallmatmatT!(q_el_old, interp_op, q_el_new)

    @simd for j=1:mesh_new.numNodesPerElement
      @simd for k=1:mesh_new.numDofPerNode
        q_new[mesh_new.dofs[k, j, i]] = q_el_new[k, j]
      end
    end

  end  # end i

  return nothing
end


text:
using Router

route("/") do
  Router.serve_static_file("/welcome.html")
end

max_stars_repo_name: anirudhmurali/hello-julia-genie

max_issues_repo_name: anirudhmurali/hello-julia-genie

max_forks_repo_name: anirudhmurali/hello-julia-genie


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://discourse.julialang.org/t/empty-container-reduce-error-diff-eq/64956

# %%
using DifferentialEquations
using LinearAlgebra
using Plots

# %%
function plot_2ndorder(sol; l1=:bottomright, l2=:bottomright, size=(720, 300), kwargs...)
    t = range(sol.prob.tspan...; length=400)
    n = Base.size(sol, 1) ÷ 2
    v = vcat((t -> sol(t)[1:n]').(t)...)
    d = vcat((t -> sol(t)[n+1:end]').(t)...)
    P = plot(t, d; label=permutedims(["\$d_{$i}\$" for i in 1:n]), legend=l1)
    Q = plot(t, v; label=permutedims(["\$v_{$i}\$" for i in 1:n]), legend=l2)
    plot(P, Q; size, kwargs...)
end

# %%
function g_naive!(dv, v, d, p, t)
    m, f, c, K = p
    # The following is equivalent to dv .= M \ (f - (C * v + K * d)),
    # where K = Matrix, C = Diagonal(c), M = Diagonal(m)
    dv .= m .\ (f .- (c .* v .+ K * d))
    return
end

# %%
m = Float64[1, 2, 3, 4]
f = Float64[-10, -10, -10, -10]
c = Float64[1, 1, 1, 1]
K = Float64[
     2 -1  0  0
    -1  2 -1  0
     0  1  2 -1
     0  0 -1  2
]
p = (m, f, c, K)
v0 = Float64[0, 0, 0, 0]
d0 = Float64[3, 1, -1, -3]
tspan = (0.0, 10.0)

# %%
prob = SecondOrderODEProblem(g_naive!, v0, d0, tspan, p)
sol = solve(prob)
plot_2ndorder(sol)

# %%
function g!(dv, v, d, p, t)
    m, f, c, K = p
    # The following is equivalent to dv .= M \ (f - (C * v + K * d)),
    # where K = Matrix, C = Diagonal(c), M = Diagonal(m)
    mul!(dv, K, d)
    @. dv = m \ (f - (c * v + dv))
    return
end

# %%
prob = SecondOrderODEProblem(g!, v0, d0, tspan, p)
sol = solve(prob)
plot_2ndorder(sol)

# %%
using BenchmarkTools

@btime g_naive!($(similar(v0)), $v0, $d0, $p, 0.0)
@btime g!($(similar(v0)), $v0, $d0, $p, 0.0)

# %%

max_stars_repo_path: 0012/SecondOrderODEProblem example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0012/SecondOrderODEProblem example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0012/SecondOrderODEProblem example.jl

max_forks_repo_name: genkuroki/public


text:
function palindrome_permutation(str :: ASCIIString)
  """
  Convert a character to a number according
  to the corresponding shift from 'a'
  """
  function char_to_int(ch :: Char)
    a :: Int = int('a')
    ch :: Int = int(ch)
    if a <= ch <= int('z')
      return ch - a
    else
      return -1
    end
  end

  """
  Toggle the bit at a specific index in a bit vector
  """
  function toggle_bit(bit_vector :: Int, idx :: Int)
    if idx <= 0
      return bit_vector
    end
    return bit_vector $= 1 << idx
  end

  """
  Create a bit vecotor for a given string
  """
  function create_bit_vector(str :: ASCIIString)
    bit_vector :: Int = 0
    for ch in str
      bit_vector = toggle_bit(bit_vector, char_to_int(ch))
    end
    return bit_vector
  end

  """
  Check whether a bit vecotr has only one bit
  """
  function has_one_bit(bit_vector :: Int)
    ((bit_vector - 1) & bit_vector) == 0
  end

  bv :: Int  = create_bit_vector(lowercase(str))

  return (bv == 0) || has_one_bit(bv) # satisfies either all even or only one odd
end

# Time  : O(n)
# Space : O(1)

using Base: Test

@test palindrome_permutation("Tact Coa") == true
@test palindrome_permutation("Tact Coad") == false

max_stars_repo_name: Kiandr/crackingcodinginterview

max_issues_repo_name: Kiandr/crackingcodinginterview

max_forks_repo_name: Kiandr/crackingcodinginterview


text:
import .StatsFuns

for lfn in propertynames(StatsFuns)
  m = match(r"^(.*)log(pdf|cdf|ccdf)", string(lfn))
  if m !== nothing
    fn = Symbol(m[1], m[2])
    @eval StatsFuns.$(fn)(::Type{E}, args...; opts...) where {E<:AnyLog} = exp(E, StatsFuns.$(lfn)(args...; opts...))
  end
end


text:
function create_instances_table(db)
    createtable_query = """
    CREATE TABLE IF NOT EXISTS instances(
        id INTEGER NOT NULL,
        name TEXT NOT NULL,
        scenario INTEGER NOT NULL,
        source_path TEXT NOT NULL,
        source_type TEXT NOT NULL,
        pfn_path TEXT,
        graph_path TEXT,
        mat_path TEXT,
        ctr_path TEXT,
        matpowerm_path TEXT,
        date TEXT NOT NULL,

        nbus INTEGER, nbranch INTEGER, nbranch_unique INTEGER, ngen INTEGER,

        nb_edge INTEGER,
        nb_vertex INTEGER,
        degree_max INTEGER, degree_min INTEGER, degree_mean REAL, degree_median REAL, degree_var REAL,
        global_clustering_coefficient REAL,
        density REAL,
        diameter INTEGER,
        radius INTEGER,

        PD_max REAL, PD_min REAL, PD_mean REAL, PD_median REAL, PD_var REAL,
        QD_max REAL, QD_min REAL, QD_mean REAL, QD_median REAL, QD_var REAL,
        GS_max REAL, GS_min REAL, GS_mean REAL, GS_median REAL, GS_var REAL,
        BS_max REAL, BS_min REAL, BS_mean REAL, BS_median REAL, BS_var REAL,
        VM_max REAL, VM_min REAL, VM_mean REAL, VM_median REAL, VM_var REAL,
        VA_max REAL, VA_min REAL, VA_mean REAL, VA_median REAL, VA_var REAL,
        VMAX_max REAL, VMAX_min REAL, VMAX_mean REAL, VMAX_median REAL, VMAX_var REAL,
        VMIN_max REAL, VMIN_min REAL, VMIN_mean REAL, VMIN_median REAL, VMIN_var REAL,

        PG_max REAL, PG_min REAL, PG_mean REAL, PG_median REAL, PG_var REAL,
        QG_max REAL, QG_min REAL, QG_mean REAL, QG_median REAL, QG_var REAL,
        QMAX_max REAL, QMAX_min REAL, QMAX_mean REAL, QMAX_median REAL, QMAX_var REAL,
        QMIN_max REAL, QMIN_min REAL, QMIN_mean REAL, QMIN_median REAL, QMIN_var REAL,
        PMAX_max REAL, PMAX_min REAL, PMAX_mean REAL, PMAX_median REAL, PMAX_var REAL,
        PMIN_max REAL, PMIN_min REAL, PMIN_mean REAL, PMIN_median REAL, PMIN_var REAL,

        BR_R_max REAL, BR_R_min REAL, BR_R_mean REAL, BR_R_median REAL, BR_R_var REAL,
        BR_X_max REAL, BR_X_min REAL, BR_X_mean REAL, BR_X_median REAL, BR_X_var REAL,
        BR_B_max REAL, BR_B_min REAL, BR_B_mean REAL, BR_B_median REAL, BR_B_var REAL,

        UNIQUE(name, scenario),
        PRIMARY KEY(id)
    )
    """
    SQLite.execute(db, createtable_query)
end

function create_decompositions_table(db)
    createtable_query = """
    CREATE TABLE IF NOT EXISTS decompositions(
        id INTEGER NOT NULL,
        uuid TEXT NOT NULL,
        origin_id INTEGER NOT NULL,
        origin_name TEXT NOT NULL,
        origin_scenario INTEGER NOT NULL,

        extension_alg TEXT,
        preprocess_path TEXT,
        date TEXT,

        clique_path TEXT,
        cliquetree_path TEXT,
        graph_path TEXT,

        nb_added_edge_dec INTEGER,
        nb_edge INTEGER,
        nb_vertex INTEGER,
        degree_max INTEGER, degree_min INTEGER, degree_mean REAL, degree_median REAL, degree_var REAL,
        global_clustering_coefficient REAL,
        density REAL,
        diameter INTEGER,
        radius INTEGER,

        nb_clique INTEGER,
        clique_size_max INTEGER, clique_size_min INTEGER, clique_size_mean REAL, clique_size_median REAL, clique_size_var REAL,
        nb_lc INTEGER,

        date_solving TEXT,
        solve_log_path TEXT,
        objective REAL, nb_iter INTEGER, m REAL,
        solving_time REAL,

        PRIMARY KEY(id),
        FOREIGN KEY(origin_name, origin_scenario) REFERENCES instances(name, scenario)
    )
    """
    SQLite.execute(db, createtable_query)
end

function create_merges_table(db)
    createtable_query = """
    CREATE TABLE IF NOT EXISTS merges(
        in_id INTEGER NOT NULL,
        out_id INTEGER NOT NULL,
        heuristics TEXT NOT NULL,
        treshold_name TEXT NOT NULL,
        treshold_percent REAL NOT NULL,

        nb_added_edge INT NOT NULL,

        PRIMARY KEY(in_id, out_id),
        FOREIGN KEY(in_id) REFERENCES decompositions(id),
        FOREIGN KEY(out_id) REFERENCES decompositions(id)
    )
    """
    SQLite.execute(db, createtable_query)
end

function create_combinations_table(db)
    createtable_query = """
    CREATE TABLE IF NOT EXISTS combinations(
        in_id1 INTEGER NOT NULL,
        in_id2 INTEGER NOT NULL,
        out_id INTEGER NOT NULL,
        how TEXT NOT NULL,
        extension_alg TEXT NOT NULL,

        PRIMARY KEY(in_id1, in_id2, out_id),
        FOREIGN KEY(in_id1) REFERENCES decompositions(id),
        FOREIGN KEY(in_id2) REFERENCES decompositions(id),
        FOREIGN KEY(out_id) REFERENCES decompositions(id)
    )
    """
    SQLite.execute(db, createtable_query)
end

function setup_db(name; delete_if_exists=false)
    db = SQLite.DB(name)
    if delete_if_exists
        delete_table(db, table) = DBInterface.execute(db, "DROP TABLE IF EXISTS $table")
        delete_table(db, "instances")
        delete_table(db, "decompositions")
        delete_table(db, "merges")
        delete_table(db, "combinations")
    end
    create_instances_table(db)
    create_decompositions_table(db)
    create_merges_table(db)
    create_combinations_table(db)
    return db
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
function SimpleRatioAmbiguityAcceptance(arg0::jdouble)
    return SimpleRatioAmbiguityAcceptance((jdouble,), arg0)
end

function accept(obj::SimpleRatioAmbiguityAcceptance, arg0::Vector{IntegerLeastSquareSolution})
    return jcall(obj, "accept", IntegerLeastSquareSolution, (Vector{IntegerLeastSquareSolution},), arg0)
end

function number_of_candidates(obj::SimpleRatioAmbiguityAcceptance)
    return jcall(obj, "numberOfCandidates", jint, ())
end


text:
"A package for some basic plane geometry computations."
module PlaneGeometry

export
    # Geometric objects
    GeoObject, Point, Edge, Triangle, Circle,
    vertices, edges,

    # distances
    distance, squaredist, midpoint, concurrent,

    # triangles
    equipoints, outer_equitri, outer_equitriangles, isequilateral, median,

    # circles
    circumcircle,

    # utilities
    does_thmhold, @code_md

include("elementary.jl")
include("distance.jl")
include("triangle.jl")
include("circle.jl")
include("simpy.jl")
include("plots.jl")
include("util.jl")
include("theorem/theorem.jl")

using .Theorems.Napoleon
export
    napoleon_proof, napoleon_draw, napoleon_rand, napoleon_tri

using .Theorems.Centroid
export
    centroid_proof, centroid_draw, centroid_rand, centroid

using .GeoPlots
export
    plot, plot!, shape

using .Theorems.Ceva
export
    ceva_proof, ceva_draw, ceva_rand, ceva, ceva_check

function __init__()
    global A, B, C, O, tri

    A = Point(0, 0)
    B = Point(1, 3)
    C = Point(4, 2)
    O = Point(2, 2)
    tri = Triangle(A, B, C)
end

end # module


text:
export stationary_polynomial, stationary_mean, stationary_variance, stationary_covariance_ellipsoid, max_entropy_measure

"""
	stationary_pop(MP::MarkovProcess, v::APL, d::Int, solver)

returns SOS program of degree `d` for compuitation of a **lower** bound
on the expecation of a polynomial observable ``v(x)`` at steady state of the
Markov process `MP`.
"""
function stationary_pop(MP::MarkovProcess, p::APL, order::Int, solver, P::Partition = trivial_partition(MP.X))
    model = SOSModel(solver)
    @variable(model, s)
    w = Dict(v => (props(P.graph, v)[:cell] isa Singleton ?
                   @variable(model) :
                   @variable(model, [1], Poly(monomials(MP.x, 0:order)))[1]) for v in vertices(P.graph))

    for v in vertices(P.graph)
        add_stationarity_constraints!(model, MP, v, P, props(P.graph, v)[:cell], w, s-p)
    end

    for e in edges(P.graph)
        add_coupling_constraints!(model, MP, e, P, w)
    end
    @objective(model, Max, s)
    return model, w
end

"""
	stationary_polynomial(MP::MarkovProcess, v::APL, d::Int, solver)

returns a **lower** bound on the expecation of a polynomial observables ``v(x)``
at steady state of the Markov process `MP`. The bound is computed based on an
SOS program over a polynomial of degree at most `d`; the bounds can be
tightened by increasing `d`. The program is solved with `solver`.
"""
function stationary_polynomial(MP::MarkovProcess, v::APL, d::Int, solver, P::Partition = trivial_partition(MP.X))
    model, w = stationary_pop(MP, v, d, solver, P)
    optimize!(model)
    return Bound(objective_value(model), model, P, Dict(key => value(w[key]) for key in keys(w)))
end

stationary_polynomial(MP::MarkovProcess, v::Num, d::Int, solver, P::Partition = trivial_partition(MP.X)) = stationary_polynomial(MP, polynomialize_expr(v, MP.poly_vars), d, solver, P)

"""
	stationary_mean(MP::MarkovProcess, v::APL, d::Int, solver)

returns **lower** and **upper** bound on the observable ``v(x)`` at steady state of
the Markov process `MP`. Bounds are computed based on SOS programs over a
polynomial of degree at most `d`; the bounds can be tightened by increasing
`d`. The program is solved with `solver`.
"""
function stationary_mean(MP::MarkovProcess, v::APL, d::Int, solver, P::Partition = trivial_partition(MP.X))
    lb = stationary_polynomial(MP, v, d, solver, P)
    ub = stationary_polynomial(MP, -v, d, solver, P)
	ub.value *= -1
    return lb, ub
end

stationary_mean(RP::ReactionProcess, S, d::Int, solver, P::Partition = trivial_partition(RP.JumpProcess.X)) = stationary_mean(RP.JumpProcess, RP.species_to_state[S], d, solver, P)
stationary_mean(LP::LangevinProcess, S, d::Int, solver, P::Partition = trivial_partition(LP.DiffusionProcess.X)) = stationary_mean(LP.DiffusionProcess, LP.species_to_state[S], d, solver, P)
stationary_mean(MP::MarkovProcess, v::Num, d::Int, solver, P::Partition = trivial_partition(MP.X)) = stationary_mean(MP, polynomialize_expr(v, MP.poly_vars), d, solver, P)

"""
	stationary_mean(rn::ReactionSystem, S0::Dict, S, d::Int, solver,
			scales = Dict(s => 1 for s in species(rn));
			auto_scaling = false)

returns **lower** and **upper** bound on the mean of species `S` of the reaction
network `rn` with initial condition `S0` (for all species!). The bound is based
on an SOS program of order `d` solved via `solver`; the bounds can be tightened
by increasing `d`.

For numerical stability, it is recommended to provide scales of the expected
magnitude of molecular counts for the different species at steady state.
If the system is **closed** it is also possible to enable auto_scaling which will
find the maximum molecular counts for each species under stoichiometry
constraints (via LP).

If the initial condition of the reaction network under investigation is
unknown or irrelevant, simply call

	stationary_mean(rn::ReactionSystem, S, d::Int, solver,
			scales = Dict(s => 1 for s in species(rn))).
"""
function stationary_mean(rn::ReactionSystem, S0::Dict, S, d::Int, solver,
	 					 scales = Dict(s => 1 for s in species(rn));
						 params::Dict = Dict(), auto_scaling = false)
	RP, S0 = reaction_process_setup(rn, S0, scales = scales, auto_scaling = auto_scaling, solver = solver, params = params)
 	return stationary_mean(RP.JumpProcess, RP.species_to_state[S], d, solver)
end

function stationary_mean(rn::ReactionSystem, S, d::Int, solver,
					     scales = Dict(s => 1 for s in species(rn));
						 params::Dict = Dict())
	RP = reaction_process_setup(rn, scales = scales, params = params)
	return stationary_mean(RP.JumpProcess, RP.species_to_state[S], d, solver)
end

"""
	stationary_variance(MP::MarkovProcess, v::APL, d::Int, solver)

returns SOS program of degree `d` for computation of an **upper** bound on the
variance of a polynomial observables `v` at steady state of the Markov process
`MP`.
"""
function stationary_variance(MP::MarkovProcess, p::APL, d::Int, solver, P::Partition = trivial_partition(MP.X))
    model = SOSModel(solver)
	w = Dict(v => (props(P.graph, v)[:cell] isa Singleton ?
                   @variable(model) :
                   @variable(model, [1], Poly(monomials(MP.x, 0:d)))[1]) for v in vertices(P.graph))
    @variable(model, s)
    @variable(model, S[1:2])
    @constraint(model, [1 S[1]; S[1] S[2]] in PSDCone())

	for v in vertices(P.graph)
        add_stationarity_constraints!(model, MP, v, P, props(P.graph, v)[:cell], w, s + p^2 + 2*S[1]*p)
    end

	for e in edges(P.graph)
        add_coupling_constraints!(model, MP, e, P, w)
    end

    @objective(model, Max, s-S[2])
    optimize!(model)
    return Bound(-objective_value(model), model, P, Dict(key => value(w[key]) for key in keys(w)))
end

stationary_variance(RP::ReactionProcess, x, d::Int, solver, P::Partition = trivial_partition(RP.JumpProcess.X)) = stationary_variance(RP.JumpProcess, RP.species_to_state[x], d, solver, P)
stationary_variance(MP::MarkovProcess, x::Num, d::Int, solver, P::Partition = trivial_partition(MP.X)) = stationary_variance(MP, polynomialize_expr(x, MP.poly_vars), d, solver, P)

"""
	stationary_variance(rn::ReactionSystem, S0, x, d::Int, solver,
			    scales = Dict(s => 1 for s in species(rn));
			    auto_scaling = false)

returns **upper** bound on the variance of species `S` of the reaction
network rn with initial condition `S0` (for all species!). The bound is based
on an SOS program of degree `d` solved via `solver`; the bound can be tightened
by increasing `d`.

For numerical stability, it is recommended to provide scales of the expected
magnitude of molecular counts for the different species at steady state. If
the system is **closed** it is also possible to enable `auto_scaling` which will
find the maximum molecular counts for each species under stoichiometry
constraints (via LP).

If the initial condition of the reaction network under investigation is
unknown or irrelevant, simply call

	stationary_variance(rn::ReactionSystem, S, d::Int, solver,
			    scales = Dict(s => 1 for s in species(rn)))
"""
function stationary_variance(rn::ReactionSystem, S0, S, d::Int, solver,
	 						 scales = Dict(s => 1 for s in species(rn));
							 auto_scaling = false, params::Dict = Dict())
	RP, S0 = reaction_process_setup(rn, S0, scales = scales, auto_scaling = auto_scaling, solver = solver, params = params)
 	return stationary_variance(RP.JumpProcess, RP.species_to_state[S], d, solver)
end

function stationary_variance(rn::ReactionSystem, S, d::Int, solver,
							 scales = Dict(s => 1 for s in species(rn));
							 params::Dict = Dict())
	RP = reaction_process_setup(rn, scales = scales, params = params)
	return stationary_variance(RP.JumpProcess, RP.species_to_state[S], d, solver)
end

@doc raw"""
	stationary_covariance_ellipsoid(MP::MarkovProcess, v::Vector{<:APL}, d::Int, solver)

returns an **upper** on the volume of the covariance ellipsoid of a vector of
polynomial observables ``v(x)``, i.e., ``\text{det}(\mathbb{E} [v(x)v(x)^\top] - \mathbb{E}[v(x)] \mathbb{E}[v(x)]^\top)``, at steady state of the
Markov process `MP`.
The bounds are computed via an SOS program of degree `d`, hence can be tightened
by increasing `d`. This computation requires a `solver` that can handle
exponential cone constraints.
"""
function stationary_covariance_ellipsoid(MP::MarkovProcess, v::Vector{<:APL}, d::Int, solver, P::Partition = trivial_partition(MP.X))
    n = length(v)
    model = SOSModel(solver)
	w = Dict(v => (props(P.graph, v)[:cell] isa Singleton ?
                   @variable(model) :
                   @variable(model, [1], Poly(monomials(MP.x, 0:d)))[1]) for v in vertices(P.graph))
    @variable(model, s)
    @variable(model, S[1:n+1, 1:n+1], PSD)
    @variable(model, U[1:2*n, 1:2*n], PSD)
    @variable(model, r[1:n])
    @variable(model, q[1:n])

	@constraint(model, S[1:n,1:n] .== U[1:n, 1:n])
    @constraint(model, [i in 1:n, j in 1:i], -2*U[n+i,j] - (i == j ? U[n+i,n+i] - r[i] : 0) == 0)
    @constraint(model, [i in 1:n], [-1, q[i], r[i]] in MOI.DualExponentialCone())

	for v in vertices(P.graph)
        flag = (v_target == v)
        add_stationarity_constraints!(model, MP, v, P, props(P.graph, v)[:cell], w, s + (v'*S[1:n, 1:n]*v + 2*S[end,1:n]'*v))
    end

	for e in edges(P.graph)
        add_coupling_constraints!(model, MP, e, P, w)
    end

    @objective(model, Max, s-S[n+1,n+1]-sum(q))
    optimize!(model)
    return Bound(exp(-objective_value(model)), model, P, Dict(key => value(w[key]) for key in keys(w)))
end

stationary_covariance_ellipsoid(RP::ReactionProcess, v::Vector, d::Int, solver, P::Partition = trivial_partition(RP.JumpProcess.X)) =
								stationary_variance(RP.JumpProcess, [RP.species_to_state[x] for x in v], d, solver, P)
stationary_covariance_ellipsoid(MP::MarkovProcess, v::Vector{Num}, d::Int, solver, P::Partition = trivial_partition(MP.X)) =
								stationary_variance(MP, polynomialize_expr(v, MP.poly_vars), d, solver, P)

@doc raw"""
	stationary_covariance_ellipsoid(rn::ReactionSystem, S0::Dict, S::AbstractVector, d::Int, solver,
					scales = Dict(s => 1 for s in species(rn));
					auto_scaling = false)

returns an **upper** on the volume of the covariance ellipsoid of any subset `S`
of the chemical species in the reaction network `rn`, i.e., ``\text{det}(\mathbb{E}[SS^\top] - \mathbb{E}[S] \mathbb{E}[S]^\top)``,
at steady state of the associated jump process. The reaction network is assumed
to have the deterministic initial state `S0` (all species must be included here!).
The bounds are computed via an SOS program of degree `d`, hence can be tightened
by increasing `d`. This computation requires a solver that can deal with
exponential cone constraints.

For numerical stability, it is recommended to provide scales of the expected
magnitude of molecular counts for the different species at steady state. If
the system is **closed** it is also possible to enable auto_scaling which will
find the maximum molecular counts for each species under stoichiometry
constraints (via LP).

If the initial condition of the reaction network under investigation is
unknown or irrelevant, simply call

	stationary_covariance_ellipsoid(rn::ReactionSystem, S, d::Int, solver,
					scales = Dict(s => 1 for s in species(rn)))
"""
function stationary_covariance_ellipsoid(rn::ReactionSystem, S0::Dict, S::AbstractVector, d::Int, solver,
	 									 scales = Dict(s => 1 for s in species(rn));
										 auto_scaling = false, params::Dict = Dict())
	RP, x0 = reaction_process_setup(rn, S0, scales = scales, auto_scaling = auto_scaling, solver = solver, params = params)
 	return stationary_covariance_ellipsoid(RP.JumpProcess, [RP.species_to_state[x] for x in S], d, solver)
end

function stationary_covariance_ellipsoid(rn::ReactionSystem, S::Vector, d::Int, solver,
	 									 scales = Dict(s => 1 for s in species(rn));
										 params::Dict = Dict())
	RP = reaction_process_setup(rn, scales = scales, params = params)
	return stationary_covariance_ellipsoid(RP.JumpProcess, [RP.species_to_state[x] for x in S], d, solver)
end

"""
	stationary_probability_mass(MP::MarkovProcess, X::BasicSemialgebraicSet, d::Int,
							solver)

returns **lower** and **upper** bounds on the probability mass associated with the set `X`.
`d` refers to the order of the relaxation used, again the bounds will tighten
monotonically with increasing order. solver refers to the optimizer used to solve
the semidefinite programs which optimal values furnish the bounds. This is the
weakest formulation that can be used to compute bounds on the probabiltiy mass
associated with a Basic semialgebraic set. For sensible results the set `X` should have non-
empty interior. In order to improve the bounds the user should supply a carefully defined
partition of the state space. In that case it is most sensible to choose `X` as a subset of the element
of said partition. Then, one should call

    stationary_probability_mass(MP::MarkovProcess, v::Int, order::Int, solver,
                                 P::Partition)
where `v` refers to the vertex of the partition graph that corresponds to the set `X`.
"""
function stationary_probability_mass(MP::MarkovProcess, X::BasicSemialgebraicSet, d::Int, solver)
	P = split_state_space(MP, X)
 	return stationary_probability_mass(MP, 1, d, solver, P)
end

function stationary_probability_mass(MP::MarkovProcess, v::Int, order::Int, solver, P::Partition)
	model, w = stationary_indicator(MP, v, order, P, solver; sense = 1) # Max
	optimize!(model)
	ub = Bound(-objective_value(model), model, P, Dict(v => value(w[v]) for v in vertices(P.graph)))
	model, w = stationary_indicator(MP, v, order, P, solver; sense = -1) # Min
	optimize!(model)
	lb = Bound(objective_value(model), model, P, Dict(v => value(w[v]) for v in vertices(P.graph)))
 	return lb, ub
end

function stationary_indicator(MP::MarkovProcess, v_target::Int, order::Int, P::Partition, solver; sense = 1)
    model = SOSModel(solver)
    @variable(model, s)
    w = Dict(v => (props(P.graph, v)[:cell] isa Singleton ?
                   @variable(model) :
                   @variable(model, [1], Poly(monomials(MP.x, 0:order)))[1]) for v in vertices(P.graph))
    cons = Dict()
    for v in vertices(P.graph)
        flag = (v_target == v)
        cons[v] = add_stationarity_constraints!(model, MP, v, P, props(P.graph, v)[:cell], w, s + flag*sense)
    end

    for e in edges(P.graph)
        add_coupling_constraints!(model, MP, e, P, w)
    end

    @objective(model, Max, s)
    return model, w
end

"""
    approximate_stationary_measure(MP::MarkovProcess, v::APL, order::Int, solver, P::Partition;
                                   side_infos::BasicSemialgebraicSet)

returns approximate values for the stationary measure on the partition `P`.
`v` is a polynomial observable whose expectation is minimized when determining the approximation.
The choice of `v` shall be understood as means to regularize the problem.
"""
function approximate_stationary_measure(MP::MarkovProcess, v::APL, order::Int, solver, P::Partition,
                                        side_infos = FullSpace())
    model = SOSModel(solver)
    @variable(model, s)
    w = Dict(v => (props(P.graph, v)[:cell] isa Singleton ?
                    @variable(model) :
                    @variable(model, [1], Poly(monomials(MP.x, 0:order)))[1]) for v in vertices(P.graph))
    ineqs = inequalities(side_infos)
    eqs = equalities(side_infos)
    if !isempty(ineqs)
        @variable(model, μ[ineqs] >= 0)
        ineq_slack = sum(μ[ineq]*(ineq - ineq.a[end]) for ineq in ineqs)
        ineq_obj =  sum(μ[ineq]*ineq.a[end] for ineq in ineqs)
    else
        ineq_slack = 0
        ineq_obj = 0
    end
    if !isempty(eqs)
        @variable(model, λ[eqs])
        eq_slack = sum(λ[eq]*(eq - eq.a[end]) for eq in eqs)
        eq_obj = sum(λ[eq]*eq.a[end] for eq in eqs)
    else
        eq_slack = 0
        eq_obj = 0
    end
    slack = - v + s - eq_slack - ineq_slack
    cons = Dict()
    for v in vertices(P.graph)
        cons[v] = add_stationarity_constraints!(model, MP, v, P, props(P.graph, v)[:cell], w, slack)
    end

    for e in edges(P.graph)
        add_coupling_constraints!(model, MP, e, P, w)
    end

    @objective(model, Max, s - eq_slack - ineq_slack)
    optimize!(model)
    dist = []
    for v in vertices(P.graph)
        if props(P.graph, v)[:cell] isa Singleton
            push!(dist, dual(cons[v]))
        elseif props(P.graph, v)[:cell] isa Vector{BasicSemialgebraicSet}
            push!(dist, sum(dual.(cons[v])).a[end])
        elseif props(P.graph, v)[:cell] isa BasicSemialgebraicSet
            push!(dist, dual(cons[v]).a[end])
        end
    end
    return Bound(objective_value(model), model, P, Dict(key => value(w[key]) for key in keys(w))), dist
end

"""
    max_entropy_measure(MP::MarkovProcess, order::Int, solver, P::Partition;
                                   side_infos::BasicSemialgebraicSet)

returns approximate values for the stationary measure which maximizes the entropy on the partition `P`.
"""
function max_entropy_measure(MP::MarkovProcess, order::Int, solver, P::Partition,
                             side_infos = FullSpace())
    model = SOSModel(solver)
    w = Dict(v => (props(P.graph, v)[:cell] isa Singleton ?
                   @variable(model) :
                   @variable(model, [1], Poly(monomials(MP.x, 0:order)))[1]) for v in vertices(P.graph))
    @variable(model, u[vertices(P.graph)])
    @variable(model, q[vertices(P.graph)])
    @variable(model, s)
    ineqs = inequalities(side_infos)
    eqs = equalities(side_infos)
    if !isempty(ineqs)
        @variable(model, μ[ineqs] >= 0)
        ineq_slack = sum(μ[ineq]*(ineq - ineq.a[end]) for ineq in ineqs)
        ineq_obj =  sum(μ[ineq]*ineq.a[end] for ineq in ineqs)
    else
        ineq_slack = 0
        ineq_obj = 0
    end
    if !isempty(eqs)
        @variable(model, λ[eqs])
        eq_slack = sum(λ[eq]*(eq - eq.a[end]) for eq in eqs)
        eq_obj = sum(λ[eq]*eq.a[end] for eq in eqs)
    else
        eq_slack = 0
        eq_obj = 0
    end
    cons = Dict()
    for v in vertices(P.graph)
        cons[v] = add_stationarity_constraints!(model, MP, v, P, props(P.graph, v)[:cell], w, s + q[v] - eq_slack - ineq_slack )
    end

    for e in edges(P.graph)
        add_coupling_constraints!(model, MP, e, P, w)
    end

    @constraint(model, [v in vertices(P.graph)], [-1, q[v], u[v]] in MOI.DualExponentialCone())

    @objective(model, Max, -sum(u) + s)
    optimize!(model)
    dist = []
    for v in vertices(P.graph)
        if props(P.graph, v)[:cell] isa Singleton
            if dual(cons[v]) isa Number
                push!(dist, dual(cons[v]))
            else
                push!(dist, dual(cons[v]).a[end])
            end
        elseif props(P.graph, v)[:cell] isa Vector{BasicSemialgebraicSet}
            push!(dist, sum(dual.(cons[v])).a[end])
        elseif props(P.graph, v)[:cell] isa BasicSemialgebraicSet
            push!(dist, dual(cons[v]).a[end])
        end
    end
    return Bound(-objective_value(model), model, P, Dict(key => value(w[key]) for key in keys(w))), dist
end


text:
export smRecovery

###################################################################
# System Matrix recovery using DCT-sparsity and low-rank constraing
###################################################################
function smRecovery(y::Matrix{T}, samplingIdx::Array{Int64}, params::Dict) where T
  shape = params[:shape]

  # sampling operator
  P = SamplingOp(samplingIdx, shape)

  # sparsifying transformation and regularization
  if !haskey(params, :sparseTrafo)
    params[:sparseTrafo] = DCTOp(ComplexF64,params[:shape],2)
  end
  reg1Name = get(params, :reg1, "L1")
  reg1 = Regularization(reg1Name, params[:λ_l1]; params...)

  # setup low rank regularization
  λ_lr = get(params, :λ_lr, 0.0)
  ρ_lr = get(params, :ρ_lr, 0.0)
  if λ_lr==0 || ρ_lr==0
    # no LR regularization
    reg = [reg1]
    ρ = [params[:ρ_l1]]
  elseif length(shape)==2
    # 2d case (use nuclear norm)
    regNN = regularizationNN(λ_lr, shape)
    reg = [reg1, regNN]
    ρ = [params[:ρ_l1], ρ_lr]
  else
    # 3d case (use tensor nuclear norm)
    regTNN1 = regularizationTNN(λ_lr, shape, 1)
    regTNN2 = regularizationTNN(λ_lr, shape, 2)
    regTNN3 = regularizationTNN(λ_lr, shape, 3)
    reg = [reg1, regTNN1, regTNN2, regTNN3]
    ρ = [params[:ρ_l1], ρ_lr, ρ_lr, ρ_lr]
  end

  # precalculate inverse for split Bregman iteration
  precon = invPrecon(samplingIdx, prod(shape); params...)
  params[:precon] = precon
  params[:iterationsCG] = 1

  # normalized measurement
  y2, y_norm = getNormalizedMeasurement(y)

  solver = [SplitBregman(P; reg=deepcopy(reg), ρ=ρ, params...) for i=1:Threads.nthreads()]

  sfMat = zeros(ComplexF64,prod(shape),size(y,2))
  @time Threads.@threads for k=1:size(y,2)
    t = Threads.threadid()
    sfMat[:,k] .= solve(solver[t], y2[k])
    # undo normalization
    sfMat[:,k] *= y_norm[k]
  end

  return sfMat
end

#########################
# low rank regularization
#########################
function regularizationTNN(λ::AbstractFloat, shape::NTuple{3,Int64}, mode::Int64)
  proxMap = (x,y;kargs...)->proxTNN!(x,y,shape,mode)
  reg = Regularization(prox! = proxMap, λ=λ, params=Dict{Symbol,Any}())
end

function proxTNN!(x, λ::Real, shape::NTuple{3,Int64},mode::Int64; kargs...)
  nx,ny,nz = shape
  # flattened tensor
  x = reshape(x,shape)
  if mode==1
    x_flat = reshape(x,nx,ny*nz)
  elseif mode==2
    x_flat = reshape(permutedims(x,[2,3,1]),ny,nx*nz)
  else
    if mode!=3
      @error "mode in proxTNN! should be 1,2 or 3"
    end
    x_flat = reshape(permutedims(x,[3,1,2]),nz,nx*ny)
  end
  # perform singular value thresholding
  usv = svd!(x_flat)
  proxL1!(usv.S,λ)
  x_flat[:] .= vec(usv.U*Diagonal(usv.S)*usv.Vt)
  # undo flattening operation
  if mode==1
    x[:] .= vec(x_flat)
  elseif mode==2
    x[:] .= vec(permutedims(reshape(x_flat,ny,nz,nx),[3,1,2]))
  else
    x[:] .= vec(permutedims(reshape(x_flat,nz,nx,ny),[2,3,1]))
  end
end

# low rank regularization for 2d (using SVD)
function regularizationNN(λ::Real, shape::NTuple{2,Int64} )
  proxMap = (x,y;kargs...)->proxNN!(x,y,shape;kargs...)
  return Regularization(prox! = proxMap, λ=λ)
end

# singular-value thresholding
function proxNN!(x, λ::Real, shape::NTuple{2,Int64}; kargs...)
  U,S,V = svd(reshape(x,shape))
  proxL1!(S,λ)
  x[:] .= vec(U*Matrix(Diagonal(S))*V')
end

# fixed rank constraing for 3d (using HOSVD)
function FRTruncation!(x::Array{T,3}, rank::NTuple{3,Int64}) where T
  guvw_r = hosvd(real.(x),rank)
  guvw_i = hosvd(imag.(x),rank)
  x[:] .= compose(guvw_r)+1im*compose(guvw_i)
end

# fixed rank constraing for 2d (using SVD)
function FRTruncation!(x, rank::Int64)
  U,S,V = svd(reshape(x, shape))
  S[rank+1:end] .= 0.0
  x[:] .= vec(U*Matrix(Diagonal(S))*V')
end

##########################
# inverse of normal matrix
##########################
mutable struct diagPrecon{T}
  diag::Vector{T}
end

function ldiv!(y::Vector{T},P::diagPrecon{T},x::Vector{T}) where T
  y[:] .= P.diag .* x
end

function ldiv!(P::diagPrecon{T},x::Vector{T}) where T
  x[:] .= P.diag .* x
end

\(P::diagPrecon{T}, x::Vector{T}) where T = P.diag .* x

function invPrecon(samplingIdx, ncol; ρ::Vector{Float64}=[1.0], kargs...)
  normalP = normalPDiag(samplingIdx, ncol)
  diag = [normalP[i]+sum(ρ) for i=1:length(normalP)]
  diag = diag.^-1
  return diagPrecon(diag)
end

function normalPDiag(samplingIdx, ncol; T=ComplexF64)
  normalP = zeros(T,ncol)
  for i=1:length(samplingIdx)
    normalP[samplingIdx[i]] = 1.0
  end
  return normalP
end

function getNormalizedMeasurement(y::Matrix{T}) where T
  y_norm = [norm(y[:,k]) for k=1:size(y,2)]
  y2 = [ComplexF64.(y[:,k]/y_norm[k]) for k=1:size(y,2)]
  return y2, y_norm
end


text:
"""
  `EventNotIdle <: Exception`

Exception thrown when an event is scheduled ([`schedule`](@ref)) that has already been scheduled or is being processed.

Only used internally.
"""

type EventNotIdle <: Exception end

"""
  `EventProcessing <: Exception`

Exception thrown:

- when a callback is added to an event ([`append_callback`](@ref)) that is being processed or
- when an event is scheduled ([`schedule!`](@ref)) that is being processed.

Only used internally.
"""
type EventProcessed <: Exception end

"""
  `StopSimulation <: Exception`

Exception that stops the simulation. A return value can be set.

**Fields**:

- `value :: Any`

**Constructor**:

`StopSimulation(value::Any=nothing)`
"""
type StopSimulation <: Exception
  value :: Any
  function StopSimulation(value::Any=nothing)
    new(value)
  end
end

type EmptySchedule <: Exception end

type InterruptException <: Exception
  cause :: Any
  function InterruptException(cause::Any)
    inter = new()
    inter.cause = cause
    return inter
  end
end


text:
@testset "Rosenbrock" begin

    function test_result(result::Vector, fitness::Float64, N::Int, tol::Float64)
        @test length(result) == N
    	@test ≈(result, ones(N), atol=tol)
        @test ≈(fitness, 0.0, atol=tol)
    end

	# Objective function
	rosenbrock(x::Vector{Float64}) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2

	# Parameters
	N = 2
	terminate(σ) = σ < 1e-10

	# Testing: (15/5+100)-ES
	# with isotropic mutation operator y' := y + σ(N_1(0, 1), ..., N_N(0, 1))
	result, fitness, cnt = es(rosenbrock, N;
        initStrategy = strategy(σ = 1.0),
        recombination = average, mutation = isotropic,
        μ = 15, ρ = 5, λ = 100, iterations = 1000)
	println("(15/5+100)-ES => F: $(fitness), C: $(cnt), OBJ: $(result)")
	test_result(result, fitness, N, 1e-1)

	# Testing: (15/15+100)-σ-Self-Adaptation-ES
	# with isotropic mutation operator y' := y + σ(N_1(0, 1), ..., N_N(0, 1))
	result, fitness, cnt = es(rosenbrock, N;
		initPopulation = [.5, .5],
        initStrategy = strategy(σ = 1.0, τ = 1/sqrt(2*N)),
		recombination = average, srecombination = averageSigma1,
		mutation = isotropic, smutation = isotropicSigma,
		μ = 15, λ = 100, iterations = 1000)
	println("(15/15+100)-σ-SA-ES-IS => F: $(fitness), C: $(cnt), OBJ: $(result)")
	test_result(result, fitness, N, 1e-1)

	# Testing: (15/15+100)-σ-Self-Adaptation-ES
	# with non-isotropic mutation operator y' := y + (σ_1 N_1(0, 1), ..., σ_N N_N(0, 1))
	result, fitness, cnt = es(rosenbrock, N;
		initPopulation = rand(N,25),
        initStrategy = strategy(σ = .5ones(N), τ = 1/sqrt(2*N), τ0 = 1/sqrt(N)),
		recombination = average, srecombination = averageSigmaN,
		mutation = anisotropic, smutation = anisotropicSigma,
		μ = 15, λ = 100, iterations = 1000)
	println("(15/15+100)-σ-SA-ES-AS => F: $(fitness), C: $(cnt), OBJ: $(result)")
	test_result(result, fitness, N, 1e-1)

	# Testing: CMA-ES
	result, fitness, cnt = cmaes(rosenbrock, N;
    	μ = 3, λ = 12, iterations = 100_000, tol = 1e-3)
    println("(3/3,12)-CMA-ES => F: $(fitness), C: $(cnt), OBJ: $(result)")
    test_result(result, fitness, N, 1e-1)

    # Testing: GA
    result, fitness, cnt = ga(rosenbrock, N;
    	initPopulation = (n -> rand(n)),
    	populationSize = 100,
    	ɛ = 0.1,
        selection = sus,
        crossover = intermediate(0.25),
        mutation = domainrange(fill(0.5,N)))
    println("GA(p=50,x=0.8,μ=0.1) => F: $(fitness), C: $(cnt), OBJ: $(result)")
    test_result(result, fitness, N, 3e-1)

end


text:
const TAIL_MIN_ASSETS = 0.
const TAIL_NUM_YEARS = 5

function performancebyein(data::NCCSData; tailminassets = TAIL_MIN_ASSETS,
  numyears = TAIL_NUM_YEARS, lastyear = 2015,
  outputpath::String = TABLE_PATH,
  identifierfields::Vector{Symbol} = [:ein, :name, :fisyr, :datefiscal],
  return1y::Symbol = :lreturn,
  wealthfield::Symbol = :adjnetassets,
  sumfields::Vector{Symbol} = [:monthscovered, :lt30ret, :lsp500ret,
    :bmff3, :bmff5, :bmsp500],
  annualizefields::Bool = true)::DataFrame

  #make some assertions on what must be present in the field lists
  @assert :monthscovered ∈ sumfields
  @assert :fisyr ∈ identifierfields

  local firstyear::Int = lastyear - numyears
  local fields2annualize::Vector{Symbol}

  #partition off the part we want
  validfirmsdf = view(data.df, (data.df[wealthfield] .> tailminassets) .&
    (data.df[:fisyr] .≤ lastyear) .& (data.df[:fisyr] .> firstyear),
    [identifierfields; return1y; wealthfield; sumfields;]
    )

  #form the consolidated df

  taildf::DataFrame = DataFrame(ein=validfirmsdf[:ein])
  N::Int = size(taildf, 1)

  #pre-allocate the columns
  for f ∈ setdiff(names(validfirmsdf), [:ein])
    local newtype::Type
    local origtype::Type = eltype(data.df[f])

    if Float64 <: origtype
      newtype = MFloat64
    elseif Int <: origtype
      newtype = MInt
    elseif Date <: origtype
      newtype = MDate
    elseif String <: origtype
      newtype = MString
    else
      newtype = MSymbol
    end

    println("field $f: oldtype= $(origtype) newtype = $newtype")

      taildf[f] = Vector{Union{Missing, newtype}}(missing, N)
  end
  taildf[:N] = Vector{MInt}(undef, N)

  #create the same set of fields but annualized
  #NOTE: Annualized fields are not logged
  if annualizefields
    for f ∈ [return1y; sumfields]
      f_ann::Symbol = Symbol(f, "_ann")
      taildf[f_ann] = similar(taildf[f])
    end

    fields2annualize = setdiff([return1y; sumfields;], [:monthscovered])
  end

  taildfindex::Dict = Dict(taildf[i,:ein] => taildf[i,:] for i::Int  ∈ 1:N)

  #consolidate the eprformacne data
  for subdf ∈ groupby(validfirmsdf, :ein)
    tailrow = taildfindex[subdf[1,:ein]]

    #identifier fields receive the value from the first row
    #NOTE: No verification that identifier fields are the same for a given org
    for f ∈ identifierfields
      tailrow[f] = subdf[1,f]
    end

    tailrow[return1y] = sum(subdf[return1y])

    #subtotal the sum fields
    for f ∈ sumfields
      tailrow[f] = sum(subdf[f])
    end

    #if we want annualized copies
    if annualizefields
      yearscoveredinv::MFloat64 = 1.0 / (tailrow[:monthscovered] / 12.0 )

      for f ∈ fields2annualize
        tailrow[Symbol(f,"_ann")] = exp(tailrow[f] * yearscoveredinv) - 1.0
      end
    end

    #get average assets
    tailrow[wealthfield] = mean(subdf[wealthfield]) ./ 1_000_000

    tailrow[:N] = size(subdf,1)
  end

  taildf = taildf[completecases(taildf[[identifierfields; return1y; wealthfield]]),:]
  taildf = taildf[taildf[:N] .== numyears,:]

  return taildf
end

#writes out the data to the tails
function summarizesimpletails(data::NCCSData;
  numyears = TAIL_NUM_YEARS, lastyear = 2015,
  return1y::Symbol = :lreturn,
  numberpertail::Int = 200,
  outputpath::String = TABLE_PATH,
  datasuffix::String = DATA_SUFFIX)::Nothing

  firstyear::Int = lastyear - numyears
  outpath = "$outputpath\\tails_$(firstyear)_$(lastyear)_$return1y"

  #get a frame with cumulative returns by year
  taildf::DataFrame = performancebyein(data, numyears=numyears, lastyear=lastyear,
    return1y=return1y)

  #organize the return dataframe
  sort!(taildf, [return1y])
  N::Int = size(taildf, 1)

  taildfbottom::SubDataFrame = view(taildf, 1:numberpertail,:)
  taildftop::SubDataFrame = view(taildf, (N-numberpertail + 1):N,:)

  #finally write out all the tail data
  taildftop |> CSV.write("$(outpath)_top_$datasuffix.csv")
  taildfbottom |> CSV.write("$(outpath)_bottom_$datasuffix.csv")

  return nothing
end

#extract the tails
function extractpersistencetails(df::DataFrame, period::Int, endyear::Int, lbenchmark::NSymbol=nothing;
    outputpath::String = WORKING_PATH, threshold::Float64 = TAIL_THRESHOLD,
    minentriesperfile::Int = TAIL_MIN_ENTRIES_PER_FILE,
    consecutive::Bool = false, datasuffix::String=DATA_SUFFIX)::Nothing

  local subdf::SubDataFrame
  local topdf::SubDataFrame
  local bottomdf::SubDataFrame
  local topdfout::DataFrame
  local bottomdfout::DataFrame

  #get the field symbols given the period and benchmark
  local (percentilesym, Lpercentilesym) = buildpercentilesym(period, lbenchmark)
  local returnsym = Symbol("lreturn$(period)yr")

  #make the windows non-overlapping
  subdf = view(df, df[:fisyr] .≤ endyear, :)
  subdf = view(subdf, ((endyear .- subdf[:fisyr]) .% period) .== 0, :)

  #no missing data
  subdf = view(subdf, (!ismissing).(subdf[percentilesym]), :)

  #split out views for the perisstent tails
  topdf = view(subdf, (subdf[percentilesym] .≥ (1.0 - threshold)), :)
  bottomdf = view(subdf, (subdf[percentilesym] .≤ threshold), :)

  #use a different function for consecutive tails

  topdfout = aggregatetails(topdf, period, returnsym, percentilesym,
    minentriesperfile=minentriesperfile, consecutive=consecutive)
  bottomdfout = aggregatetails(bottomdf, period, returnsym, percentilesym,
    minentriesperfile=minentriesperfile, consecutive=consecutive)

  #sort to enable easy ranking
  topdfout = sort(topdfout,
    (order(:numappearances, rev=true),
    order(:avgreturnannual, rev=true)))

  bottomdfout = sort(bottomdfout,
    (order(:numappearances, rev=true),
    order(:avgreturnannual, rev=false)))

  #build the file outname
  outname::String = "perstails_per$(period)_minrec$(minentriesperfile)_thresh$threshold"
  consecutive && (outname = "$(outname)_consec")

  topdfout |> CSV.write("$outputpath\\$(outname)_top_$datasuffix.csv")
  bottomdfout |> CSV.write("$outputpath\\$(outname)_bottom_$datasuffix.csv")

  return nothing
end

#constructs the output dataframe
function aggregatetails(df::AbstractDataFrame, period::Int,
  returnsym::Symbol, percentilesym::Symbol;
  minentriesperfile::Int = TAIL_MIN_ENTRIES_PER_FILE,
  consecutive::Bool = false)

  local outdf::DataFrame
  local aggregatefunction!::Function
  local validrows::Vector{Bool}
  local appearancescutoff::Int

  #we want info at the non-profit level, so create a way to that that info
  outdf = DataFrame(ein = unique(df[:ein]))
  outN::Int = size(outdf,1)

  #this contains the information we want to look at
  outdf[:name] = Vector{MString}(undef, outN)
  outdf[:categorynamefiltered] = Vector{MSymbol}(undef, outN)
  outdf[:numappearances] = Vector{MInt}(undef, outN)
  outdf[:avgreturn] = Vector{MFloat64}(undef, outN)
  outdf[:avgreturnannual] = Vector{MFloat64}(undef, outN)
  outdf[:avgquantile] = Vector{MFloat64}(undef, outN)
  outdf[:appearances] = Vector{MString}(undef, outN)

  outindex::Dict = Dict(outdf[i,:ein] => outdf[i,:] for i ∈ 1:size(outdf,1))

  #NOTE: the kernel function called to aggreagte for each ein is set here
  if consecutive
    aggregatefunction! = ((subdf::AbstractDataFrame, r::DataFrameRow)->
      aggregateconsecutivetails!(subdf,r,period, returnsym, percentilesym)::Nothing)
  else
    aggregatefunction! = ((subdf::AbstractDataFrame, r::DataFrameRow)->
      aggregatetotaltails!(subdf,r,period, returnsym, percentilesym)::Nothing)
  end

  for ssubdf::SubDataFrame ∈ groupby(df, :ein)
    outrow::DataFrameRow = outindex[ssubdf[1,:ein]]
    aggregatefunction!(ssubdf, outrow)
  end

  validrows = falses(size(outdf, 1))
  appearancescutoff = maximum(outdf[:numappearances])

  while (sum(validrows) < minentriesperfile) && (appearancescutoff ≥ 1)
    appearancescutoff -= 1
    validrows .= ((i::MInt) -> (!ismissing(i)) && (i > appearancescutoff)).(outdf[:numappearances])
  end

   #outdf = DataFrame(view(outdf, validrows, :))

  return outdf[validrows, :]
end

#does the processing of the dataframe in the event we want a consecutive sequence
function aggregateconsecutivetails!(df::AbstractDataFrame,
    outrow::DataFrameRow, period::Int, returnsym::Symbol,
    percentilesym::Symbol)::Nothing

  N::Int = size(df, 1)

  local currentctr::Int = 0
  local maxconsecutive = 0

  local runningquantiletotal::Float64 = 0.0
  local runningreturntotal::Float64 = 0.0
  local runningstring::String = ""

  for i ∈ 1:N
    #first check if its a consecutive appearance
    currentctr += 1
    year::Int = df[i,:fisyr]

    #currentctr==1 is a special case
    if (i == 1) || (df[i-1,:fisyr] == year - period)
      runningreturntotal += df[i, returnsym]
      runningquantiletotal += df[i, percentilesym]
      runningstring = "$(runningstring)$(year), "

      if (maxconsecutive < currentctr)
        #record data on the max consecutive tail run
        maxconsecutive = currentctr
        outrow[:avgreturn] = runningreturntotal
        outrow[:avgquantile] = runningquantiletotal
        outrow[:appearances] = runningstring
      end
    else #otherwise reset the counter and the running totals
      runningreturntotal = df[i, returnsym]
      runningquantiletotal = df[i, percentilesym]
      runningstring = "$year, "
      currentctr = 1
    end
    #println("year: $year currentctr: $currentctr maxconsecutive: $maxconsecutive")
  end

  #=if maxconsecutive* 6 ≠ length(outrow[:appearances])
    println("output: $outrow")
    println(df[[:ein, :fisyr, returnsym, percentilesym]])
    println("period: $period, maxconsecutive: $maxconsecutive")
  end=#
  #@assert maxconsecutive* 6 == length(outrow[:appearances])

  outrow[:avgreturn] *= 1/maxconsecutive
  outrow[:avgreturnannual] = outrow[:avgreturn] / period
  outrow[:avgquantile] *= 1/maxconsecutive
  outrow[:numappearances] = maxconsecutive
  outrow[:name] = df[1,:name]
  outrow[:categorynamefiltered] = df[1,:categorynamefiltered]

  return nothing
end

#in this simpler case, we care only about tail appearances not consecutive appearances
function aggregatetotaltails!(df::AbstractDataFrame,
  outrow::DataFrameRow, period::Int, returnsym::Symbol, percentilesym::Symbol)::Nothing

  outrow[:numappearances] = size(df,1)
  outrow[:avgreturn] = mean(df[returnsym])
  outrow[:avgreturnannual] = outrow[:avgreturn] / period
  outrow[:avgquantile] = mean(df[percentilesym])
  outrow[:appearances] = join((string).(df[:fisyr]), ", ")
  outrow[:name] = df[1,:name]
  outrow[:categorynamefiltered] = df[1,:categorynamefiltered]

  return nothing
end

function extractpersistencetails(data::NCCSData;
  periods::Vector{Int} = PERSISTENCE_PERIODS,
  persistenceyears::Vector{Int} = PERSISTENCE_YEARS
  )::Nothing

  for p ∈ periods, y ∈ persistenceyears
    extractpersistencetails(data.df, p, y, consecutive = false)
    extractpersistencetails(data.df, p, y, consecutive = true)
  end

  return nothing
end

function summarizetails(data::NCCSData)
  #NOTE: Uncomment the below when we want simple tails
  #summarizesimpletails(data)
  extractpersistencetails(data)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
const G0 = GRAVITY_ACCEL

"""
    atmosphere_isa(height)

Calculate temperature, pressure, density and sound velocity for the
given geopotential height according to International Standard Atmosphere 1976.

# References

- [1] U.S. Standard Atmosphere, 1976, U.S. Government Printing Office,
        Washington, D.C., 1976

From: https://en.wikipedia.org/wiki/U.S._Standard_Atmosphere

| Layer | h (m) | p (Pa)  | T (K)  | ``α`` (K/m) |
|-------|-------|---------|--------|-------------|
| 0     | 0     | 101325  | 288.15 | -0.0065     |
| 1     | 11000 | 22632.1 | 216.65 | 0           |
| 2     | 20000 | 5474.89 | 216.65 | 0.001       |
| 3     | 32000 | 868.019 | 228.65 | 0.0028      |
| 4     | 47000 | 110.906 | 270.65 | 0           |
| 5     | 51000 | 66.9389 | 270.65 | -0.0028     |
| 6     | 71000 | 3.95642 | 214.65 | -0.002      |
"""
function atmosphere_isa(height)

    if 0 <= height < 11000  # Troposphere
        alpha = -0.0065  # K/m
        T0 = 288.15      # K
        p0 = 101325.0    # Pa

        T = T0 + alpha * height
        p = p0 * (T0 / T) ^ (G0 / (R_AIR * alpha))

    elseif 11000 <= height < 20000  # Tropopause
        T = 216.65    # K
        p0 = 22632.1  # Pa
        h0 = 11000    # m

        p = p0 * exp(-G0 * (height - h0) / (R_AIR * T))

    elseif 20000 <= height < 32000  # Stratosphere 1
        alpha = 0.001    # K/m
        T0 = 216.65      # K
        p0 = 5474.89     # Pa
        h0 = 20000       # m

        T = T0 + alpha * (height - h0)
        p = p0 * (T0 / T) ^ (G0 / (R_AIR * alpha))

    elseif 32000 <= height < 47000  # Stratosphere 2
        alpha = 0.0028   # K/m
        T0 = 228.65      # K
        p0 = 868.019     # Pa
        h0 = 32000       # m

        T = T0 + alpha * (height - h0)
        p = p0 * (T0 / T) ^ (G0 / (R_AIR * alpha))

    elseif 47000 <= height < 51000  # Stratopause
        T = 270.65    # K
        p0 = 110.906  # Pa
        h0 = 47000    # m
        h0 = 47000    # m

        p = p0 * exp(-G0 * (height - h0) / (R_AIR * T))

    elseif 51000 <= height < 71000  # Mesosphere 1
        alpha = -0.0028  # K/m
        T0 = 270.65      # K
        p0 = 66.9389     # Pa
        h0 = 51000       # m

        T = T0 + alpha * (height - h0)
        p = p0 * (T0 / T) ^ (G0 / (R_AIR * alpha))

    elseif 71000 <= height <= 84500  # Mesosphere 2
        alpha = -0.002   # K/m
        T0 = 214.65      # K
        p0 = 3.95642     # Pa
        h0 = 71000       # m

        T = T0 + alpha * (height - h0)
        p = p0 * (T0 / T) ^ (G0 / (R_AIR * alpha))

    else
        throw(DomainError(height, "height must be between 0 m and 84500 m"))
    end

    rho = p / (R_AIR * T)
    a   = sqrt(GAMMA_AIR * R_AIR * T)

    return [T, p, rho, a]
end

# TODO: document difference in results
"""
    atmosphere_f16(height)

Standard Atmosphere implemented in Stevens F16 model.

# Arguments
- alt: altitude [m].

# Returns
- T: temperature [K].
- P: pressure [Pa].
- RHO: density [kg/m³].
- a: Speed of sound [m/s].

# Notes
- This model is implemented just for F16 complete model validation purposes.
- For ISA 1978 atmospheric model, use atmosphere_isa.
- Note limited precision in tests against `atmosphere_isa`.

# References
- [1] Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control
 and simulation: dynamics, controls design, and autonomous systems. John Wiley
 & Sons. (page 715)
"""
function atmosphere_f16(alt)

    # The caluclations are prepared for input in ft
    ALT = alt * M2FT

    RO = 0.002377  # sluf/ft³
    TFAC = 1.0 - ALT*0.703E-05
    T = 519.0*TFAC  # ºR

    if ALT >= 35000.0 T=390.0 end
    RHO = RO*(TFAC^4.14)
    # will be calculated in aerostate
    # RMACH = VT/sqrt(1.4*1716.3*T)
    # QBAR = 0.5*RHO*VT*VT
    PS = 1715.0*RHO*T  # psf
    a = sqrt(1.4*1716.3*T)
    return [T*RANK2KEL, PS*PSF2PA, RHO*SLUGFT3_2_KGM3, a*FT2M]
end


text:
"""
$(TYPEDEF)

A system of difference equations.

# Fields
$(FIELDS)

# Example

```
using ModelingToolkit

@parameters σ=28.0 ρ=10.0 β=8/3 δt=0.1
@variables t x(t)=1.0 y(t)=0.0 z(t)=0.0
D = Difference(t; dt=δt)

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

@named de = DiscreteSystem(eqs,t,[x,y,z],[σ,ρ,β]) # or
@named de = DiscreteSystem(eqs)
```
"""
struct DiscreteSystem <: AbstractTimeDependentSystem
    """The differential equations defining the discrete system."""
    eqs::Vector{Equation}
    """Independent variable."""
    iv::Sym
    """Dependent (state) variables. Must not contain the independent variable."""
    states::Vector
    """Parameter variables. Must not contain the independent variable."""
    ps::Vector
    """Array variables."""
    var_to_name
    """Control parameters (some subset of `ps`)."""
    ctrls::Vector
    """Observed states."""
    observed::Vector{Equation}
    """
    Name: the name of the system
    """
    name::Symbol
    """
    systems: The internal systems. These are required to have unique names.
    """
    systems::Vector{DiscreteSystem}
    """
    defaults: The default values to use when initial conditions and/or
    parameters are not supplied in `DiscreteProblem`.
    """
    defaults::Dict
    """
    type: type of the system
    """
    connection_type::Any
    function DiscreteSystem(discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed, name, systems, defaults, connection_type; checks::Bool = true)
        if checks
            check_variables(dvs, iv)
            check_parameters(ps, iv)
            all_dimensionless([dvs;ps;iv;ctrls]) ||check_units(discreteEqs)
        end
        new(discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed, name, systems, defaults, connection_type)
    end
end

"""
    $(TYPEDSIGNATURES)

Constructs a DiscreteSystem.
"""
function DiscreteSystem(
                   eqs::AbstractVector{<:Equation}, iv, dvs, ps;
                   controls = Num[],
                   observed = Num[],
                   systems = DiscreteSystem[],
                   name=nothing,
                   default_u0=Dict(),
                   default_p=Dict(),
                   defaults=_merge(Dict(default_u0), Dict(default_p)),
                   connection_type=nothing,
                   kwargs...,
                  )
    name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro"))
    eqs = collect(eqs)
    iv′ = value(iv)
    dvs′ = value.(dvs)
    ps′ = value.(ps)
    ctrl′ = value.(controls)

    if !(isempty(default_u0) && isempty(default_p))
        Base.depwarn("`default_u0` and `default_p` are deprecated. Use `defaults` instead.", :DiscreteSystem, force=true)
    end
    defaults = todict(defaults)
    defaults = Dict(value(k) => value(v) for (k, v) in pairs(defaults))

    var_to_name = Dict()
    process_variables!(var_to_name, defaults, dvs′)
    process_variables!(var_to_name, defaults, ps′)

    sysnames = nameof.(systems)
    if length(unique(sysnames)) != length(sysnames)
        throw(ArgumentError("System names must be unique."))
    end
    DiscreteSystem(eqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, name, systems, defaults, connection_type, kwargs...)
end

function DiscreteSystem(eqs, iv=nothing; kwargs...)
    eqs = collect(eqs)
    # NOTE: this assumes that the order of algebric equations doesn't matter
    diffvars = OrderedSet()
    allstates = OrderedSet()
    ps = OrderedSet()
    # reorder equations such that it is in the form of `diffeq, algeeq`
    diffeq = Equation[]
    algeeq = Equation[]
    # initial loop for finding `iv`
    if iv === nothing
        for eq in eqs
            if !(eq.lhs isa Number) # assume eq.lhs is either Differential or Number
                iv = iv_from_nested_difference(eq.lhs)
                break
            end
        end
    end
    iv = value(iv)
    iv === nothing && throw(ArgumentError("Please pass in independent variables."))
    for eq in eqs
        collect_vars_difference!(allstates, ps, eq.lhs, iv)
        collect_vars_difference!(allstates, ps, eq.rhs, iv)
        if isdifferenceeq(eq)
            diffvar, _ = var_from_nested_difference(eq.lhs)
            isequal(iv, iv_from_nested_difference(eq.lhs)) || throw(ArgumentError("A DiscreteSystem can only have one independent variable."))
            diffvar in diffvars && throw(ArgumentError("The difference variable $diffvar is not unique in the system of equations."))
            push!(diffvars, diffvar)
            push!(diffeq, eq)
        else
            push!(algeeq, eq)
        end
    end
    algevars = setdiff(allstates, diffvars)
    # the orders here are very important!
    return DiscreteSystem(append!(diffeq, algeeq), iv, vcat(collect(diffvars), collect(algevars)), ps; kwargs...)
end

"""
    $(TYPEDSIGNATURES)

Generates an DiscreteProblem from an DiscreteSystem.
"""
function DiffEqBase.DiscreteProblem(sys::DiscreteSystem,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();
                                    eval_module = @__MODULE__,
                                    eval_expression = true,
                                    kwargs...)
    dvs = states(sys)
    ps = parameters(sys)
    eqs = equations(sys)
    eqs = linearize_eqs(sys, eqs)
    defs = defaults(sys)
    iv = get_iv(sys)

    if parammap isa Dict
        u0defs = merge(parammap, defs)
    elseif eltype(parammap) <: Pair
        u0defs = merge(Dict(parammap), defs)
    elseif eltype(parammap) <: Number
        u0defs = merge(Dict(zip(ps, parammap)), defs)
    else
        u0defs = defs
    end
    if u0map isa Dict
        pdefs = merge(u0map, defs)
    elseif eltype(u0map) <: Pair
        pdefs = merge(Dict(u0map), defs)
    elseif eltype(u0map) <: Number
        pdefs = merge(Dict(zip(dvs, u0map)), defs)
    else
        pdefs = defs
    end

    u0 = varmap_to_vars(u0map,dvs; defaults=u0defs)

    rhss = [eq.rhs for eq in eqs]
    u = dvs
    p = varmap_to_vars(parammap,ps; defaults=pdefs)

    f_gen = generate_function(sys; expression=Val{eval_expression}, expression_module=eval_module)
    f_oop, _ = (@RuntimeGeneratedFunction(eval_module, ex) for ex in f_gen)
    f(u,p,iv) = f_oop(u,p,iv)
    fd = DiscreteFunction(f, syms=Symbol.(dvs))
    DiscreteProblem(fd,u0,tspan,p;kwargs...)
end

function linearize_eqs(sys, eqs=get_eqs(sys); return_max_delay=false)
    unique_states = unique(operation.(states(sys)))
    max_delay = Dict(v=>0.0 for v in unique_states)

    r = @rule ~t::(t -> istree(t) && any(isequal(operation(t)), operation.(states(sys))) && is_delay_var(get_iv(sys), t)) => begin
        delay = get_delay_val(get_iv(sys), first(arguments(~t)))
        if delay > max_delay[operation(~t)]
            max_delay[operation(~t)] = delay
        end
        nothing
    end
    SymbolicUtils.Postwalk(r).(rhss(eqs))

    if any(values(max_delay) .> 0)

        dts = Dict(v=>Any[] for v in unique_states)
        state_ops = Dict(v=>Any[] for v in unique_states)
        for v in unique_states
            for eq in eqs
                if isdifferenceeq(eq) && istree(arguments(eq.lhs)[1]) && isequal(v, operation(arguments(eq.lhs)[1]))
                    append!(dts[v], [operation(eq.lhs).dt])
                    append!(state_ops[v], [operation(eq.lhs)])
                end
            end
        end

        all(length.(unique.(values(state_ops))) .<= 1) || error("Each state should be used with single difference operator.")

        dts_gcd = Dict()
        for v in keys(dts)
            dts_gcd[v] = (length(dts[v]) > 0) ? first(dts[v]) : nothing
        end

        lin_eqs = [
            v(get_iv(sys) - (t)) ~ v(get_iv(sys) - (t-dts_gcd[v]))
            for v in unique_states if max_delay[v] > 0 && dts_gcd[v]!==nothing for t in collect(max_delay[v]:(-dts_gcd[v]):0)[1:end-1]
        ]
        eqs = vcat(eqs, lin_eqs)
    end
    if return_max_delay return eqs, max_delay end
    eqs
end

function get_delay_val(iv, x)
    delay = x - iv
    isequal(delay > 0, true) && error("Forward delay not permitted")
    return -delay
end

check_difference_variables(eq) = check_operator_variables(eq, Difference)

function generate_function(
        sys::DiscreteSystem, dvs = states(sys), ps = parameters(sys);
        kwargs...
    )
    eqs = equations(sys)
    foreach(check_difference_variables, eqs)
    # substitute x(t) by just x
    rhss = [eq.rhs for eq in eqs]

    u = map(x->time_varying_as_func(value(x), sys), dvs)
    p = map(x->time_varying_as_func(value(x), sys), ps)
    t = get_iv(sys)

    build_function(rhss, u, p, t; kwargs...)
end


text:
#
# Test signal
#

## Li & Stoica fig1
spectral_lines_f = [0.04, 0.05, 0.08, 0.09,    0.16:0.005:0.35]
spectral_lines_p = vec(repmat([pi/4], length(spectral_lines_f)))
spectral_lines_a = [10, 10, 10, 2,   vec(repmat([0.5], length(spectral_lines_f)-4, 1))]

## Simulation properties
N = 128
noise_mean = 0
noise_vari = 1

## Create example signal
function create_example(N, f, a, p)
    z = vec(zeros(N, 1))
    n = 1.0:N
    for i in 1:length(f)
        z = z + a[i] * sin(2 * pi * f[i] * n + p[i])
    end
    return z
end

# For passing in known cov matrix
function covariance_matrix{T <: FloatingPoint}(s::Array{T}; M::Int=int(length(s)/2))
    N = length(s)
    L = N - M + 1
    CovMatrix = zeros(M, M)
    for l = 0:N-M
        CovMatrix += s[l+1 : 1 : l+M-1+1] * s[l+1 : 1 : l+M-1+1]'
    end
    CovMatrix ./ L
end

# For reproducible results
srand(1)

#
# Test example from Li & Stoica
#

signal = create_example(N, spectral_lines_f, spectral_lines_a, spectral_lines_p)
noise  = sqrt(noise_vari) * randn(size(signal))
sn     = signal + noise
MSE    = mean((abs(apes(sn, spectral_lines_f)) .- spectral_lines_a).^2)

@test_approx_eq abs(apes(sn, spectral_lines_f)[1:4])  [10.150160442095107, 9.871766839666593, 10.02406973575788, 2.131756870868149]
@test_approx_eq MSE  0.10794065421715185

#
# Test at lower SNR
#

signal  = create_example(N, spectral_lines_f, spectral_lines_a, spectral_lines_p)
noise   = sqrt(noise_vari) * randn(size(signal))
sn      = signal + noise * 50
MSE     = mean((abs(apes(sn, spectral_lines_f)) .- spectral_lines_a).^2)
MSE_w_Q = mean((abs(apes(sn, spectral_lines_f, Q = covariance_matrix(noise))) .- spectral_lines_a).^2)

@test_approx_eq MSE     92.84644647200092
@test_approx_eq MSE_w_Q 9.936179958764319