text:
using CUDA, Knet, JLD2, FileIO, Test

@testset "jld2convert" begin
    f = joinpath(tempdir(), "test.jld2")
    struct CuStruct; a::CuArray; end
    struct KnetStruct; a::KnetArray; end

    for (atype, stype) in ((KnetArray, KnetStruct), (CuArray, CuStruct))

        # Single Array
        x = atype(rand(4))
        save(f, "x", x)
        @test x == load(f, "x")

        # In Array container
        a = [x, x, x]
        save(f, "a", a)
        @test a == load(f, "a")

        # In Tuple container
        t = (x, x, x)
        save(f, "t", t)
        @test t == load(f, "t")

        # In Dict container
        d = Dict(1=>x)
        save(f, "d", d)
        @test d == load(f, "d")

        # In Struct container
        s = stype(x)
        save(f, "s", s)
        @test x == load(f, "s").a
    end
end


text:
using ..GLMakie: enable_SSAO, transparency_weight_scale

@enum Shape CIRCLE RECTANGLE ROUNDED_RECTANGLE DISTANCEFIELD TRIANGLE
@enum CubeSides TOP BOTTOM FRONT BACK RIGHT LEFT

struct Grid{N,T <: AbstractRange}
    dims::NTuple{N,T}
end
Base.ndims(::Grid{N,T}) where {N,T} = N

Grid(ranges::AbstractRange...) = Grid(ranges)
function Grid(a::Array{T,N}) where {N,T}
    s = Vec{N,Float32}(size(a))
    smax = maximum(s)
    s = s ./ smax
    Grid(ntuple(Val{N}) do i
        range(0, stop=s[i], length=size(a, i))
    end)
end

Grid(a::AbstractArray, ranges...) = Grid(a, ranges)

"""
This constructor constructs a grid from ranges given as a tuple.
Due to the approach, the tuple `ranges` can consist of NTuple(2, T)
and all kind of range types. The constructor will make sure that all ranges match
the size of the dimension of the array `a`.
"""
function Grid(a::AbstractArray{T,N}, ranges::Tuple) where {T,N}
    length(ranges) = ! N && throw(ArgumentError(
        "You need to supply a range for every dimension of the array. Given: $ranges
        given Array: $(typeof(a))"
    ))
    Grid(ntuple(Val(N)) do i
        range(first(ranges[i]), stop=last(ranges[i]), length=size(a, i))
    end)
end

Base.length(p::Grid) = prod(size(p))
Base.size(p::Grid) = map(length, p.dims)
function Base.getindex(p::Grid{N,T}, i) where {N,T}
    inds = ind2sub(size(p), i)
    return Point{N,eltype(T)}(ntuple(Val(N)) do i
        p.dims[i][inds[i]]
    end)
end

Base.iterate(g::Grid, i=1) = i <= length(g) ? (g[i], i + 1) : nothing

GLAbstraction.isa_gl_struct(x::Grid) = true
GLAbstraction.toglsltype_string(t::Grid{N,T}) where {N,T} = "uniform Grid$(N)D"
function GLAbstraction.gl_convert_struct(g::Grid{N,T}, uniform_name::Symbol) where {N,T}
    return Dict{Symbol,Any}(
        Symbol("$uniform_name.start") => Vec{N,Float32}(minimum.(g.dims)),
        Symbol("$uniform_name.stop") => Vec{N,Float32}(maximum.(g.dims)),
        Symbol("$uniform_name.lendiv") => Vec{N,Cint}(length.(g.dims) .- 1),
        Symbol("$uniform_name.dims") => Vec{N,Cint}(map(length, g.dims))
    )
end
function GLAbstraction.gl_convert_struct(g::Grid{1,T}, uniform_name::Symbol) where T
    x = g.dims[1]
    return Dict{Symbol,Any}(
        Symbol("$uniform_name.start") => Float32(minimum(x)),
        Symbol("$uniform_name.stop") => Float32(maximum(x)),
        Symbol("$uniform_name.lendiv") => Cint(length(x) - 1),
        Symbol("$uniform_name.dims") => Cint(length(x))
    )
end

struct GLVisualizeShader <: AbstractLazyShader
    paths::Tuple
    kw_args::Dict{Symbol,Any}
    function GLVisualizeShader(paths::String...; view=Dict{String,String}(), kw_args...)
        # TODO properly check what extensions are available
        @static if !Sys.isapple()
            view["GLSL_EXTENSIONS"] = "#extension GL_ARB_conservative_depth: enable"
            view["SUPPORTED_EXTENSIONS"] = "#define DEPTH_LAYOUT"
        end
        args = Dict{Symbol, Any}(kw_args)
        args[:view] = view
        args[:fragdatalocation] = [(0, "fragment_color"), (1, "fragment_groupid")]
        new(map(x -> loadshader(x), paths), args)
    end
end

function assemble_robj(data, program, bb, primitive, pre_fun, post_fun)
    transp = get(data, :transparency, Observable(false))
    overdraw = get(data, :overdraw, Observable(false))
    pre = if !isnothing(pre_fun)
        _pre_fun = GLAbstraction.StandardPrerender(transp, overdraw)
        ()->(_pre_fun(); pre_fun())
    else
        GLAbstraction.StandardPrerender(transp, overdraw)
    end
    robj = RenderObject(data, program, pre, nothing, bb, nothing)
    post = if haskey(data, :instances)
        GLAbstraction.StandardPostrenderInstanced(data[:instances], robj.vertexarray, primitive)
    else
        GLAbstraction.StandardPostrender(robj.vertexarray, primitive)
    end
    robj.postrenderfunction = if !isnothing(post_fun)
        () -> (post(); post_fun())
    else
        post
    end
    robj
end

function assemble_shader(data)
    shader = data[:shader]
    delete!(data, :shader)
    glp = get(data, :gl_primitive, GL_TRIANGLES)
    return assemble_robj(
        data, shader, Rect3f(), glp,
        get(data, :prerender, nothing),
        get(data, :postrender, nothing)
    )
end

"""
Converts index arrays to the OpenGL equivalent.
"""
to_index_buffer(x::GLBuffer) = x
to_index_buffer(x::TOrSignal{Int}) = x
to_index_buffer(x::VecOrSignal{UnitRange{Int}}) = x
to_index_buffer(x::TOrSignal{UnitRange{Int}}) = x
"""
For integers, we transform it to 0 based indices
"""
to_index_buffer(x::AbstractVector{I}) where {I <: Integer} = indexbuffer(Cuint.(x .- 1))
function to_index_buffer(x::Observable{<: AbstractVector{I}}) where I <: Integer
    indexbuffer(lift(x -> Cuint.(x .- 1), x))
end

"""
If already GLuint, we assume its 0 based (bad heuristic, should better be solved with some Index type)
"""
function to_index_buffer(x::VectorTypes{I}) where I <: Union{GLuint,LineFace{GLIndex}}
    indexbuffer(x)
end

to_index_buffer(x) = error(
    "Not a valid index type: $(typeof(x)).
    Please choose from Int, Vector{UnitRange{Int}}, Vector{Int} or a signal of either of them"
)

function visualize(@nospecialize(main), @nospecialize(s), @nospecialize(data))
    data = _default(main, s, copy(data))
    return assemble_shader(data)
end

function output_buffers(transparency = false)
    if transparency
        """
        layout(location=2) out float coverage;
        """
    elseif enable_SSAO[]
        """
        layout(location=2) out vec3 fragment_position;
        layout(location=3) out vec3 fragment_normal_occlusion;
        """
    else
        ""
    end
end

function output_buffer_writes(transparency = false)
    if transparency
        scale = transparency_weight_scale[]
        """
        float weight = color.a * max(0.01, $scale * pow((1 - gl_FragCoord.z), 3));
        coverage = 1.0 - clamp(color.a, 0.0, 1.0);
        fragment_color.rgb = weight * color.rgb;
        fragment_color.a = weight;
        """
    elseif enable_SSAO[]
        """
        fragment_color = color;
        fragment_position = o_view_pos;
        fragment_normal_occlusion.xyz = o_normal;
        """
    else
        "fragment_color = color;"
    end
end

ext: jl


text:
function prs_predict(x::Array{Float64, 2}, y::Array{Float64, 2},
  xx::Array{Float64, 2}; )

# Some radial basis functions
rbf_thin_plate(r::Float64, c::Float64 = 1.0) = (r == 0.0) ? 0.0 : (r^2 * log(c*r))

rbf_gaussian(r::Float64, c::Float64 = 1.0) = exp(-c*r^2)

rbf_cubic(r::Float64, c::Float64 = 1.0) = r^3

rbf_multiquadric(r::Float64, c::Float64 = 1.0) = sqrt(r^2 + c^2)

# Calculate the pairwise kernel matrix for the points in x taken to be the
# individual columns.
function calc_kernel_matrix(x::Array{Float64, 2}, centers::Array{Float64, 2},
  phifunc = rbf_gaussian, distfunc = norm, c = 1.0)
  num_points = size(x, 2)
  num_centers = size(centers, 2)
  k_matrix = zeros(Float64, num_points, num_centers)
  for i in 1:num_points
    xi = x[:,i]
    for j in 1:num_centers
      k_matrix[i, j] = phifunc(norm(xi - centers[:,j]), c)
    end
  end
  k_matrix
end

