text:
module SSCant

using PyCall
using ProgressMeter
import ForwardDiff
import ODE
import DSP
import Dierckx
@pyimport lockin

"Fits a straight line through a set of points, `y = a₁ + a₂ * x`

From CurveFit.jl."
linear_fit(x, y) = hcat(ones(x), x) \ y

"""Zero function"""
u(t) = 0

"""
get_specified(x, y, x_specified)
For sorted x, x_specified, get y data points corresponding to x_specified.
"""
function get_specified(x, y, x_specified)
    y_specified = Vector{eltype(y)}()
    j = 1
    for i = eachindex(x)
        if x[i] == x_specified[j]
            push!(y_specified, y[i])
            j += 1
        end
    end
    return y_specified
end

function chunk_range(t0, tf, Tchunk)
    tchunks = collect(t0:Tchunk:tf)
    if tchunks[end] != tf
        push!(tchunks, tf)
    end
    return tchunks
end

function make_spline(f, t0, tf, dt; k=3)
    T = tf - t0
    N = round(Int, T/dt) + 1
    t = linspace(t0, tf, N)
    ft = map(f, t)
    F_spline = Dierckx.Spline1D(t, ft, k=k)
    F(t) = Dierckx.evaluate(F_spline, t)
    return F
end

function make_spline(f, t0, tf, dt, prev; k=3)
    T = tf - t0
    N = round(Int, T/dt) + 1
    t = linspace(t0, tf, N)

    ft = Vector{typeof(prev)}(length(t))
    ft[1] = prev

    for i = 2:length(t)
        ft[i] = f(t[i])
    end

    F_spline = Dierckx.Spline1D(t, ft, k=k)
    F(t) = Dierckx.evaluate(F_spline, t)

    return F
end

"""
Params
------
x: position
v: velocity
omega_c: angular frequency (optional)
f_c: frequency (optional)
"""
function amplitude(x, v; omega_c=nothing, f_c=nothing)
    if omega_c == nothing
        omega_c = 2 * pi * f_c
    end
    (x.^2 + (v ./ omega_c).^2).^0.5
end

function phase(x, v; omega_c=nothing, f_c=nothing)
    if omega_c == nothing
        omega_c = 2 * pi * f_c
    end
    DSP.unwrap(atan2(-v./omega_c, x))
end

function fit_phase(t, phi)
    mb = linear_fit(t, phi)
    dphi = phi - (t * mb[2] + mb[1])
    df = zeros(size(dphi))
    df[2:end] = diff(dphi)

    return dphi, df, mb[2] / (2*pi)
end

"""
- f_c: cantilever frequency (0.05 MHz)
- Q:  quality factor (1000)
- k: spring constant (1 µN/µm)
- C_z: first derivative of tip capacitance (1e-6 pF / µm)
- C_zz: second derivative of tip capacitance (-1e-5 pF / µm)
- C_b: cantilever body capacitance (0.1 pF)
- R_t: Cantilever tip capacitance (1e-4 Mohm)
- C_1: C1 capacitance (600 pF)
- R_1: R1 resistance (3e-6 Mohm)
- tau_s: sample response time (15e3 us)
"""
immutable CantileverParams
    f_c::Float64
    Q::Float64
    k::Float64
    C_z::Float64
    C_zz::Float64
    C_b::Float64
    R_t::Float64
    C_1::Float64
    R_1::Float64
    tau_s::Float64
    T::Float64
    f_det::Float64
end

function CantileverParams(;f_c=0.05, Q=1000., k=1.,
                 C_z=1.6e-4, C_zz=-1.6e-5, C_b=0.1,
                 R_t=1e-4, C_1=600., R_1=3e-6,
                 tau_s=15e3, T=293., f_det=0.2)
    CantileverParams(f_c, Q, k,
                 C_z, C_zz, C_b,
                 R_t, C_1, R_1,
                 tau_s, T, f_det)
end

"""
Params
------

- f_c: cantilever frequency (0.05 MHz)
- Q:  quality factor (1000)
- k: spring constant (1 µN/µm)
- C_z: first derivative of tip capacitance (1e-6 pF / µm)
- C_zz: second derivative of tip capacitance (-1e-5 pF / µm)
- C_b: cantilever body capacitance (0.1 pF)
- R_t: Cantilever tip capacitance (1e-4 Mohm)
- C_1: C1 capacitance (600 pF)
- R_1: R1 resistance (3e-6 Mohm)
- tau_s: sample response time (15e3 us)
- F: Force as a function of time (µN)
- V: Voltage as a function of time (V)
- U: percentange change in capacitance (unitless)

Returns
-------

"""
function ss_cant_1tc(f_c=0.05, Q=1000., k=1,
                 C_z=-1e-6, C_zz=1e-5, C_b=0.1,
                 R_t=1e-4, C_1=600., R_1=3e-6,
                 tau_s=15e3, F=u, V=u, U=u)

    omega_c = 2*pi*f_c
    omega_r = omega_c / Q
    m = k / omega_c ^ 2
    mC_b2 = m * C_b^2

    omega_11 = (C_1 * R_1)^-1
    omega_t1 = (R_t * C_1)^-1

    omega_tb = (R_t * C_b)^-1

    F_x(x) = [x[2], (
            -omega_c^2 * x[1] - omega_r * x[2]
            + C_z/(2mC_b2) * (1 + x[5]) * x[3]^2
            + (C_zz - 2*C_z^2/C_b)/(2mC_b2) * (1 + x[5]) * x[1]*x[3]^2),

            (-omega_tb * x[3]
            + C_z/(R_t * C_b^2) * (1 + x[5]) * x[1] * x[3]
            + omega_t1 * x[4]),

            (omega_tb * x[3]
            - C_z/(R_t * C_b^2) * (1 + x[5]) * x[1] * x[3]
            - (omega_11 + omega_t1) * x[4]),

            -x[5] / tau_s
            ]

    function F_tx(t, x)
        return F_x(x) + [0., F(t) / m, 0,
                         V(t) / R_1, U(t) / tau_s]
    end

    F_x, F_tx
end

ss_cant_1tc(cant::CantileverParams, F=u, V=u, U=u) = ss_cant_1tc(cant.f_c,
    cant.Q, cant.k, cant.C_z, cant.C_zz, cant.C_b, cant.R_t, cant.C_1,
    cant.R_1, cant.tau_s, F, V, U)

"""
Params
------

- f_c: cantilever frequency (0.05 MHz)
- Q:  quality factor (1000)
- k: spring constant (1 µN/µm)
- C_z: first derivative of tip capacitance (1e-6 pF / µm)
- C_zz: second derivative of tip capacitance (-1e-5 pF / µm)
- C_b: cantilever body capacitance (0.1 pF)
- R_t: Cantilever tip capacitance (1e-4 Mohm)
- C_1: C1 capacitance (600 pF)
- R_1: R1 resistance (3e-6 Mohm)
- tau_s: sample response time (15e3 us)
- F: Force as a function of time (µN)
- V: Voltage as a function of time (V)
- U: percentange change in capacitance (unitless)
- Pn: Surface potential input (V/s)

Returns
-------

"""
function ss_cant_1tc_volt_noise(f_c=0.05, Q=1000., k=1,
                 C_z=1e-6, C_zz=-1e-5, C_b=0.1,
                 R_t=1e-4, C_1=600., R_1=3e-6,
                 tau_s=15e3, F=u, V=u, U=u, Pn=u)

    omega_c = 2*pi*f_c
    omega_r = omega_c / Q
    m = k / omega_c ^ 2
    mC_b2 = m * C_b^2

    omega_11 = (C_1 * R_1)^-1
    omega_t1 = (R_t * C_1)^-1

    omega_tb = (R_t * C_b)^-1

    F_x(x) = [x[2], (
            -omega_c^2 * x[1] - omega_r * x[2]
            - C_z/(2mC_b2) * (1 + x[5]) * x[3]^2
            - C_zz/(2mC_b2) * (1 + x[5]) * x[1]*x[3]^2),

            (-omega_tb * x[3]
            + C_z/(R_t * C_b^2) * (1 + x[5]) * x[1] * x[3]
            + omega_t1 * x[4]),

            (omega_tb * x[3]
            - C_z/(R_t * C_b^2) * (1 + x[5]) * x[1] * x[3]
            - (omega_11 + omega_t1) * x[4]
            - x[6] / R_1
            ),

            -x[5] / tau_s,

            0.0
            ]

    function F_tx(t, x)
        return F_x(x) + [0., F(t) / m, 0,
                         V(t) / R_1, U(t) / tau_s, Pn(t)]
    end

    F_x, F_tx
end

ss_cant_1tc_volt_noise(
    cant::CantileverParams, F=u, V=u, U=u, Pn=u
    ) = ss_cant_1tc_volt_noise(cant.f_c,
    cant.Q, cant.k, cant.C_z, cant.C_zz, cant.C_b, cant.R_t, cant.C_1,
    cant.R_1, cant.tau_s, F, V, U, Pn)

function init_steady_state_1tc(cp::SSCant.CantileverParams, A0, phi0, t0, V, U)
    [A0*cos(phi0), -2*pi*cp.f_c*A0*sin(phi0),
     cp.C_b*V(t0), cp.C_1*V(t0), U(t0)]
end

function init_steady_state_1tc_volt_noise(cp::SSCant.CantileverParams, A0, phi0, t0, V, U, Pn)
    [A0*cos(phi0), -2*pi*cp.f_c*A0*sin(phi0),
     cp.C_b*V(t0), cp.C_1*V(t0), U(t0), Pn(t0)]
end

function ss_cant_det(;f_c=0.05, k=1, Q=10000,
                 C_z=1e-6, C_zz=-1e-5, C_b=0.1,
                 R_t=1e-4, C_1=600., R_1=3e-6,
                 tau_s=15e3, fdet=0.2, F=u, V=u, U=u, P=u)

    omega_det = 2*pi*fdet
    omega_c = 2*pi*f_c
    omega_r = omega_c / Q
    m = k / omega_c ^ 2
    mC_b2 = m * C_b^2

    omega_11 = (C_1 * R_1)^-1
    omega_t1 = (R_t * C_1)^-1

    omega_tb = (R_t * C_b)^-1

    Ftx(t, x) = [x[2],
        (
        -omega_c^2 * x[1] - omega_r  * x[2]
            - C_z/(2mC_b2) * (1 + x[5]) * x[3]^2
            - (C_zz - 2*C_z^2/C_b)/(2mC_b2) * (1 + x[5]) * x[1]*x[3]^2

        + F(t) / m
        ),

            (-omega_tb * x[3]
            + C_z/(R_t * C_b^2) * (1 + x[5]) * x[1] * x[3]
            + omega_t1 * x[4]),

            (omega_tb * x[3]
            - C_z/(R_t * C_b^2) * (1 + x[5]) * x[1] * x[3]
            - (omega_11 + omega_t1) * x[4]
            - x[6] / R_1

            + V(t) / R_1
            ),

            -x[5] / tau_s + U(t) / tau_s,

            P(t),
            omega_det * (x[1] - x[7])
            ]

    Ftx
end

ss_cant_det(cant::SSCant.CantileverParams,
F=u, V=u, U=u, P=u) = ss_cant_det(f_c = cant.f_c, k=cant.k, Q=cant.Q, C_z=cant.C_z,
C_zz=cant.C_zz, C_b=cant.C_b, R_t=cant.R_t, C_1=cant.C_1, R_1=cant.R_1,
tau_s=cant.tau_s, F=F, V=V, U=U, P=P, fdet=cant.f_det)

function init_ss_cant_det(cant::SSCant.CantileverParams, A0, phi0, t0; V=u, U=u, P=u)
    resp = 1 / (1 + cant.f_c * im / cant.f_det)
    A_x = A0 / abs(resp)
    phi_x = phi0 - angle(resp)
    [A_x*cos(phi_x), -2*pi*cant.f_c*A_x*sin(phi_x), cant.C_b*V(t0),
     cant.C_1*V(t0), U(t0), P(t0), A0*cos(phi0)]
end

function ss_cant_xv(f_c=0.05, Q=1000., k=1., F=u)

    omega_c = 2*pi*f_c
    omega_r = omega_c / Q
    m = k / omega_c ^ 2

    F_x(x) = [x[2], -omega_c^2 * x[1] - omega_r * x[2]]

    function F_tx(t, x)
        return F_x(x) + [0., F(t) / m]
    end

    F_x, F_tx
end

ss_cant_xv(cant::CantileverParams, F=u, V=u, U=u) = ss_cant_xv(
    cant.f_c, cant.Q, cant.k, F)

type TimeSeries
    cant::CantileverParams
    y0
    t
    x
    v
    q_t
    q_1
    x_5
    A
    phi
    dphi
    df
    f0
end

function TimeSeries(cant::CantileverParams, y0, t, y)
    x = y[:, 1]
    v = y[:, 2]
    q_t = y[:, 3]
    q_1 = y[:, 4]
    x_5 = y[:, 5]
    A = amplitude(x, v, f_c=cant.f_c)
    phi = phase(x, v, f_c=cant.f_c)
    dphi, df, f0 = fit_phase(t, phi)

    TimeSeries(cant, y0, t, x, v, q_t, q_1, x_5,
               A, phi, dphi, df, f0)
end

type ResampledTimeSeries
    cant::CantileverParams
    ts::TimeSeries
    dt
    m
    t
    x
    A
    phi
    dphi
    df
    f0
end

type RawTimeSeries
    cant::CantileverParams
    ts::TimeSeries
end

immutable DownsampledTimeSeries
    cant::CantileverParams
    y0::Vector{Float64}
    dt0::Float64
    fs::Float64
    f0::Float64
    phi0::Float64
    t::Vector{Float64}
    A::Vector{Float64}
    phi::Vector{Float64}
    df::Vector{Float64}
end

immutable DownsampledTimeSeriesRaw
    cant::CantileverParams
    y0::Vector{Float64}
    t_raw::Vector{Float64}
    y_raw::Array{Float64,2}
    dt0::Float64
    fs::Float64
    f0::Float64
    phi0::Float64
    t::Vector{Float64}
    A::Vector{Float64}
    phi::Vector{Float64}
    df::Vector{Float64}
end

function squeeze_y(y)
    y_out = zeros(Float64, length(y), length(y[1]))
    for i in 1:length(y)
        y_out[i, :] = y[i]
    end
    y_out
end

function jacobian(dx)
    jac_t(t, x) = ForwardDiff.jacobian(dx, x)

    return jac_t
end

"""
Params
------

- cant: CantileverParams
- y0: array containing [x, v, q_t, q_1, x_5] at time zero
      The units are µm, µm / µs, pC, pC, unitless (percentage change)
"""
function resample_sscant_sim(cant::CantileverParams, y0::Vector{Float64},
                             T_pre::Float64,
                             T::Float64, dt::Float64;
                             F=u, V=u, U=u, coeff_ratio=8., args...)
    fs = 1/dt

    F_tx, jac_t = make_Ftx_jac(ss_cant_1tc, cant, F, V, U)

    t, _y = ODE.ode23s(F_tx, y0, [-T_pre, T], jacobian=jac_t; maxstep=dt*0.5, args...)
    y = squeeze_y(_y)

    ts = TimeSeries(cant, y0, t, y)

    spline = Dierckx.Spline1D(ts.t, ts.x, k=3)
    t = collect(-T_pre:dt:T)
    x = Dierckx.evaluate(spline, t)

    li = lockin.LockIn(t, x, fs)

    li[:lock2](fp_ratio=0.1, fc_ratio=0.4,
               coeff_ratio=coeff_ratio, window="blackman")
    li[:phase](tf=0.)

    ResampledTimeSeries(cant, ts, dt, li[:m], t, x, li[:A],
                        li[:phi], li[:dphi], li[:df], li[:f0])

end

function raw_sscant_sim(cant::CantileverParams, y0::Vector{Float64},
                             T_pre::Float64,
                             T::Float64, dt::Float64;
                             F=u, V=u, U=u, coeff_ratio=8., args...)
    fs = 1/dt

    F_tx, jac_t = make_Ftx_jac(ss_cant_1tc, cant, F, V, U)

    t, _y = ODE.ode23s(F_tx, y0, [-T_pre, T], jacobian=jac_t; maxstep=dt*0.5, args...)
    y = squeeze_y(_y)

    ts = TimeSeries(cant, y0, t, y)

    RawTimeSeries(cant, ts)

end

function make_Ftx_jac(ss_cant, cant::CantileverParams, x...)

    Fx, F_tx = ss_cant(cant, x...)
    jac = jacobian(Fx)

    return F_tx, jac
end

function chunk_sscant_sim(cant::CantileverParams, y0::Vector{Float64},
                             T_pre::Float64,
                             T::Float64, dt::Float64;
                             F=u, V=u, U=u, coeff_ratio=8.,
                             Tchunk=1000., f_dec=nothing, args...)
    fs = 1/dt

    if f_dec == nothing
        f_dec = 2*cant.f_c
    end

    dec = trunc(Int, fs / f_dec)

    F_tx, jac = make_Ftx_jac(ss_cant_1tc, cant, F, V, U)

    Tendpts = linspace(-T_pre, T, ceil(Int, (T_pre+T)/Tchunk))
    for i = 1:length(Tendpts)-1
        t, _y = ODE.ode23s(F_tx, y0, [Tendpts[i], Tendpts[i+1]],
                           jacobian=jac; maxstep=dt*0.5)
        y = squeeze_y(_y)
        y0 = _y[end]
        spline = Dierckx.Spline1D(t, y[:, 1], k=3)
        t_equal = collect(Tendpts[i]:dt:Tendpts[i+1])
        x = Dierckx.evaluate(spline, t_equal)

        if i == 1
            global li = lockin.LockIn(t_equal, x, fs)
            li[:lock2](fp_ratio=0.1, fc_ratio=0.4,
                       coeff_ratio=coeff_ratio, window="blackman")
            li[:phase](tf=0.)
            global slock = lockin.FIRStateLock(li[:fir], dec, li[:f0],
                                        -li[:mb][end], t0=-T_pre, fs=fs)
            slock[:filt](x)
        else
            slock[:filt](x[2:end])
        end

    end

    tout = slock[:get_t]()
    A = abs(slock[:output])
    phi = DSP.unwrap(angle(slock[:z]))
    df = zeros(phi)
    df[2:end] = diff(phi)

    DownsampledTimeSeries(cant, y0, dt, fs / dec,
                          slock[:f0], slock[:phi0][end], tout,
                          A, phi, df)
end

function downsample_sscant_sim(cant::CantileverParams, y0::Vector{Float64},
                             T_pre::Float64,
                             T::Float64, dt::Float64;
                             F=u, V=u, U=u, coeff_ratio=8.,
                             f_resampled=nothing, args...)
    if f_resampled == nothing
        f_resampled = cant.f_c
    end

    fs = 1/dt

    F_tx, jac_t = make_Ftx_jac(ss_cant_1tc, cant, F, V, U)

    t, _y = ODE.ode23s(F_tx, y0, [-T_pre, T], jacobian=jac_t; maxstep=dt*0.5, args...)
    y = squeeze_y(_y)

    ts = TimeSeries(cant, y0, t, y)

    spline = Dierckx.Spline1D(ts.t, ts.x, k=3)
    t = collect(-T_pre:dt:T)
    x = Dierckx.evaluate(spline, t)

    li = lockin.LockIn(t, x, fs)

    li[:lock2](fp_ratio=0.1, fc_ratio=0.4, coeff_ratio=coeff_ratio, window="blackman")
    li[:phase](tf=0.)
    factor = trunc(fs / f_resampled)
    li[:decimate](factor)
    phi0 = li[:mb][2]

    DownsampledTimeSeries(cant, y0, dt, li[:dec_fs],
                          li[:dec_f0], li[:mb][2], li[:dec_t],
                          li[:dec_A], li[:dec_phi], li[:dec_df])
end

# Since a discrete delta function has a magnitude $\delta[0] = 1/\Delta t$,
# we can model the thermal force on the cantilever as a random,
# Gaussian force with,

# $$\sigma = \sqrt{\frac{k_c k_B T}{\pi Q f_c \Delta t}}$$
# (see doi:10.1088/0034-4885/29/1/306).

function resample_sscant_noise(cant::CantileverParams, y0::Vector{Float64},
                             T_pre::Float64,
                             T::Float64, dt::Float64;
                             V=u, U=u, coeff_ratio=8., seed=-1, args...)
    fs = 1/dt
    k_B = 1.38064852e-11

    dt_force = 0.05/cant.f_c
    sigma = sqrt(cant.k*k_B*cant.T/(pi*cant.Q*cant.f_c*dt_force))
    if seed > 0
        srand(seed)
    end
    N_force = round(Int, T/dt_force)+1

    F_spline = Dierckx.Spline1D(linspace(0, T, N_force), randn(N_force) * sigma,
                              k=1)

    function F(t)
        Dierckx.evaluate(F_spline, t)
    end

    F_tx, jac_t = make_Ftx_jac(ss_cant_1tc, cant, F, V, U)

    t, _y = ODE.ode23s(F_tx, y0, [-T_pre, T], jacobian=jac_t; maxstep=dt*0.5, args...)
    y = squeeze_y(_y)

    ts = TimeSeries(cant, y0, t, y)

    spline = Dierckx.Spline1D(ts.t, ts.x, k=3)
    t = collect(-T_pre:dt:T)
    x = Dierckx.evaluate(spline, t)

    li = lockin.LockIn(t, x, fs)

    li[:lock2](fp_ratio=0.1, fc_ratio=0.4,
               coeff_ratio=coeff_ratio, window="blackman")
    li[:autophase](tf=0.)

    ResampledTimeSeries(cant, ts, dt, li[:m], t, x, li[:A],
                        li[:phi], li[:dphi], li[:df], li[:f0])

end

# Another idea would be to "remember" the last value the

function re_down_sscant_noise3(cant::CantileverParams, y0::Vector{Float64},
                             T_pre::Float64,
                             T::Float64, dt::Float64;
                             V=u, U=u, coeff_ratio=8., seed=-1, args...)
    fs = 1/dt
    k_B = 1.38064852e-11

    dt_force = 0.05/cant.f_c
    sigma = sqrt(cant.k*k_B*cant.T/(pi*cant.Q*cant.f_c*dt_force))
    if seed > 0
        srand(seed)
    end
    N_force = round(Int, T/dt_force)+1
    N = round(Int, (T+T_pre)/dt_force)+1

    force(t) = t >= 0. ? randn() * sigma : 0.
    times = linspace(-T_pre, T, N)
    F_spline = Dierckx.Spline1D(times, map(force, times),
                              k=3)

    function F(t)
        Dierckx.evaluate(F_spline, t)
    end

    F_tx, jac_t = make_Ftx_jac(ss_cant_1tc, cant, F, V, U)

    t, _y = ODE.ode23s(F_tx, y0, [-T_pre, T], jacobian=jac_t; maxstep=dt*0.5, args...)
    y = squeeze_y(_y)

    ts = TimeSeries(cant, y0, t, y)

    spline = Dierckx.Spline1D(ts.t, ts.x, k=3)
    t = collect(-T_pre:dt:T)
    x = Dierckx.evaluate(spline, t)

    li = lockin.LockIn(t, x, fs)

    li[:lock2](fp_ratio=0.1, fc_ratio=0.4,
               coeff_ratio=coeff_ratio, window="blackman")
    li[:phase](tf=0.)
    factor = trunc(fs/cant.f_c)
    li[:decimate](factor)

    r = ResampledTimeSeries(cant, ts, dt, li[:m], t, x, li[:A],
                        li[:phi], li[:dphi], li[:df], li[:f0])

    down = DownsampledTimeSeries(cant, y0, dt, li[:dec_fs],
                          li[:dec_f0], li[:mb][2], li[:dec_t],
                          li[:dec_A], li[:dec_phi], li[:dec_df])

    return r, down

end

function down_sscant_noise3(cant::CantileverParams, y0::Vector{Float64},
                             T_pre::Float64,
                             T::Float64, dt::Float64;
                             V=u, U=u, coeff_ratio=8., seed=-1, args...)
    fs = 1/dt
    k_B = 1.38064852e-11

    dt_force = 0.05/cant.f_c
    sigma = sqrt(cant.k*k_B*cant.T/(pi*cant.Q*cant.f_c*dt_force))
    if seed > 0
        srand(seed)
    end
    N_force = round(Int, T/dt_force)+1
    N = round(Int, (T+T_pre)/dt_force)+1

    force(t) = randn() * sigma
    times = linspace(-T_pre, T, N)
    F_spline = Dierckx.Spline1D(times, map(force, times),
                              k=3)

    function F(t)
        Dierckx.evaluate(F_spline, t)
    end

    F_tx, jac_t = make_Ftx_jac(ss_cant_1tc, cant, F, V, U)

    t, _y = ODE.ode23s(F_tx, y0, [-T_pre, T], jacobian=jac_t; maxstep=dt*0.5, args...)
    y = squeeze_y(_y)

    ts = TimeSeries(cant, y0, t, y)

    spline = Dierckx.Spline1D(ts.t, ts.x, k=3)
    t = collect(-T_pre:dt:T)
    x = Dierckx.evaluate(spline, t)

    li = lockin.LockIn(t, x, fs)

    li[:lock2](fp_ratio=0.1, fc_ratio=0.4,
               coeff_ratio=coeff_ratio, window="blackman")
    li[:phase](tf=0.)
    factor = trunc(fs/cant.f_c)
    li[:decimate](factor)

    down = DownsampledTimeSeries(cant, y0, dt, li[:dec_fs],
                          li[:dec_f0], li[:mb][2], li[:dec_t],
                          li[:dec_A], li[:dec_phi], li[:dec_df])

    return down

end

# Eliminate Spline interpolation of x. Use tspan argument instead.
function down_sscant_noise4(cant::CantileverParams, y0::Vector{Float64},
                             T_pre::Float64,
                             T::Float64, dt::Float64;
                             V=u, U=u, coeff_ratio=8., maxstep=dt*0.5,
                             seed=-1, fs_force=20*cant.f_c, args...)
    fs = 1/dt
    k_B = 1.38064852e-11
    t0 = -T_pre
    tf = t0 + T

    if seed > 0
        srand(seed)
    end

    # Cut off force noise at 20•fc by default
    dt_force = 1/fs_force
    sigma = sqrt(cant.k*k_B*cant.T/(pi*cant.Q*cant.f_c*dt_force))

    force(t) = randn() * sigma

    times = t0:dt:tf

    F = make_spline(force, times[1], times[end], dt_force)

    F_tx, jac_t = make_Ftx_jac(ss_cant_1tc, cant, F, V, U)
    t, _y = ODE.ode23s(F_tx, y0, times, jacobian=jac_t; points=:specified,
                       maxstep=maxstep, args...)

    y_specified = squeeze_y(_y)

    return TimeSeries(cant, y0, times, y_specified)

end