function rbf_predict(x::Array{Float64, 2}, y::Array{Float64, 2},
  xx::Array{Float64, 2}; c = nothing, rbf = rbf_gaussian)

  # Set c to ok values unless it was specified
  if c == nothing && rbf == rbf_gaussian
    D = size(x, 1)
    c = 1/(D*D)
  elseif c == nothing
    c = 1.0
  end

  # Scale x per row (i.e. covariate)
  x_row_means = mean(x, 2)
  x_row_std = std(x, 2)
  xt = broadcast(-, x, x_row_means)
  xsc = broadcast(/, xt, x_row_std)

  # Scale y per its mean and std
  ymean = mean(y)
  ystd = std(y)
  ysc = (y - ymean) / ystd

  # Dor LSR to get lambdas
  A = calc_kernel_matrix(xsc, xsc, rbf, norm, c)
  lambdas = A\(ysc')

  # Scale xx
  xxt = broadcast(-, xx, x_row_means)
  xxsc = broadcast(/, xxt, x_row_std)

  # Predict for the test points
  B = calc_kernel_matrix(xxsc, xsc, rbf, norm, c)
  ypredsc = (B * lambdas)'

  # Transform back to original range
  ypred = ymean + ystd * ypredsc

  return ypred

end

function r2(ypred, y)
  SSR = sum( (y .- ypred).^2  )
  SST = sum( (y - mean(y)).^2 )
  1.0 - SSR / SST
end

# Test function
tf(x) = sum(x.^4)

N = 100
D = 10
xtraining = randn(D, N)
ytraining = mapslices(tf, xtraining, 1)
xtest = randn(D, 5*D)
y = mapslices(tf, xtest, 1)

ypred_tp = rbf_predict(xtraining, ytraining, xtest; rbf = rbf_thin_plate)
r2_tp = r2(ypred_tp, y)

ypred_multiquadric = rbf_predict(xtraining, ytraining, xtest; rbf = rbf_multiquadric)
r2_multiquadric = r2(ypred_multiquadric, y)

ypred_cubic = rbf_predict(xtraining, ytraining, xtest; rbf = rbf_cubic)
r2_cubic = r2(ypred_cubic, y)

ypred_gaussiand = rbf_predict(xtraining, ytraining, xtest; rbf = rbf_gaussian, c = 1/D)
r2_gaussiand = r2(ypred_gaussiand, y)

ypred_gaussiand2 = rbf_predict(xtraining, ytraining, xtest; rbf = rbf_gaussian, c = 1/(D*D))
r2_gaussiand2 = r2(ypred_gaussiand2, y)


text:
using Documenter, TimeAxisArrays

makedocs(
    modules = [TimeAxisArrays]
)

deploydocs(
    deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo   = "github.com/GordStephen/TimeAxisArrays.jl.git",
    julia  = "release",
    osname = "linux"
)

max_stars_repo_name: GordStephen/TimeAxisArrays

max_issues_repo_name: GordStephen/TimeAxisArrays

max_forks_repo_name: GordStephen/TimeAxisArrays


text:
export iauApci13
"""
For a terrestrial observer, prepare star-independent astrometry
parameters for transformations between ICRS and geocentric CIRS
coordinates.  The caller supplies the date, and SOFA models are used
to predict the Earth ephemeris and CIP/CIO.

The parameters produced by this function are required in the
parallax, light deflection, aberration, and bias-precession-nutation
parts of the astrometric transformation chain.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   date1  double      TDB as a 2-part...
   date2  double      ...Julian Date (Note 1)

Returned:
   astrom iauASTROM*  star-independent astrometry parameters:
      pmt    double       PM time interval (SSB, Julian years)
      eb     double[3]    SSB to observer (vector, au)
      eh     double[3]    Sun to observer (unit vector)
      em     double       distance from Sun to observer (au)
      v      double[3]    barycentric observer velocity (vector, c)
      bm1    double       sqrt(1-|v|^2): reciprocal of Lorenz factor
      bpn    double[3][3] bias-precession-nutation matrix
      along  double       unchanged
      xpl    double       unchanged
      ypl    double       unchanged
      sphi   double       unchanged
      cphi   double       unchanged
      diurab double       unchanged
      eral   double       unchanged
      refa   double       unchanged
      refb   double       unchanged
   eo     double*     equation of the origins (ERA-GST)

Notes:

   1. The TDB date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TDB)=2450123.7 could be expressed in any of these ways, among
   others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in cases
   where the loss of several decimal digits of resolution is
   acceptable.  The J2000 method is best matched to the way the
   argument is handled internally and will deliver the optimum
   resolution.  The MJD method and the date & time methods are both
   good compromises between resolution and convenience.  For most
   applications of this function the choice will not be at all
   critical.

   TT can be used instead of TDB without any significant impact on
   accuracy.

   2. All the vectors are with respect to BCRS axes.

   3. In cases where the caller wishes to supply his own Earth
   ephemeris and CIP/CIO, the function iauApci can be used instead
   of the present function.

   4. This is one of several functions that inserts into the astrom
   structure star-independent parameters needed for the chain of
   astrometric transformations ICRS <-> GCRS <-> CIRS <-> observed.

   The various functions support different classes of observer and
   portions of the transformation chain:

         functions         observer        transformation

      iauApcg iauApcg13    geocentric      ICRS <-> GCRS
      iauApci iauApci13    terrestrial     ICRS <-> CIRS
      iauApco iauApco13    terrestrial     ICRS <-> observed
      iauApcs iauApcs13    space           ICRS <-> GCRS
      iauAper iauAper13    terrestrial     update Earth rotation
      iauApio iauApio13    terrestrial     CIRS <-> observed

   Those with names ending in "13" use contemporary SOFA models to
   compute the various ephemerides.  The others accept ephemerides
   supplied by the caller.

   The transformation from ICRS to GCRS covers space motion,
   parallax, light deflection, and aberration.  From GCRS to CIRS
   comprises frame bias and precession-nutation.  From CIRS to
   observed takes account of Earth rotation, polar motion, diurnal
   aberration and parallax (unless subsumed into the ICRS <-> GCRS
   transformation), and atmospheric refraction.

   5. The context structure astrom produced by this function is used by
   iauAtciq* and iauAticq*.

Called:
   iauEpv00     Earth position and velocity
   iauPnm06a    classical NPB matrix, IAU 2006/2000A
   iauBpn2xy    extract CIP X,Y coordinates from NPB matrix
   iauS06       the CIO locator s, given X,Y, IAU 2006
   iauApci      astrometry parameters, ICRS-CIRS
   iauEors      equation of the origins, given NPB matrix and s

This revision:   2013 October 9

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
function iauApci13(date1::Real, date2::Real)
   # Allocate return value
   ref_astrom = Ref{iauASTROM}(iauASTROM())
   ref_eo     = Ref{Float64}(0.0)

   ccall((:iauApci13, libsofa_c), Cvoid,
           (Cdouble, Cdouble, Ref{iauASTROM}, Ref{Cdouble}),
           convert(Float64, date1),
           convert(Float64, date2),
           ref_astrom, ref_eo)

   return ref_astrom[], ref_eo[]
end


text:
cipher(c) = isdigit(c) ? c : 'a' + ('z' - c)

decode(input) = map(cipher, replace(lowercase(input), r"[\s,.]" => ""))

function encode(input)
    # r"(\w{5})(?=\w)" => s"\1 " adds space after every 5 letters
    replace(decode(input), r"(\w{5})(?=\w)" => s"\1 ")
end

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism


text:
using MultidimensionalSpectroscopy, PyPlot, QuantumOptics, LinearAlgebra, FFTW, Colors,
        Printf, DelimitedFiles

## good source of information
# chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Book%3A_Time_Dependent_Quantum_Mechanics_and_Spectroscopy_(Tokmakoff)/13%3A_Coupling_of_Electronic_and_Nuclear_Motion/13.01%3A_The_Displaced_Harmonic_Oscillator_Model
# www.scm.com/doc/ADF/Input/Vibrationally_resolved_el_spectra.html

# make sure to set script directory as pwd()
cd(@__DIR__)

# let pyplot make popup figures
pygui(true)

# custom colormap (neg-0-pos : purple-blue-white-green-yellow-red)
cmp = create_colormap("bright");

# do you want to calculate 2D spectra ?
calc_2d = true

# Huang-Rhys factor (excited state; w.r.t. ground state)
D_HR = .3
# displacement
d = sqrt(D_HR)

# Huang-Rhys factor (second excited state for ESA)
D_HRf = .01
df = sqrt(D_HRf)

# create Hilbert spaces
b_3ls = NLevelBasis(3)   # Hilbert-space of system  (3 electronic levels)                Basis: {|ui⟩}
b_vib = FockBasis(1)     # Hilbert-space of oscillator (1, 2, ... vibrational states)    Basis: {|vi⟩}
b = b_3ls ⊗ b_vib       # combined basis

# operators of 3 level system
j21 = transition(b_3ls,2,1)                     # |e⟩⟨g|
j12 = dagger(j21)                               # |g⟩⟨e|
j32 = transition(b_3ls,3,2)                     # |e⟩⟨f|
j23 = dagger(j32)                               # |f⟩⟨e|

# operators of vibrational mode
at = create(b_vib)                              # ...
a  = dagger(at)
D  = displace(b_vib,d)                          # displacement operator for excited state
Df = displace(b_vib,df)                         # displacement operator for second excited state

# ground states Hamiltonian
Eg = 0                                          # energy of electronic level
ωg = 0.15                                       # energy of vibration
Hg_el = Eg * j12 * j21                          # electronic part of Hamiltonian        == HG = g * Eg ⊗ dagger(g)
#Hg_vib = ωg/2 * (at*a + one(b_vib))
Hg_vib = ωg * (number(b_vib) + identityoperator(b_vib) * 1/2)   # vibrational part of H

#Hg = Hg_el ⊗ one(b_vib) + one(b_3ls) ⊗ Hg_vib
Hg = Hg_el ⊗ one(b_vib) + j12*j21 ⊗ Hg_vib    # combined H for ground state -- H = |0⟩⟨0| + #TODO fill in
#Hg = Hg_el ⊗ Hg_vib

# #TODO Check: HAVE TO created excited states w/o displaced potential and the use displacement
# operator to act on transition dipole moment ... or calculate FC factors and use
# this as transition dipole operator ...

## excited state Hamiltonian
Ee = 1.65                                       # energy of excited electronic level
ωe = 0.15                                       # energy of vibration
He_el = Ee * j21 * j12                          # electronic part of H
#He_vib = dagger(D) * Hg_vib * D # issue with D and dagger(D)
#He_vib = D * Hg_vib * dagger(D) # issue with D and dagger(D)
He_vib = ωe * (number(b_vib) + identityoperator(b_vib) * 1/2)    #vibrational part of H
#He = He_el ⊗ one(b_vib) + one(b_3ls) ⊗ He_vib
He = He_el ⊗ one(b_vib) + j21*j12 ⊗ He_vib
#He = He_el ⊗ He_vib

## 2nd excited state Hamiltonian
Ef = 2 * Ee + 0.3 #4.45                                                       # energy of second exc. el. level
ωf = 0.15                                                               # energy of vibration
Hf_el = Ef * j32 * j23                                                  # el. part of H
Hf_vib = ωf * (number(b_vib) + identityoperator(b_vib) * 1/2)           # vib. part of H
Hf = Hf_el ⊗ one(b_vib) + j32*j23 ⊗ Hf_vib

## construct system Hamiltonian
H = Hg + He + Hf
#println("Display Hamiltonian")
#display(dense(H).data[1:10,1:4])

## diagonalize system Hamiltonian
states   = eigvecs(dense(H).data)
energies = eigvals(dense(H).data)

# same for ground state
energies_g, states_g = eigenstates(dense(Hg_vib))
# and excited state
energies_e, states_e = eigenstates(dense(He_vib))

## number of eigentstates
N = length(energies)

## plot energy levels
figure(figsize=(6,4))
subplot(121)
# visualize ground state
scatter(zeros(Int8(ceil(N/3))), energies[1:Int8(ceil(N/3))])
# and arbitrary potential energy surface (harmonic)
x = [-1:0.01:2*d;]
m = 100                                         # arbitrary
f = 1/2 * m * ωg^2
plot(x,f*x.^2,"b")

#TODO use transformation into position basis instead
#xx = PositionBasis(0,10,200)
#xxx = position(xx)
#### go to position basis ?? ?
#xpoints = samplepoints(b)

# visualize excited state
m = 100                                          # arbitrary
f = 1/2 * m * ωe^2
scatter(zeros(Int8(floor(N/3))) .+ d, energies[Int8(floor(N/3))+1:end-Int8(floor(N/3))])
plot(x,f * (x.-d).^2 .+ Ee,"y")

# visualize second exc. state
m = 100                                         # arbitrary
f = 1/2 * m * ωf^2
scatter(zeros(Int8(floor(N/3))) .+ df, energies[end-Int8(floor(N/3))+1:end])
plot(x,f * (x.-df).^2 .+ Ef,"g"); ylim(0, 6)

## ????
#for i=1:length(states_g)
#    plot(xpoints, abs2.(states_g[i].data).*40 .+ energies_g[i],color="k",linewidth=1)
#    plot(xpoints, abs2.(states_e[i].data).*40 .+ energies_e[i],color="g",linewidth=1)
#end

## calculate Franck-Condon factors between ground and excited state
FC = complex(zeros(length(energies_g),length(energies_e)))
for ii = 1:length(energies_g)
    for jj = 1:length(energies_e)
        FC[ii,jj] = abs(conj(states_g[1]' * (D*states_e[jj])) * states_g[1]' * (D*states_e[ii]))
        #FC[ii,jj] = abs(conj(states_g[1]' * states_e[jj]) * (states_g[1]' * states_e[ii]))                 # this with D somewhere else ?
    end
end

## ?????? #TODO CHECK!!! the previous algorith corresponds to calculating the FC-factors
l = complex(zeros(length(energies_g),length(energies_e)))
for ii = 1:length(energies_g)
    for jj = 1:length(energies_e)
        l[ii,jj] = exp(-d^2) * d^(2*jj) / factorial(jj)
    end
end

## ??????
ln = zeros(5,)
for n = 1:5
    ln[n] = exp(-d^2) * d^(2*n) / factorial(n)
end

## # TODO what's going on here ?
D.data = FC

## TDM operator between ground and excited state in NLevelBasis only
μ12 = j12 + j21
μ12 = μ12 ⊗ D

## TDM operator between excited and second excited state in NLevelBasis only
μ23 = j23 + j32
μ23 = μ23 ⊗ Df

## initialize system in the ground state
Psi0 = nlevelstate(b_3ls,1) ⊗ fockstate(b_vib,0)
rho0 = Psi0 ⊗ dagger(Psi0)                                                                  # alt.: rho0 = dm(Psi0)

## TODO: used thermally populated ground state instead
#T = 0.01                                                               # Temperature T in what unit ?
#rho0 = thermalstate(H,T)

## Lindblad operators with rates Γ
Γ = [sqrt(0.05), sqrt(0.065)]
L = Γ .* [one(b_3ls) ⊗ (a), (j12) ⊗ one(b_vib)]

## time dimension
tlist = [0:0.12:20;]*2*π

## calculate correlation function for ground/excited state transition
corr = timecorrelations.correlation(tlist, rho0, H, L, μ12, μ12)

## could do the same in the following way
#tout, rhot = timeevolution.master(tlist,μ12*rho0,H,L)
#corr = expect(μ23,rhot)

## plot correlation function
subplot(222)
plot(tlist,real.(corr))

## zeropad to smooth data
zp = 11
corr    = zeropad(corr,zp)
tnew, ~ = interpt(tlist,zp)

#TODO : correctly get absorption and emission! (with correct frequency sign)
## get absorption spectrum
ω_abs, spec_abs = timecorrelations.correlation2spectrum(tnew, corr;
                                                        normalize_spec=true)

## get emission spectrum (just sign switch so far ... calc relaxation in excited state ?)
corr = []
corr = timecorrelations.correlation(tlist, rho0, H, L, μ12, μ12)
ω_em, spec_em = timecorrelations.correlation2spectrum(tlist, corr;
                                                        normalize_spec=true)

## plot absorption spectrum
subplot(224)
plot(ω_abs,spec_abs)

## plot stick spectrum of Franck-Condon factors
temp = FC
temp = temp ./ maximum(real(temp))
for jj = 1:length(energies_g)
    plot(-([Ee-Eg, Ee-Eg] .+ 2*(jj-1)*ωg), [0, real(temp[jj,jj])], color="r",linewidth=1)
end

tight_layout()

## calculate 2D spectra
if calc_2d
    F = L                                           # works with method "lindblad"
    zp = 10                                         # zeropad data up to 10^zp values

    ## calculate 2d spectra at
    T = [0.]

    ## initialize output array/structure
    spectra2d = Array{out2d}(undef, length(T))

    ## use multithreading, if you can:
    Threads.@threads for i = 1:length(T)
    # ... or else:
    #for i = 1:length(T)
        spectra2d[i] = make2Dspectra(tlist,rho0,H,F,μ12,μ23,T[i],
                                            "lindblad";debug=true,use_sub=false,
                                                zp=zp, t2coh=false);
    end

    ## crop 2D data and increase dw (step>1) if necessary
    spectra2d = [crop2d(spectra2d[i],1.4;w_max=2.2,step=1) for i = 1:length(T)]

end

if calc_2d

    ## plot 2D spectra for each(?) T
    # decide what to plot
    rep  = "absorptive"                             # "absorptive", "dispersive", "absolute"
    scal = "lin"                                    # "lin", "tan", "cube", "asinh"

    ## make subplot layout
    nplots = length(T);                             # total number of 2D spectra
    ncols  = Int32(ceil(sqrt(nplots)));
    nrows  = Int32(ceil(nplots / ncols));

    ## determine maximum value in dataset out2d[:].full2d[:,:]
    maxi = maximum([maximum(real(spectra2d[i].full2d)) for i in 1:length(spectra2d)])

    ## plot 2D spectra
    fig, ax = subplots(nrows,ncols,sharex=true,sharey=true,figsize=(ncols*3.2,nrows*3))
    if nplots == 1
        ax = [ax]
    end
    fig.suptitle(rep * " 2D spectrum (" * scal * ". scaling)")
    k = 0
    for i = 1:ncols
        for j = 1:nrows
            global k += 1
            if k > nplots
                continue
            end
            sca(ax[i,j])
            ax[i,j].set_aspect="equal"
            plot2d(spectra2d[k].ω,round.(spectra2d[k].full2d,digits=1);repr=rep,scaling=scal,norm=maxi)
            title("2D spectrum at $(T[k]) fs")
        end
    end
    tight_layout()
    subplots_adjust(top=0.85)

    ## #TODO: remove the following of make more universal
    ω = spectra2d[1].ω

    #plot([Ee-Eg, Ee-Eg] .+ 0*ωg,[ω[1], ω[end]],"k--",linewidth=0.5)
    #plot([Ee-Eg, Ee-Eg] .+ 2*ωg,[ω[1], ω[end]],"k--",linewidth=0.5)
    #plot([Ee-Eg, Ee-Eg] .+ 4*ωg,[ω[1], ω[end]],"k--",linewidth=0.5)
    #plot([ω[1], ω[end]],[Ee-Eg, Ee-Eg] .+ 0*ωg,"k--",linewidth=0.5)
    #plot([ω[1], ω[end]],[Ee-Eg, Ee-Eg] .+ 2*ωg,"k--",linewidth=0.5)
    #plot([ω[1], ω[end]],[Ee-Eg, Ee-Eg] .+ 4*ωg,"k--",linewidth=0.5)

end

max_stars_repo_name: lbmnky/CMDS

max_issues_repo_name: lbmnky/CMDS

max_forks_repo_name: lbmnky/CMDS


text:
function gen_exponential_multistep_starting_values(psi::WaveFunction, dt::Number, N::Int;
          final_iteration=false, combine_first::Bool=false)
    acc = wave_function(psi.m)
    psi0 = wave_function(psi.m)
    rhs_back = WaveFunction[wave_function(psi.m) for j=1:N]
    gen_rhs!(rhs_back[1], psi)
    t0 = get_time(psi)
    copy!(psi0, psi)
    for K=2:N
        for R=1:(final_iteration&&K==N?2:1)
            copy!(psi, psi0)
            for J=2:K
                V  = Matrix{Float64}(inv(Rational{Int}[k^m//factorial(m) for k=-J+2:K-J, m=0:K-2]))
                if combine_first
                    copy!(acc, rhs_back[J-1])
                    add_apply_A!(psi, acc, 1.0)
                    add_phi_A!(acc, psi, dt, 1, dt)
                    set_time!(psi, get_time(psi)+dt)
                else
                    propagate_A!(psi, dt)
                end
                for m=(combine_first?2:1):K-1
                    set!(acc, 0.0)
                    for k=1:K-1
                        axpy!(acc, rhs_back[k], V[m,k])
                    end
                    add_phi_A!(acc, psi, dt, m, dt)
                end
                gen_rhs!(rhs_back[J], psi)
            end
        end
    end
    rhs_back
end

type ExponentialMultistepEquidistantTimeStepperIterator
    psi::WaveFunction
    t0::Real
    tend::Real
    dt::Real
    steps::Int
    V1::Matrix{Float64}
    V2::Matrix{Float64}
    N::Int
    N1::Int
    N2::Int
    iter::Int
    ptr::Int
    rhs_back::Vector{WaveFunction} # storage for previous solution values
    psi0::WaveFunction
    acc::WaveFunction
    combine_first::Bool
end

function exponential_multistep_equidistant_time_stepper(psi::WaveFunction,
         t0::Real, tend::Real, dt::Real, N1::Int;
         iter::Int=0, N2::Int=(iter>0?N1+1:0), steps::Int=-1,
         psi_back=nothing,  final_iteration=false, combine_first::Bool=true)
    m = psi.m
    set_propagate_time_together_with_A!(m, true)
    set_time!(psi, t0)
    N = max(N1, N2)
    if psi_back!=nothing
        rhs_back = psi_back
        @assert length(psi_back)==N
        for k=1:N
            gen_rhs!(rhs, rhs_back[k])
            copy!(rhs_back[k], rhs)
        end
    else
        rhs_back = gen_exponential_multistep_starting_values(psi, dt, N,
                       final_iteration=final_iteration, combine_first=combine_first)
    end
    V1 = Matrix{Float64}(inv(Rational{Int}[n^m//factorial(m) for n=-N1+1:0, m=0:N1-1]))
    V2 = Matrix{Float64}(inv(Rational{Int}[n^m//factorial(m) for n=-N2+2:1, m=0:N2-1]))

    ptr = N
    acc = wave_function(m)
    psi0 = wave_function(m)
    ExponentialMultistepEquidistantTimeStepperIterator(psi, t0, tend, dt, steps, V1, V2,
                 N, N1, N2, iter, ptr, rhs_back, psi0, acc, combine_first)
end

Base.start(tsi::ExponentialMultistepEquidistantTimeStepperIterator) = (tsi.steps<0 ?
    tsi.t0 + (tsi.N-1)*tsi.dt : tsi.N-1)

Base.done(tsi::ExponentialMultistepEquidistantTimeStepperIterator, t) = (tsi.steps<0 ? t >= tsi.tend : t>=tsi.steps )

function Base.next(tsi::ExponentialMultistepEquidistantTimeStepperIterator, t)
    if tsi.iter>=1
        copy!(tsi.psi0, tsi.psi)
    end

    #predictor
    if tsi.combine_first
        k1 = mod(tsi.ptr-1, tsi.N)+1
        copy!(tsi.acc, tsi.rhs_back[k1])
        add_apply_A!(tsi.psi, tsi.acc, 1.0)
        add_phi_A!(tsi.acc, tsi.psi, tsi.dt, 1, tsi.dt)
        set_time!(tsi.psi, get_time(tsi.psi)+tsi.dt)
    else
        propagate_A!(tsi.psi, tsi.dt)
    end
    for m=(tsi.combine_first?2:1):tsi.N1
        set!(tsi.acc, 0.0)
        for k=1:tsi.N1
            k1 = mod(k-tsi.N1+tsi.ptr-1, tsi.N)+1
            axpy!(tsi.acc, tsi.rhs_back[k1], tsi.V1[m,k])
        end
        add_phi_A!(tsi.acc, tsi.psi, tsi.dt, m, tsi.dt)
    end

    tsi.ptr = mod(tsi.ptr, tsi.N) + 1
    gen_rhs!(tsi.rhs_back[tsi.ptr], tsi.psi)

    for iter=1:tsi.iter
        copy!(tsi.psi, tsi.psi0)

        #corrector
        if tsi.combine_first
            k1 = mod(tsi.ptr-2, tsi.N)+1
            copy!(tsi.acc, tsi.rhs_back[k1])
            add_apply_A!(tsi.psi, tsi.acc, 1.0)
            add_phi_A!(tsi.acc, tsi.psi, tsi.dt, 1, tsi.dt)
            set_time!(tsi.psi, get_time(tsi.psi)+tsi.dt)
        else
            propagate_A!(tsi.psi, tsi.dt)
        end
        for m=(tsi.combine_first?2:1):tsi.N2
            set!(tsi.acc, 0.0)
            for k=1:tsi.N2
                k1 = mod(k-tsi.N2+tsi.ptr-1, tsi.N)+1
                axpy!(tsi.acc, tsi.rhs_back[k1], tsi.V2[m,k])
            end
            add_phi_A!(tsi.acc, tsi.psi, tsi.dt, m, tsi.dt)
        end

        gen_rhs!(tsi.rhs_back[tsi.ptr], tsi.psi)
    end

    if tsi.steps<0
        t1 = t + tsi.dt < tsi.tend ? t + tsi.dt : tsi.tend
    else
        t1 = t+1
    end
    return t1, t1
end

function global_orders(psi::WaveFunction, reference_solution::WaveFunction,
                       t0::Real, tend::Real, dt::Real, N1::Int;
                       iter::Int=0, N2::Int=(iter>0?N1+1:0),
                       rows=8, final_iteration=false, combine_first::Bool=true, version::Int=1)
    @assert psi.m==reference_solution.m
    tab = Array(Float64, rows, 3)

    wf_save_initial_value = clone(psi)
    copy!(wf_save_initial_value, psi)

    steps = Int(floor((tend-t0)/dt))
    dt1 = dt
    err_old = 0.0
    println("             dt         err      p")
    println("-----------------------------------")
    for row=1:rows
        if version==2
            for t in exponential_multistep2_equidistant_time_stepper(psi, t0, tend, dt1, N1, N2=N2, steps=steps,
                         iter=iter, final_iteration=final_iteration)
            end
        else
            for t in exponential_multistep_equidistant_time_stepper(psi, t0, tend, dt1, N1, N2=N2, steps=steps,
                         iter=iter, final_iteration=final_iteration, combine_first=combine_first)
            end
        end
        err = distance(psi, reference_solution)
        if (row==1) then
            @printf("%3i%12.3e%12.3e\n", row, Float64(dt1), Float64(err))
            tab[row,1] = dt1
            tab[row,2] = err
            tab[row,3] = 0
        else
            p = log(err_old/err)/log(2.0);
            @printf("%3i%12.3e%12.3e%7.2f\n", row, Float64(dt1), Float64(err), Float64(p))
            tab[row,1] = dt1
            tab[row,2] = err
            tab[row,3] = p
        end
        err_old = err
        dt1 = 0.5*dt1
        steps = 2*steps
        copy!(psi,wf_save_initial_value)
    end
end

function gen_exponential_multistep2_starting_values(psi::WaveFunction, dt::Number, N::Int;
          final_iteration=false)
    psi0 = wave_function(psi.m)
    rhs_back = WaveFunction[wave_function(psi.m) for j=1:N]
    gen_rhs!(rhs_back[1], psi)
    t0 = get_time(psi)
    copy!(psi0, psi)
    for K=2:N
        for R=1:(final_iteration&&K==N?2:1)
            copy!(psi, psi0)
            for k=1:K-1
                propagate_A!(rhs_back[k], -(K-2)*dt)
            end
            for J=2:K
                C = Vector{Float64}(Rational{Int}[n^m for m=0:K-2, n=-J+2:K-J]\Rational{Int}[1//(m+1) for m=0:K-2])
                for k=1:K-1
                   axpy!(psi, rhs_back[k], C[k]*dt)
                end
                propagate_A!(psi, dt)
                gen_rhs!(rhs_back[J], psi)
                propagate_A!(rhs_back[J], -dt)
                for k=1:K
                    propagate_A!(rhs_back[k], dt)
                end
            end

        end
    end
    rhs_back
end

type ExponentialMultistep2EquidistantTimeStepperIterator
    psi::WaveFunction
    t0::Real
    tend::Real
    dt::Real
    steps::Int
    C1::Vector{Float64}
    C2::Vector{Float64}
    N::Int
    N1::Int
    N2::Int
    iter::Int
    ptr::Int
    rhs_back::Vector{WaveFunction} # storage for previous solution values
    psi0::WaveFunction
end

function exponential_multistep2_equidistant_time_stepper(psi::WaveFunction,
         t0::Real, tend::Real, dt::Real, N1::Int;
         iter::Int=0, N2::Int=(iter>0?N1+1:0), steps::Int=-1,
         psi_back=nothing,  final_iteration=false)
    m = psi.m
    set_propagate_time_together_with_A!(m, true)
    set_time!(psi, t0)
    N = max(N1, N2)
    if psi_back!=nothing
        rhs_back = psi_back
        @assert length(psi_back)==N
        for k=1:N
            gen_rhs!(rhs, rhs_back[k])
            copy!(rhs_back[k], rhs)
        end
        for k=0:N-1
            propagate_A!(rhs_back[N-k], dt*k)
        end
    else
        rhs_back = gen_exponential_multistep2_starting_values(psi, dt, N,
                       final_iteration=final_iteration)
        #rhs_back = gen_exponential_multistep_starting_values(psi, dt, N,
        #               final_iteration=final_iteration)
        #for k=0:N-1
        #    propagate_A!(rhs_back[N-k], dt*k)
        #end
    end
    C1 = Vector{Float64}(Rational{Int}[n^m for m=0:N1-1, n=-N1+1:0]\Rational{Int}[1//(m+1) for m=0:N1-1])
    C2 = Vector{Float64}(Rational{Int}[n^m for m=0:N2-1, n=-N2+2:1]\Rational{Int}[1//(m+1) for m=0:N2-1])
    ptr = N
    psi0 = wave_function(m)
    ExponentialMultistep2EquidistantTimeStepperIterator(psi, t0, tend, dt, steps, C1, C2,
                 N, N1, N2, iter, ptr, rhs_back, psi0)
end

Base.start(tsi::ExponentialMultistep2EquidistantTimeStepperIterator) = (tsi.steps<0 ?
    tsi.t0 + (tsi.N-1)*tsi.dt : tsi.N-1)

Base.done(tsi::ExponentialMultistep2EquidistantTimeStepperIterator, t) = (tsi.steps<0 ? t >= tsi.tend : t>=tsi.steps )

function Base.next(tsi::ExponentialMultistep2EquidistantTimeStepperIterator, t)
    if tsi.iter>=1
        copy!(tsi.psi0, tsi.psi)
    end

    #predictor
    for k=1:tsi.N1
        k1 = mod(k-tsi.N1+tsi.ptr-1, tsi.N)+1
        axpy!(tsi.psi, tsi.rhs_back[k1], tsi.C1[k]*tsi.dt)
    end
    propagate_A!(tsi.psi, tsi.dt)

    tsi.ptr = mod(tsi.ptr, tsi.N) + 1
    gen_rhs!(tsi.rhs_back[tsi.ptr], tsi.psi)
    propagate_A!(tsi.rhs_back[tsi.ptr], -tsi.dt)

    for iter=1:tsi.iter
        copy!(tsi.psi, tsi.psi0)

        #corrector
        for k=1:tsi.N2
            k1 = mod(k-tsi.N2+tsi.ptr-1, tsi.N)+1
            axpy!(tsi.psi, tsi.rhs_back[k1], tsi.C2[k]*tsi.dt)
        end
        propagate_A!(tsi.psi, tsi.dt)

        gen_rhs!(tsi.rhs_back[tsi.ptr], tsi.psi)
        propagate_A!(tsi.rhs_back[tsi.ptr], -tsi.dt)
    end

    for k=1:tsi.N
        propagate_A!(tsi.rhs_back[k], tsi.dt)
    end

    if tsi.steps<0
        t1 = t + tsi.dt < tsi.tend ? t + tsi.dt : tsi.tend
    else
        t1 = t+1
    end
    return t1, t1
end


text:
"""
    parse_file(file; import_all)

Parses a Matpower .m `file` or PTI (PSS(R)E-v33) .raw `file` into a
PowerModels data structure. All fields from PTI files will be imported if
`import_all` is true (Default: false).
"""
function parse_file(file::String; import_all=false, validate=true)
    if endswith(file, ".m")
        pm_data =  parse_matpower(file, validate=validate)
    elseif endswith(lowercase(file), ".raw")
        @info("The PSS(R)E parser currently supports buses, loads, shunts, generators, branches, transformers, and dc lines")
        pm_data =  parse_psse(file; import_all=import_all, validate=validate)
    else
        pm_data = parse_json(file, validate=validate)
    end

    # TODO:  not sure if this relevant for all three file types, or only .m, JJS 3/7/19
    move_genfuel_and_gentype!(pm_data)

    return pm_data
end

""
function parse_json(file_string::String; kwargs...)
    open(file_string) do f
        pm_data = parse_json(f, kwargs...)
    end
    return pm_data
end

""
function parse_json(io::IO; validate=true)
    data_string = read(io, String)
    pm_data = JSON.parse(data_string)
    if validate
        check_network_data(pm_data)
    end
    return pm_data
end

"""
Runs various data quality checks on a PowerModels data dictionary.
Applies modifications in some cases.  Reports modified component ids.
"""
function check_network_data(data::Dict{String,Any})
    mod_bus = Dict{Symbol,Set{Int}}()
    mod_gen = Dict{Symbol,Set{Int}}()
    mod_branch = Dict{Symbol,Set{Int}}()
    mod_dcline = Dict{Symbol,Set{Int}}()

    check_conductors(data)
    make_per_unit(data)
    check_connectivity(data)

    mod_branch[:xfer_fix] = check_transformer_parameters(data)
    mod_branch[:vad_bounds] = check_voltage_angle_differences(data)
    mod_branch[:mva_zero] = check_thermal_limits(data)
    mod_branch[:orientation] = check_branch_directions(data)
    check_branch_loops(data)

    mod_dcline[:losses] = check_dcline_limits(data)

    mod_bus[:type] = check_bus_types(data)
    check_voltage_setpoints(data)

    check_storage_parameters(data)

    gen, dcline = check_cost_functions(data)
    mod_gen[:cost_pwl] = gen
    mod_dcline[:cost_pwl] = dcline

    simplify_cost_terms(data)

    return Dict(
        "bus" => mod_bus,
        "gen" => mod_gen,
        "branch" => mod_branch,
        "dcline" => mod_dcline
    )
end

#=
function row_to_typed_dict(row_data, columns)
    @info("call to depreciated function  row_to_typed_dict, use row_to_typed_dict")
    return row_to_typed_dict(row_data, columns)
end

function row_to_dict(row_data, columns)
    @info("call to depreciated function  row_to_dict, use row_to_dict")
    return row_to_dict(row_data, columns)
end
=#


text:
module TestInsert
using Test
using StaticArrays
using Accessors
using Accessors: insert

@testset "test insert" begin
    @testset "function" begin
        @test @inferred(insert( (b=2, c=3), @optic(_.a), 1 )) == (b=2, c=3, a=1)
        let A = [1, 2]
            @test insert(A, @optic(_[2]), 3) == [1, 3, 2]
            @test_throws BoundsError insert(A, @optic(_[4]), 3)
            @test_throws Exception insert(A, @optic(_[1, 3]), 3)
            @test insert(A, first, 3) == [3, 1, 2]
            @test A == [1, 2]  # not changed
        end
        @test insert((1,2), last, 3) == (1, 2, 3)
        Base.thisminor(VERSION) == v"1.6" && @inferred(insert((1,2), last, 3))
        @test @inferred(insert(SVector(1,2), @optic(_[1]), 3)) == SVector(3, 1, 2)
        @test @inferred(insert(SVector(1,2), last, 3)) == SVector(1, 2, 3)
        let D = Dict(:a => 1)
            @test insert(D, @optic(_[:b]), 2) == Dict(:a => 1, :b => 2)
            @test D == Dict(:a => 1)  # not changed
        end
        @test insert((a=1, b=(2, 3)), @optic(_.b[2]), "xxx") === (a=1, b=(2, "xxx", 3))
        @test_broken begin
            @inferred insert((a=1, b=(2, 3)), @optic(_.b[2]), "xxx")
            true
        end
    end

    @testset "macro" begin
        x = (b=2, c=3)
        @test @insert(x.a = 1) == (b=2, c=3, a=1)
        x = [1, 2]
        @test @insert(x[3] = 3) == [1, 2, 3]
        x = (a=(b=(1, 2),), c=1)
        @test @insert(x.a.b[1] = 0) == (a=(b=(0, 1, 2),), c=1)
    end
end

end


text:
## Substitutions and numeric conversions

## subs
"""

`subs` is used to subsitute a value in an expression with another
value.

Examples:

```
x,y = symbols("x,y")
ex = (x-y)*(x+2y)
subs(ex, (y, y^2))
subs(ex, (x,1), (y,2))
subs(ex, (x,y^3), (y,2))
subs(ex, y, 3)
```

There is a curried form of `subs` to use with the chaining `|>` operator

```
ex |> subs(x,e)
```

Since version 0.4, the use of pairs gives a convenient alternative:

```
subs(ex, x=>1, y=>2)
ex |> subs(x=>1, y=>2)
```

There were some older convenience forms, but these will be deprecated, as they don't work as expected when a variable has assumptions. These are:

```
subs(ex, :y, pi)
subs(ex, x=1, y=pi)
## or
ex |> subs(:x, e)
ex |> subs(x=e)
## or as of version 0.4:
ex(x=2, y=3)
```

The `replace` function is related, but not identical to subs.

"""
subs{T <: SymbolicObject}(ex::T, y::@compat(Tuple{SymbolicTypes, Any})) =
    object_meth(ex, :subs, Sym(y[1]), convert(Sym,y[2]))
subs{T <: SymbolicObject}(ex::T, y::@compat(Tuple{SymbolicTypes, Any}), args...) = subs(subs(ex, y), args...)
subs{T <: SymbolicObject, S<:SymbolicObject}(ex::T, y::S, val) = subs(ex, (y,val))
subs{T <: SymbolicObject}(ex::T, dict::Dict) = subs(ex, dict...)
subs{T <: SymbolicObject}(ex::T, d::Pair...) = subs(ex, [(p.first, p.second) for p in d]...)

function subs{T <: SymbolicObject, S <: Symbol}(ex::T, y::S, val)
    warn("Calling subs with a symbol and not a symbolic variable is deprecated")
    subs(ex, (y,val))
end

## curried versions to use with |>
subs(x::SymbolicObject, y) = ex -> subs(ex, x, y)
subs(;kwargs...) = ex -> subs(ex; kwargs...)
subs(dict::Dict) = ex -> subs(ex, dict...)
subs(d::Pair...) = ex -> subs(ex, [(p.first, p.second) for p in d]...)

function subs(x::Symbol, y)
    warn("Calling `subs` with a symbol and not a symbolic variable is deprecated")
    ex -> subs(ex, Sym(x), y)
end

## Convenience method for keyword arguments
## Will deprecate this
function subs{T <: SymbolicObject}(ex::T; kwargs...)
    warn("""
Calling `subs` with keyword arguments will be deprecated. From v0.4 onward, the use of pairs, as in
`subs(ex, var1=>val1, var2=>val2)` is suggested.
""")
    subs(ex, kwargs...)
end

# """
#
# As of version 0.4, expressions are callable. The following works
# using pairs, a dictionary or, one can not specify the variables and
# use the order of the free symbols;

# ```
# ex = (x-3)*(y^2
# ex(x=>1, y=>2)
# ex(Dict(x=>1, y=>2))
# ex(1,2)  ## uses order of free_symbols
# ```
#     """

## different conversions

## helper, as :is_rational will find 1.2 rational...
function _is_rational(ex::Sym)
    ex.x[:is_rational] == nothing && return false
    ex.x[:is_rational] && denom(ex).x[:is_integer]
end

## evalf, n, N
## We want to convert to numeric:
## evalf -> keep as sympy expression
## N bring back into julia using convert
## In SymPy N and evalf are the same, so this could be confusing!!!

"""

Convert a `Sym` value to a numeric Julian value.

The `N` function of SymPy is an alias for `evalf`. Within SymPy, either may be used to
find numeric values from symbolic values.

For example, symbolic roots can be computed numerically, even if not
available symbolically, by calling `N` on the values.

Using `SymPy` within `Julia` makes having two such functions useful:

* one to do the equivalent of SymPy's `evalf` function
* one to *also* convert these expressions back into `Julia` objects.

We use `N` to return a `Julia` object and `evalf` to return a symbolic
object.

Examples:
```
x = Sym("x")
p = subs(x, x, pi)
N(p)                            # float version of pi
evalf(p, 60)                    # 60 digits of pi, as a symbolic value
N(p, 60)                        # when a precision is given, "Big" values are returned
r = subs(x,x,1.2)
N(r)                            # float
q = subs(x, x, 1//2)
N(q)                            # 1//2
z = solve(x^2 + 1)[1]           # -ⅈ
N(z)                            # 0 - 1im
evalf(z)

rts = solve(x^5 - x + 1)
[N(r) for r in rts]             # numeric solutions to quintic
```

The `evalf` function is similar, though it leaves the expression as a symbolic object.
This breaks the similarity of N and evalf for sympy users.

Throws a `DomainError` if no conversion is possible, such as when the expression still has symbolic values.

`N` is type unstable.

"""
function N(ex::Sym)
    ## more work than I'd like
    if ex.x[:is_integer] == nothing
        return(N(ex[:evalf]()))
    end
    if ex.x[:is_integer]
        for T in [Int, BigInt]
            try (return(convert(T, ex))) catch e end
        end
    elseif _is_rational(ex)
        return N(numer(ex)) // N(denom(ex))
        ## `convert(Rational, ex)))` fails on `Sym(4//3)`
    elseif ex.x[:is_real]
        for T in [Irrational, Float64] ## BigFloat???
              try (return(convert(T, ex))) catch e end
        end
    elseif ex.x[:is_complex]
        try
            r, i = ex[:re](), ex[:im]()
            r, i = promote(N(r), N(i))
            return(Complex(r, i))
        catch e
        end
    end
    throw(DomainError())
end
N(x::Number) = x
N(m::Array{Sym}) = map(N, m)
"""
`N` can take a precision argument.

When given as an integer greater than 16, we try to match the digits of accuracy using `BigFloat` precision on conversions to floating point.

"""
function N(x::Sym, digits::Int)
    ## check
    digits <= 16 && return(N(x))
    if x.x[:is_integer] == nothing
        out = x[:evalf](digits)
        return( N(out, digits) )
    end

    ex = evalf(x, digits)
    if x.x[:is_integer]
        return(convert(BigInt, x))
    elseif _is_rational(x)
        return N(numer(x)) / N(denom(x))
    elseif x.x[:is_real]
        p = round(Int,log2(10)*digits)

        out = setprecision(p) do
            convert(BigFloat, ex)
        end
        return(out)
    elseif x.x[:is_complex]
        r, i = ex[:re](), ex[:im]()
        u, v = promote(N(r, digits), N(i, digits))
        return(Complex(u, v))
    end

    throw(DomainError())
end

"""

The `evalf` function has keyword possibilities, such as `subs` being given directly through a `Dict(Sym,Any)`.
Unlike `N`, `evalf` returns an object of type `Sym`.

Examples
```
x = Sym("x")
evalf(x, subs=Dict([(x,1/2)]))
```
"""

function evalf(x::Sym, args...; kwargs...)
    x[:evalf](args...; kwargs...)
end


text:
#################### Mathematical Operators ####################
function cummean(x::AbstractArray)
    return mapslices(cummean, x, dims = 1)
end

function cummean(x::AbstractVector)
    y = similar(x, Float64)
    xs = 0.0
    fill!(y, xs)

    c = 0
    for i in 1:length(x)
        if !ismissing(x[i])
            c += 1
            xs += x[i]
        end
        y[i] = xs / c
    end
    return y
end

## Csorgo S and Faraway JJ. The exact and asymptotic distributions of the
## Cramer-von Mises statistic. Journal of the Royal Statistical Society,
## Series B, 58: 221-234, 1996.
function pcramer(q::Real)
    p = 0.0
    for k in 0:3
        c1 = 4.0 * k + 1.0
        c2 = c1^2 / (16.0 * q)
        p += gamma(k + 0.5) / factorial(k) * sqrt(c1) * exp(-c2) * besselk(0.25, c2)
    end
    return p / (pi^1.5 * sqrt(q))
end

# This sorting stack was sourced from https://rosettacode.org/wiki/Natural_sorting#Julia
splitbynum(x) = split(x, r"(?<=\D)(?=\d)|(?<=\d)(?=\D)")
numstringtonum(arr) = [(n = tryparse(Float32, e)) != nothing ? n : e for e in arr]
function natural(x, y)
    xarr = numstringtonum(splitbynum(x))
    yarr = numstringtonum(splitbynum(y))
    for i in 1:min(length(xarr), length(yarr))
        if typeof(xarr[i]) != typeof(yarr[i])
            a = string(xarr[i]); b = string(yarr[i])
        else
             a = xarr[i]; b = yarr[i]
        end
        if a == b
            continue
        else
            return a < b
        end
    end
    return length(xarr) < length(yarr)
end

function _dict2namedtuple(d::Dict)
    t_keys = ntuple(x -> Symbol(collect(keys(d))[x]), length(keys(d)))
    t_vals = ntuple(x -> collect(values(d))[x], length(values(d)))
    return NamedTuple{t_keys}(t_vals)
end

function _namedtuple2dict(n::NamedTuple)
    return Dict(pairs(n))
end

# adapted from https://github.com/JuliaLang/julia/blob/3fdfb6734bb6ed7fc805a484183077dd7924b7c0/base/namedtuple.jl#L222
function merge_union(a::NamedTuple{an}, b::NamedTuple{bn}) where {an,bn}
    if @generated
        names = Base.merge_names(an, bn)
        types = merge_union_types(names, a, b)

        values = map(names) do n
            if Base.sym_in(n, an)
                if Base.sym_in(n, bn)
                    :(union(getfield(a, $(QuoteNode(n))), getfield(b, $(QuoteNode(n)))))
                else
                    :(getfield(a, $(QuoteNode(n))))
                end
            else
                :(getfield(b, $(QuoteNode(n))))
            end
        end

        return :(NamedTuple{$names,$types}(($(values...),)))
    else
        names = Base.merge_names(an, bn)
        types = merge_union_types(names, typeof(a), typeof(b))

        values = map(names) do n
            if Base.sym_in(n, an)
                if Base.sym_in(n, bn)
                    union(getfield(a, n), getfield(b, n))
                else
                    getfield(a, n)
                end
            else
                getfield(b, n)
            end
        end

        return NamedTuple{names,types}(values)
    end
end

# adapted from https://github.com/JuliaLang/julia/blob/3fdfb6734bb6ed7fc805a484183077dd7924b7c0/base/namedtuple.jl#L192
Base.@pure function merge_union_types(names::Tuple{Vararg{Symbol}}, a::Type{<:NamedTuple}, b::Type{<:NamedTuple})
    an = Base._nt_names(a)
    bn = Base._nt_names(b)

    types = map(names) do n
        if Base.sym_in(n, an)
            if Base.sym_in(n, bn)
                promote_type(fieldtype(a, n), fieldtype(b, n))
            else
                fieldtype(a, n)
            end
        else
            fieldtype(b, n)
        end
    end

    return Tuple{types...}
end

# promote element types of the tail of a NamedTuple
function promote_eltype_namedtuple_tail(::NamedTuple{k,v}) where {k,v}
    return promote_eltype_tuple_type(Base.tuple_type_tail(v))
end

# promote element types of a tuple
promote_eltype_tuple_type(::Type{Tuple{}}) = Any
promote_eltype_tuple_type(::Type{Tuple{T}}) where T = T
function promote_eltype_tuple_type(t::Type{<:Tuple})
    Base.promote_eltype(Base.tuple_type_head(t), promote_eltype_tuple_type(Base.tuple_type_tail(t)))
end

"""
    cskip(x)

Wrapper for `collect(skipmissing(x))`.
"""
cskip(x) = collect(skipmissing(x))

# check if `T` and recursive `eltype` of `T` are concrete types
function isconcretetype_recursive(T)
    return isconcretetype(T) && (eltype(T) === T || isconcretetype_recursive(eltype(T)))
end

"""
    concretize(x)

Return a container whose elements are as concrete as possible by, e.g., replacing
element types of `Real` or `Union{Missing,Float64}` with `Float64`, if appropriate.

Note that this function is only type stable if the containers are already concretely
typed. In this case the input `x` is returned, and hence input and output share the
same underlying data. In general, however, the input and the output do not share the same
data.
"""
concretize(x) = x
function concretize(x::AbstractArray)
    if isconcretetype_recursive(typeof(x))
        return x
    else
        xnew = map(concretize, x)
        T = mapreduce(typeof, promote_type, xnew)
        if T <: eltype(xnew)
            return convert(AbstractArray{T}, xnew)
        else
            return xnew
        end
    end
end

function concretize(x::Chains)
    value = x.value
    if isconcretetype_recursive(value)
        return x
    else
        return Chains(concretize(value), x.logevidence, x.name_map, x.info)
    end
end


text:
# BED
# ===
#
# Reader and writer of the BED file format.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

immutable BED <: FileFormat end

"""Metadata for BED interval records"""
type BEDMetadata
    used_fields::Int # how many of the first n fields are used
    name::StringField
    score::Int
    thick_first::Int
    thick_last::Int
    item_rgb::RGB{Float32}
    block_count::Int
    block_sizes::Vector{Int}
    block_firsts::Vector{Int}
end

function BEDMetadata()
    return BEDMetadata(0, StringField(), 0, 0, 0, RGB{Float32}(0.0, 0.0, 0.0),
                       0, Int[], Int[])
end

function Base.copy(metadata::BEDMetadata)
    return BEDMetadata(
        metadata.used_fields, copy(metadata.name),
        metadata.score, metadata.thick_first, metadata.thick_last,
        metadata.item_rgb, metadata.block_count,
        metadata.block_sizes[1:metadata.block_count],
        metadata.block_firsts[1:metadata.block_count])
end

@compat function Base.:(==)(a::BEDMetadata, b::BEDMetadata)
    if a.used_fields != b.used_fields
        return false
    end

    n = a.used_fields
    ans = (n < 1 || a.name == b.name) &&
          (n < 2 || a.score == b.score) &&
          (n < 4 || a.thick_first == b.thick_first) &&
          (n < 5 || a.thick_last == b.thick_last) &&
          (n < 6 || a.item_rgb == b.item_rgb) &&
          (n < 7 || a.block_count == b.block_count)
    if !ans
        return false
    end

    if n >= 8
        for i in 1:a.block_count
            if a.block_sizes[i] != b.block_sizes[i]
                return false
            end
        end
    end

    if n >= 9
        for i in 1:a.block_count
            if a.block_sizes[i] != b.block_sizes[i]
                return false
            end
        end
    end

    return true
end

"An `Interval` with associated metadata from a BED file"
typealias BEDInterval Interval{BEDMetadata}

function Base.open(filepath::AbstractString, mode::AbstractString, ::Type{BED};
                   n_fields::Integer=-1)
    io = open(filepath, mode)
    if mode[1] == 'r'
        return open(BufferedInputStream(io), BED)
    elseif mode[1] ∈ ('w', 'o')
        return BEDWriter(io, n_fields)
    end
    error("invalid open mode")
end

# Parser
# ------

type BEDParser <: AbstractParser
    state::Ragel.State

    # intermediate values used during parsing
    red::Float32
    green::Float32
    blue::Float32
    block_size_idx::Int
    block_first_idx::Int

    function BEDParser(input::BufferedInputStream)
        return new(Ragel.State(bedparser_start, input), 0, 0, 0, 1, 1)
    end
end

function Intervals.metadatatype(::BEDParser)
    return BEDMetadata
end

function Base.eltype(::Type{BEDParser})
    return BEDInterval
end

function Base.eof(parser::BEDParser)
    return eof(parser.state.stream)
end

function Base.open(input::BufferedInputStream, ::Type{BED})
    return BEDParser(input)
end

function IntervalCollection(interval_stream::BEDParser)
    intervals = collect(BEDInterval, interval_stream)
    return IntervalCollection{BEDMetadata}(intervals, true)
end

include("bed-parser.jl")

# Writer
# ------

type BEDWriter{T<:IO} <: AbstractWriter
    output::T

    # number of written additional fields (-1: inferred from the first interval)
    n_fields::Int
end

function Base.flush(writer::BEDWriter)
    # TODO: This can be removed on Julia v0.5
    # (because flush will be defined for IOBuffer).
    if applicable(flush, writer.output)
        flush(writer.output)
    end
end

Base.close(writer::BEDWriter) = close(writer.output)

function Base.write(writer::BEDWriter, interval::BEDInterval)
    if writer.n_fields == -1
        writer.n_fields = interval.metadata.used_fields
    end

    if interval.metadata.used_fields < writer.n_fields
        error("interval doesn't have enough additional fields")
    end

    output = writer.output
    n_fields = writer.n_fields
    n = 0
    n += write(
        output,
        interval.seqname, '\t',
        dec(interval.first - 1), '\t',
        dec(interval.last))

    if n_fields ≥ 1
        n += write(output, '\t', interval.metadata.name)
    end

    if n_fields ≥ 2
        n += write(output, '\t', dec(interval.metadata.score))
    end

    if n_fields ≥ 3
        n += write(output, '\t', Char(interval.strand))
    end

    if n_fields ≥ 4
        n += write(output, '\t', dec(interval.metadata.thick_first - 1))
    end

    if n_fields ≥ 5
        n += write(output, '\t', dec(interval.metadata.thick_last))
    end

    if n_fields ≥ 6
        rgb = interval.metadata.item_rgb
        n += write(
            output, '\t',
            dec(round(Int, 255 * rgb.r)), ',',
            dec(round(Int, 255 * rgb.g)), ',',
            dec(round(Int, 255 * rgb.b)))
    end

    if n_fields ≥ 7
        n += write(output, '\t', dec(interval.metadata.block_count))
    end

    if n_fields ≥ 8
        block_sizes = interval.metadata.block_sizes
        if !isempty(block_sizes)
            n += write(output, '\t', dec(block_sizes[1]))
            for i in 2:length(block_sizes)
                n += write(output, ',', dec(block_sizes[i]))
            end
        end
    end

    if n_fields ≥ 9
        block_firsts = interval.metadata.block_firsts
        if !isempty(block_firsts)
            n += write(output, '\t', dec(block_firsts[1] - 1))
            for i in 2:length(block_firsts)
                n += write(output, ',', dec(block_firsts[i] - 1))
            end
        end
    end

    n += write(output, '\n')
    return n
end

function Base.write(out::IO, interval::BEDInterval)
    print(out, interval.seqname, '\t', interval.first - 1, '\t', interval.last)
    write_optional_fields(out, interval)
    write(out, '\n')
end

function write_optional_fields(out::IO, interval::BEDInterval, leadingtab::Bool=true)
    if interval.metadata.used_fields >= 1
        if leadingtab
            write(out, '\t')
        end
        print(out, interval.metadata.name)
    else return end

    if interval.metadata.used_fields >= 2
        print(out, '\t', interval.metadata.score)
    else return end

    if interval.metadata.used_fields >= 3
        print(out, '\t', interval.strand)
    else return end

    if interval.metadata.used_fields >= 4
        print(out, '\t', interval.metadata.thick_first - 1)
    else return end

    if interval.metadata.used_fields >= 5
        print(out, '\t', interval.metadata.thick_last)
    else return end

    if interval.metadata.used_fields >= 6
        item_rgb = interval.metadata.item_rgb
        print(out, '\t',
              round(Int, 255 * item_rgb.r), ',',
              round(Int, 255 * item_rgb.g), ',',
              round(Int, 255 * item_rgb.b))
    else return end

    if interval.metadata.used_fields >= 7
        print(out, '\t', interval.metadata.block_count)
    else return end

    if interval.metadata.used_fields >= 8
        block_sizes = interval.metadata.block_sizes
        if !isempty(block_sizes)
            print(out, '\t', block_sizes[1])
            for i in 2:length(block_sizes)
                print(out, ',', block_sizes[i])
            end
        end
    else return end

    if interval.metadata.used_fields >= 9
        block_firsts = interval.metadata.block_firsts
        if !isempty(block_firsts)
            print(out, '\t', block_firsts[1] - 1)
            for i in 2:length(block_firsts)
                print(out, ',', block_firsts[i] - 1)
            end
        end
    end
end


text:
"""
Sparse transform (SPRIGHT) decoding main file. Logic flow:

1. Generate a signal from input_signal.jl
2. Subsample from query.jl
3. Peel using reconstruct.jl
"""
module SparseTransforms
    export all_methods, spright, method_test, method_report
    using ProgressMeter
    include("reconstruct.jl")
    export fwht, bin_to_dec, dec_to_bin, binary_ints, sign_spright, expected_bin
    export Signal, TestSignal, InputSignal, LazySignal, get_subsignal, get_random_sparse_signal, get_random_delta_sparse_signal
    export get_D, get_b, get_Ms, subsample_indices, compute_delayed_subtransform
    export singleton_detection, bin_cardinality
    export decode_with

    all_methods = Dict(
        "query" => [:simple, :random],
        "delays" => [:identity_like, :random, :nso],
        "reconstruct" => [:noiseless, :mle, :nso],
        "code" => [:none]
    )

    function spright(signal::Signal, methods_list::Array{Symbol,1}; verbose::Bool=false, report::Bool=false)
        return transform(signal, methods_list, fwht; verbose=verbose, report=report)
    end

    """
    Full SPRIGHT encoding and decoding. Implements Algorithms 1 and 2 from [2].
    (numbers) in the comments indicate equation numbers in [2].

    Arguments
    ---------
    signal : TestSignal object.
    The signal to be transformed / compared to.

    methods_list : Array{Symbol,1}
    The three symbols [query_method, delays_method, reconstruct_method].
    All implemented methods are available in `all_methods`: if you add a new method, make sure to update `all_methods`.

        query_method : Symbol
        The method to generate the sparsity coefficient and the subsampling matrices.
        Currently implemented methods:
            :simple : choose some predetermined matrices based on problem size.

        delays_method : Symbol
        The method to generate the matrix of delays.
        Currently implemented methods:
            :identity_like : return a zero row and the identity matrix vertically stacked together.
            :random : make a random set of delays.
            :nso : set delays according to the NSO-SPRIGHT algorithm.

        reconstruct_method : Symbol
        The method to detect singletons.
        Currently implemented methods:
            :noiseless : decode according to [2], section 4.2, with the assumption the signal is noiseless.
            :mle : naive noisy decoding; decode by taking the maximum-likelihood singleton that could be at that bin.
            :nso : reconstruct according to the NSO-SPRIGHT algorithm.

    transform : Function
    The base transform: either `fwht` or `fft`.

    verbose : Bool
    Whether to print intermediate steps.

    report : Bool
    Whether to report on algorithm performance.

    Returns
    -------
    wht : Array{Float64,1}
    The WHT constructed by subsampling and peeling.
    """
    function transform(signal::Signal, methods_list::Array{Symbol,1}, transform::Function; verbose::Bool=false, report::Bool=false)
        for (method_type, method_name) in zip(["query", "delays", "reconstruct", "code"], methods_list)
            impl_methods = all_methods[method_type]
            @assert method_name in impl_methods "$method_type method must be one of $impl_methods"
        end
        query_method, delays_method, reconstruct_method, code = methods_list
        @assert (reconstruct_method != :so) || (code != :none)
        # check the condition for p_failure > eps
        # upper bound on the number of peeling rounds, error out after that point

        num_peeling = 0
        locs = BitArray[]
        strengths = Float64[]
        wht = Dict{Int64,Float64}()
        b = get_b(signal; method=query_method)
        Ms = get_Ms(signal.n, b; method=query_method)
        peeling_max = 2^b
        N, B = 2^signal.n, 2^b

        Us = Array{Float64}[]
        if report
            used = Set{Int}()
        end

        if delays_method == :nso
            num_delays = signal.n * Int64(ceil(log2(signal.n))) # idk
        elseif delays_method == :so
            linearrithmic = signal.n * Int64(ceil(log2(signal.n)))
            num_delays = [linearrithmic,    # P1
                          signal.n,         # P2
                          linearrithmic]    # P3
        else
            num_delays = signal.n + 1
        end
        D = get_D(signal.n; method=delays_method, num_delays=num_delays, code=code)
        if reconstruct_method == :mle
            K = binary_ints(signal.n)
            S = (-1) .^(D * K)
        end

        # subsample, make the observation [U] matrices
        for M in Ms
            U, used_i = compute_delayed_subtransform(signal, M, D, transform)
            push!(Us, U)
            if report
                used = union(used, used_i)
            end
        end

        cutoff = 4 * signal.noise_sd ^ 2 * (2 ^ (signal.n - b)) * sum(num_delays) # noise threshold

        # K is the binary representation of all integers from 0 to 2 ^ n - 1.
        select_froms = Array{Int64}[]
        for M in Ms
            selects = [0]
            if reconstruct_method == :mle
                selects = @pipe M' * K |> transpose |> Bool.(_) |> eachrow |> bin_to_dec.(_)
            end
            push!(select_froms, selects)
        end
        # `select_froms` is the collection of 'j' values and associated indices
        # so that we can quickly choose from the coefficient locations such that M.T * k = j as in (20)
        # example: ball j goes to bin at "select_froms[i][j]"" in stage i

        # begin peeling
        # index convention for peeling: 'i' goes over all M/U values
        # i.e. it refers to the index of the subsampling group (zero-indexed - off by one from the paper).
        # 'j' goes over all columns of the WHT subsample matrix, going from 0 to 2 ** b - 1.
        # e.g. (i, j) = (0, 2) refers to subsampling group 0, and aliased bin 2 (10 in binary)
        # which in the example of section 3.2 is the multiton X[0110] + X[1010] + W1[10]

        # a multiton will just store the (i, j)s in a list
        # a singleton will map from the (i, j)s to the true (binary) values k.
        # e.g. the singleton (0, 0), which in the example of section 3.2 is X[0100] + W1[00]
        # would be stored as the dictionary entry (0, 0): array([0, 1, 0, 0]).

        multitons_found = true
        iters = 0
        max_iters = 2 ^ (b + 1)
        while multitons_found && (num_peeling < peeling_max) && (iters < max_iters)

            # first step: find all the singletons and multitons.
            singletons = Dict{Tuple{Int64,Int64},Tuple{BitArray,Float64}}() # dictionary from (i, j) values to the true index of the singleton, k.
            multitons = Tuple{Int64,Int64}[] # list of (i, j) values indicating where multitons are.

            for (i, (M, U, select_from)) in enumerate(zip(Ms, Us, select_froms))
                col_gen = U |> eachrow |> enumerate
                for (j, col) in col_gen
                    if col⋅col > cutoff
                        selection = findall(==(j-1), select_from) # pick all the k such that M.T @ k = j
                        slice = 0
                        if reconstruct_method == :mle
                            slice = S[:, selection]
                        end
                        k, ρ, s_k = singleton_detection(
                            col;
                            method=reconstruct_method,
                            selection=selection,
                            S_slice=slice,
                            n=signal.n,
                            num_delays=num_delays,
                            D=D
                        ) # find the best fit singleton
                        residual = col - ρ * s_k
                        if (expected_bin(k, M) != j) || (residual ⋅ residual > cutoff)
                            push!(multitons, (i, j))
                        else # declare as singleton
                            singletons[(i, j)] = (k, ρ)
                        end # if residual norm > cutoff
                    end # if col norm > cutoff
                end # for col
            end # for U, select_from

            # all singletons and multitons are discovered
            if verbose
                println("Singletons:")
                for (ston_key, ston_value) in singletons
                    println(ston_key, bin_to_dec(ston_value[1]))
                end

                println("Multitons : $multitons")
            end

            # WARNING: this is not a correct thing to do
            # in the last iteration of peeling, everything will be singletons and there
            # will be no multitons
            if length(multitons) == 0 # no more multitons, and can construct final transform
                multitons_found = false
            end

            # balls to peel
            balls_to_peel = Set{Int64}()
            ball_values = Dict{Int64,Float64}()
            for (_, (k, ρ)) in singletons
                ball = bin_to_dec(k)
                push!(balls_to_peel, ball)
                ball_values[ball] = ρ
            end

            if verbose
                println("Balls to be peeled")
                println(balls_to_peel)
            end
            # peel
            iters += 1
            for ball in balls_to_peel
                num_peeling += 1
                k = dec_to_bin(ball, signal.n)

                push!(locs, k)
                push!(strengths, ball_values[ball])

                for (l, M) in enumerate(Ms)
                    peel = @pipe M' * k |> _ .% 2 |> Bool.(_) |> bin_to_dec |> _ + 1
                    signature_in_stage = (-1) .^ (D * k)
                    to_subtract = ball_values[ball] * signature_in_stage
                    U = Us[l]
                    U_slice = U[peel, :]
                    Us[l][peel,:] = U_slice - to_subtract
                    if verbose
                        println("Peeled ball $(bin_to_dec(k)) off bin ($l, $peel)")
                    end
                end # for peel
            end # for ball
        end # while

        loc = Set{Int64}()
        for (k, value) in zip(locs, strengths) # iterating over (i, j)s
            idx = bin_to_dec(k) # converting 'k's of singletons to decimals
            push!(loc, idx)
            if !haskey(wht, idx)
                wht[idx] = value / (2 ^ b)
            end
        end

        wht = filter(f -> abs(last(f)) ≥ signal.noise_sd, wht)
        if report
            return wht, length(used)
        else
            return wht
        end
    end
end # module


text:
import REPL.LineEdit

function LineEdit.refresh_line(s)
    LineEdit.refresh_multi_line(s)
    OhMyREPL.Prompt.rewrite_with_ANSI(s)
end


text:
using LightGraphs
using MonteCarloX
using Random
using Distributions
import Distributions.pdf
import Distributions.cdf
using StatsBase
import StatsBase.kldivergence

include("utils.jl")

""" Testing reweighting on 2D Ising model"""
function test_ising_reweighting(;verbose = false)
    log_P(E, beta) = -beta * E
    P(E, beta) = exp(-beta * E)

    list_beta = [0.0,0.3,0.7,1.0,1.5,2.0]
    pass = true

    rng = MersenneTwister(1000)
    system = constructIsing([8,8], rng)
    nearest_neighbors(index) = outneighbors(system.lattice, index)
    N = length(system.spins)

    for beta in list_beta
        samples = 1000
        list_energy = zeros(samples)
        # thermalization 100 sweeps
        for sweep in 1:samples + 100
            # Metropolis:
            for i in 1:N
                update_spin_flip(system, beta, rng)
            end
            if sweep > 100
                list_energy[sweep - 100] = E(system)
            end
        end
        H_meas = fit(Histogram, list_energy, minimum(list_energy):4:maximum(list_energy) + 4, closed = :left)
        P_meas = normalize(H_meas, mode = :probability)
        # reweighting both with timeseries (list_energy) and histrogram
        beta_source = beta
        beta_target = beta + 0.05
        log_P_source(E) = log_P(E, beta_source)
        log_P_target(E) = log_P(E, beta_target)
        P_source(E) = P(E, beta_source)
        P_target(E) = P(E, beta_target)
        E_ref = analytic_expectation_value_E(beta_target, log_dos_beale_8x8)
        P_ref_source = BoltzmannDistribution(beta_source, log_dos_beale_8x8).pdf
        P_ref_target = BoltzmannDistribution(beta_target, log_dos_beale_8x8).pdf
        kld_ref = kldivergence(P_meas, x->P_ref_source[x])
        if verbose
            println("result analytic = $(E_ref)")
        end
        ###########################################################################
        E_reweight_list1 = MonteCarloX.expectation_value_from_timeseries(log_P_target, log_P_source, list_energy, list_energy)
        E_reweight_list1_error = abs((E_reweight_list1 - E_ref) / E_ref)
        if verbose
            println("result expectation_value_from_timeseries_log: <E> = $(E_reweight_list1); difference from exact = $(E_reweight_list1_error)")
        end
        pass &= E_reweight_list1_error < 0.1
        ###########################################################################
        P_reweight_list = MonteCarloX.distribution_from_timeseries(log_P_target, log_P_source, list_energy, minimum(list_energy):4:maximum(list_energy) + 4)
        kld_source = kldivergence(P_reweight_list, x->P_ref_source[x])
        kld_target = kldivergence(P_reweight_list, x->P_ref_target[x])
        if verbose
            println("result distribution_from_timeseries_log: kld_target ($(kld_target)) !< kld_source ($(kld_source))")
        end
        # only compare reweighted distribution to target and source. There result is comparable.
        pass &= kld_target < kld_source
        ###########################################################################
        E_reweight_hist1 = MonteCarloX.expectation_value_from_histogram(E->E, log_P_target, log_P_source, H_meas)
        E_reweight_hist1_error = abs((E_reweight_hist1 - E_ref) / E_ref)
        if verbose
            println("result expectation_value_from_histogram_log-1: <E> = $(E_reweight_hist1); difference from exact = $(E_reweight_hist1_error)")
            println("result expectation_value_from_histogram_log-1: <E> = $(E_reweight_hist1); difference from timeseries = $(E_reweight_hist1 - E_reweight_list1)")
        end
        pass &= E_reweight_hist1_error < 0.1
        ###########################################################################
        hist_obs = zero(H_meas)
        for E in list_energy
            index = StatsBase.binindex(hist_obs, E)
            hist_obs.weights[index] += E
        end
        E_reweight_hist2 = MonteCarloX.expectation_value_from_histogram(log_P_target, log_P_source, H_meas, hist_obs)
        E_reweight_hist2_error = abs(E_reweight_hist2 - E_reweight_list1)
        if verbose
            println("result expectation_value_from_histogram_log: <E> = $(E_reweight_hist2); difference from timeseries = $(E_reweight_hist2_error)")
        end
        pass &= E_reweight_hist2_error < 0.1
        ###########################################################################
        # TODO: canonical reweighting
    end

    # TODO: multi-histogram reweighting

    return pass
end

""" Testing Metropolis accept on 2D Ising model"""
function test_ising_metropolis(;verbose = false)
    list_beta = [0.0,0.3,0.7,1.0,1.5,2.0]
    pass = true

    rng = MersenneTwister(1000)
    system = constructIsing([8,8], rng)
    nearest_neighbors(index) = outneighbors(system.lattice, index)
    N = length(system.spins)

    for beta in list_beta
        samples = 1000
        list_energy = zeros(samples)
        # thermalization 100 sweeps
        for sweep in 1:samples + 100
            # Metropolis:
            for i in 1:N
                update_spin_flip(system, beta, rng)
            end
            if sweep > 100
                list_energy[sweep - 100] = E(system)
            end
        end
        P_meas = fit(Histogram, list_energy, minimum(list_energy):4:maximum(list_energy) + 4, closed = :left)
        P_meas = normalize(P_meas, mode = :probability)
        P_true = BoltzmannDistribution(beta, log_dos_beale_8x8).pdf
        # KL divergence sum P(E)logP(E)/Q(E)
        # P=P_meas, Q=P_true s.t P(E)=0 simply ignored
        kld = abs(kldivergence(P_meas, x->P_true[x]))
        if verbose
            println(beta, " ", kld, " ", length(keys(P_meas)))
        end
        pass &= kld < 0.1
    end

    return pass
end

""" Testing Cluster update on 2D Ising model"""
function test_ising_cluster(;verbose = false)
    # avoid beta_c = 0.44
    # why does it consistently fail around beta=0.6
    list_beta = [0,0.3,0.5,1,1.5,2.0]
    pass = true

    rng = MersenneTwister(1000)
    system = constructIsing([8,8], rng)
    nearest_neighbors(index) = outneighbors(system.lattice, index)
    N = length(system.spins)

    cluster_algorithm = MonteCarloX.ClusterWolff()

    for beta in list_beta
        samples = Int(1e4)
        list_energy = zeros(samples)

        # thermalization 100 sweeps
        for sweep in 1:samples + 100
            for i in 1:5
                MonteCarloX.update(cluster_algorithm, rng, system.spins, nearest_neighbors, beta)
            end
            if sweep > 100
                list_energy[sweep - 100] = E(system)
            end
        end

        #
        P_meas = fit(Histogram, list_energy, minimum(list_energy):4:maximum(list_energy) + 4, closed = :left)
        P_meas = normalize(P_meas, mode = :probability)
        P_true = BoltzmannDistribution(beta, log_dos_beale_8x8).pdf

        kld = abs(kldivergence(P_meas, x->P_true[x]))
        if verbose
            println("result kldivergence P_meas to P_true = $(kld)")
        end
        pass &= kld < 0.02
    end

    return pass
end

###############################################################################
###############################################################################
### Details functions

struct IsingSystem{F}
    dims::Vector{Int}
    lattice::SimpleGraph
    spins::Vector{Int}
    nearest_neighbors::F
end

function constructIsing(dims::Vector{Int}, rng::AbstractRNG)::IsingSystem
    lattice = LightGraphs.SimpleGraphs.grid(dims, periodic = true)
    spins   = rand(rng, [-1,1], nv(lattice))
    nearest_neighbors(i) = outneighbors(lattice, i)
    return IsingSystem(dims, lattice, spins, nearest_neighbors)
end

function M(system::IsingSystem)::Int
    return abs(sum(system.spins))
end

function analytic_expectation_value_E(beta::Float64, log_dos::Vector{Tuple{Int,Float64}})::Float64
    mean_E  = 0.0
    norm    = 0.0
    for (E, log_d) in log_dos
        mean_E  +=     E * exp(log_d - beta * E)
        norm    +=         exp(log_d - beta * E)
    end
    mean_E /= norm
    return mean_E
end

# J=1
# E=-J*sum s_i s_j
function E(system::IsingSystem)::Int
    E = 0
    for i in 1:length(system.spins)
        E += E_local(system, i)
    end
    return E / 2
end

function E_local(system::IsingSystem, index::Int)::Int
    e = 0
    for j in system.nearest_neighbors(index)
        e -= system.spins[index] * system.spins[j]
    end
    return e
end

function update_spin_flip(system::IsingSystem, beta::Float64, rng::AbstractRNG)::Int
    index = rand(rng, 1:length(system.spins))
    dE        = -2 * E_local(system, index)
    if accept(Metropolis(), rng, beta, dE)
        system.spins[index] *= -1
    else
        dE = 0
    end
    return dE
end

# """ Definition of the Boltzmann distribution for hypothesis testing."""
struct BoltzmannDistribution <: ContinuousUnivariateDistribution
    pdf::Dict{Any,Float64}
    #cdf::Dict{Any,Float64}

    BoltzmannDistribution(beta, log_dos) = new(initialize_BoltzmannDistribution(beta, log_dos))
    # BoltzmannDistribution(beta, log_dos, log_Z) = new(Float64(beta), Dict{Any,Float64}(log_dos), Float64(log_Z))
end

# Probability distribution (for both precompute in constructor...)
pdf(d::BoltzmannDistribution, E::Real) = d.pdf[E]

# cumulated probability distribution
#cdf(d::BoltzmannDistribution, E::Real) = d.cdf[E]

function initialize_BoltzmannDistribution(beta, log_dos)
    # partition sum
    log_Z = -Inf
    for (E, log_d) in log_dos
        log_Z = log_sum(log_Z, log_d - beta * E)
    end

    pdf = Dict{Int64,Float64}()
    cdf = Dict{Int64,Float64}()

    log_cdf = 0.0
    for (E, log_d) in log_dos
        log_pdf = log_d - beta * E - log_Z
        log_cdf = log_sum(log_cdf, log_pdf)
        pdf[E] = exp(log_pdf)
        cdf[E] = exp(log_cdf)
    end
    #return pdf, cdf
    return pdf
end

log_dos_beale_8x8 = [ (-128, 0.6931471805599453), (-124, 0.0), (-120, 4.852030263919617), (-116, 5.545177444479562), (-112, 8.449342524508063), (-108, 9.793672686528922), (-104, 11.887298863200714), (-100, 13.477180596840947), (-96, 15.268195474147658), (-92, 16.912371686315282), (-88, 18.59085846191256), (-84, 20.230089202801466), (-80, 21.870810400320693), (-76, 23.498562234123614), (-72, 25.114602234581373), (-68, 26.70699035290573), (-64, 28.266152815389898), (-60, 29.780704423363996), (-56, 31.241053997806176), (-52, 32.63856452513369), (-48, 33.96613536105969), (-44, 35.217576663643314), (-40, 36.3873411250109), (-36, 37.47007844691906), (-32, 38.46041522581422), (-28, 39.35282710786369), (-24, 40.141667825183845), (-20, 40.82130289691285), (-16, 41.38631975325592), (-12, 41.831753810069756), (-8, 42.153328313883975), (-4, 42.34770636939425), (0, 42.41274640460084), (4, 42.34770636939425), (8, 42.153328313883975), (12, 41.831753810069756), (16, 41.38631975325592), (20, 40.82130289691285), (24, 40.141667825183845), (28, 39.35282710786369), (32, 38.46041522581422), (36, 37.47007844691906), (40, 36.3873411250109), (44, 35.217576663643314), (48, 33.96613536105969), (52, 32.63856452513369), (56, 31.241053997806176), (60, 29.780704423363996), (64, 28.266152815389898), (68, 26.70699035290573), (72, 25.114602234581373), (76, 23.498562234123614), (80, 21.870810400320693), (84, 20.230089202801466), (88, 18.59085846191256), (92, 16.912371686315282), (96, 15.268195474147658), (100, 13.477180596840947), (104, 11.887298863200714), (108, 9.793672686528922), (112, 8.449342524508063), (116, 5.545177444479562), (120, 4.852030263919617), (124, 0.0), (128, 0.6931471805599453) ];

"""
Kullback-Leibler divergence between two distributions

distributions are here considered to be dictionaries

valid for n-dimensional dictionaries (args=Tuple or larger?)
"""
function kldivergence(P::Dict, Q::Dict)
        ##KL divergence sum P(args)logP(args)/Q(args)
        ## P=P_meas, Q=P_true s.t P(args)=0 simply ignored
        #
        kld = 0.0
        for (args, p) in P
            q = Q[args]
            kld += p * log(p / q)
        end
        return kld
end


text:
export PropHook

immutable PropHook
    _hook::AbstractString
    val
end

function diff(a::PropHook, b::PropHook)
    val = diff(a.val, b.val)
    val === nothing ? nothing : PropHook(b._hook, val)
end


text:
using Roots

function opt2τσ(γ,cFrac)
    c = cFrac*γ^2
    lhs(τ) = exp(2τ/γ)*(γ^2-c)*(γ-2)+2exp(τ/γ)*γ*(γ-τ*(γ-1))-γ*(γ^2-c)
    find_zero(lhs,(0,20))
end

function optReward(γ,cFrac)
    A₁(τ,γ) = exp(2τ/γ)*((τ-1)*γ^3-(3τ-2)*γ^2 + 2τ*γ) -2exp(τ/γ)*γ^2*(1-γ)^2 + 2γ^4 + (τ-3)*γ^3 - τ*γ^2
    A₂(τ,γ) = (γ-1)*(exp(2τ/γ)*(γ-2) + γ)
    A₃(τ,γ) = γ^3+(τ-2)*γ^2-τ*γ-exp(2τ/γ)*(γ^3-(τ+2)*γ^2 + 3τ*γ-2τ)
    g(τ,γ,c) = (A₁(τ,γ) - c*A₂(τ,γ))/A₃(τ,γ)

    c = cFrac*γ^2
    τ̂ = opt2τσ(γ,cFrac)
    g(τ̂,γ,c)
end

# opt2τσ(0.4,0.7)

# optReward(0.4,0.5)


text:
# Faster B-spline space

@doc raw"""
B-spline space for lower polynomial degree
```math
\mathcal{P}[p,k]
```
This type `FastBSplineSpace` is faster than `BSplineSpace`, but the degree must be equal or less than `MAX_DEGREE`.
"""
struct FastBSplineSpace{p} <: AbstractBSplineSpace
    knots::Knots
    function FastBSplineSpace(p::Integer, knots::Knots)
        if p < 0
            throw(DomainError(p, "degree of polynominal must be non-negative"))
        elseif p > MAX_DEGREE
            throw(DomainError(p, "FastBSpline supports only degree 0,...,$(MAX_DEGREE)"))
        end
        new{p}(knots)
    end
end
function FastBSplineSpace{q}(p::Integer, knots::Knots) where {q}
    FastBSplineSpace(p,knots)
end

"""
convert AbstractBSplineSpace to FastBSplineSpace
"""
function FastBSplineSpace(P::AbstractBSplineSpace)
    return FastBSplineSpace(degree(P), knots(P))
end

# """
# convert AbstractBSplineSpace to FastBSplineSpace
# """
# function Base.convert(::Type{<:FastBSplineSpace}, P::AbstractBSplineSpace)
#     return FastBSplineSpace(degree(P), knots(P))
# end

function degree(::FastBSplineSpace{p}) where {p}
    return p
end

function knots(P::FastBSplineSpace)
    return P.knots
end


text:
include("../../ApproximateCOSMO.jl/src/COSMO.jl")
using Main.COSMO, Random
include("../src/closest_correlation.jl"); include("../utils/plot_utils.jl")
using RandomCorrelationMatrices

n = 300
C = RandomCorrelationMatrices.randcormatrix(n, 0.1) + 0.1*Symmetric(randn(n, n))
all(eigvals(C) .>= 0) && @warn "The perturbed correlation matrix is still positive definite!"

settings = COSMO.Settings(verbose=true, verbose_timing=true,
    eps_abs=1e-7, eps_rel=1e-7, max_iter=5000,
    rho=.1, adaptive_rho=false,
    check_termination=40, check_infeasibility=1000)
    # kkt_solver=COSMO.with_options(COSMO.IndirectReducedKKTSolver, solver=:CG)
model = create_closest_correlation_model(C, settings; psdcone_type=COSMO.PsdConeTriangleLanczos)
COSMO.optimize!(model)


text:
doc = """Geocoding NYC sales data by cross-referencing lots
with geospatial shapefiles.

Usage:
  NYC_dataprep.jl [--sales=<sales_xls>...] [--parcels=<parcels_shp>...] [--condos=<condos_dbf>]
"""
import CSV
using DataFrames
using Printf
import Dates
import Shapefile
import DBFTables
import ExcelFiles
import ExcelReaders
import Statistics
import Proj4
import GeoInterface
using GeoInterface: coordinates
import Proj4: transform
using DocOpt  # import docopt function

function main()
    salesfiles, parcelfiles, condosfile = parse_docopt()

    # concatenate dataframe
    sales_df = vcat(read_sales_df.(salesfiles)...)
    # load parcel centroids and store information in dictionary
    sbl_to_taxinfo_map = Dict(union(sbl_to_taxinfo_pairs_from_shapefile.(parcelfiles)...))

    condo_bbl_map = DBFTables.Table(condosfile) |> map_condos_to_tax_lot
    sales_taxinfo = extract_sales_taxinfo(sales_df, sbl_to_taxinfo_map, condo_bbl_map)

    # add additional columns to sales dataframe
    sales_df[!,:XCoord] = [info.x for info in sales_taxinfo]
    sales_df[!,:YCoord] = [info.y for info in sales_taxinfo]
    sales_df[!,:SchDistr] = [info.schcode for info in sales_taxinfo]
    summary_of_missing(sales_df)

    CSV.write(Base.stdout, sales_df; dateformat="Y-mm-dd")
end

function parse_docopt()
    args = docopt(doc, version=v"2.0.0")

    salesfiles = args["--sales"]
    parcelfiles = args["--parcels"]
    condosfile = args["--condos"]

    print(stderr, "salesfiles: \n")
    show(stderr, salesfiles)
    print(stderr, "\n")
    print(stderr, "parcelfiles: \n")
    show(stderr, parcelfiles)
    print(stderr, "\n")
    print(stderr, "condosfile: \n")
    show(stderr, condosfile)
    print(stderr, "\n")
    return salesfiles, parcelfiles, condosfile
end

function summary_of_missing(sales_df)
    # proportion of lots that could not be geocoded
    print(stderr, "fraction of sales not geocoded: \n")
    show(stderr, Statistics.mean(ismissing.(sales_df[!,:XCoord])))
    print(stderr, "\n")
    # proportion of lots that could not be geocoded amongst 1- 2- 3-family dwellings
    dwellings = occursin.("DWELLINGS", sales_df.BUILDING_CLASS_CATEGORY)
    print(stderr, "fraction of dwellings not geocoded: \n")
    show(stderr, Statistics.mean(ismissing.(sales_df[dwellings,:XCoord])))
    print(stderr, "\n")
end

fix_alphanumeric_type(x::Real) = (x % 1.0 ≈ 0.0) ? string(Int(x)) : string(x)
fix_alphanumeric_type(x::String) = string(strip(x))
fix_alphanumeric_type(x) = string(x) # handle some dates in Apartment Number column
function read_sales_df(filepath)
    file = ExcelReaders.openxl(filepath)
    sheet1 = file.workbook.sheet_names()[1]
    # first 4 rows contain descriptions of the data
    sales_df = ExcelFiles.load(filepath, sheet1; skipstartrows=4) |> DataFrame
    # change spaces to underscores
    rename!(name -> CSV.normalizename(String(name)), sales_df)
    # clean up types
    for intkey in (:BOROUGH, :BLOCK, :LOT, :ZIP_CODE, :RESIDENTIAL_UNITS,
                   :COMMERCIAL_UNITS, :TOTAL_UNITS, :YEAR_BUILT)
        sales_df[!,intkey] = convert.(Int, sales_df[!,intkey])
    end
    for zeronullkey in (:LAND_SQUARE_FEET, :GROSS_SQUARE_FEET, :COMMERCIAL_UNITS,
                        :TOTAL_UNITS, :RESIDENTIAL_UNITS, :SALE_PRICE, :YEAR_BUILT)
        sales_df[!,zeronullkey] = [x==0 ? missing : x for x in sales_df[!,zeronullkey]]
    end
    for alphanumkey in (:TAX_CLASS_AT_PRESENT, :TAX_CLASS_AT_TIME_OF_SALE, :APARTMENT_NUMBER)
        sales_df[!,alphanumkey] = fix_alphanumeric_type.(sales_df[!,alphanumkey])
    end

    for (key,dtype) in zip(names(sales_df), eltype.(eachcol(sales_df)))
        if dtype <: String
            sales_df[!,key] = String.(strip.(sales_df[!, key]))
        end
    end
    return sales_df
end

struct TaxInfo
    x::Union{Float64, Missing}
    y::Union{Float64, Missing}
    schcode::Union{String, Missing}
    address::Union{String, Missing}
end

function transform(src::Proj4.Projection, dest::Proj4.Projection, multicoords::AbstractVector{V} where V<:AbstractVector)
    return transform.(Ref(src), Ref(dest), multicoords)
end
Shapefile.Point(xy::Vector{Float64}) = Shapefile.Point(xy...)
function transform(src, dist, geom::G) where G<:GeoInterface.AbstractGeometry
    G(transform(src, dist, coordinates(geom)))
end
function transform(src, dist, p::Missing) missing end
function transform_epsg(coords; epsg_from::Int, epsg_to::Int)
    proj_from = Proj4.Projection(Proj4.epsg[epsg_from])
    proj_to = Proj4.Projection(Proj4.epsg[epsg_to])
    return transform(proj_from, proj_to, coords)
end

function sbl_to_taxinfo_pairs_from_shapefile(filepath)
    # read shapefile
    table = Shapefile.Table(filepath)
    # obtain the geometries (points)
    points = Shapefile.shapes(table)
    # school district code
    sch_code = table.SCH_CODE
    # and address
    address = [ismissing(st_nbr) ? missing : st_nbr*" "*street*", "*zipcode
               for (st_nbr, street, zipcode)
               in zip(table.LOC_ST_NBR, table.LOC_STREET, table.LOC_ZIP)]
    projected_points = transform_epsg.(points; epsg_from=26918, epsg_to=2263)
    info = [TaxInfo(ismissing(p) ? missing : p.x,
                  ismissing(p) ? missing : p.y,
                  s, addr)
          for (p, s, addr)
          in zip(projected_points, sch_code, address)]
    # create a vector of pairs of SBL codes and points
    sbl_to_taxinfo_pairs = zip(table.SBL, info)
    return sbl_to_taxinfo_pairs
end

function map_condos_to_tax_lot(dtm_condo)
    condo_bbl_map = Dict{String,String}()
    for row in dtm_condo
        if ismissing(row.UNIT_BBL)
            continue
        elseif ismissing(row.CONDO_BA02)
            continue
        end
        condo_bbl_map[row.UNIT_BBL] = row.CONDO_BA02
    end
    return condo_bbl_map
end

function extract_sales_taxinfo(sales_df, sbl_to_taxinfo_map, condo_bbl_map)
    sales_sbl = sbl_from_components.(sales_df.BOROUGH, sales_df.BLOCK, sales_df.LOT)
    missing_info = TaxInfo(missing,missing,missing,missing)
    sales_taxinfo = TaxInfo[
        # if the SBL code is directly mapped to a TaxInfo record, use that
        if haskey(sbl_to_taxinfo_map, sbl)
            sbl_to_taxinfo_map[sbl]
        elseif haskey(condo_bbl_map, sbl)
            # otherwise it might be a condo, so
            # map the condo lot to a tax lot and
            # try again
            tax_sbl = condo_bbl_map[sbl]
            if haskey(sbl_to_taxinfo_map, tax_sbl)
                sbl_to_taxinfo_map[tax_sbl]
            else
                missing_info
            end
        else
            missing_info
        end
        for sbl in sales_sbl
    ]
    return sales_taxinfo
end

function sbl_from_components(borough::Int, block::Int, lot::Int)
    @sprintf("%01d%05d%04d", borough, block, lot)
end

main()


text:
#Found using a recursive function to be the Fastest. Take a look in the O2-Workbench.jl for other ways
function traverse_recursive_max(node, start_value)
    highest_value = start_value
    if (node.value > highest_value)
        highest_value = node.value
    end
    if node.next == nothing
        return highest_value
    end
    traverse_recursive_max(node.next,highest_value)
end

traversemax(node) = traverse_recursive_max(node,node.value)


text:
immutable DNormal{T <: Real} <: AbstractNormal
    μ::T
    σ::T

    DNormal(μ, σ) = new(μ, σ)
end

# DNormal constructors
DNormal{T <: Real}(μ::T, σ::T) = DNormal{T}(μ, σ)
DNormal(μ::Real, σ::Real) = DNormal(promote(μ, σ)...)

# Forwarding from Normal when one argument is differentiated
Normal(μ::FDN, σ::FDN) = DNormal(μ, σ)
Normal(μ::FDN, σ::Real) = DNormal(μ, σ)
Normal(μ::Real, σ::FDN) = DNormal(μ, σ)


text:
module {{{PKG}}}

export question1,
       question2,
       question3,
       question4,
       question5

include("assignment.jl")

end


text:
@testset "$(rpad("Variational Partitioned Runge-Kutta Tableaus",80))" begin

    @test typeof(TableauVPGLRK(1)) <: TableauVPRK
    @test typeof(TableauVPGLRK(2)) <: TableauVPRK
    @test typeof(TableauVPGLRK(3)) <: TableauVPRK
    @test typeof(TableauVPGLRK(4)) <: TableauVPRK

    @test typeof(TableauVPSRK3()) <: TableauVPRK

    @test typeof(TableauVPLobattoIIIA(2)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIA(3)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIA(4)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIA(5)) <: TableauVPRK

    @test typeof(TableauVPLobattoIIIB(2)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIB(3)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIB(4)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIB(5)) <: TableauVPRK

    @test typeof(TableauVPLobattoIIIC(2)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIC(3)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIC(4)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIC(5)) <: TableauVPRK

    @test typeof(TableauVPLobattoIIIC̄(2)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIC̄(3)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIC̄(4)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIC̄(5)) <: TableauVPRK

    @test typeof(TableauVPLobattoIIID(2)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIID(3)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIID(4)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIID(5)) <: TableauVPRK

    @test typeof(TableauVPLobattoIIIE(2)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIE(3)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIE(4)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIE(5)) <: TableauVPRK

    @test typeof(TableauVPLobattoIIIF(2)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIF(3)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIF(4)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIF(5)) <: TableauVPRK

    @test typeof(TableauVPLobattoIIIG(2)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIG(3)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIG(4)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIG(5)) <: TableauVPRK

    @test typeof(TableauVPLobattoIIIAIIIA(2)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIAIIIA(3)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIAIIIA(4)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIAIIIA(5)) <: TableauVPRK

    @test typeof(TableauVPLobattoIIIBIIIB(2)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIBIIIB(3)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIBIIIB(4)) <: TableauVPRK
    @test typeof(TableauVPLobattoIIIBIIIB(5)) <: TableauVPRK

    @test typeof(TableauVPRadauIIAIIA(2)) <: TableauVPRK
    @test typeof(TableauVPRadauIIAIIA(3)) <: TableauVPRK
    @test typeof(TableauVPRadauIIAIIA(4)) <: TableauVPRK
    @test typeof(TableauVPRadauIIAIIA(5)) <: TableauVPRK

    @test TableauVPGLRK(1) == TableauVPGLRK(Float64, 1)
    @test TableauVPSRK3()  == TableauVPSRK3(Float64)

    @test TableauVPLobattoIIIA(2) == TableauVPLobattoIIIA(Float64, 2)
    @test TableauVPLobattoIIIB(2) == TableauVPLobattoIIIB(Float64, 2)
    @test TableauVPLobattoIIIC(2) == TableauVPLobattoIIIC(Float64, 2)
    @test TableauVPLobattoIIIC̄(2) == TableauVPLobattoIIIC̄(Float64, 2)
    @test TableauVPLobattoIIID(2) == TableauVPLobattoIIID(Float64, 2)
    @test TableauVPLobattoIIIE(2) == TableauVPLobattoIIIE(Float64, 2)
    @test TableauVPLobattoIIIF(2) == TableauVPLobattoIIIF(Float64, 2)
    @test TableauVPLobattoIIIF(2) == TableauVPLobattoIIIF(Float64, 2)
    @test TableauVPLobattoIIIG(2) == TableauVPLobattoIIIG(Float64, 2)
    @test TableauVPRadauIIAIIA(2) == TableauVPRadauIIAIIA(Float64, 2)

end


text:
#!/usr/bin/env julia

using Luxor, Colors

using Test

using Random
Random.seed!(42)

function blend_mode_test(fname)
    Drawing(600, 600, fname)
    origin()
    # transparent, no background
    fontsize(15)
    setline(1)
    tiles = Tiler(600, 600, 4, 5, margin=30)
    modes = length(Luxor.blendingmodes)
    setgrey(0.15)
    for (pos, n) in tiles
        n > modes && break
        gsave()
        translate(pos)
        box(O, tiles.tilewidth-10, tiles.tileheight-10, :clip)

        # diagonal
        diag = (Point(-tiles.tilewidth/2, -tiles.tileheight/2),
                Point(tiles.tilewidth/2,  tiles.tileheight/2))
        upper = between(diag, 0.4)
        lower = between(diag, 0.6)

        # first red shape uses default blend operator
        setcolor(0.7, 0, 0, .8)
        circle(upper, tiles.tilewidth/4, :fill)

        # second blue shape shows results of blend operator
        setcolor((0, 0, 0.9, 0.4))
        blendingmode = Luxor.blendingmodes[mod1(n, modes)]
        setmode(blendingmode)
        circle(lower, tiles.tilewidth/4, :fill)

        clipreset()
        grestore()

        gsave()
        translate(pos)
        text(Luxor.blendingmodes[mod1(n, modes)], O.x, O.y + tiles.tilewidth/2, halign=:center)
        grestore()
    end
    @test finish() == true
end

fname = "blend_mode_test.png"
blend_mode_test(fname)
println("...finished test: output in $(fname)")


text:
using PkgBake
using Test

@testset "options probe" begin
    @test PkgBake.have_trace_compile() == true

end