# Eliminate Spline interpolation of x. Use tspan argument instead.
function chunk_sscant_noise4(cant::CantileverParams, y0::Vector{Float64},
                             T_pre::Float64,
                             T::Float64, dt::Float64, Tchunk;
                             V=u, U=u, coeff_ratio=8., maxstep=dt*0.5, seed=-1,
                             fs_force=20*cant.f_c, args...)
    fs = 1/dt
    k_B = 1.38064852e-11
    t0 = -T_pre
    tf = t0 + T

    if seed > 0
        srand(seed)
    end

    # Cut off force noise at 20•fc by default
    dt_force = 1/fs_force
    sigma = sqrt(cant.k*k_B*cant.T/(pi*cant.Q*cant.f_c*dt_force))
    force(t) = randn() * sigma

    times = t0:dt:tf
    N_total = length(times)
    N_chunk_pts = Int(Tchunk/dt)
    i_endpts = chunk_range(1, N_total, N_chunk_pts)

    y_specified = Matrix{Float64}(N_total, length(y0))
    y_specified[1, :] = y0

    y1 = y0
    F_prev_endpoint = 0.

    for i in eachindex(i_endpts[1 : end-1])
        j0 = i_endpts[i]
        jf = i_endpts[i+1]
        curr_times = times[j0 : jf]

        if i == 1
            F = make_spline(force, curr_times[1], curr_times[end], dt_force)
        else
            F = make_spline(force, curr_times[1], curr_times[end], dt_force,
                            F_prev_endpoint)
        end

        F_prev_endpoint = F(curr_times[end])

        F_tx, jac_t = make_Ftx_jac(ss_cant_1tc, cant, F, V, U)

        t, _y = ODE.ode23s(F_tx, y1, curr_times, jacobian=jac_t;
                           points=:specified,
                           maxstep=maxstep, args...)

        y1 = _y[end]

        _y_specified = squeeze_y(_y)
        y_specified[j0+1 : jf, :] = _y_specified[2 : end, :]
    end

    return TimeSeries(cant, y0, times, y_specified)

end

function N_chunks(T, Tchunk, dt, t0=0.)
    times = t0:dt:(t0+T)
    N_total = length(times)
    N_chunk_pts = Int(Tchunk/dt)
    i_endpts = chunk_range(1, N_total, N_chunk_pts)
    return length(i_endpts)-1
end

# Add voltage noise term.
function chunk_sscant_noise_volt_noise(cant::CantileverParams, y0::Vector{Float64},
                             T_pre::Float64,
                             T::Float64, dt::Float64, Tchunk;
                             V=u, U=u, maxstep=dt*0.5, seed=-1,
                             fs_force=20*cant.f_c, sigma_potential=1e-5,
                             dt_potential=10., progressbar=nothing,args...)
    fs = 1/dt
    k_B = 1.38064852e-11
    t0 = -T_pre
    tf = t0 + T

    if seed > 0
        srand(seed)
    end

    # Cut off force noise at 20•fc by default
    dt_force = 1/fs_force
    sigma = sqrt(cant.k * k_B * cant.T / (pi * cant.Q * cant.f_c * dt_force))
    force(t) = randn() * sigma
    potential_noise(t) = randn() * sigma_potential

    times = t0:dt:tf
    N_total = length(times)
    N_chunk_pts = Int(Tchunk/dt)
    i_endpts = chunk_range(1, N_total, N_chunk_pts)

    y_specified = Matrix{Float64}(N_total, length(y0))
    y_specified[1, :] = y0

    y1 = y0

    F_prev =0.
    Pn_prev = 0.

    if progressbar == true
        progressbar = Progress(length(i_endpts[1 : end-1]))
    end

    for i in eachindex(i_endpts[1 : end-1])
        j0 = i_endpts[i]
        jf = i_endpts[i+1]
        curr_times = times[j0 : jf]

        # Might be nice to properly handle endpoints, by using
        # F_prev_endpoint as first function value. Then evaluate the function
        # at t[2:end]
        if i == 1
            F = make_spline(force, curr_times[1], curr_times[end], dt_force)
            Pn = make_spline(potential_noise, curr_times[1], curr_times[end],
                         dt_potential)
        else
            F = make_spline(force, curr_times[1], curr_times[end], dt_force,
                            F_prev)
            Pn = make_spline(potential_noise, curr_times[1], curr_times[end],
                         dt_potential, Pn_prev)
        end
        F_prev = F(curr_times[end])
        Pn_prev = Pn(curr_times[end])

        F_tx, jac_t = make_Ftx_jac(ss_cant_1tc_volt_noise, cant, F, V, U, Pn)

        t, _y = ODE.ode23s(F_tx, y1, curr_times, jacobian=jac_t;
                           points=:specified,
                           maxstep=maxstep, args...)

        y1 = _y[end]

        _y_specified = squeeze_y(_y)
        y_specified[j0+1 : jf, :] = _y_specified[2 : end, :]

        if progressbar != nothing
            next!(progressbar)
        end
    end

    return times, y_specified

end

function make_jacobian(f)
    jac(t, x) = ForwardDiff.jacobian(x_ -> f(t, x_), x)
    return jac
end

# Add voltage noise term.
function chunk_sscant_noise_volt_noise_det(cant::CantileverParams, y0::Vector{Float64},
                             T_pre::Float64,
                             T::Float64, dt::Float64, Tchunk;
                             V=u, U=u, maxstep=dt*0.5, seed=-1,
                             fs_force=20*cant.f_c, sigma_potential=1e-5,
                             dt_potential=10., progressbar=nothing,args...)
    fs = 1/dt
    k_B = 1.38064852e-11
    t0 = -T_pre
    tf = t0 + T

    if seed > 0
        srand(seed)
    end

    # Cut off force noise at 20•fc by default
    dt_force = 1/fs_force
    sigma = sqrt(cant.k * k_B * cant.T / (pi * cant.Q * cant.f_c * dt_force))
    force(t) = randn() * sigma
    potential_noise(t) = randn() * sigma_potential

    times = t0:dt:tf
    N_total = length(times)
    N_chunk_pts = Int(Tchunk/dt)
    i_endpts = chunk_range(1, N_total, N_chunk_pts)

    y_specified = Matrix{Float64}(N_total, length(y0))
    y_specified[1, :] = y0

    y1 = y0

    F_prev =0.
    Pn_prev = 0.

    for i in eachindex(i_endpts[1 : end-1])
        j0 = i_endpts[i]
        jf = i_endpts[i+1]
        curr_times = times[j0 : jf]

        # Might be nice to properly handle endpoints, by using
        # F_prev_endpoint as first function value. Then evaluate the function
        # at t[2:end]
        if i == 1
            F = make_spline(force, curr_times[1], curr_times[end], dt_force)
            Pn = make_spline(potential_noise, curr_times[1], curr_times[end],
                         dt_potential)
        else
            F = make_spline(force, curr_times[1], curr_times[end], dt_force,
                            F_prev)
            Pn = make_spline(potential_noise, curr_times[1], curr_times[end],
                         dt_potential, Pn_prev)
        end
        F_prev = F(curr_times[end])
        Pn_prev = Pn(curr_times[end])

        F_tx = ss_cant_det(cant, u, V, U, Pn)
        jac_t = make_jacobian(F_tx)

        t, _y = ODE.ode23s(F_tx, y1, curr_times, jacobian=jac_t;
                           points=:specified,
                           maxstep=maxstep, args...)

        y1 = _y[end]

        _y_specified = squeeze_y(_y)
        y_specified[j0+1 : jf, :] = _y_specified[2 : end, :]

        if progressbar != nothing
            next!(progressbar)
        end
    end

    return times, y_specified

end

# # Add voltage noise term.
# function chunk_sscant_noise_volt_noise_sep_pulse(cant::CantileverParams, y0::Vector{Float64},
#                              T_pre::Float64,
#                              T::Float64, dt::Float64, Tchunk, tp;
#                              V=u, U=u, maxstep=dt*0.5, seed=-1,
#                              fs_force=20*cant.f_c, sigma_potential=1e-5,
#                              dt_potential=10., progressbar=nothing,args...)
#     fs = 1/dt
#     k_B = 1.38064852e-11
#     t0 = -T_pre
#     tf = t0 + T

#     if seed > 0
#         srand(seed)
#     end

#     # Cut off force noise at 20•fc by default
#     dt_force = 1/fs_force
#     sigma = sqrt(cant.k * k_B * cant.T / (pi * cant.Q * cant.f_c * dt_force))
#     force(t) = randn() * sigma
#     potential_noise(t) = randn() * sigma_potential

#     times = t0:dt:tf

#     N_total = length(times)
#     N_chunk_pts = Int(Tchunk/dt)
#     i_endpts = chunk_range(1, N_total, N_chunk_pts)

#     y_specified = Matrix{Float64}(N_total, length(y0))
#     y_specified[1, :] = y0

#     y1 = y0

#     F_prev =0.
#     Pn_prev = 0.

#     i=1
#     pulse_

#     while i < length(i_endpts[1 : end-1])
#         j0 = i_endpts[i]
#         jf = i_endpts[i+1]
#         curr_times = times[j0 : jf]
#         curr_ti = curr_times[1]
#         curr_tf = curr_times[end]

#         # Might be nice to properly handle endpoints, by using
#         # F_prev_endpoint as first function value. Then evaluate the function
#         # at t[2:end]
#         if i == 1
#             F = make_spline(force, curr_times[1], curr_times[end], dt_force)
#             Pn = make_spline(potential_noise, curr_times[1], curr_times[end],
#                          dt_potential)
#         else
#             F = make_spline(force, curr_times[1], curr_times[end], dt_force,
#                             F_prev)
#             Pn = make_spline(potential_noise, curr_times[1], curr_times[end],
#                          dt_potential, Pn_prev)
#         end
#         F_prev = F(curr_times[end])
#         Pn_prev = Pn(curr_times[end])

#         F_tx, jac_t = make_Ftx_jac(ss_cant_1tc_volt_noise, cant, F, V, U, Pn)

#         t, _y = ODE.ode23s(F_tx, y1, curr_times, jacobian=jac_t;
#                            points=:specified,
#                            maxstep=maxstep, args...)

#         y1 = _y[end]

#         _y_specified = squeeze_y(_y)
#         y_specified[j0+1 : jf, :] = _y_specified[2 : end, :]

#         if progressbar != nothing
#             next!(progressbar)
#         end
#     end

#     return times, y_specified

# end

function sscant_noise3(cant::CantileverParams, y0::Vector{Float64},
                             T_pre::Float64,
                             T::Float64, dt::Float64;
                             V=u, U=u, coeff_ratio=8., seed=-1, args...)
    fs = 1/dt
    k_B = 1.38064852e-11

    dt_force = 0.05/cant.f_c
    sigma = sqrt(cant.k*k_B*cant.T/(pi*cant.Q*cant.f_c*dt_force))
    if seed > 0
        srand(seed)
    end
    N_force = round(Int, T/dt_force)+1
    N = round(Int, (T+T_pre)/dt_force)+1

    force(t) = t >= 0. ? randn() * sigma : 0.
    times = linspace(-T_pre, T, N)
    F_spline = Dierckx.Spline1D(times, map(force, times),
                              k=3)

    function F(t)
        Dierckx.evaluate(F_spline, t)
    end

    Fx, F_tx = ss_cant_1tc(cant, F, V, U)
    jac_t = jacobian(Fx)

    t, _y = ODE.ode23s(F_tx, y0, [-T_pre, T], jacobian=jac_t; maxstep=dt*0.5, args...)
    y = squeeze_y(_y)

    ts = TimeSeries(cant, y0, t, y)

    spline = Dierckx.Spline1D(ts.t, ts.x, k=3)
    t = collect(-T_pre:dt:T)
    x = Dierckx.evaluate(spline, t)

    li = lockin.LockIn(t, x, fs)

    li[:lock2](fp_ratio=0.1, fc_ratio=0.4,
               coeff_ratio=coeff_ratio, window="blackman")
    li[:phase](tf=0.)
    factor = trunc(fs/cant.f_c)
    li[:decimate](factor)

    down = DownsampledTimeSeries(cant, y0, dt, li[:dec_fs],
                          li[:dec_f0], li[:mb][2], li[:dec_t],
                          li[:dec_A], li[:dec_phi], li[:dec_df])

    return down

end

function chunk_sscant_sim_noise(cant::CantileverParams, y0::Vector{Float64},
                             T_pre::Float64,
                             T::Float64, dt::Float64;
                              V=u, U=u, coeff_ratio=8.,
                             Tchunk=1000., f_dec=nothing,
                             seed=-1, args...)
    fs = 1/dt
    k_B = 1.38064852e-11

    dt_force = 0.05/cant.f_c
    sigma = sqrt(cant.k*k_B*cant.T/(pi*cant.Q*cant.f_c*dt_force))

    if seed > 0
        srand(seed)
    end

    N_force = round(Int, T/dt_force)+1
    N = round(Int, (T+T_pre)/dt_force)+1

    force(t) = t >= 0. ? randn() * sigma : 0.
    times = linspace(-T_pre, T, N)
    F_spline = Dierckx.Spline1D(times, map(force, times),
                              k=3)

    function F(t)
        Dierckx.evaluate(F_spline, t)
    end

    if f_dec == nothing
        f_dec = 2*cant.f_c
    end

    dec = trunc(Int, fs / f_dec)

    F_tx, jac = make_Ftx_jac(ss_cant_1tc, cant, F, V, U)

    Tendpts = linspace(-T_pre, T, ceil(Int, (T_pre+T)/Tchunk))
    for i = 1:length(Tendpts)-1
        t, _y = ODE.ode23s(F_tx, y0, [Tendpts[i], Tendpts[i+1]],
                           jacobian=jac; maxstep=dt*0.5)
        y = squeeze_y(_y)

        y0 = _y[end]
        spline = Dierckx.Spline1D(t, y[:, 1], k=3)
        t_equal = collect(Tendpts[i]:dt:Tendpts[i+1])
        x = Dierckx.evaluate(spline, t_equal)

        if i == 1
            global li = lockin.LockIn(t_equal, x, fs)
            li[:lock2](fp_ratio=0.1, fc_ratio=0.4,
                       coeff_ratio=coeff_ratio, window="blackman")
            li[:phase](tf=0.)
            global slock = lockin.FIRStateLock(li[:fir], dec, li[:f0],
                                        -li[:mb][end], t0=-T_pre, fs=fs)
            slock[:filt](x)
        else
            slock[:filt](x[2:end])
        end

    end

    tout = slock[:get_t]()
    A = abs(slock[:z_out])
    phi = DSP.unwrap(angle(slock[:z_out]))
    df = zeros(phi)
    df[2:end] = diff(phi) / (2pi)
    df[1] = df[2]

    DownsampledTimeSeries(cant, y0, dt, fs / dec,
                          slock[:f0], slock[:phi0], tout,
                          A, phi, df)
end

function chunk_all_sscant_sim_noise_outall(cant::CantileverParams, y0::Vector{Float64},
                             T_pre::Float64,
                             T::Float64, dt::Float64; xdot=ss_cant_1tc,
                              V=u, U=u, coeff_ratio=8., V_append=0.01, U_append=0.005,
                             Tchunk=1000., f_dec::Float64=2*cant.f_c,
                             seed=-1, force_ratio=20.,
                             fp=nothing,
                             fc=nothing,
                             progressbar=nothing)
    fs = 1/dt
    const k_B = 1.38064852e-11

    dt_force = 1./(cant.f_c * force_ratio)

    sigma = sqrt(cant.k*k_B*cant.T/(pi*cant.Q*cant.f_c*dt_force))

    if seed > 0
        srand(seed)
    end

    force(t) = randn() * sigma

    Ndec = trunc(Int, fs / f_dec)

    Tendpts = collect(-T_pre:Tchunk:T)
    if Tendpts[end] != T
        push!(Tendpts, T)
    end

    t_save = Vector{Float64}()
    y_save = Vector{Vector{Float64}}()
    delta_y_save = [cant.C_b*V_append, cant.C_1*V_append, U_append]

    for i = 1:length(Tendpts)-1
        N = round(Int, (Tendpts[i+1] - Tendpts[i])/dt_force) + 1
        t_force = linspace(Tendpts[i], Tendpts[i+1], N)
        F_spline = Dierckx.Spline1D(t_force, map(force, t_force),
                              k=3)

        F(t) = Dierckx.evaluate(F_spline, t)

        F_tx, jac = make_Ftx_jac(xdot, cant, F, V, U)

        t, _y = ODE.ode23s(F_tx, y0, [Tendpts[i], Tendpts[i+1]],
                           jacobian=jac; maxstep=dt*0.5)
        y = squeeze_y(_y)

        y0 = _y[end]

        spline = Dierckx.Spline1D(t, y[:, 1], k=3)

        # may need to revert these changes for some unknown reason.

        t_equal = collect(Tendpts[i]:dt:Tendpts[i+1])
        x = Dierckx.evaluate(spline, t_equal)

        if i == 1
            push!(t_save, t[1])
            push!(y_save, _y[1])
            global li = lockin.LockIn(t_equal, x, fs)
            if (fp != nothing) && (fc != nothing)
                li[:lock2](fp=fp, fc=fc,
                       coeff_ratio=coeff_ratio,
                       window="blackman")
            else
                li[:lock2](fp_ratio=0.1, fc_ratio=0.4,
                       coeff_ratio=coeff_ratio,
                       window="blackman")
            end
            li[:phase](tf=0.)

            global slock = lockin.FIRStateLock(li[:fir], Ndec, li[:f0corr],
                                        -li[:mb][end], t0=-T_pre, fs=fs)

            slock[:filt](x)
        else
            slock[:filt](x[2:end])
        end

        for i = eachindex(_y)
            if any(abs(_y[i][3:end] - y_save[end][3:end]) .> delta_y_save)
                push!(t_save, t[i])
                push!(y_save, _y[i])
            end
        end

        if progressbar != nothing
            next!(progressbar)
        end

    end

    y_save_sq = squeeze_y(y_save)

    tout = slock[:get_t]()
    A = abs(slock[:z_out])
    phi = DSP.unwrap(angle(slock[:z_out]))
    df = zeros(phi)
    df[2:end] = diff(phi) / (2pi)
    df[1] = df[2]

    DownsampledTimeSeriesRaw(cant, y0, t_save, y_save_sq, dt, fs / Ndec,
                          slock[:f0], slock[:phi0], tout,
                          A, phi, df)
end

function chunk_all_sscant_sim_noise(cant::CantileverParams, y0::Vector{Float64},
                             T_pre::Float64,
                             T::Float64, dt::Float64; xdot=ss_cant_1tc,
                              V=u, U=u, coeff_ratio=8.,
                             Tchunk=1000., f_dec=nothing,
                             seed=-1, force_ratio=20., progressbar=nothing)
    fs = 1/dt
    k_B = 1.38064852e-11

    dt_force = 1./(cant.f_c * force_ratio)

    sigma = sqrt(cant.k*k_B*cant.T/(pi*cant.Q*cant.f_c*dt_force))

    if seed > 0
        srand(seed)
    end

    force(t) = t >= 0. ? randn() * sigma : 0.

    if f_dec == nothing
        f_dec = fs
    end

    dec = trunc(Int, fs / f_dec)

    Tendpts = collect(-T_pre:Tchunk:T)
    if Tendpts[end] != T
        push!(Tendpts, T)
    end

    N_Tendpts = length(Tendpts)-1

    if progressbar == true
        progressbar = Progress(N_Tendpts)
    end

    for i = 1:N_Tendpts
        N = round(Int, (Tendpts[i+1] - Tendpts[i])/dt_force) + 1
        t_force = linspace(Tendpts[i], Tendpts[i+1], N)
        F_spline = Dierckx.Spline1D(t_force, map(force, t_force),
                              k=3)
        function F(t)
            Dierckx.evaluate(F_spline, t)
        end
        F_tx, jac = make_Ftx_jac(xdot, cant, F, V, U)

        t, _y = ODE.ode23s(F_tx, y0, [Tendpts[i], Tendpts[i+1]],
                           jacobian=jac, maxstep=dt*0.5)
        y = squeeze_y(_y)

        y0 = _y[end]

        spline = Dierckx.Spline1D(t, y[:, 1], k=3)
        t_equal = collect(Tendpts[i]:dt:Tendpts[i+1])

        x = Dierckx.evaluate(spline, t_equal)

        if i == 1
            li = lockin.LockIn(t_equal, x, fs)
            li[:lock2](fp_ratio=0.1, fc_ratio=0.4,
                       coeff_ratio=coeff_ratio, window="blackman")

            li[:phase]()
            global slock = lockin.FIRStateLock(li[:fir], dec, li[:f0],
                                        -li[:mb][end], t0=-T_pre, fs=fs)
            slock[:filt](x)
        else
            slock[:filt](x[2:end])
        end

        if progressbar != nothing
            next!(progressbar)
        end

    end

    tout = slock[:get_t]()
    A = abs(slock[:z_out])
    phi = DSP.unwrap(angle(slock[:z_out]))
    df = zeros(phi)
    df[2:end] = diff(phi) / (2pi)
    df[1] = df[2]

    DownsampledTimeSeries(cant, y0, dt, fs / dec,
                          slock[:f0], slock[:phi0], tout,
                          A, phi, df)
end

function chunk_all_sscant_sim_brownian(cant::CantileverParams, y0::Vector{Float64},
                             T_pre::Float64,
                             T::Float64, dt::Float64, dt_out::Float64;
                             xdot=ss_cant_1tc,
                             V=u, U=u, coeff_ratio=8.,
                             Tchunk=1000.,
                             seed=-1, force_ratio=20., progressbar=nothing)
    fs = 1/dt
    k_B = 1.38064852e-11

    dt_force = 1./(cant.f_c * force_ratio)

    sigma = sqrt(cant.k*k_B*cant.T/(pi*cant.Q*cant.f_c*dt_force))

    if seed > 0
        srand(seed)
    end

    force(t) = t >= 0. ? randn() * sigma : 0.

    Tendpts = collect(-T_pre:Tchunk:T)
    if Tendpts[end] != T
        push!(Tendpts, T)
    end

    N_Tendpts = length(Tendpts)-1

    if progressbar == true
        progressbar = Progress(N_Tendpts)
    end

    t_collected = Vector{Float64}()
    x_collected = Vector{Float64}()

    Npts_expected = round(Int, (T+t_pre)/dt)

    sizehint!(t_collected, Npts_expected)
    sizehint!(x_collected, Npts_expected)

    for i = 1:N_Tendpts

        N = round(Int, (Tendpts[i+1] - Tendpts[i])/dt_force) + 1
        t_force = linspace(Tendpts[i], Tendpts[i+1], N)

        F_spline = Dierckx.Spline1D(t_force, map(force, t_force),
                              k=3)

        F(t) = Dierckx.evaluate(F_spline, t)

        F_tx, jac = make_Ftx_jac(xdot, cant, F, V, U)

        t, _y = ODE.ode23s(F_tx, y0, [Tendpts[i], Tendpts[i+1]],
                           jacobian=jac; maxstep=dt*0.5)
        y = squeeze_y(_y)

        y0 = _y[end]

        spline = Dierckx.Spline1D(t, y[:, 1], k=3)
        t_equal = collect(Tendpts[i]:dt_out:Tendpts[i+1])
        x = Dierckx.evaluate(spline, t_equal)

        append!(t_collected, t_equal[2:end])
        append!(x_collected, x[2:end])

        if progressbar != nothing
            next!(progressbar)
        end

    end

    t_collected, x_collected, y0
end

end # module


text:
module DIIS
using LinearAlgebra
import Base.push!
import Base.length

export DIISManager

"""
    DIISManager{T1<:AbstractFloat,
                T2<:AbstractFloat}

Structure holding arrays necessary for extrapolation with the `DIIS` method.

# Fields
| Name | Description |
|:-----|:------------|
|`vecs`    |  Vectors to extrapolate from |
|`errs`    |  Error vectors to build B matrix from |
|`max_vec` |  Max number of vectors to hold |
"""
struct DIISManager{T1<:AbstractFloat,T2 <: AbstractFloat}
    vecs::Array{AbstractArray{T1},1}
    errs::Array{AbstractArray{T2},1}
    max_vec::Int64
end

function DIISManager{T1,T2}(;size=6) where { T1 <: AbstractFloat,
                                             T2 <: AbstractFloat }
    vecs = Array{AbstractArray{T1}}(undef,0)
    errs = Array{AbstractArray{T2}}(undef,0)
    DIISManager{T1,T2}(vecs,errs,size)
end

function length(M::DIISManager)
    length(M.vecs)
end

function push!(M::DIISManager{T1,T2}, V::AbstractArray, E::AbstractArray) where { T1 <: AbstractFloat,
                                                                  T2 <: AbstractFloat }
    # If the number of vectors stored has reach its maximum
    # the vector associated with the biggest error is replaced with the new one
    if length(M)+1 > M.max_vec
        norms = norm.(M.errs)
        _, idx = findmax(norms)
        deleteat!(M.vecs,idx)
        deleteat!(M.errs,idx)
    end
    push!(M.vecs,deepcopy(V))
    push!(M.errs,deepcopy(E))
end

"""
    Fermi.DIIS.extrapolate(M::DIISManager{T1,T2}; add_res=false) where { T1 <: AbstractFloat,
                                                              T2 <: AbstractFloat }

Produces a new guess vector using *direct inversion in the iterative subspace* given
the information in `M` which is a `DIISManager`.
"""
function extrapolate(M::DIISManager{T1,T2}; add_res=false) where {T1 <: AbstractFloat, T2 <: AbstractFloat}

    if length(M) < 2
        throw(BoundsError("DIIS cannot extrapolate from one vector"))
    end

    # Solves the equation for the new vector
    N = length(M) + 1

    # Create the B matrix
    B = ones(T1, N, N)
    B[end,end] = zero(T1)
    for (n1, e1) in enumerate(M.errs[1:end])
        for (n2, e2) in enumerate(M.errs[1:end])
            # Bij = (ei⋅ej)
            B[n1,n2] = sum(e1 .* e2)
        end
    end

    # Create the residual vector
    resid = zeros(T1, N)
    resid[end] = one(T1)

    # Get coefficients
    ci = B\resid
    out = zero(M.vecs[1])

    # Compute the new vector as Pnew = ∑ ci * Pi
    add_res ? outw = zeros(T1,size(M.vecs[1])) : nothing
    out = sum(ci[1:end-1] .* M.vecs)    # Need to pop the last element as it is λ not ci

    # Add this extra stuff if requested
    if add_res
        out += sum(ci[1:end-1] .* M.errs)
    end

    return out
end
end #module

ext: jl

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
"""
    Phases
"""
module Phases

abstract type AbstractFittingPhase end
abstract type AbstractTrainingPhase <: AbstractFittingPhase end

struct TrainingPhase <: AbstractTrainingPhase end

struct ValidationPhase <: AbstractFittingPhase end
struct TestPhase <: AbstractFittingPhase end
struct InitializationPhase <: AbstractFittingPhase end
struct CleanupPhase <: AbstractFittingPhase end

export
    AbstractFittingPhase,
    AbstractTrainingPhase,
    TrainingPhase,
    ValidationPhase,
    TestPhase,
    InitializationPhase,
    CleanupPhase
end # module

using .Phases


text:
# Context management

export
    CuContext, destroy!, CuCurrentContext, activate,
    synchronize, device

## construction and destruction

"""
    CuContext(dev::CuDevice, flags::CUctx_flags=CTX_SCHED_AUTO)
    CuContext(f::Function, ...)

Create a CUDA context for device. A context on the GPU is analogous to a process on the CPU,
with its own distinct address space and allocated resources. When a context is destroyed,
the system cleans up the resources allocated to it.

Contexts are unique instances which need to be `destroy`ed after use. For automatic
management, prefer the `do` block syntax, which implicitly calls `destroy`.
"""
mutable struct CuContext
    handle::CUcontext
    owned::Bool
    valid::Bool

    """
    The `owned` argument indicates whether the caller owns this context. If so, the context
    should get destroyed when it goes out of scope. If not, it is up to the caller to do so.
    """
    function CuContext(handle::CUcontext, owned=true)
        handle == C_NULL && return new(C_NULL, false, true)

        # we need unique context instances for garbage collection reasons
        #
        # refcounting the context handle doesn't work, because having multiple instances can
        # cause the first instance (if without any use) getting gc'd before a new instance
        # (eg. through getting the current context, _with_ actual uses) is created
        #
        # instead, we force unique instances, and keep a reference alive in a global dict.
        # this prevents contexts from getting collected, requiring the user to destroy it.
        ctx = get!(context_instances, handle) do
            obj = new(handle, owned, true)
            finalizer(unsafe_destroy!, obj)
            return obj
        end

        if owned && !ctx.owned
            # trying to get a non-owned handle on an already-owned context is common,
            # eg. construct a context, and call a function doing CuCurrentContext()
            #
            # the inverse isn't true: constructing an owning-object on using a handle
            # which was the result of a non-owning API call really shouldn't happen
            warn("Ownership conflict on context $ctx")
        end

        ctx
    end
end
const context_instances = Dict{CUcontext,CuContext}()

isvalid(ctx::CuContext) = ctx.valid
function invalidate!(ctx::CuContext)
    ctx.valid = false
    nothing
end

function unsafe_destroy!(ctx::CuContext)
    # finalizers do not respect _any_ oder during process teardown
    # (ie. it doesn't respect active instances carefully set-up in `gc.jl`)
    # TODO: can we check this only happens during teardown?
    if ctx.owned && isvalid(ctx)
        cuCtxDestroy(ctx)
        invalidate!(ctx)
    end
end

Base.unsafe_convert(::Type{CUcontext}, ctx::CuContext) = ctx.handle

Base.:(==)(a::CuContext, b::CuContext) = a.handle == b.handle
Base.hash(ctx::CuContext, h::UInt) = hash(ctx.handle, h)

"""
    destroy!(ctx::CuContext)

Mark a context for destruction.

This does not immediately destroy the context, as there might still be dependent resources
which have not been collected yet. The context will get freed as soon as all outstanding
instances have been finalized.
"""
function destroy!(ctx::CuContext)
    delete!(context_instances, ctx.handle)
    return
end

Base.deepcopy_internal(::CuContext, ::IdDict) =
    error("CuContext cannot be copied")

@enum_without_prefix CUctx_flags CU_

function CuContext(dev::CuDevice, flags::CUctx_flags=CTX_SCHED_AUTO)
    handle_ref = Ref{CUcontext}()
    cuCtxCreate(handle_ref, flags, dev)
    CuContext(handle_ref[])
end

"""
    CuCurrentContext()

Return the current context, or `nothing` if there is no active context.
"""
function CuCurrentContext()
    handle_ref = Ref{CUcontext}()
    cuCtxGetCurrent(handle_ref)
    if handle_ref[] == C_NULL
        return nothing
    else
        return CuContext(handle_ref[], false)
    end
end

"""
    push!(CuContext, ctx::CuContext)

Pushes a context on the current CPU thread.
"""
Base.push!(::Type{CuContext}, ctx::CuContext) =
    cuCtxPushCurrent(ctx)

"""
    pop!(CuContext)

Pops the current CUDA context from the current CPU thread, and returns that context.
"""
function Base.pop!(::Type{CuContext})
    handle_ref = Ref{CUcontext}()
    cuCtxPopCurrent(handle_ref)
    CuContext(handle_ref[], false)
end

"""
    activate(ctx::CuContext)

Binds the specified CUDA context to the calling CPU thread.
"""
activate(ctx::CuContext) = cuCtxSetCurrent(ctx)

function CuContext(f::Function, args...)
    ctx = CuContext(args...)    # implicitly pushes
    try
        f(ctx)
    finally
        @assert pop!(CuContext) == ctx
        destroy!(ctx)
    end
end

## context properties

"""
    device()
    device(ctx::Cucontext)

Returns the device for a context.
"""
function device(ctx::CuContext)
    push!(CuContext, ctx)
    device_ref = Ref{CUdevice}()
    cuCtxGetDevice(device_ref)
    pop!(CuContext)
    return CuDevice(Bool, device_ref[])
end
function device()
    device_ref = Ref{CUdevice}()
    cuCtxGetDevice(device_ref)
    return CuDevice(Bool, device_ref[])
end

"""
    synchronize()

Block for the current context's tasks to complete.
"""
synchronize() = cuCtxSynchronize()

## cache config

export cache_config, cache_config!

@enum_without_prefix CUfunc_cache CU_

function cache_config()
    config = Ref{CUfunc_cache}()
    cuCtxGetCacheConfig(config)
    return config[]
end

function cache_config!(config::CUfunc_cache)
    cuCtxSetCacheConfig(config)
end

## shared memory config

export shmem_config, shmem_config!

@enum_without_prefix CUsharedconfig CU_

function shmem_config()
    config = Ref{CUsharedconfig}()
    cuCtxGetSharedMemConfig(config)
    return config[]
end

function shmem_config!(config::CUsharedconfig)
    cuCtxSetSharedMemConfig(config)
end

## limits

export limit, limit!

@enum_without_prefix CUlimit CU_

function limit(lim::CUlimit)
    val = Ref{Csize_t}()
    cuCtxGetLimit(val, lim)
    return Int(val[])
end

limit!(lim::CUlimit, val) = cuCtxSetLimit(lim, val)


text:
# Type definition
struct FourierSeries{T,S<:Real} <: FunctionSpace{T}
    coeffs::Vector{T}
    period::S
    function FourierSeries(coeffs::Vector{T}, period=1) where {T}
        isodd(length(coeffs)) ||
            throw(ArgumentError("expect odd length for coefficient vector"))
        S = typeof(period)
        new{T,S}(coeffs, period)
    end
end

# Basic properties
period(f::FourierSeries) = f.period
domain(f::FourierSeries) = (zero(period(f)), period(f))

nummodes(F::FourierSeries) = (length(F.coeffs)-1) >> 1
coefficients(F::FourierSeries) = F.coeffs

function Base.:(==)(F1::FourierSeries, F2::FourierSeries)
    F1.period == F2.period || return false
    K = max(nummodes(F1), nummodes(F2))
    for k = -K:K
        F1[k] == F2[k] || return false
    end
    return true
end

# Indexing, getting and setting coefficients
Base.eachindex(F::FourierSeries) = UnitRange(-nummodes(F), nummodes(F))
function Base.getindex(F::FourierSeries, k)
    if abs(k) > nummodes(F)
        return zero(F.coeffs[1])
    else
        j = ifelse(k == 0, 1, ifelse(k < 0, -2*k+1, 2*k))
        @inbounds getindex(F.coeffs, j)
    end
end
function Base.setindex!(F::FourierSeries, v, k)
    if abs(k) > nummodes(F)
        setnummodes!(F, abs(k))
    end
    j = ifelse(k == 0, 1, ifelse(k < 0, -2*k+1, 2*k))
    @inbounds setindex!(F.coeffs, v, j)
end

# Use as function
function (F::FourierSeries)(x)
    xred = x/period(F)
    x = mod1(xred, one(xred))
    K = nummodes(F)
    factor = 2pi
    f = zero(F[0]) + F[0]*(cos(0*factor)+im*sin(0*factor)) # for type stability
    for k = 1:K
        f += (F[k]+F[-k])*cos(factor*k*x) + im*(F[k]-F[-k])*sin(factor*k*x)
    end
    return f
end

# Change number of coefficients
function truncate!(F::FourierSeries; Kmax::Integer = nummodes(F), tol::Real = 0)
    Ktol = findlast(x->norm(x)>=tol, F.coeffs)
    Kmax = min(Kmax, Ktol === nothing ? 0 : Ktol>>1)
    if Kmax < nummodes(F)
        resize!(F.coeffs, 2*Kmax+1)
    end
    return F
end

function setnummodes!(F::FourierSeries, K::Int)
    while nummodes(F) < K
        push!(F.coeffs, zero(F.coeffs[1]), zero(F.coeffs[1]))
    end
    if nummodes(F) > K
        resize!(F.coeffs, 2*K+1)
    end
    return F
end

# special purpose constructor
Base.similar(F::FourierSeries, ::Type{T} = scalartype(F)) where {T} =
    FourierSeries([zero(T)*F[0]], period(F))

Base.zero(F::FourierSeries) = FourierSeries([zero(F[0])], period(F))
Base.one(F::FourierSeries) = FourierSeries([one(F[0])], period(F))

# Arithmetic (out of place)
Base.copy(F::FourierSeries) = FourierSeries(copy.(coefficients(F)), period(F))

Base.:-(F::FourierSeries) = FourierSeries(.-coefficients(F), period(F))
Base.:+(F::FourierSeries) = FourierSeries(.+coefficients(F), period(F))

const Const = Union{Number,AbstractArray}
Base.:*(F::FourierSeries, a::Const) =
    FourierSeries([f*a for f in coefficients(F)], period(F))
Base.:*(a::Const, F::FourierSeries) =
    FourierSeries([a*f for f in coefficients(F)], period(F))
Base.:/(F::FourierSeries, a::Const) =
    FourierSeries([f/a for f in coefficients(F)], period(F))
Base.:\(a::Const, F::FourierSeries) =
    FourierSeries([a\f for f in coefficients(F)], period(F))

function Base.:+(F1::FourierSeries, F2::FourierSeries)
    domain(F1) == domain(F2) || throw(DomainMismatch())
    K1 = nummodes(F1)
    K2 = nummodes(F2)
    K = max(K1, K2)
    F = FourierSeries(sizehint!([F1[0]+F2[0]], 2*K+1), period(F1))
    for k = 1:K
        if k <= K1 && k <= K2
            F[k] = F1[k] + F2[k]
            F[-k] = F1[-k] + F2[-k]
        elseif k <= K1
            F[k] = copy(F1[k])
            F[-k] = copy(F1[-k])
        else
            F[k] = copy(F2[k])
            F[-k] = copy(F2[-k])
        end
    end
    return F
end

function Base.:-(F1::FourierSeries, F2::FourierSeries)
    domain(F1) == domain(F2) || throw(DomainMismatch())
    K1 = nummodes(F1)
    K2 = nummodes(F2)
    K = max(K1, K2)
    F = FourierSeries(sizehint!([F1[0] - F2[0]], 2*K+1), period(F1))
    for k = 1:K
        if k <= K1 && k <= K2
            F[k] = F1[k] - F2[k]
            F[-k] = F1[-k] - F2[-k]
        elseif k <= K1
            F[k] = copy(F1[k])
            F[-k] = copy(F1[-k])
        else
            F[k] = -(F2[k])
            F[-k] = -(F2[-k])
        end
    end
    return F
end

Base.:*(F1::FourierSeries, F2::FourierSeries) = truncmul(F1, F2)

function truncmul(F1::FourierSeries, F2::FourierSeries;
                    Kmax::Integer = nummodes(F1) + nummodes(F2), tol::Real = 0)
    domain(F1) == domain(F2) || throw(DomainMismatch())
    F = FourierSeries(sizehint!([zero(F1[0])*zero(F2[0])], 2*Kmax+1), period(F1))
    truncmul!(F, F1, F2, true, false; Kmax = Kmax, tol = tol)
    return F
end

function Base.conj(F::FourierSeries)
    Fc = FourierSeries(map(conj, coefficients(F)), period(F))
    K = nummodes(F)
    for k = 1:K
        Fc[k], Fc[-k] = Fc[-k], Fc[k]
    end
    return Fc
end
function Base.adjoint(F::FourierSeries)
    Fc = FourierSeries(map(adjoint, coefficients(F)), period(F))
    K = nummodes(F)
    for k = 1:K
        Fc[k], Fc[-k] = Fc[-k], Fc[k]
    end
    return Fc
end
Base.transpose(F::FourierSeries) = FourierSeries(map(transpose, coefficients(F)), period(F))

LinearAlgebra.tr(F::FourierSeries) = FourierSeries(map(tr, coefficients(F)), period(F))

Base.real(F::FourierSeries) = rmul!(F+conj(F), 1//2)
Base.imag(F::FourierSeries) = rmul!(F-conj(F), 1//(2*im))

# Arithmetic (in place / mutating methods)
function Base.copy!(Fdst::FourierSeries, Fsrc::FourierSeries)
    domain(Fdst) == domain(Fsrc) || throw(DomainMismatch())
    setnummodes!(Fdst, nummodes(Fsrc))
    if eltype(Fdst) <: Number
        copy!(coefficients(Fdst), coefficients(Fsrc))
    else
        for k in eachindex(Fdst)
            copy!(Fdst[k], Fsrc[k])
        end
    end
    return Fdst
end
function LinearAlgebra.rmul!(F::FourierSeries, α::Number)
    if eltype(F) <: Number
        rmul!(coefficients(F), α)
    else
        for k in eachindex(F)
            rmul!(F[k], α)
        end
    end
    return F
end
function LinearAlgebra.lmul!(α::Number, F::FourierSeries)
    if eltype(F) <: Number
        lmul!(α, coefficients(F))
    else
        for k in eachindex(F)
            lmul!(α, F[k])
        end
    end
    return F
end

function LinearAlgebra.axpy!(α::Number, Fx::FourierSeries, Fy::FourierSeries)
    domain(Fx) == domain(Fy) || throw(DomainMismatch())
    Kx = nummodes(Fx)
    setnummodes!(Fy, max(Kx, nummodes(Fy)))
    Ky = nummodes(Fy)
    if eltype(Fy) <: Number
        if Ky > Kx
            LinearAlgebra.axpy!(α, coefficients(Fx), view(coefficients(Fy), 1:2Kx+1))
        else
            LinearAlgebra.axpy!(α, coefficients(Fx), coefficients(Fy))
        end
    else
        K = nummodes(Fx)
        Threads.@threads for k in -K:K
            axpy!(α, Fx[k], Fy[k])
        end
    end
    return Fy
end
function LinearAlgebra.axpby!(α::Number, Fx::FourierSeries, β::Number, Fy::FourierSeries)
    domain(Fx) == domain(Fy) || throw(DomainMismatch())
    Kx = nummodes(Fx)
    setnummodes!(Fy, max(Kx, nummodes(Fy)))
    Ky = nummodes(Fy)
    if eltype(Fy) <: Number
        if Ky > Kx
            LinearAlgebra.axpby!(α, coefficients(Fx), β, view(coefficients(Fy), 1:2Kx+1))
            lmul!(β, view(coefficients(Fy), (2Kx+2):(2Ky+1)))
        else
            LinearAlgebra.axpby!(α, coefficients(Fx), β, coefficients(Fy))
        end
    else
        K = nummodes(Fy)
        Threads.@threads for k in -K:K
            if abs(k) <= Kx
                axpby!(α, Fx[k], β, Fy[k])
            else
                lmul!(β, Fy[k])
            end
        end
    end
    return Fy
end

LinearAlgebra.mul!(Fdst::FourierSeries, α::Number, Fsrc::FourierSeries) =
    truncmul!(Fdst, α, Fsrc)
LinearAlgebra.mul!(Fdst::FourierSeries, Fsrc::FourierSeries, α::Number) =
    truncmul!(Fdst, Fsrc, α)
LinearAlgebra.mul!(F::FourierSeries, F1::FourierSeries, F2::FourierSeries,
                    α = true, β = false) = truncmul!(F, F1, F2, α, β)

function truncmul!(Fdst::FourierSeries, α₁::Number, Fsrc::FourierSeries,
                    α₂ = true, β = false;
                    Kmax::Integer = nummodes(Fsrc), tol::Real = 0)
    α = α₁ * α₂
    domain(Fdst) == domain(Fsrc) || throw(DomainMismatch())
    setnummodes!(Fdst, Kmax)
    if eltype(Fdst) <: Number
        mul!(coefficients(Fdst), α, coefficients(Fsrc))
    else
        for k in eachindex(Fdst)
            mul!(Fdst[k], α, Fsrc[k])
        end
    end
    return Fdst
end
function truncmul!(Fdst::FourierSeries, Fsrc::FourierSeries, α₁::Number,
                    α₂ = true, β = false;
                    Kmax::Integer = nummodes(Fsrc), tol::Real = 0)
    α = α₁ * α₂
    domain(Fdst) == domain(Fsrc) || throw(DomainMismatch())
    setnummodes!(Fdst, Kmax)
    if eltype(Fdst) <: Number
        mul!(coefficients(Fdst), coefficients(Fsrc), α)
    else
        for k in eachindex(Fdst)
            mul!(Fdst[k], Fsrc[k], α)
        end
    end
    return Fdst
end
function truncmul!(F::FourierSeries, F1::FourierSeries, F2::FourierSeries,
                    α = true, β = false;
                    Kmax::Integer = max(nummodes(F), nummodes(F1)+nummodes(F2)), tol::Real = 0)
    domain(F) == domain(F1) == domain(F2) || throw(DomainMismatch())
    K1 = nummodes(F1)
    K2 = nummodes(F2)
    K = min(Kmax, max(iszero(β) ? 0 : nummodes(F), K1+K2))
    setnummodes!(F, K)
    Threads.@threads for k = -K:K
        fk = F[k]
        if fk isa AbstractArray
            T = eltype(fk)
            if β == 0
                fill!(fk, T(β))
            elseif β != 1
                rmul!(fk, T(β))
            end
            for k1 = max(-K1, k-K2):min(K1, k+K2)
                k2 = k - k1
                fk1 = F1[k1]
                fk2 = F2[k2]
                if fk1 isa AbstractArray && fk2 isa AbstractArray
                    mul!(fk, fk1, fk2, T(α), T(true))
                elseif fk1 isa AbstractArray
                    axpy!(T(fk2*α), fk1, fk)
                elseif fk2 isa AbstractArray
                    axpy!(T(fk1*α), fk2, fk)
                else
                    @warn "unexpected branch"
                    fk .= fk + (fk1*fk2*α)
                end
            end
        else
            fk *= β
            for k1 = max(-K1, k-K2):min(K1, k+K2)
                k2 = k - k1
                fk1 = F1[k1]
                fk2 = F2[k2]
                fk += fk1*fk2*α
            end
            F[k] = fk
        end
    end
    if tol != 0
        truncate!(F; tol = tol)
    end
    return F
end

# Inner product and norm
function localdot(F1::FourierSeries, F2::FourierSeries)
    domain(F1) == domain(F2) || throw(DomainMismatch())
    K1 = nummodes(F1)
    K2 = nummodes(F2)
    K = K1+K2
    F = FourierSeries(sizehint!([zero(dot(F1[0], F2[0]))], 2*K+1), period(F1))
    for k = -K:K
        fk = F[k]
        for k1 = max(-K1, -K2-k):min(K1, K2-k)
            k2 = k + k1
            fk += dot(F1[k1], F2[k2])
        end
        F[k] = fk
    end
    return F
end

function LinearAlgebra.dot(F1::FourierSeries, F2::FourierSeries)
    domain(F1) == domain(F2) || throw(DomainMismatch())
    K = min(nummodes(F1), nummodes(F2))
    return sum(k->dot(F1[k], F2[k]), -K:K)
end
LinearAlgebra.norm(F::FourierSeries) = norm(norm(Fk) for Fk in coefficients(F))

function LinearAlgebra.isapprox(x::FourierSeries, y::FourierSeries;
                                atol::Real=0,
                                rtol::Real=defaulttol(x[0]))
    return norm(x-y) <= max(atol, rtol*max(norm(x), norm(y)))
end

function differentiate(F::FourierSeries)
    dF = 0im * F
    ω = (2*pi)/period(F)
    K = nummodes(F)
    for k = 1:K
        f = im*ω*k
        dF[-k] = -f*F[-k]
        dF[+k] = +f*F[+k]
    end
    return dF
end
function integrate(F::FourierSeries, (a,b) = domain(F))
    p = period(F)
    if b-a == p
        return F[0]*p
    else
        error("not yet implemented")
    end
end

# Simple Fourier transform: typically not needed for large number of points, but should
# work with matrix valued functions etc
fit(f, ::Type{FourierSeries}, (a,b)::Tuple{Real,Real}; kwargs...) =
    fit(f, FourierSeries, b-a; kwargs...)

function fit(f, ::Type{FourierSeries}, period = 1; Kmax = 10, tol = 1e-12)
    K = Kmax
    x = (0:2*K) * (period / (2K+1))
    fx = map(f, x)
    ω = 2*pi*0/period
    integrand = exp.((-im*ω) .* x) .* fx
    F = FourierSeries([sum(integrand)/(2K+1)], period)
    isrealf = all(isreal, fx)
    for k = 1:Kmax
        ω = 2*pi*k/period
        integrand .= exp.((-im*ω) .* x) .* fx
        F[k] = sum(integrand)/(2K+1)
        if isrealf
            F[-k] = conj(F[k])
        else
            integrand .= exp.((+im*ω) .* x) .* fx
            F[-k] = sum(integrand)/(2K+1)
        end
    end
    return F
end

# function Base.inv(F::FourierSeries; tol=1e-12, Kmax = 2*nummodes(F))
#     K = nummodes(F)
#     while true
#         Fx = map(F, (0:2*K)/(2*K+1))
#         Gx = map(inv, Fx)
#         krange = vcat([0],[(-1)^iseven(l)*((l+1)>>1) for l=1:2*K])
#         factor = -im*2*pi/(2*K+1)
#         coeffG = let factor=factor, Gx=Gx, K = K
#             map(krange) do k
#                 sum(exp(factor*k*(l-1))*Gx[l] for l = 1:(2*K+1))/(2*K+1)
#             end
#         end
#         G = FourierSeries(coeffG)
#         if norm(G*F-one(G))<tol || K == Kmax
#             return G
#         else
#             K = min(2*K, Kmax)
#         end
#     end
# end
#
# function inv2(F::FourierSeries; tol = 1e-12, Kmax = 2*nummodes(F))
#     Fi = inv(F; tol = tol, Kmax = Kmax)
#     Fi2 = linsolve(x->setnummodes!(x*F, nummodes(x)), one(F), Fi; tol = tol)
#
#     K = nummodes(F)
#     while true
#         Fx = map(F, (0:2*K)/(2*K+1))
#         Gx = map(inv, Fx)
#         krange = vcat([0],[(-1)^iseven(l)*((l+1)>>1) for l=1:2*K])
#         factor = -im*2*pi/(2*K+1)
#         coeffG = let factor=factor, Gx=Gx, K = K
#             map(krange) do k
#                 sum(exp(factor*k*(l-1))*Gx[l] for l = 1:(2*K+1))/(2*K+1)
#             end
#         end
#         G = FourierSeries(coeffG)
#         if norm(G*F-one(G))<tol || K == Kmax
#             return G
#         else
#             K = min(2*K, Kmax)
#         end
#     end
#     return Fi2
# end
#
# function Base.sqrt(F::FourierSeries; tol=1e-12, Kmax = 200)
#     K = nummodes(F)
#     while true
#         Fx = map(F, (0:2*K)/(2*K+1))
#         Gx = map(sqrt, Fx)
#         krange = vcat([0],[(-1)^iseven(l)*((l+1)>>1) for l=1:2*K])
#         factor = -im*2*pi/(2*K+1)
#         coeffG = let factor=factor, Gx=Gx
#             map(krange) do k
#                 sum(exp(factor*k*(l-1))*Gx[l] for l = 1:(2*K+1))/(2*K+1)
#             end
#         end
#         G = FourierSeries(coeffG)
#         if norm(G*G-F)<tol || K == Kmax
#             return G
#         else
#             K = min(2*K, Kmax)
#         end
#     end
# end


text:
include("../CNSE.jl")
include("ThreeDimensionalCompressibleNavierStokesEquations.jl")

function Config(
    name,
    resolution,
    domain,
    params;
    numerical_flux_first_order = RoeNumericalFlux(),
    Nover = 0,
    boundary = (1, 1),
    boundary_conditons = (FluidBC(Impenetrable(FreeSlip()), Insulating()),),
)
    mpicomm = MPI.COMM_WORLD
    ArrayType = ClimateMachine.array_type()

    println(string(resolution.Nᶻ) * " elems in the vertical")

    vert_range =
        grid1d(domain.min_height, domain.max_height, nelem = resolution.Nᶻ)

    println(
        string(resolution.Nʰ) * "x" * string(resolution.Nʰ) * " elems per face",
    )

    topology = StackedCubedSphereTopology(
        mpicomm,
        resolution.Nʰ,
        vert_range;
        boundary = boundary,
    )

    println("poly order is " * string(resolution.N))
    println("OI order is " * string(Nover))

    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = resolution.N + Nover,
        meshwarp = equiangular_cubed_sphere_warp,
    )

    model = CNSE3D{FT}(
        nothing,
        (domain.min_height, domain.max_height),
        ClimateMachine.Orientations.SphericalOrientation(),
        NonLinearAdvectionTerm(),
        ConstantViscosity{FT}(μ = params.μ, ν = params.ν, κ = params.κ),
        nothing,
        nothing,
        boundary_conditons;
        cₛ = params.cₛ,
        ρₒ = params.ρₒ,
    )

    dg = DGModel(
        model,
        grid,
        numerical_flux_first_order,
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    return Config(name, dg, Nover, mpicomm, ArrayType)
end

function cnse_init_aux!(::CNSE3D, aux, geom)
    @inbounds begin
        aux.x = geom.coord[1]
        aux.y = geom.coord[2]
        aux.z = geom.coord[3]
    end

    return nothing
end


text:
# SVD Algorithms
# --------------

# Truncated SVD using implicitly restarted Lonczos iterations.
function tsvd(A::Matrix{T}, k::Integer) where T
    m, n = size(A)
    F, = svds(A, nsv=k)
    U, Σ, V = F.U, F.S, F.V
    @assert size(U, 2) == length(Σ) == size(V, 2) == k
    return U, Σ, V
end

# H. Li, G. C. Linderman, et al. "Algorithm 971: An Implementation of a
# Randomized Algorithm for Principal Component Analysis", ACM Transactions on
# Mathematical Software (TOMS), 2017
# DOI: https://doi.org/10.1145/3004053
function rsvd(A::Matrix{T}, k::Integer; its=3, l=k+5) where T
    m, n = size(A)
    l = min(l, m, n)
    @assert 0 < k ≤ l ≤ min(m, n)
    Q::Matrix{T} = rand(eltype(T), n, l) .- T(0.5)
    Y = A*Q
    F = lu!(Y)
    for i in 1:its
        Y = A'F.L
        F = lu!(Y)
        Y = A*F.L
        if i < its
            F = lu!(Y)
        else
            F = qr!(Y)
        end
    end
    Q = Matrix(F.Q)
    B = Q'A
    W, Σ, V = svd(B)
    U = Q*W
    return U[:,1:k], Σ[1:k], V[:,1:k]
end


text:
# Bar geometry summarizes data as vertical bars.
immutable BarGeometry <: Gadfly.GeometryElement
    # How bars should be positioned if they are grouped by color.
    # Valid options are:
    #   :stack -> place bars on top of each other (default)
    #   :dodge -> place bar next to each other
    position::Symbol

    # :vertical or :horizontal
    orientation::Symbol

    default_statistic::Gadfly.StatisticElement

    function BarGeometry(default_statistic=Gadfly.Stat.identity();
                         position::Symbol=:stack, orientation::Symbol=:vertical)
        new(position, orientation, default_statistic)
    end
end

const bar = BarGeometry

function histogram(; position=:stack, bincount=nothing,
                   minbincount=3, maxbincount=150, orientation::Symbol=:vertical)
    BarGeometry(Gadfly.Stat.histogram(bincount=bincount,
                                      minbincount=minbincount,
                                      maxbincount=maxbincount,
                                      orientation=orientation),
                position=position,
                orientation=orientation)
end

function element_aesthetics(::BarGeometry)
    [:x, :y, :color]
end

function default_statistic(geom::BarGeometry)
    geom.default_statistic
end

# Render a single color bar chart
function render_colorless_bar(geom::BarGeometry,
                              theme::Gadfly.Theme,
                              aes::Gadfly.Aesthetics,
                              orientation::Symbol)
    bar_form = empty_form
    if orientation == :horizontal
        for (y_min, y_max, x) in zip(aes.ymin, aes.ymax, aes.x)
            geometry_id = Gadfly.unique_svg_id()
            bar_form = combine(bar_form,
                compose(rectangle(0.0, y_min*cy + theme.bar_spacing/2,
                                  x, (y_max - y_min)*cy - theme.bar_spacing),
                        svgid(geometry_id), svgclass("geometry")))
        end
    else
        for (x_min, x_max, y) in zip(aes.xmin, aes.xmax, aes.y)
            geometry_id = Gadfly.unique_svg_id()
            bar_form = combine(bar_form,
                compose(rectangle(x_min*cx + theme.bar_spacing/2, 0.0,
                                  (x_max - x_min)*cx - theme.bar_spacing, y),
                        svgid(geometry_id), svgclass("geometry")))
        end
    end
    bar_form
end

# Render a bar chart grouped by discrete colors and stacked.
function render_colorful_stacked_bar(geom::BarGeometry,
                                     theme::Gadfly.Theme,
                                     aes::Gadfly.Aesthetics,
                                     orientation::Symbol)
    bar_form = empty_form

    if orientation == :horizontal
        minvar = :ymin
    else
        minvar = :xmin
    end

    stack_height = Dict()
    for x_min in getfield(aes, minvar)
        stack_height[x_min] = 0.0
    end

    # preserve factor orders of pooled data arrays
    if isa(aes.color, PooledDataArray)
        idxs = sortperm(aes.color.refs, rev=true)
    else
        idxs = 1:length(aes.color)
    end

    if orientation == :horizontal
        for i in idxs
            y_min, y_max, x, c = aes.ymin[i], aes.ymax[i], aes.x[i], aes.color[i]
            geometry_id = Gadfly.unique_svg_id()
            bar_form = combine(bar_form,
                compose(rectangle(stack_height[y_min],
                                  y_min*cy + theme.bar_spacing/2,
                                  x, (y_max - y_min)*cy - theme.bar_spacing),
                        fill(c),
                        svgid(geometry_id),
                        svgclass("geometry")))
            stack_height[y_min] += x
        end
    else
        for i in idxs
            x_min, x_max, y, c = aes.xmin[i], aes.xmax[i], aes.y[i], aes.color[i]
            geometry_id = Gadfly.unique_svg_id()
            bar_form = combine(bar_form,
                compose(rectangle(x_min*cx + theme.bar_spacing/2,
                                  stack_height[x_min],
                                  (x_max - x_min)*cx - theme.bar_spacing, y),
                        fill(c),
                        svgid(geometry_id),
                        svgclass("geometry")))
            stack_height[x_min] += y
        end
    end
    bar_form
end

# Render a bar chart grouped by discrete colors and stuck next to each other.
function render_colorful_dodged_bar(geom::BarGeometry,
                                    theme::Gadfly.Theme,
                                    aes::Gadfly.Aesthetics,
                                    orientation::Symbol)
    if orientation == :horizontal
        minvar = :ymin
    else
        minvar = :xmin
    end
    bar_form = empty_form
    dodge_width = Dict()
    dodge_count = Dict()
    for x_min in getfield(aes, minvar)
        dodge_width[x_min] = 0.0cx
        dodge_count[x_min] = get(dodge_count, x_min, 0) + 1
    end

    # preserve factor orders of pooled data arrays
    if isa(aes.color, PooledDataArray)
        idxs = sortperm(aes.color.refs, rev=true)
    else
        idxs = 1:length(aes.color)
    end

    if orientation == :horizontal
        for i in idxs
            y_min, y_max, x, c = aes.ymin[i], aes.ymax[i], aes.x[i], aes.color[i]
            geometry_id = Gadfly.unique_svg_id()
            barwidth = ((y_max - y_min)*cy - theme.bar_spacing) / dodge_count[y_min]
            bar_form = combine(bar_form,
                compose(rectangle(0.0,
                                  y_min*cy + dodge_width[y_min] + theme.bar_spacing/2,
                                  0.0, barwidth),
                        fill(c),
                        svgid(geometry_id),
                        svgclass("geometry")))
            dodge_width[y_min] += barwidth
        end
    else
        for i in idxs
            x_min, x_max, y, c = aes.xmin[i], aes.xmax[i], aes.y[i], aes.color[i]
            geometry_id = Gadfly.unique_svg_id()
            barwidth = ((x_max - x_min)*cx - theme.bar_spacing) / dodge_count[x_min]
            bar_form = combine(bar_form,
                compose(rectangle(x_min*cx + dodge_width[x_min] + theme.bar_spacing/2,
                                  0.0, barwidth, y),
                        fill(c),
                        svgid(geometry_id),
                        svgclass("geometry")))
            dodge_width[x_min] += barwidth
        end
    end
    bar_form
end

# Render bar geometry.
#
# Args:
#   geom: bar geometry
#   theme: the plot's theme
#   aes: some aesthetics
#
# Returns
#   A compose form.
#
function render(geom::BarGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    if geom.orientation == :horizontal
        Gadfly.assert_aesthetics_defined("Geom.bar", aes, :y)
        if (is(aes.xmin, nothing) || is(aes.xmax, nothing)) && is(aes.x, nothing)
            error("Geom.bar required \"x\" to be bound or both \"x_min\" and \"x_max\".")
        end
        var = :y
        minvar = :ymin
        maxvar = :ymax
    else
        Gadfly.assert_aesthetics_defined("Geom.bar", aes, :x)
        if (is(aes.ymin, nothing) || is(aes.ymax, nothing)) && is(aes.y, nothing)
            error("Geom.bar required \"y\" to be bound or both \"y_min\" and \"y_max\".")
        end
        var = :x
        minvar = :xmin
        maxvar = :xmax
    end

    if getfield(aes, minvar) === nothing
        aes2 = Gadfly.Aesthetics()
        values = getfield(aes, var)
        T = eltype(values)

        span = zero(T)
        unique_count = length(set(values))
        if unique_count > 1
            span = (maximum(values) - minimum(values)) / (unique_count - 1)
        end

        if span == zero(T)
            span = one(T)
        end

        T = promote_type(eltype(values), typeof(span/2))
        setfield!(aes2, minvar, Array(T, length(values)))
        setfield!(aes2, maxvar, Array(T, length(values)))

        for (i, x) in enumerate(values)
            getfield(aes2, minvar)[i] = x - span/2
            getfield(aes2, maxvar)[i] = x + span/2
        end
        aes = inherit(aes, aes2)
    end

    if aes.color === nothing
        form = render_colorless_bar(geom, theme, aes, geom.orientation)
    elseif geom.position == :stack
        form = render_colorful_stacked_bar(geom, theme, aes, geom.orientation)
    elseif geom.position == :dodge
        form = render_colorful_dodged_bar(geom, theme, aes, geom.orientation)
    else
        error("$(geom.position) is not a valid position for the bar geometry.")
    end

    compose(canvas(units_inherited=true),
            form,
            svgattribute("shape-rendering", "crispEdges"),
            fill(theme.default_color),
            stroke(nothing))
end


text:
# This file was generated, do not modify it. # hide
select!(power, [
    :utc_timestamp,
    :DE_solar_generation_actual,
    :DE_wind_generation_actual]);


text:
@testset "ArraySlice" begin
    a = [1,2,3,4]
    b = ArraySlice(a, 2:3)
    @test b[1] == a[2]
    b[1] = 100
    @test a[2] == 100
    @test length(b) == 2
    b = b >> 1
    @test length(b) == 3
    @test b[3] == 4
    b = b << 1
    @test length(b) == 2
    @test b[1] == 3
    b = b >> -1
    @test length(b) == 1
    b = b << -2
    @test length(b) == 3
    @test b[1] == 1
end


text:
function sparsetest(lib::PolyhedraLibrary)
    # Sparse H-representation, fixandeliminate and ⊆
    # 3x_1 + 2x_2        ≤  2
    # 2x_1        + 3x_3 ≥  1
    #  x_1               ≥  0
    #         x_2        ≥ -1
    #                x_3 ≤  1
    h = HalfSpace(sparsevec([1, 2], [3, 2], 3), 2) ∩
        HalfSpace(sparsevec([1, 3], [-2, -3]), -1) ∩
        HalfSpace(sparsevec([1], [-1], 3),      0) ∩
        HalfSpace(sparsevec([2], [-1], 3),      1) ∩
        HalfSpace(sparsevec([3], [1]),          1)
    v = convexhull([0,     1,  1//3],
                   [4//3, -1, -5//9],
                   [0,    -1,  1//3],
                   [0,     1,  1],
                   [4//3, -1,  1],
                   [0,    -1,  1])
    p = polyhedron(h, lib)
    @test !vrepiscomputed(p)
    @test p ⊆ HalfSpace([3, 2, 0], 2)
    @test !vrepiscomputed(p)
    @test !(p ⊆ HalfSpace([3, 2, 0], 1))
    @test !vrepiscomputed(p)
    inequality_fulltest(p, h)
    generator_fulltest(p, v)
    q = polyhedron(h, lib)
    @test !vrepiscomputed(q)
    q = fixandeliminate(q, [1], [1])
    # If lib supported sparse, hvectortype(q) should be sparse hence hvectortype(p) should also be sparse
    @test Polyhedra.hvectortype(typeof(p)) == Polyhedra.hvectortype(typeof(q))
    @test !vrepiscomputed(q)
    @test q ⊆ HalfSpace([2, -3], 0) # This also checks that the solver was ot dropped during fixandeliminate
    @test !vrepiscomputed(q)
    @test !(q ⊆ HalfSpace([2, -3], -1))
    @test !vrepiscomputed(q)
    hfix = HalfSpace([ 2,  0], -1) ∩
           HalfSpace([ 0, -3],  1) ∩
           HalfSpace([-1,  0],  1) ∩
           HalfSpace([ 0,  1],  1)
    vfix = convexhull([-1//2, -1//3],
                      [-1,    -1//3],
                      [-1//2,  1],
                      [-1,     1])
    inequality_fulltest(q, hfix)
    generator_fulltest(q, vfix)
end


text:
using WaveProp
using WaveProp.Geometry
using WaveProp.IO
using WaveProp.Mesh
using WriteVTK

@testset "VTK export" begin
    # This test should simply not throw an error
    Geometry.clear!()
    Ω, M = WaveProp.IO.gmsh_sphere()
    vtk_mesh_file(M, joinpath(@__DIR__,"ball")) |> vtk_save
    rm(joinpath(@__DIR__,"ball.vtu"))
    vtk_mesh_file(M, Ω, joinpath(@__DIR__,"ball")) |> vtk_save
    rm(joinpath(@__DIR__,"ball.vtu"))
    vtk_mesh_file(M, external_boundary(Ω), joinpath(@__DIR__,"sphere")) |> vtk_save
    rm(joinpath(@__DIR__,"sphere.vtu"))
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
"""
Ledgers

This package provides support for general financial ledgers

See README.md for the full documentation

Copyright 2019-2020, Eric Forgy, Scott P. Jones and other contributors

Licensed under MIT License, see LICENSE.md
"""
module Ledgers

using UUIDs, StructArrays, AbstractTrees
import Instruments
using Instruments: Position
# Although importing Instruments is enough to extend the methods below, we still import
# them explicitly so we can use them like `Ledgers.position` instead of `Instruments.position`.
# Unfortunately, `Base` also exports `position`.
import Instruments: instrument, currency, symbol, position, amount

export Identifier, AccountId, AccountType, AccountWrapper, AccountNumber
export LedgerAccount, Ledger, Entry, Account, AccountGroup, ledgeraccount, add_account
export id, balance, credit!, debit!, post!
# Although the methods below are exported by Instruments, we export them explicitly
# because we do not expect users to use `Instruments` directly.
export instrument, currency, symbol, position, amount

abstract type Identifier end

struct AccountId <: Identifier
    value::UUID
end

AccountId() = AccountId(uuid4())

struct AccountNumber
    value::String
end

AccountNumber() = AccountNumber("")

Base.convert(::Type{AccountNumber}, value::String) = AccountNumber(value)

abstract type AccountType{P <: Position} end

abstract type AccountWrapper{P <: Position} <: AccountType{P} end

account(acc::AccountType) = acc

account(acc::AccountWrapper) = acc.account

mutable struct LedgerAccount{P <: Position} <: AccountType{P}
    id::AccountId
    balance::P
end

function LedgerAccount(balance::P; ledger=nothing) where {P <: Position}
    ledger === nothing && return LedgerAccount{P}(AccountId(), balance)
    acc = LedgerAccount{P}(AccountId(), balance)
    push!(ledger.accounts, acc)
    acc
end

LedgerAccount(::Type{P}; ledger=nothing) where {P <: Position} =
    LedgerAccount(P(0), ledger=ledger)

ledgeraccount(acc::LedgerAccount) = acc

ledgeraccount(acc::AccountType) = account(acc).account

struct LedgerId <: Identifier
    value::UUID
end

LedgerId() = LedgerId(uuid4())

struct Account{P <: Position} <: AccountType{P}
    account::LedgerAccount{P}
    number::AccountNumber
    name::String
    isdebit::Bool
end

struct Ledger{P <: Position}
    id::LedgerId
    indexes::Dict{AccountId,Int}
    accounts::StructArray{LedgerAccount{P}}
end

function Ledger(accounts::Vector{LedgerAccount{P}}; id=LedgerId()) where {P <: Position}
    indexes = Dict{AccountId,Int}()
    for (index, account) in enumerate(accounts)
        indexes[account.id] = index
    end
    Ledger{P}(id, indexes, StructArray(ledgeraccount.(accounts)))
end

Ledger(::Type{P}) where {P <: Position} = Ledger(Vector{LedgerAccount{P}}())

function add_account(ledger::Ledger{P}, acc::LedgerAccount{P}) where {P <: Position}
    push!(ledger.accounts, ledgeraccount(acc))
    ledger.indexes[id(acc)] = length(ledger.accounts)
    ledger
end

struct AccountGroup{A <: AccountType}
    id::AccountId
    number::AccountNumber
    name::String
    isdebit::Bool
    indexes::Dict{AccountId,Int}
    accounts::StructArray{A}
    subgroups::StructArray{AccountGroup{A}}
end

function AccountGroup(
        ::Type{A};
        name::String="$(symbol(A)) Accounts",
        number="",
        isdebit=true,
        id=AccountId(),
        parent::Union{Nothing,AccountGroup{A}}=nothing
    ) where {A <: AccountType}
    indexes = Dict{AccountId,Int}()
    accounts = StructArray(Vector{A}())
    subgroups = StructArray(Vector{AccountGroup{A}}())
    if parent === nothing
        return AccountGroup{A}(id, number, name, isdebit, indexes, accounts, subgroups)
    else
        acc = AccountGroup{A}(id, number, name, isdebit, indexes, accounts, subgroups)
        push!(parent.subgroups, acc)
        return acc
    end
end

function Account(
        ::Type{P};
        name::String="$(symbol(P)) Account",
        number="",
        isdebit=true,
        parent::Union{Nothing,<:AccountGroup}=nothing,
        ledger::Union{Nothing,<:Ledger}=nothing
    ) where {P <: Position}
    account = LedgerAccount(P, ledger=ledger)
    parent === nothing &&
        return Account{P}(account, number, name, isdebit)
    acc = Account{P}(account, number, name, isdebit)
    add_account(parent, acc)
    acc
end

account(acc::AccountGroup) = acc

function add_account(grp::AccountGroup{A}, acc::A) where {A <: AccountType}
    push!(grp.accounts, acc)
    grp.indexes[id(acc)] = length(grp.accounts)
    grp
end

function add_account(grp::AccountGroup{A}, acc::AccountGroup{A}) where {A <: AccountType}
    push!(grp.subgroups, acc)
    grp
end

struct Entry{A1 <: AccountType, A2 <: AccountType}
    debit::A1
    credit::A2
end

id(acc::AccountType) = ledgeraccount(acc).id

id(acc::AccountGroup) = acc.id

balance(acc::AccountType) = ledgeraccount(acc).balance

function balance(group::AccountGroup{<:AccountType{P}}) where {P <: Position}
    btot = zero(P)
    for acc in group.accounts.account
        btot += balance(acc)
    end
    for grp in group.subgroups
        btot += balance(grp)
    end
    btot
end

Instruments.symbol(::Type{A}) where {P, A<:AccountType{P}} = symbol(P)

Instruments.currency(::Type{A}) where {P, A<:AccountType{P}} = currency(P)

Instruments.instrument(::Type{A}) where {P, A<:AccountType{P}} = instrument(P)

Instruments.position(::Type{A}) where {P, A<:AccountType{P}} = P

Instruments.amount(acc::AccountType) = amount(balance(acc))

debit!(acc::LedgerAccount, amt::Position) = (acc.balance += amt)

credit!(acc::LedgerAccount, amt::Position) = (acc.balance -= amt)

function post!(entry::Entry, amt::Position)
    debit!(ledgeraccount(entry.debit), amt)
    credit!(ledgeraccount(entry.credit), amt)
    entry
end

Base.getindex(ledger::Ledger, ix) = ledger.accounts[ix]

Base.getindex(ledger::Ledger, id::AccountId) =
    ledger.accounts[ledger.indexes[id]]

Base.getindex(ledger::Ledger, array::AbstractVector{<:AccountId}) =
    ledger.accounts[broadcast(id -> ledger.indexes[id], array)]

Base.getindex(grp::AccountGroup, ix) = grp.accounts.accounts[ix]

Base.getindex(grp::AccountGroup, id::AccountId) =
    grp.accounts.accounts[grp.indexes[id]]

Base.getindex(grp::AccountGroup, array::AbstractVector{<:AccountId}) =
    grp.accounts[broadcast(id -> grp.indexes[id], array)]

# struct EntityId <: Identifier
#     value::UUID
# end

# EntityId() = EntityId(uuid4())

# struct Entity
#     id::EntityId
#     name::String
#     ledgers::Dict{Type{<:Position},Ledger{<:AccountType}}
# end

# const chartofaccounts = Dict{String,AccountGroup{<:Cash}}()

# iscontra(a::AccountGroup) = !isequal(a,a.parent) && !isequal(a.parent,getledger(a)) && !isequal(a.parent.isdebit,a.isdebit)

# function loadchart(ledgername,ledgercode,csvfile)
#     data, headers = readdlm(csvfile,',',String,header=true)
#     nrow,ncol = size(data)

#     ledger = add(AccountGroup(ledgername,ledgercode))
#     for i = 1:nrow
#         row = data[i,:]
#         number = row[1]
#         name = row[2]
#         parent = chartofaccounts[row[3]]
#         isdebit = isequal(row[4],"Debit")
#         add(AccountGroup(parent,name,number,isdebit))
#     end
#     return ledger
# end

# function trim(a::AccountGroup,newparent::AccountGroup=AccountGroup(a.parent.name,a.parent.number))
#     newaccount = isequal(a,a.parent) ? newparent : AccountGroup(newparent,a.name,a.number,a.isdebit,a.balance)
#     for subaccount in a.accounts
#         balance(subaccount).amount > 0. && trim(subaccount,newaccount)
#     end
#     return newaccount
# end

function _print_account(io::IO, acc)
    iobuff = IOBuffer()
    isempty(account(acc).number.value) || print(iobuff,"[$(account(acc).number)] ")
    print(iobuff,account(acc).name,": ")
    account(acc).isdebit ? print(iobuff,balance(acc)) : print(iobuff,-balance(acc))
    print(io,String(take!(iobuff)))
end

Base.show(io::IO, id::Identifier) = print(io, id.value)

Base.show(io::IO, number::AccountNumber) = print(io, number.value)

Base.show(io::IO, acc::LedgerAccount) = print(io, "$(string(id(acc))): $(balance(acc))")

Base.show(io::IO, acc::AccountType) = _print_account(io, acc)

Base.show(io::IO, acc::AccountGroup) = print_tree(io, acc)

Base.show(io::IO, entry::Entry) = print_tree(io, entry)

Base.show(io::IO, ledger::Ledger) = print_tree(io, ledger)

Base.show(io::IO, a::Vector{<:AccountType}) = print_tree(io, a)

AbstractTrees.printnode(io::IO, acc::Ledger{P}) where {P <: Position} =
    print(io, "$(symbol(P)) Ledger: [$(acc.id)]")

AbstractTrees.children(acc::Ledger) = Vector(acc.accounts)

AbstractTrees.printnode(io::IO, acc::AccountGroup) = _print_account(io, acc)

AbstractTrees.children(acc::AccountGroup) = vcat(Vector(acc.subgroups), Vector(acc.accounts))

AbstractTrees.printnode(io::IO, b::Vector{<:AccountType}) =
    isempty(b) ? print(io, "Accounts: None") : print(io, "Accounts:")

AbstractTrees.children(b::Vector{<:AccountType}) = b

AbstractTrees.children(entry::Entry) = [entry.debit, entry.credit]

AbstractTrees.printnode(io::IO, ::Entry) = print(io, "Entry:")

end # module Ledgers


text:
using FinancialMonteCarlo, Distributions, Test;
@show "Test Parameters"
S0 = 100.0;
r = 0.02;
d = 0.01;
T = 1.0;
Nsim = 100;
Nstep = 30;

sigma = 0.2;

McConfig = MonteCarloConfiguration(Nsim, Nstep);
mc1 = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.AntitheticMC());
rfCurve = ZeroRate(r);
rfCurve2 = ZeroRate([0.00, 0.02], T);

@show "Test Subordinated Brownian Motion Parameters"
drift = 0.0;
@test_throws(ErrorException, SubordinatedBrownianMotion(-sigma, drift, InverseGaussian(1.0, 1.0)))
@test_throws(ErrorException, SubordinatedBrownianMotion(-sigma, rfCurve2.r, InverseGaussian(1.0, 1.0)))


text:
VERSION >= v"0.4" && __precompile__()

module Vega

    using JSON, ColorBrewer, Compat, KernelDensity, NoveltyColors
    import Base.print

    export VegaAxis, VegaData, VegaMark, VegaPadding, VegaScale,
           VegaTransform, VegaVisualization, VegaFormat, VegaScene,
           #VegaJSON, VegaCSV, VegaTSV,
           VegaDataRef, VegaValueRef,
           VegaMarkPropertySet, VegaMarkProperties, VegaLegend, VegaStream, VegaMarkFrom, VegaPredicate

    export tojs, print

    export barplot, choropleth, lineplot, scatterplot, areaplot, heatmap, piechart,
           histogram, popchart, waterfall, wordcloud, groupedbar, rugplot, boxplot,
           bubblechart, streamplot, stemleaf, asterplot, jointplot, dotplot, ribbonplot,
           horizon, punchcard

    export xlab!, ylab!, xlim!, ylim!, title!, legend!, text!

    export default_scales!, default_axes!, add_data!
    export colorscheme!, coord_flip!, stroke!, hline!, vline!, hover!, jitter!, layer

    #Create base color library
    colorpalettes = merge(ColorBrewer.colorSchemes, NoveltyColors.ColorDict)

    #Import helper code
    include("primitives/definetypes.jl")
    include("primitives/constructors.jl")
    include("primitives/tojs.jl")

    include("render.jl")

    # Intermediates
    include("intermediates/modifiers.jl")
    include("intermediates/defaults.jl")
    include("intermediates/add_data.jl")
    include("intermediates/xlab_ylab.jl")
    include("intermediates/xlim_ylim.jl")
    include("intermediates/title_text.jl")
    include("intermediates/hline_vline.jl")
    include("intermediates/layer.jl")

    # Higher-level API
    include("derived/choropleth.jl")
    include("derived/barplot.jl")
    include("derived/bubblechart.jl")
    include("derived/boxplot.jl")
    include("derived/groupedbar.jl")
    include("derived/lineplot.jl")
    include("derived/scatterplot.jl")
    include("derived/areaplot.jl")
    include("derived/heatmap.jl")
    include("derived/piechart.jl")
    include("derived/histogram.jl")
    include("derived/popchart.jl")
    include("derived/ribbonplot.jl")
    include("derived/rugplot.jl")
    include("derived/stemleaf.jl")
    include("derived/streamplot.jl")
    include("derived/waterfall.jl")
    include("derived/wordcloud.jl")
    include("derived/asterplot.jl")
    include("derived/jointplot.jl")
    include("derived/dotplot.jl")
    include("derived/horizon.jl")
    include("derived/punchcard.jl")

    if Pkg.installed("Escher") != nothing
        include("escher_integration.jl")
    end

    if Pkg.installed("Blink") != nothing
        include("blink_integration.jl")
        BLINK_INSTALLED = true
    else
        BLINK_INSTALLED = false
    end

    print(x::VegaVisualization) = print(JSON.json(tojs(x)))

end

max_stars_repo_name: jpthor/vega