text:
using StrongFieldApproximation
using ElectricFields
using Unitful
using UnitfulAtomic
using FFTW

using PyPlot
using Jagot.plotting
plot_style("ggplot")

using Jagot

function savedocfig(name,dir="figures")
    fig = gcf()
    filename = joinpath(@__DIR__, "src", dir, "$(name).svg")
    savefig(filename,
            transparent=false,
            facecolor=fig.get_facecolor())
    # close(fig)
    if isfile(filename)
        println("Saved $(name) to $(filename)")
    else
        @warn "Saving $(name) to $(filename) failed"
    end
end

function savedfigure(fun::Function, name, args...; kwargs...)
    cfigure(name, args...; kwargs...) do
        fun()
    end
    savedocfig(name)
end

spectrum(d::AbstractVector{<:Number}) = fftshift(fft(d))

function spectrum(d::AbstractVector{<:AbstractVector{T}}) where {T}
    D = Matrix{complex(T)}(undef, length(d), 3)
    for (i,dd) in enumerate(d)
        D[i,:] .= dd
    end
    fft!(D, 1)
    fftshift(D, 1)
end

function hhg_example()
    @field(F) do # Linear polarization
        λ = 800.0u"nm"
        I₀ = 3e14u"W/cm^2"
        τ = 6.2u"fs"
        toff = 10.0u"fs"
        tmax = 13.0u"fs"
        env = :trunc_gauss
    end;

    @field(F2) do # Slightly elliptical polarization
        λ = 800.0u"nm"
        I₀ = 3e14u"W/cm^2"
        τ = 6.2u"fs"
        toff = 10.0u"fs"
        tmax = 13.0u"fs"
        env = :trunc_gauss
        ξ = 0.2
    end;

    ndt = 300 # Steps per cycle
    Iₚ = 0.5 # Hydrogen

    # d will be a vector of scalars; by limiting the "memory" of the
    # integrals, we can include only the short trajectory.
    d = induced_dipole(Iₚ, F, ndt, memory=floor(Int, 0.65ndt));

    # d2 will be a vector of 3d vectors
    d2 = induced_dipole(Iₚ, F2, ndt, memory=floor(Int, 0.65ndt));

    t = timeaxis(F, ndt)
    tplot = 24.2e-3t
    q = fftshift(fftfreq(length(t), ndt))
    qsel = ind(q,0):ind(q,100)
    D = spectrum(d)
    D2 = spectrum(d2)
    cutoff = 3.17austrip(ponderomotive_potential(F)) + Iₚ

    savedfigure("hhg_example", figsize=(8,8)) do
        csubplot(311) do
            plot(tplot, field_amplitude(F, t))
            plot(tplot, field_amplitude(F2, t), "--")
            legend(["1d", "3d x", "3d y", "3d z"])
            axes_labels_opposite(:x)
            xlabel(L"$t$ [fs]")
            ylabel(L"$F(t)$ [au]")
        end
        csubplot(312,nox=true) do
            plot(tplot, d)
            plot(tplot, d2, "--")
            legend(["1d", "3d x", "3d y", "3d z"])
            ylabel(L"\langle\mathbf{r}\rangle(t)")
        end
        csubplot(313) do
            semilogy(q[qsel], abs2.(D[qsel]))
            semilogy(q[qsel], abs2.(D2[qsel,:]), "--")
            legend(["1d", "3d x", "3d y", "3d z"])
            axvline(cutoff/photon_energy(F), linestyle=":", color="black")
            xlabel(L"Harmonic order of 800 nm [$q$]")
            ylabel(L"|\mathbf{r}(q)|^2")
        end
    end
end

function streaking_example()
    @field(pump) do
        ħω = 80u"eV"
        I₀ = 1e8u"W/cm^2"
        τ = 200u"as"
        toff = 400u"as"
        tmax = 500u"as"
        env = :trunc_gauss
    end;

    @field(probe) do
        λ = 800.0u"nm"
        I₀ = 1e12u"W/cm^2"
        τ = 2.66u"fs"
        toff = 4.0u"fs"
        tmax = 5.0u"fs"
        env = :trunc_gauss
    end;

    F = pump + delay(probe, 1.0u"fs");

    ndt = 100

    nk = 100
    nθ = 2 # Forward and backward spectra only
    k,kmag,θ = momentum_grid(50u"eV", 80u"eV", nk, nθ,
                             spacing=:energy);

    Iₚ = 14u"eV" # "Krypton"
    c = photoelectron_spectrum(k, Iₚ, F, ndt);

    x = 27.211kmag.^2/2

    savedfigure("streaking_single_delay") do
        plot(x, abs2.(c))
        legend(["Forward", "Backward"])
        xlabel(L"$W_k$ [eV]")
    end

    delays = range(-4u"fs", stop=4u"fs", length=100)

    forward = zeros(ComplexF64, nk, length(delays))
    backward = zeros(ComplexF64, nk, length(delays))

    for (j,d) in enumerate(delays)
        println(j,"/",length(delays))
        cc = photoelectron_spectrum(k, Iₚ, pump + delay(probe, d), ndt, verbosity=0)
        forward[:,j] = cc[:,1]
        backward[:,j] = cc[:,2]
    end

    savedfigure("streaking_spectrogram", figsize=(8,10)) do
        csubplot(311) do
            plot(ustrip(delays), vector_potential(probe, austrip.(delays)))
            axes_labels_opposite(:x)
            xlabel(L"$t$ [fs]")
            ylabel(L"$A_{\mathrm{IR}}(t)$ [au]")
            margins(0,0.1)
        end
        csubplot(312, nox=true) do
            plot_map(ustrip(delays), x, abs2.(forward))
            ylabel(L"$W_k$ forward [eV]")
        end
        csubplot(313) do
            plot_map(ustrip(delays), -x, abs2.(backward))
            xlabel("Pump–probe delay [fs]")
            ylabel(L"$W_k$ backward [eV]")
        end
    end

end

function ppt_example()
    I = 10 .^ range(10, stop=15, length=100) * u"W/cm^2";
    Iau = ElectricFields.Iaustrip.(I);
    λ = [2000u"nm", 800u"nm", 500u"nm"]
    ω = austrip.(2π*1u"c" ./ λ)

    Iₚ = 0.5 # Hydrogen

    rates = StrongFieldApproximation.IonizationRates.PPT.(Iₚ, Iau, ω', 0, 0);

    savedfigure("ppt_example") do
        loglog(ustrip.(I), rates)
        xlabel(L"$I$ [W/cm$^2$]")
        ylabel(L"Ionization rates [s$^{-1}$]")
        legend(string.(λ))
    end
end

macro echo(expr)
    println(expr)
    :(@time $expr)
end

@info "Documentation plots"
mkpath("docs/src/figures")
@echo hhg_example()
@echo streaking_example()
@echo ppt_example()


text:
export AnalogInput10V, AnalogOutput10V

mutable struct AnalogInput10V <: AbstractDevice
    i::Int32
    stream::LabStream
    AnalogInput10V(i::Int32) = new(i)
end
mutable struct AnalogOutput10V <: AbstractDevice
    i::Int32
    stream::LabStream
    AnalogOutput10V(i::Int32) = new(i)
end
AnalogInput10V(i::Int64) = AnalogInput10V(convert(Int32, i))
AnalogOutput10V(i::Int64) = AnalogOutput10V(convert(Int32, i))

initialize(::AnalogInput10V) = nothing
initialize(::AnalogOutput10V) = nothing

close(::AnalogInput10V) = nothing
close(input::AnalogOutput10V) = ccall((:comedi_write_zero, comedipath), Void, (Int32, Int32, Int32), Int32(0), Int32(1), input.i)

getwritecommand(stream::LabStream, input::AnalogInput10V, val) = error("Can't write to device $input")
getreadcommand(stream::LabStream, output::AnalogOutput10V, val) = error("Can't read from device $output")

function getwritecommand(stream::BeagleBoneStream, input::AnalogOutput10V, val::Float64)
    abs(val) <= 10 || error("Volatage $val not in range [-10,10]")
    return ("analogin10", input.i, val)
end
function getreadcommand(stream::BeagleBoneStream, input::AnalogInput10V)
    return ("analogin10", input.i, val)
end

function getwritecommand(stream::ComediStream, input::AnalogOutput10V, val::Float64)
    abs(val) <= 10 || error("Volatage $val not in range [-10,10]")
    return (Int32(0),Int32(1),input.i,val)
end
function getreadcommand(stream::ComediStream, input::AnalogInput10V)
    return (Int32(0),Int32(0),input.i)
end

ext: jl


text:
"""
Charged Particle in a singular magnetic field of the form
``B(x,y,z) = (x^2 + y^2)^{-3/2} e_z``.
"""
module SingularField

    import ElectromagneticFields.Singular

    export charged_particle_3d_iode, hamiltonian, angular_momentum

    Singular.@code() # inject magnetic field code

    const q₀ = [1., 0., 0., 0., -1., 0.]

    include("charged_particle_3d.jl")

end


text:
using Test, ModelObjectsLH,  ModelParams

## Parameter vector test
import ModelParams.get_pvector

mdl = ModelParams;

## Test basic operations
function pvectorTest()
    @testset "Pvector" begin
        pv = ParamVector(ObjectId(:pv1));
        # println(pv);
        @test length(pv) == 0

        p = Param(:p1, "param1", "sym1", [1.2 3.4; 4.5 5.6])
        ModelParams.append!(pv, p)
        @test length(pv) == 1
        nParam, nElem = mdl.n_calibrated_params(pv; isCalibrated = is_calibrated(p));
        @test nParam == 1
        @test nElem == 4

        pOut = retrieve(pv, :p11)
        @test isnothing(pOut);
        @test !param_exists(pv, :p11)
        pOut = retrieve(pv, :p1)
        @test pOut.name == :p1
        @test pv[1].name == :p1;
        @test param_exists(pv, :p1)

        p = Param(:p2, "param2", "sym2", [2.2 4.4; 4.5 5.6])
        ModelParams.append!(pv, p)
        # println(pv)
        ModelParams.remove!(pv, :p1)
        @test length(pv) == 1
        @test param_exists(pv, :p2)
        @test !param_exists(pv, :p1)

        # Test replace
        p2 = Param(:p2, "param2", "sym2", 1.0);
        ModelParams.replace!(pv, p2);
        p22 = retrieve(pv, :p2);
        @test ModelParams.value(p22) == 1

        # Retrieve non-existing
        p3 = retrieve(pv, :notThere);
        @test isnothing(p3);

        # Parameter value
        pValue = param_value(pv, :p2);
        @test pValue == p2.value
        pValue = param_value(pv, :notThere);
        @test isnothing(pValue)
    end
end

function set_status_test()
    @testset "Set status" begin
        n = 7;
        pvec = mdl.make_test_pvector(n);
        set_calibration_status_all_params!(pvec, false);
        @test n_calibrated_params(pvec; isCalibrated = true) == (0, 0);
        set_calibration_status_all_params!(pvec, true);
        @test n_calibrated_params(pvec; isCalibrated = true)[1] == n;

        set_default_values_all_params!(pvec);
        for p in pvec
            @test value(p) ≈ mdl.default_value(p);
        end
    end
end

function iter_test()
    @testset "Iteration" begin
        n = 5;
        pv = mdl.make_test_pvector(n);
        j = 0;
        for p in pv
            j += 1;
            @test isa(p, AbstractParam);
        end
        @test j == n
	end
end

## Vector to Dict and back
function pvectorDictTest()
    @testset "Pvector Dict" begin
        isCalibrated = true;
        pv = mdl.make_test_pvector(4);
        d = make_dict(pv; isCalibrated = true, valueType = :value);
        for pName in keys(d)
            p1 = retrieve(pv, pName);
            @test d[pName] == mdl.value(p1);
        end
        # p1 = retrieve(pv, :p1);
        # p3 = retrieve(pv, :p3);
        # @test d[:p1] == p1.value
        # @test d[:p3] == p3.value
        @test length(d) == n_calibrated_params(pv; isCalibrated = isCalibrated)[1];

        # Make vector and its inverse (make Dict from vector)
        # Vector contains transformed parameters
        vv = mdl.make_value_vector(pv, 1);
        @test isa(get_values(pv, vv),  Vector{Float64});
        pValueV = Vector{Float64}();
        for pName in keys(d)
            p1 = retrieve(pv, pName);
            pValue = transform_param(pv.pTransform, p1);
            append!(pValueV, pValue);
        end
        # p1Value = transform_param(pv.pTransform, p1);
        # p3Value = transform_param(pv.pTransform, p3);
        @test get_values(pv, vv) == pValueV
        # @test all(ModelParams.lb(vv) .≈ pv.pTransform.lb)

        # pDict, _ = vector_to_dict(pv, vv, isCalibrated);
        # @test length(pDict) == 2
        # @test pDict[:p1] == p1.value
        # @test pDict[:p3] == p3.value
    end
end

function set_values_test()
    @testset "Set values" begin
        isCalibrated = true;
        pv = mdl.make_test_pvector(7);
        d = make_dict(pv; isCalibrated, valueType = :value);
        # Change values in `pv`
        pv2 = mdl.make_test_pvector(7; offset = 1.0);
        # d2 = make_dict(pv2; isCalibrated, useValues = true);
        ModelParams.set_own_values_from_pvec!(pv, pv2, true);
        # Check that values are now higher using a `Dict`
        d3 = make_dict(pv; isCalibrated, valueType = :value);
        @test isequal(keys(d), keys(d3));
        for key in keys(d)
            @test all(d3[key] .> d[key])
        end

        p = pv[1];
        @test p isa AbstractParam;
        if default_value(p) isa Real
            lbnd = -10.0;
            ubnd = 20.0;
        else
            sz = size(default_value(p));
            lbnd = fill(-10.0, sz);
            ubnd = fill(20.0, sz);
        end
        set_bounds!(pv, p.name; lb = lbnd, ub = ubnd);
        @test mdl.lb(p) == lbnd;
        @test mdl.ub(p) == ubnd;
	end
end

function compare_test()
    @testset "Compare" begin
        pvec1 = mdl.make_test_pvector(7);
        pvec2 = mdl.make_test_pvector(7);
        pMiss1, pMiss2, pDiff = compare_params(pvec1, pvec2);
        @test isempty(pMiss1);
        @test isempty(pMiss2);
        @test isempty(pDiff);
    end
end

## Finding ParamVector in object
mutable struct M1 <: ModelObject
    objId :: ObjectId
    pvec :: ParamVector
end

mutable struct M2 <: ModelObject
    objId :: ObjectId
    x :: ParamVector
    y :: Float64
end

get_pvector(m :: M2) = m.x;

mutable struct M3 <: ModelObject
    objId :: ObjectId
    y :: Float64
end

function get_pvector_test()
    @testset "get pvector" begin
        m1 = M1(ObjectId(:m1),  ParamVector(ObjectId(:m1)));
        @test ModelParams.get_pvector(m1) == m1.pvec

        m2 = M2(ObjectId(:m2),  ParamVector(ObjectId(:m2)), 1.2);
        @test ModelParams.get_pvector(m2) == m2.x

        m3 = M3(ObjectId(:m3),  1.2);
        pvec3 = ModelParams.get_pvector(m3);
        @test length(pvec3) == 0
    end
end

function report_test()
    @testset "Pvector reporting" begin
        pv = mdl.make_test_pvector(9);

        pList = ModelParams.calibrated_params(pv; isCalibrated = true);
        for p in pList
            @test is_calibrated(p)
        end

        pvEmpty = ParamVector(ObjectId(:pvEmpty));
        @test isempty(pvEmpty)
        @test isempty(ModelParams.calibrated_params(pvEmpty; isCalibrated = false))

        println("-----  Calibrated parameters")
        ModelParams.report_params(pv, true)
        nParam, nElem = ModelParams.n_calibrated_params(pv; isCalibrated = true);
        @test nParam == 5
        @test nElem > 10

        println("-----  Fixed parameters")
        ModelParams.report_params(pv, false)
        nParam, nElem = ModelParams.n_calibrated_params(pv; isCalibrated = false);
        @test nParam == 4
        @test nElem > 8

        println("-----  Close to bounds");
        report_params(pv, true; closeToBounds = true);
    end
end

@testset "ParamVector" begin
    pvectorTest();
    set_status_test();
    iter_test()
    get_pvector_test()
    pvectorDictTest()
    set_values_test()
    report_test()
end

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


text:
using Base

"""
"""
function calculate!(cascade::Cascade{Data}, args...; kwargs...)
    data = collect_data(cascade)
    set_value!(data, calculate(data, args...; kwargs...))
    return cascade
end

"""
"""
calculate(v::AbstractMatrix, fun::Function, args...) = matrix(fun.(vectorize(v), args...))
calculate(data::Vector{Data}, args...; kwargs...) = calculate(get_value(data), args...; kwargs...)
calculate(cascade, args...; kwargs...) = calculate!(copy(cascade), args...; kwargs...)

"Check that there are two dimensions, MAXIMUM, and that one of these dimensions is ONE."
isvector(mat::AbstractMatrix) = length(size(mat)) == 2 && (1 in size(mat))
isvector(vec::AbstractVector) = true

"Drops zero values."
dropzero(mat::Matrix) = mat[all.(eachrow(mat.>=1E-10)),:]
dropzero(vec::Vector) = vec[abs.(vec).>=1E-10]

"""
This function calculates the graphical width of each bar based on the number of cascade
steps, ``N_{step}``:
```math
w_{step} = \\dfrac{WIDTH - \\left(N_{step}+1\\right) SEP}{N_{step}}
```
"""
width(steps::Integer; space=SEP, margin=SEP/2) = (WIDTH - 2*margin - space*steps)/steps

"""
"""
function cumulative_v(v::AbstractArray; shift=0.0, kwargs...)
    N = size(v,1)
    L = matrix(LinearAlgebra.UnitLowerTriangular, N; value=1.0) + shift*I
    return L * v
end

cumulative_v(x::Cascade{Data}; kwargs...) = cumulative_v(collect_value(x); kwargs...)

"""
    correlate(v::AbstractVector, A::AbstractMatrix)
This function applies a correlation matrix `A` to the vector `v`.

```math
\\vec{x}' = \\prod_{i=N}^1 \\left(S_{i,i} A+I\\right) \\vec{x}
```

where ``S_{i,i}``, defined using [`pick`](@ref), is a sparse matrix `[i,i]=1`,
``A`` is a random rotation matrix produced by [`random_rotation`](@ref).
``S_{i,i} A`` selects the ``i^\\text{th}`` row of ``A``, which defines how steps `1:i-1`
impact step `i`. Once correlations have been applied to the `i`th step,
this new value will be used when applying correlations between steps `i` and `i+1`.
Taking this approach applies correlations to ``A`` sequentially to propagate nonlinearities
to subsequent steps.
"""
function correlate(v::AbstractArray, A::AbstractMatrix, order::AbstractVector{Int})
    N = size(v,1)

    # Save a list of the A1 = A[ii,:] matrix, where each row ii is individually selected and
    # consistent with the natural order. Then, reorder this list to match the permutation order.
    lst_A = pick_from(A; dims=2)
    lst_A = lst_A[order]

    # Apply the interaction defined in each row individually and sequentially.
    # If A is a 3x3 matrix, ordered [2,3,1], this would look like.
    #   v[1] = A[2,:] * v
    #   v[2] = A[3,:] * A[2,:] * v
    #   v[3] = A[1,:] * A[3,,:] * A[2,:] * v
    # Note how the interactions defined in A are applied in the permuted order,
    # but the order in which the VALUES are stored remains unchanged.
    prod_A = [*(lst_A[ii:-1:1]...) for ii in 1:N]
    lst_v = [*(lst_A[ii:-1:1]...) * v for ii in 1:N]

    # Now, populate the final, correlated value with the iith ROW (the current investment)
    # of the iith list element (the impact of ALL investments, IN THE PERMUTED ORDER)
    # on ALL values thus far.
    vout = copy.(v)
    [vout[ii,:] .= lst[ii,:] for (ii,lst) in zip(order, lst_v)]
    return update_stop!(convert(Matrix, vout))
end

"""
    correlate!(x::Cascade)
This function applies the correlation matrix defined in `x.correlation` to
the values in `x.steps` using [`correlate`](@ref)

# Returns
- `x::Cascade` with the following updates:
    1. `x.steps`, with correlations applied;
    1. `x.stop`, to the new cumulative sum after correlations have been applied; and
    1. `x.correlation_applied=true`, to avoid the risk of applying the correlation matrix
        multiple times.
"""
function correlate!(cascade::Cascade{Data}; kwargs...)
    if !cascade.iscorrelated
        data = collect_data(cascade)

        v = get_value(data)
        A = collect_correlation(cascade)
        order = collect_permutation(cascade.permutation)

        set_value!(data, correlate(v, A, order))
        cascade.iscorrelated = true
    end

    return cascade
end

"""
"""
function Base.permute!(cascade::Cascade; kwargs...)
    if !cascade.ispermuted
        cascade.steps = cascade.steps[get_permutation(cascade.permutation)]
        cascade.ispermuted = true
    end
    return cascade
end

"""
```jldoctest
julia> pick(3, 4)
4×4 SparseArrays.SparseMatrixCSC{Int64,Int64} with 1 stored entry:
  [3, 3]  =  1

julia> pick(1:2, 4)
4×4 SparseArrays.SparseMatrixCSC{Int64,Int64} with 2 stored entries:
  [1, 1]  =  1
  [2, 2]  =  1
```
"""
pick(idx, dim) = SparseArrays.sparse(fill([idx;],2)..., 1, fill(dim,2)...)
pick(dim) = [pick(ii,dim) for ii in 1:dim]

"""
    pick_from(A; kwargs...)
This method returns a list ``\\vec{A}`` of all correlation factors ``A_i`` calculated using
`pick_from(ii,A; kwargs...)`:

```math
\\vec{A} = \\begin{pmatrix}
A_1 & A_2 & \\dots A_{N}
\\end{pmatrix}
```

    pick_from(ii, A; kwargs)
This function returns the correlation matrix ``A_i`` with all off-diagonal elements,
with the exception of those in ROW or COLUMN `ii`, set to zero.

# Arguments
- `A <: AbstractMatrix`, a random correlation matrix ``A \\in \\mathbb{R}^{N\\times N}``,
    produced by [`random_rotation`]
- `ii::Int`, non-zero ROW OR COLUMN

# Keyword Arguments
- `dims::Int = 1, 2`: Matrix dimension from which to select.

# Returns
- `A_i <: AbstractMatrix`, defined as

```math
\\begin{aligned}
A_{i,\\circ} &= S_i A' + I
\\\\
A_{\\circ,j} &= A' S_j + I
\\end{aligned}
```

where
- ``S_i \\in \\mathbb{Z}^{N\\times N}`` is a sparse matrix produced by [`pick`](@ref)
    with `[i,i]=1`
- ``A' \\in \\mathbb{R}^{N\\times N}`` is defined:

```math
A' = \\begin{cases}
A_{i,j} & i\\neq j \\
0       & i=j
\\end{cases}
```

All diagonal elements of ``A`` are set to zero prior to performing any calculations since
selecting a row using ``S_i`` zeros all other elements on the diagonal. Explicitly adding
``I`` ensures that applying correlations to the `i`th step will not impact any other steps.
"""
function pick_from(ii::Int, A; dims)
    N = size(A,1)
    # Ensure that, if A has ones on the diagonal, these are made zero.
    # When one row is "picked", this will zero the other elements on the diagonal,
    # so we will add these back at each step in the iteration.
    # LinearAlgebra.tr(A)==N && (A -= I)
    A -= I

    return (dims==1 ? pick(ii,N)*A : A*pick(ii,N)) + I
end

pick_from(A; kwargs...) = [pick_from(ii,A; kwargs...) for ii in 1:size(A,1)]


text:
module Column

using DiffEqOperators: CenteredDifference
using LinearAlgebra
using RecursiveArrayTools: VectorOfArray
using NLsolve
using Dierckx
using NumericalIntegration

using ..Mesh
using ..InternalForces
using ..CrossSection
using ..MaterialModels

export Properties, SolutionControls, SectionGeometry, initialize, solve, normal_stresses, normal_strains, warping_displacements,
         deformed_shape, calculate_derivative_operators, apply_end_boundary_conditions

mutable struct Properties

   A::Array{Float64,1}
   Ix::Array{Float64,1}
   Iy::Array{Float64,1}
   J::Array{Float64,1}
   Cw::Array{Float64,1}
   xc::Array{Float64,1}
   yc::Array{Float64,1}
   xs::Array{Float64,1}
   ys::Array{Float64,1}
   xo::Array{Float64,1}
   yo::Array{Float64,1}
   Io::Array{Float64,1}
   E::Array{Float64,1}
   ν::Array{Float64,1}
   G::Array{Float64,1}
   kx::Float64
   ky::Float64
   kϕ::Float64
   hx::Float64
   hy::Float64
   uo::Array{Float64,1}
   vo::Array{Float64,1}
   ϕo::Array{Float64,1}
   P::Array{Float64,1}
   dz::Array{Float64,1}
   z::Array{Float64,1}
   dm::Array{Float64,1}
   Azz::Array{Float64,2}
   Azzzz::Array{Float64,2}
   supports::Array{Float64,1}
   xcoords::Array{Float64,1}
   ycoords::Array{Float64,1}
   Aio::Array{Float64,1}   #before plasticity begins
   Ai::Array{Float64,1}
   cxi::Array{Float64,1}
   cyi::Array{Float64,1}
   stress_strain::Dierckx.Spline1D
   E_tan_strain::Dierckx.Spline1D

end

mutable struct SolutionControls

   yield_criterion::String
   strain_limit::Float64

end

mutable struct SectionGeometry

   xcoords::Array{Float64,1}
   ycoords::Array{Float64,1}
   mesh_size::Float64

end

function calculate_derivative_operators(dz)

   NumberOfNodes=length(dz)+1

   # add extra dz on each end for padding nodes used in CenteredDifference
   dz = [dz[1]; dz; dz[end]]

   NthDerivative = 4
   DerivativeOrder = 2
   Azzzz = CenteredDifference(NthDerivative, DerivativeOrder, dz, NumberOfNodes)

   Azzzz = Array(Azzzz)   #concretization of derivative-free operator
   Azzzz = Azzzz[:,2:end-1]  #trim off ghost nodes, not needed

   NthDerivative = 2
   DerivativeOrder = 2
   Azz = CenteredDifference(NthDerivative, DerivativeOrder, dz, NumberOfNodes)
   Azz = Array(Azz)   #concretization of derivative-free operator
   Azz = Azz[:,2:end-1]  #trim off ghost nodes, not needed

   return Azzzz, Azz

end

function calculateBoundaryStencils(BCFlag, h, NthDerivative)

   #Calculate boundary conditions stencils without ghost nodes using
   #Jorge M. Souza, "Boundary Conditions in the Finite Difference Method"
   #https://cimec.org.ar/ojs/index.php/mc/article/download/2662/2607

   TaylorCoeffs =  [1; h; h^2/2; h^3/6; h^4/24]

   RHS = zeros(5)
   #row1*[A;B;C;D;E]*u2
   #row2*[A;B;C;D;E]*u2'
   #row3*[A;B;C;D;E]*u2''
   #row4*[A;B;C;D;E]*u2'''
   #row5*[A;B;C;D;E]*u2''''

   if BCFlag == 1 #simply supported end

      LHS = [1 1 1 1 0
         -1 0 1 2 0
         1 0 1 4 2
         -1 0 1 8 -6
         1 0 1 16 12]

      RHS[NthDerivative+1] = (1/TaylorCoeffs[NthDerivative+1])
      BoundaryStencil = LHS\RHS
      BoundaryStencil = ((BoundaryStencil[1:4]),(zeros(4)))  #since u''=0

   elseif BCFlag == 2 #fixed end

      LHS = [1 1 1 1 0
         -1 0 1 2 1
         1 0 1 4 -2
         -1 0 1 8 3
         1 0 1 16 -4]

      RHS[NthDerivative+1] = (1/TaylorCoeffs[NthDerivative+1])
      BoundaryStencil = LHS\RHS
      BoundaryStencil = ((BoundaryStencil[1:4]),(zeros(4)))  #since u'=0

   elseif BCFlag == 3 #free end
                #u'' u'''
      LHS = [1 1 1  0   0
           0 1 2  0   0
           0 1 4  0.5 0
           0 1 8  0   6
           0 1 16 0  0]

      RHS[NthDerivative+1] = (1/TaylorCoeffs[NthDerivative+1])
      BoundaryStencil1 = LHS\RHS  #at free end
      BoundaryStencil1 = BoundaryStencil1[1:3]   #u''=u'''=0

      # use simply supported BC to find stencil at one node in from free end
      LHS = [1 1 1 1 0
         -1 0 1 2 0
         1 0 1 4 2
         -1 0 1 8 -6
         1 0 1 16 12]

      RHS[NthDerivative+1] = (1/TaylorCoeffs[NthDerivative+1])
      BoundaryStencil2 = LHS\RHS #at one node over from free end
      BoundaryStencil2 = BoundaryStencil2[1:4]
      BoundaryStencil = ((BoundaryStencil1), (BoundaryStencil2))  #two stencils are calculated

   end

   return BoundaryStencil

end

function apply_end_boundary_conditions(A, EndBoundaryConditions, NthDerivative, dz)

   #left end
   h = dz[1]
   BCFlag = EndBoundaryConditions[1]
   BoundaryStencil = calculateBoundaryStencils(BCFlag,h,NthDerivative)

   A[1,:] .= 0.0
   A[2,:] .= 0.0

   if (BCFlag == 1) | (BCFlag == 2)   #make this cleaner, combine
      A[2,1:length(BoundaryStencil[1])] = BoundaryStencil[1]
   else
      A[1,1:length(BoundaryStencil[1])] = BoundaryStencil[1]
      A[2,1:length(BoundaryStencil[2])] = BoundaryStencil[2]
   end

   #right end
   h = dz[end]
   BCFlag = EndBoundaryConditions[2]
   BoundaryStencil = calculateBoundaryStencils(BCFlag,h,NthDerivative)

   A[end,:] .= 0.0
   A[end-1,:] .= 0.0

   if (BCFlag == 1) | (BCFlag == 2)
      A[end-1,(end-(length(BoundaryStencil[1])-1)):end] = reverse(BoundaryStencil[1])
   else
      A[end,end-(length(BoundaryStencil[1])-1):end] = reverse(BoundaryStencil[1])
      A[end-1,end-(length(BoundaryStencil[2])-1):end] = reverse(BoundaryStencil[2])
   end

   return A

end

function initialize(member_definitions, section_properties, section_geometry, material_properties, loads, springs, end_boundary_conditions, supports, imperfections, solution_controls)

   dz, z, dm = Mesh.define_line_element(member_definitions)

   NumberOfNodes=length(dz)+1

   #define property vectors
   A =  Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 1)
   Ix = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 2)
   Iy = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 3)
   J = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 4)
   Cw = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 5)
   xc = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 6)
   yc = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 7)
   xs = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 8)
   ys = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 9)

   xo = -(xc .- xs)
   yo = yc .- ys

   Io = Ix .+ Iy .+ A .* (xo.^2 + yo.^2)

   E = Mesh.create_line_element_property_array(member_definitions, dm, dz, material_properties, 4, 1)
   ν = Mesh.create_line_element_property_array(member_definitions, dm, dz, material_properties, 4, 2)
   G = E./(2 .*(1 .+ ν))

   kx = springs[1]
   ky = springs[2]
   kϕ = springs[3]

   hx = springs[4]
   hy = springs[5]

   #define imperfections
   uo = imperfections[1]
   vo = imperfections[2]
   ϕo = imperfections[3]

   #define load
   P =  loads

   #calculate derivative operators
   Azzzz,Azz = calculate_derivative_operators(dz)

   #apply left and right end boundary condition stencils to derivative operators
   NthDerivative = 4
   Azzzz = apply_end_boundary_conditions(Azzzz, end_boundary_conditions, NthDerivative, dz)

   NthDerivative = 2
   Azz = apply_end_boundary_conditions(Azz, end_boundary_conditions, NthDerivative, dz)

   #initialize cross-section geometry and triangulation, add them later if needed
   xcoords = section_geometry.xcoords
   ycoords = section_geometry.ycoords
   Ai = []
   Aio = []
   cxi = []
   cyi = []

   #initialize inelasticity
   if isempty(solution_controls.yield_criterion)

      stress_strain_curve = []

   else   #need to generalize this to any material

      σy = material_properties[1][3]
      σy1 = material_properties[1][4]
      σu = material_properties[1][5]
      σf = material_properties[1][6]
      ϵy = material_properties[1][7]
      ϵy1 = material_properties[1][8]
      ϵu = material_properties[1][9]
      ϵf = material_properties[1][10]
      n = material_properties[1][11]

      Ai, cxi, cyi = CrossSection.mesh(section_geometry.xcoords, section_geometry.ycoords, section_geometry.mesh_size)
      Aio = deepcopy(Ai)   #preserve inital cross-section cells
      steel_law  = MaterialModels.steel(σy, σy1, σu, σf, ϵy, ϵy1, ϵu, ϵf, n)

   end

   properties = Properties(A, Ix, Iy, J, Cw, xc, yc, xs, ys, xo, yo, Io, E, ν, G, kx, ky, kϕ, hx, hy, uo, vo, ϕo, P, dz, z, dm, Azz, Azzzz, supports,
   xcoords, ycoords, Aio, Ai, cxi, cyi, steel_law.stress_strain, steel_law.E_tan_strain)

   return properties

end

function equations!(properties)

   NumberOfNodes=length(properties.dz)+1

   #build identity matrix for ODE operations
   AI = Matrix(1.0I,NumberOfNodes,NumberOfNodes)

   #build operator matrix that doesn't update with load

   #start building operator matrix, LHS first
   A11 = zeros(NumberOfNodes,NumberOfNodes)
   A12 = zeros(NumberOfNodes,NumberOfNodes)
   A13 = zeros(NumberOfNodes,NumberOfNodes)
   A21 = zeros(NumberOfNodes,NumberOfNodes)
   A22 = zeros(NumberOfNodes,NumberOfNodes)
   A23 = zeros(NumberOfNodes,NumberOfNodes)
   A31 = zeros(NumberOfNodes,NumberOfNodes)
   A32 = zeros(NumberOfNodes,NumberOfNodes)
   A33 = zeros(NumberOfNodes,NumberOfNodes)

   #calculate operator quantities on LHS  AU=B
   for i = 1:NumberOfNodes
      A11[i,:] = properties.E .* properties.Iy .* properties.Azzzz[i,:] .+ properties.P .* properties.Azz[i,:] .+ properties.kx .* AI[i,:]
      A13[i,:] = properties.kx .*(properties.yo .- properties.hy) .* AI[i,:] .+ properties.P .* properties.yo .* properties.Azz[i,:]
      A22[i,:] = properties.E .* properties.Ix .* properties.Azzzz[i,:] .+ properties.P .* properties.Azz[i,:] .+ properties.ky .* AI[i,:]
      A23[i,:] = -properties.ky .* (properties.xo .-properties.hx) .* AI[i,:] - properties.P .* properties.xo .* properties.Azz[i,:]
      A31[i,:] = properties.kx.*(properties.yo .- properties.hy) .* AI[i,:] .+ properties.P .* properties.yo .* properties.Azz[i,:]
      A32[i,:] = -properties.ky.*(properties.xo .- properties.hx) .* AI[i,:] .- properties.P .* properties.xo .* properties.Azz[i,:]
      A33[i,:] = properties.E .* properties.Cw .* properties.Azzzz[i,:] .-(properties.G .* properties.J .- (properties.P .* properties.Io ./ properties.A)) .* properties.Azz[i,:]
      .+ properties.kx .* (properties.yo .- properties.hy).^2 .*AI[i,:] .+ properties.ky.*(properties.xo .- properties.hx).^2 .*AI[i,:] .+ properties.kϕ .* AI[i,:]
   end

   #calculate RHS of AU=B
   B1 = -properties.P .* properties.Azz * properties.uo - properties.P .* properties.yo .* properties.Azz * properties.ϕo
   B2 = -properties.P .* properties.Azz * properties.vo + properties.P .* properties.xo .* properties.Azz * properties.ϕo
   B3 = -properties.P .* properties.yo .* properties.Azz * properties.uo + properties.P .* properties.xo .* properties.Azz * properties.vo - properties.P .* (properties.Io ./ properties.A) .* properties.Azz * properties.ϕo

   #reduce problem to free dof
   FixedDOF = [findall(x->abs(x-properties.supports[i])<=10e-6, properties.z) for i=1:length(properties.supports)]
   FixedDOF = VectorOfArray(FixedDOF)
   FixedDOF  = convert(Array,FixedDOF)
   FreeDOF = setdiff(1:NumberOfNodes,FixedDOF)

   Am = [A11[FreeDOF,FreeDOF] A12[FreeDOF,FreeDOF] A13[FreeDOF,FreeDOF];
      A21[FreeDOF,FreeDOF] A22[FreeDOF,FreeDOF] A23[FreeDOF,FreeDOF];
      A31[FreeDOF,FreeDOF] A32[FreeDOF,FreeDOF] A33[FreeDOF,FreeDOF]]

   Bm = [B1[FreeDOF]; B2[FreeDOF]; B3[FreeDOF]]

   return Am, Bm, FreeDOF

end

function inelasticity_update_first_yield(U, properties, solution_controls, free_dof)

   num_sections = length(properties.z)

   u = zeros(num_sections)

   u[free_dof] = U[1:length(free_dof)]   #this is specific to u displacements right now, need to generalize at some point

   Myy = InternalForces.moment(properties.z, properties.dm, -u, properties.E, properties.Iy)

   for i=1:num_sections

      ϵ_axial, ϵ_flexure, ϵ_total = calculate_strain(properties.P[i], Myy[i], properties.A[i], properties.Iy[i], properties.xc[i], properties.cxi, properties.E[i])

      if maximum(abs.(ϵ_total)) > solution_controls.strain_limit

         properties.Iy[i] = 0.50 * 362.0

      end

   end

   return properties

end

#check this!!!!

function inelasticity_update!(U, properties, solution_controls, free_dof)

   num_sections = length(properties.z)

   u = zeros(num_sections)

   u[free_dof] = U[1:length(free_dof)]   #this is specific to u displacements right now, need to generalize at some point

   Myy = InternalForces.moment(properties.z, properties.dm, -u, properties.E, properties.Iy)

   Myy[1] = 0.0   #checking
   Myy[end] = 0.0

   for i = 1:num_sections
      #calculate strain
      ϵ_axial, ϵ_flexure, ϵ_total = calculate_strain(properties.P[i], Myy[i], properties.A[i], properties.Iy[i], properties.xc[i], properties.cxi, properties.E[i])

      if maximum(abs.(ϵ_total)) > solution_controls.strain_limit
         #update cross-section properties
         # properties.Ai = update_cell_area(properties.E_tan_strain, abs.(ϵ_total), properties.Aio, properties.E[i])

         Ai_temp = update_cell_area(properties.E_tan_strain, abs.(ϵ_total), properties.Aio, properties.E[i])

         if Ai_temp < properties.Ai

            properties.Ai = Ai_temp

         end

         properties.xc[i], properties.yc[i] = CrossSection.centroid_from_cells(properties.Ai, properties.cxi, properties.cyi)
         # properties.A[i] = CrossSection.area_from_cells(properties.Ai)

         properties.Iy[i] = CrossSection.moment_of_inertia_from_cells(properties.Ai, properties.cxi, properties.xc[i])

         #catch numerical issues I noticed, need to dig into this more
         if properties.A[i] <= 0.0
            # properties.A[i] = 0.0
            properties.Iy[i] = 0.0
         end

      end

   end

   # plot(properties.z, properties.Iy, seriestype = :scatter)

   return properties

end

function update_cell_area(E_tan_strain, ϵ, Aio, E)

   num_cells = length(Aio)
   Ai = zeros(Float64, num_cells)
   E_tan = zeros(Float64, num_cells)

   for i = 1:num_cells

       E_tan[i] = E_tan_strain(abs(ϵ[i]))

       Ai[i] = Aio[i] * E_tan[i]/E

   end

   return Ai

end

function residual!(R, U, K, F, free_dof, solution_controls)

   if solution_controls.yield_criterion == "first yield"

      # properties = inelasticity_update_first_yield(U, properties, solution_controls, free_dof)
      # K, F, free_dof = equations(properties)

   elseif solution_controls.yield_criterion =="follow stress-strain curve"
      # properties = inelasticity_update!(U, properties, solution_controls, free_dof)
      # K, F, free_dof = equations!(properties)

      # # println(sum(properties.Ai))

      # for i=1:length(F)

      #    R[i] = transpose(K[i,:]) * (U) - F[i]

      # end

      # return R

      # println("cdm")

   end

   # K, F, free_dof = equations!(properties)

   for i=1:length(F)

      R[i] = transpose(K[i,:]) * (U) - F[i]

   end

   return R

end

function calculate_strain(P, Myy, A, Iyy, cx, cxi, E)

   num_cells = length(cxi)
   ϵ_axial = zeros(Float64, num_cells)
   ϵ_flexure = zeros(Float64, num_cells)
   ϵ_total = zeros(Float64, num_cells)

   for i = 1:1:num_cells
         ϵ_axial[i] = -P/(A * E)
         ϵ_flexure[i] = Myy/(E * Iyy) * (cxi[i] - cx)
         ϵ_total[i] = ϵ_axial[i] + ϵ_flexure[i]

   end

   return ϵ_axial, ϵ_flexure, ϵ_total

end

function solve(properties, solution_controls)

   K, F, free_dof = equations!(properties)

   num_nodes=length(properties.z)

   u=zeros(num_nodes)
   v=zeros(num_nodes)
   ϕ=zeros(num_nodes)

   deformation_guess = K \ F    #consider revising this for large systems, it might be slow...

   solution = nlsolve((R, U) ->residual!(R, U, K, F, free_dof, solution_controls), deformation_guess)

   u[free_dof] = solution.zero[1:length(free_dof)]
   v[free_dof] = solution.zero[length(free_dof)+1:2*length(free_dof)]
   ϕ[free_dof] = solution.zero[2*length(free_dof)+1:3*length(free_dof)]

   return u, v, ϕ

end

function deformed_shape(xcoords, ycoords, zcoords, u, v, ϕ, w, scale_u, scale_v, scale_ϕ, scale_w)

   x_deformed = []
   y_deformed = []
   z_deformed = []
   num_sections = length(zcoords)

   u = scale_u .* u
   v = scale_v .* v
   ϕ = scale_ϕ .* ϕ
   w = scale_w .* w

   for i = 1:num_sections

       if i == 1
           x_deformed = xcoords .+ (u[i] .+xcoords.*cos.(ϕ[i]).-ycoords.*sin.(ϕ[i]))
           y_deformed = ycoords .+ (v[i] .+xcoords.*sin.(ϕ[i]).+ycoords.*cos.(ϕ[i]))
           z_deformed = zcoords[i] .+ w[i, :]
       else
           x_deformed = [x_deformed; xcoords .+ (u[i] .+xcoords.*cos.(ϕ[i]).-ycoords.*sin.(ϕ[i]))]
           y_deformed = [y_deformed; ycoords .+ (v[i] .+xcoords.*sin.(ϕ[i]).+ycoords.*cos.(ϕ[i]))]
           z_deformed = [z_deformed; zcoords[i] .+ w[i, :]]
       end

   end

   return x_deformed, y_deformed, z_deformed

end

function axial_stress(P, A)

   σ_axial = P / A

end

function flexural_stress(M, y, I)

   σ_flexural = M * y / I

end

function normal_stresses(xcoords, ycoords, zcoords, u, v, properties)

   num_sections = length(zcoords)
   num_cross_section_nodes = length(xcoords)

   σ_axial = zeros(Float64, (num_sections, num_cross_section_nodes))
   σ_flexural_xx = zeros(Float64, (num_sections, num_cross_section_nodes))
   σ_flexural_yy = zeros(Float64, (num_sections, num_cross_section_nodes))
   σ_total_normal = zeros(Float64, (num_sections, num_cross_section_nodes))

   Mxx = InternalForces.moment(properties.z, properties.dm, -v, properties.E, properties.Ix)
   Myy = InternalForces.moment(properties.z, properties.dm, -u, properties.E, properties.Iy)

   for i=1:num_sections

       A=properties.A[i]
       P=properties.P[i]
       Ix = properties.Ix[i]
       Iy = properties.Iy[i]

       #need a better way to include axial stresses in the deformed configuration
       # σ_axial[i, :] = axial_stress(-properties.P[i], properties.A[i]) * ones(Float64, num_cross_section_nodes)
       σ_flexural_xx[i, :] = flexural_stress.(Mxx[i], -(ycoords .- properties.yc[i]), properties.Ix[i])
       σ_flexural_yy[i, :] = flexural_stress.(Myy[i], -(xcoords .- properties.xc[i]), properties.Iy[i])
       #add warping stresses here one day...
       σ_total_normal[i,:] = σ_axial[i, :] + σ_flexural_xx[i, :] + σ_flexural_yy[i, :]

   end

   return σ_total_normal

end

function normal_stains(properties, σ_total_normal)

   ϵ_total_normal = σ_total_normal ./ properties.E

end

#integrate strains to calculate warping displacements

function warping_displacements(zcoords, ϵ_total_normal)

   num_sections = length(zcoords)
   num_cross_section_nodes = length(ϵ_total_normal[1,:])

   w = zeros(Float64, (num_sections, num_cross_section_nodes))
   w_total = zeros(Float64, (num_sections, num_cross_section_nodes))

   for i=1:num_sections

       for j = 1:num_cross_section_nodes

           w_total[i,j] = NumericalIntegration.integrate(zcoords[1:end], ϵ_total_normal[1:end,j])

           w[i, j] = w_total[i,j]/2 - NumericalIntegration.integrate(zcoords[1:i], ϵ_total_normal[1:i,j])  #not sure if this is general or not

       end

   end

   return w

end

end #module


text:
abstract type FDSpectralModel end

struct AdditiveFDModel{M1,M2} <: FDSpectralModel
    model1::M1
    model2::M2
end

freqdirspectra(FD::FDSpectralModel, ω, ϕ) = spreading(FD, ω, ϕ) * sdf(FD, ω)
freqdirspectra(aFD::AdditiveFDModel, ω, ϕ) = freqdirspectra(aFD.model1, ω, ϕ) + freqdirspectra(aFD.model2, ω, ϕ)

Base.:+(FD1::FDSpectralModel,FD2::FDSpectralModel) = AdditiveFDModel(FD1,FD2)


text:
using Jeff, Test, Measurements, Distributions, LinearAlgebra

X = range(0.01, 0.3, length=100)
Y = range(1, 1e-6, length=100)
DY = Y .* 0.001
DX = X .* 0.4

@testset "data_struct" begin
    q = zeros(Measurement, 100)
    R = zeros(Measurement, 100)
    for i in range(1, 100, step=1)
        q[i] = measurement(X[i], DX[i])
        R[i] = measurement(Y[i], DY[i])
    end
    data = Jeff.Data(q, R, "filename.dat")
    @test all(isapprox.(Measurements.value.(data.q), X, atol=1e-9))
    @test all(isapprox.(Measurements.value.(data.R), Y, atol=1e-9))
    @test all(isapprox.(Measurements.uncertainty.(data.q), DX, atol=1e-9))
    @test all(isapprox.(Measurements.uncertainty.(data.R), DY, atol=1e-9))
    @test isequal(data.filepath, "filename.dat")
end

@testset "data_four_col" begin
    data = Jeff.read_data(string(pwd(), Base.Filesystem.path_separator, "four_col.dat"))
    @test all(isapprox.(Measurements.value.(data.q), X, atol=1e-9))
    @test all(isapprox.(Measurements.value.(data.R), Y, atol=1e-9))
    @test all(isapprox.(Measurements.uncertainty.(data.q), DX, atol=1e-9))
    @test all(isapprox.(Measurements.uncertainty.(data.R), DY, atol=1e-9))
    @test isequal(data.filepath, string(pwd(), Base.Filesystem.path_separator, "four_col.dat"))
end

@testset "data_four_col_comma" begin
    data = Jeff.read_data(string(pwd(), Base.Filesystem.path_separator, "four_col.csv"); delim=',')
    @test all(isapprox.(Measurements.value.(data.q), X, atol=1e-9))
    @test all(isapprox.(Measurements.value.(data.R), Y, atol=1e-9))
    @test all(isapprox.(Measurements.uncertainty.(data.q), DX, atol=1e-9))
    @test all(isapprox.(Measurements.uncertainty.(data.R), DY, atol=1e-9))
    @test isequal(data.filepath, string(pwd(), Base.Filesystem.path_separator, "four_col.csv"))
end

@testset "data_three_col" begin
    data = Jeff.read_data(string(pwd(), Base.Filesystem.path_separator, "three_col.dat"))
    @test all(isapprox.(Measurements.value.(data.q), X, atol=1e-9))
    @test all(isapprox.(Measurements.value.(data.R), Y, atol=1e-9))
    @test all(isapprox.(Measurements.uncertainty.(data.q), X .* 0.05, atol=1e-9))
    @test all(isapprox.(Measurements.uncertainty.(data.R), DY, atol=1e-9))
    @test isequal(data.filepath, string(pwd(), Base.Filesystem.path_separator, "three_col.dat"))
end

@testset "data_two_col" begin
    data = Jeff.read_data(string(pwd(), Base.Filesystem.path_separator, "two_col.dat"))
    @test all(isapprox.(Measurements.value.(data.q), X, atol=1e-9))
    @test all(isapprox.(Measurements.value.(data.R), Y, atol=1e-9))
    @test all(isapprox.(Measurements.uncertainty.(data.q), X .* 0.05, atol=1e-9))
    @test all(isapprox.(Measurements.uncertainty.(data.R), Y .* 0.1, atol=1e-9))
    @test isequal(data.filepath, string(pwd(), Base.Filesystem.path_separator, "two_col.dat"))
end

@testset "data_transform_default" begin
    q = zeros(Measurement, 100)
    R = zeros(Measurement, 100)
    for i in range(1, 100, step=1)
        q[i] = measurement(X[i], DX[i])
        R[i] = measurement(Y[i], DY[i])
    end
    data = Jeff.Data(q, R, "test")
    @test all(isapprox.(Jeff.transform(data.R), log.(Y), atol=1e-9))
end

@testset "data_transform_rq4" begin
    q = zeros(Measurement, 100)
    R = zeros(Measurement, 100)
    for i in range(1, 100, step=1)
        q[i] = measurement(X[i], DX[i])
        R[i] = measurement(Y[i], DY[i])
    end
    data = Jeff.Data(q, R, "test")
    @test all(isapprox.(Jeff.transform(data.R, f=x->(log.(x) .* data.q .^ 4)), log.(Y) .* data.q .^ 4, atol=1e-9))
end


text:
struct KO_Centered{N} end

#= Kreiss-Oliger dissipation

These are operators of the form

  σ (-1)^(ord+3)/2 / 2^(ord+1) h^(ord+1) ∂^(ord+1) / ∂x^(ord+1)

note that, unlike the reference below, we add these directly to the state vector
(and *not* to its time-derivative)

cf: https://einsteintoolkit.org/thornguide/CactusNumerical/Dissipation/documentation.html
=#
function KO_Centered{N}(order::Int, sigma_diss::T, dx::T,
                        len::Int) where {T<:AbstractFloat,N}
    @assert (order + 1) % 2 == 0 "Only implemented for odd order."

    derivative_order    = order + 1
    approximation_order = 2
    stencil_length      = order + 2
    stencil_offset      = div(stencil_length-1,2)

    weights = calculate_weights(order + 1, order + 1, 0)

    a   = div(order + 3, 2)
    s   = (-1)^a

    stencil_coefs = s * sigma_diss / 2^(order + 1) .* weights

    PeriodicFD{T,N,typeof(stencil_coefs)}(derivative_order, approximation_order,
                                          dx, len, stencil_length,
                                          stencil_offset, stencil_coefs)
end

abstract type Filter{T,N} end

#= Kreiss-Oliger dissipation

These are operators of the form

  σ (-1)^(ord+3)/2 / 2^(ord+1) h^(ord+1) ∂^(ord+1) / ∂x^(ord+1)

note that, unlike the reference below, we add these directly to the state vector
(and *not* to its time-derivative)

cf: https://einsteintoolkit.org/thornguide/CactusNumerical/Dissipation/documentation.html
=#
function KO_kernel(order::Int, sigma_diss::T) where {T<:Real}
    @assert (order + 1) % 2 == 0 "Only implemented for odd order."

    kernel = calculate_weights(order + 1, order + 1, 0)

    a   = div(order + 3, 2)
    s   = (-1)^a
    fac = s * sigma_diss / 2^(order + 1)

    rmul!(kernel, fac)

    # adding one to the middle entry so that it becomes a low-pass filter
    mid = div(length(kernel) + 1, 2)
    kernel[mid] += 1

    kernel
end

function exp_kernel!(f::Vector{T}, γ::T, α::T) where {T<:Real}
    M  = length(f) - 1
    f .= [ exp( -α * ( (i-1)/M )^(γ*M) ) for i in eachindex(f) ]
    nothing
end
function exp_kernel(dim::Int, γ::T, α::T) where {T<:Real}
    f  = Array{T}(undef, dim)
    exp_kernel!(f, γ, α)
    f
end

struct ConvFilter{T<:Real,N,A} <: Filter{T,N}
    kernel :: Vector{T}
    _cache :: A
end

struct FftFilter{T<:Real,N,A,FT<:FFTW.r2rFFTWPlan} <: Filter{T,N}
    kernel   :: Vector{T}
    fft_plan :: FT
    _cache   :: A
end

struct KO_Filter{N} end

function KO_Filter{N}(order::Int, sigma_diss::T, Nxx...) where {T<:Real,N}
    kernel = KO_kernel(order, sigma_diss)
    nt = Threads.nthreads()
    _cache = [Array{T}(undef, Nxx...) for _ in 1:nt]
    ConvFilter{T,N,typeof(_cache)}(kernel, _cache)
end

KO_Filter(args...) = KO_Filter{1}(args...)

struct Exp_Filter{N} end

# α = -log(ϵ) where ϵ is the machine epsilon.
# for the standard choice of ϵ = 2^-52, α = 36.0437
function Exp_Filter{N}(γ::T, α::T, Nxx...) where {T<:AbstractFloat,N}
    kernel = exp_kernel(Nxx[N], γ, α)
    nt = Threads.nthreads()
    _cache = [Array{T}(undef, Nxx...) for _ in 1:nt]

    # REDFT00 is the DCT-I, see http://www.fftw.org/fftw3_doc/Real-even_002fodd-DFTs-_0028cosine_002fsine-transforms_0029.html#Real-even_002fodd-DFTs-_0028cosine_002fsine-transforms_0029

    # let's use only one thread per FFTW so that they don't trample over each other
    FFTW.set_num_threads(1)

    # we want to use the DCT-I since this basis is precisely the one we have by
    # using the Gauss-Lobatto grid points
    fft_plan = FFTW.plan_r2r(_cache[1], FFTW.REDFT00, N)

    FftFilter{T,N,typeof(_cache),typeof(fft_plan)}(kernel, fft_plan, _cache)
end
function Exp_Filter{N}(γ::T, Nxx...) where {T<:AbstractFloat,N}
    α = T(36.0437)
    Exp_Filter{N}(γ, α, Nxx...)
end

Exp_Filter(args...) = Exp_Filter{1}(args...)

function convolution!(fout::AbstractVector{T}, f::AbstractVector{T},
                      g::AbstractVector{T}) where {T}
    f_len = length(f)
    g_len = length(g)
    mid   = div(g_len + 1, 2)

    @fastmath @inbounds for i in eachindex(f)
        sum_i = zero(T)

        if mid <= i <= (f_len - mid + 1)
            @inbounds for aa in 1:g_len
                i_circ = i - (mid - aa)
                sum_i += g[aa] * f[i_circ]
            end
        else
            @inbounds for aa in 1:g_len
                # imposing periodicity
                i_circ = 1 + mod(i - (mid-aa) - 1, f_len)
                sum_i += g[aa] * f[i_circ]
            end
        end
        fout[i] = sum_i
    end

    fout
end

# convolution along N axis
function convolution!(fout::AbstractArray{T,M}, f::AbstractArray{T,M},
                      g::AbstractVector{T}, N::Integer) where {T,M}
    # make sure axis of convolution is contained in the dimensions of f
    @assert N <= M

    f_len = size(f,N)
    g_len = length(g)
    mid   = div(g_len + 1, 2)

    @fastmath @inbounds for idx in CartesianIndices(f)
        i     = idx.I[N] # convolution to be done along this direction
        sum_i = zero(T)

        if mid <= i <= (f_len - mid + 1)
            @inbounds for aa in 1:g_len
                i_circ = i - (mid - aa)
                I = Base.setindex(idx.I, i_circ, N)
                sum_i += g[aa] * f[I...]
            end
        else
            @inbounds for aa in 1:g_len
                # imposing periodicity
                i_circ = 1 + mod(i - (mid-aa) - 1, f_len)
                I = Base.setindex(idx.I, i_circ, N)
                sum_i += g[aa] * f[I...]
            end
        end
        fout[idx] = sum_i
    end

    fout
end

function (filter::ConvFilter)(f::AbstractVector)
    id = Threads.threadid()
    cache = filter._cache[id]
    @assert length(cache) == length(f)
    convolution!(cache, f, filter.kernel)
    copyto!(f, cache)
end

function (filter::ConvFilter{T,M})(f::AbstractArray{T}) where {T,M}
    id = Threads.threadid()
    cache = filter._cache[id]
    @assert length(cache) == length(f)
    convolution!(cache, f, filter.kernel, M)
    copyto!(f, cache)
end

# see, eg, D. Gottlieb, J.S. Hesthaven, "Spectral methods for hyperbolic
# problems", Journal of Computational and Applied Mathematics 128 (2001) 83–131,
# or "Idempotent filtering in spectral and spectral element methods", Journal of
# Computational Physics 220 (2006) 41-58
function (filter::FftFilter)(f::AbstractVector)
    M  = length(f) - 1
    id = Threads.threadid()
    cache = filter._cache[id]

    copyto!(cache, f)

    # compute the DCT-I of f
    mul!(f, filter.fft_plan, cache)

    # in momentum space, act with the filter kernel [eq (2.5) of "Idempotent filtering..."]
    # filter.kernel .* f
    @inbounds @simd for i in eachindex(f)
        # the division by 2*(N-1) is to invert the FFT, see below
        cache[i] = 0.5/M * filter.kernel[i] * f[i]
    end

    # now go back to position space using the DCT-I. the division by 2*(N-1)
    # above comes from the normalization used. DCT-I is its own inverse up to a
    # constant, and FFTW defines it as such, cf:
    # http://www.fftw.org/fftw3_doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html#g_t1d-Real_002deven-DFTs-_0028DCTs_0029
    mul!(f, filter.fft_plan, cache)
end

# filter along N axis
function (filter::FftFilter{T,N})(f::AbstractArray{T,M}) where {T,N,M}
    # make sure axis of filtering is contained in the dimensions of f
    @assert N <= M
    id = Threads.threadid()
    cache = filter._cache[id]
    copyto!(cache, f)

    Rpre  = CartesianIndices(size(f)[1:N-1])
    Rpost = CartesianIndices(size(f)[N+1:end])

    # use a function barrier here to separate the type-unstable portion
    _fftfilter!(f, cache, filter.kernel, filter.fft_plan, Rpre,
                size(f, N), Rpost)
end

# Multidimensional algorithm adapted from
# http://julialang.org/blog/2016/02/iteration
function _fftfilter!(fout, f, kernel, fft_plan::FFTW.Plan, Rpre, n, Rpost)
    # compute the DCT-I of f
    mul!(fout, fft_plan, f)

    @fastmath @inbounds for Ipost in Rpost
        @inbounds for i in 1:n
            @inbounds @simd for Ipre in Rpre
                f[Ipre,i,Ipost] = 0.5/(n-1) * kernel[i] * fout[Ipre,i,Ipost]
            end
        end
    end

    # go back to position space using the DCT-I
    mul!(fout, fft_plan, f)
end


text:
export generate

import Base: show

"""
    generate(initial::D, condition::C, iterator::I; scheduler::H = AsapScheduler()) where { D, C, I, H <: AbstractScheduler }

Generates an observable sequence by running a state-driven loop producing the sequence's elements, using the specified scheduler to send out observer messages.

# Arguments
- `initial`: initial state
- `condition`: condition to terminate generation (upon returning false)
- `iterator`: iteration step function
- `scheduler`: optional, scheduler-like object

# Note
`iterator` object should return objects of the same type as `initial`.

# Examples

```jldoctest
using Rocket

source = generate(1, x -> x < 3, x -> x + 1)
subscribe!(source, logger())
;

# output
[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Completed
```

See also: [`ScheduledSubscribable`](@ref), [`subscribe!`](@ref)
"""
function generate(initial::D, condition::C, iterator::I; scheduler::H = AsapScheduler()) where { D, C, I, H <: AbstractScheduler }
    return GenerateObservable{D, C, I, H}(initial, condition, iterator, scheduler)
end

@subscribable struct GenerateObservable{D, C, I, H} <: ScheduledSubscribable{D}
    initial   :: D
    condition :: C
    iterator  :: I
    scheduler :: H
end

getscheduler(observable::GenerateObservable) = observable.scheduler

function on_subscribe!(observable::GenerateObservable, actor, scheduler)
    value = observable.initial
    while observable.condition(value)
        next!(actor, value, scheduler)
        value = observable.iterator(value)
    end
    complete!(actor, scheduler)
    return voidTeardown
end

Base.show(io::IO, observable::GenerateObservable{D}) where D = print(io, "GenerateObservable($D)")


text:
#!/usr/bin/env julia

vepf = ARGS[1] #  combined.vep

function expand(locci)
    list = []
    pos = split(locci, ':')
        chr = pos[1]
        region = split(pos[2], '-')
        n = parse(Int64, region[1])
        while n <= parse(Int64, region[end])
             out = "$chr" * ":" * "$n"
            push!(list, out)
            n +=  1
        end
    return list
end

function conseq(fh=vepf)
    conseq = []
    file  = open(fh)
    lines = readlines(file)
    for line in lines
        if ismatch(r"^#", line)
            continue
        else
            cols = split(line, '\t')
            if ismatch(r",", cols[7])
                ids = split(cols[7], ',')
                append!(conseq, ids)
            else
                push!(conseq, cols[7])
            end
        end
    end
    conseq = unique(conseq)
end

function parseConseq(vepf=vepf, cons=conseq())
    file  = open(vepf)
    lines = readlines(file)
    for con in cons
        name = con * ".snp"
        open(name, "w")
        con_snp = []
        for line in lines
            if contains(line, con)
                cols = split(line)
                if contains(cols[2], "-")
                    list = expand(cols[2])
                    append!(con_snp, list)
                else
                    push!(con_snp, cols[2])
                end
            end
        end
        writedlm(name, unique(con_snp))
        println(con, '\t', length(unique(con_snp)))
    end
end

parseConseq()


text:
using SparseArrays
using LinearAlgebra
using Arpack

# Create diagonal using individual arrays from grid.
# Reshape first dimension diagonal from output into N dimensional grid using Ng values if needed
function kronSum(Ng::Array{Int64,1}, Grid::Array{Array{Float64}})
    dim = length(Ng)
    spd = Array{SparseMatrixCSC{Float64,Int64}}(undef,dim)
    spe = Array{SparseMatrixCSC{Float64,Int64}}(undef,dim)
    for ii=1:(dim)
        spd[ii] = spdiagm(0=>Grid[ii])
        spe[ii] = SparseMatrixCSC{Float64}(I, Ng[ii], Ng[ii])
    end
    for ii=2:dim
        spd[1] = kron(spd[1],spe[ii]) + kron(spe[1],spd[ii])
        #spe[1] = SparseMatrixCSC{Float64}(I, size(spd[1])[1], size(spd[1])[1])
    end
    return spd[1]
end

# 2D Stencil for k grid
function l2d(Ng::Array{Int64,1})
    dim = length(Ng)
    sp = SparseMatrixCSC{Float64}(I, Ng[1], Ng[1])
    E = sparse(2:Ng[1], 1:Ng[1]-1, 1, Ng[1], Ng[1])
    D = E+E' -2*sp;
    K = kron(D,sp) + kron(sp,D)
    return K
end

# Creates the spatial grid given maximum values and number of samples
function SpatialGrid(Ng::Array{Int64,1},GridMax::Array{Float64,1},GridMin::Array{Float64,1})
    dim = length(Ng)
    grid = Array{Array{Float64}}(undef,dim)
    for ii = 1:dim
        grid[ii] = collect(range(GridMin[ii], stop=GridMax[ii], length=Ng[ii]))
    end
    return grid
end

# Generate the potential operators along each respective dimension
function VHarm(Grid::Array{Array{Float64}}, Ng::Array{Int64,1}, Omega::Array{Float64,1}, m::Float64)
    dim = length(Ng)
    V = Array{Array{Float64}}(undef,dim)
    for ii = 1:dim
        V[ii] = zeros(Ng[ii]);
        for jj = 1:Ng[ii]
            V[ii][jj] = 0.5*m*Grid[ii][jj]^2*Omega[ii]^2;
        end
    end
    return V
end


text:
##############################################################################
##
## Fit is the only exported function
##
##############################################################################
function regife(df, m::FixedEffectModels.ModelTerm; kwargs...)
    regife(df, m.f; m.dict..., kwargs...)
end

function regife(df, f::FormulaTerm, vcov::CovarianceEstimator = Vcov.simple();
            weights::Union{Symbol, Nothing} = nothing,
            subset::Union{AbstractVector, Nothing} = nothing,
             method::Symbol = :dogleg,
             lambda::Number = 0.0,
             maxiter::Integer = 10_000,
             tol::Real = 1e-9,
             save::Union{Bool, Nothing} = nothing,
             contrasts::Dict = Dict{Symbol, Any}(),
             feformula::Union{Symbol, Expr, Nothing} = nothing,
            ifeformula::Union{Symbol, Expr, Nothing} = nothing,
            vcovformula::Union{Symbol, Expr, Nothing} = nothing,
            subsetformula::Union{Symbol, Expr, Nothing} = nothing)

    ##############################################################################
    ##
    ## Transform DataFrame -> Matrix
    ##
    ##############################################################################
    df = DataFrame(df; copycols = false)

    # to deprecate
    if vcovformula != nothing
        if (vcovformula == :simple) | (vcovformula == :(simple()))
            vcov = Vcov.Simple()
        elseif (vcovformula == :robust) | (vcovformula == :(robust()))
            vcov = Vcov.Robust()
        else
            vcov = Vcov.cluster(StatsModels.termvars(@eval(@formula(0 ~ $(vcovformula.args[2]))))...)
        end
    end
    if subsetformula != nothing
        subset = eval(evaluate_subset(df, subsetformula))
    end

    if  (ConstantTerm(0) ∉ FixedEffectModels.eachterm(f.rhs)) & (ConstantTerm(1) ∉ FixedEffectModels.eachterm(f.rhs))
        formula = FormulaTerm(f.lhs, tuple(ConstantTerm(1), FixedEffectModels.eachterm(f.rhs)...))
    end

    formula, formula_endo, formula_iv = FixedEffectModels.parse_iv(f)

    m, formula = parse_interactivefixedeffect(df, formula)
    if ifeformula != nothing # remove after depreciation
        m = OldInteractiveFixedEffectFormula(ifeformula)
    end

    ## parse formula
    if formula_iv != nothing
        error("partial_out does not support instrumental variables")
    end
    has_weights = (weights != nothing)

    ## create a dataframe without missing values & negative weightss
    vars = StatsModels.termvars(formula)
    if feformula != nothing # remove after depreciation
        vars = vcat(vars, StatsModels.termvars(@eval(@formula(0 ~ $(feformula)))))
    end
    factor_vars = [m.id, m.time]
    all_vars = unique(vcat(vars, factor_vars))
    esample = completecases(df[!, all_vars])
    esample .&= Vcov.completecases(df, vcov)
    if has_weights
        esample .&= BitArray(!ismissing(x) & (x > 0) for x in df[!, weights])
        all_vars = unique(vcat(all_vars, weights))
    end
    if subset != nothing
        if length(subset) != size(df, 1)
            error("df has $(size(df, 1)) rows but the subset vector has $(length(subset)) elements")
        end
        esample .&= BitArray(!ismissing(x) && x for x in subset)
    end
    main_vars = unique(vcat(vars, factor_vars))

    # Compute data needed for errors
    vcov_method_data = Vcov.materialize(view(df, esample,:), vcov)

     # Compute weights
     # Compute weights
     if has_weights
         weights = Weights(convert(Vector{Float64}, view(df, esample, weights)))
     else
         weights = Weights(Ones{Float64}(sum(esample)))
     end
     sqrtw = sqrt.(values(weights))
    for a in FixedEffectModels.eachterm(formula.rhs)
       if has_fe(a)
           isa(a, InteractionTerm) && error("Fixed effects cannot be interacted")
           Symbol(FixedEffectModels.fesymbol(a)) ∉ factor_vars && error("FixedEffect should correspond to id or time dimension of the factor model")
       end
    end
    fes, ids, formula = FixedEffectModels.parse_fixedeffect(df, formula)
    if feformula != nothing # remove after depreciation
        feformula = @eval(@formula(0 ~ $(feformula)))
        fes, ids = FixedEffectModels.oldparse_fixedeffect(df, feformula)
    end
    has_fes = !isempty(fes)
    has_fes_intercept = false
    ## Compute factors, an array of AbtractFixedEffects
    if has_fes
        if any([isa(fe.interaction, Ones) for fe in fes])
                formula = FormulaTerm(formula.lhs, tuple(ConstantTerm(0), (t for t in FixedEffectModels.eachterm(formula.rhs) if t!= ConstantTerm(1))...))
                has_fes_intercept = true
        end
        fes = FixedEffect[FixedEffectModels._subset(fe, esample) for fe in fes]
        feM = FixedEffectModels.AbstractFixedEffectSolver{Float64}(fes, weights, Val{:cpu})
    end

    has_intercept = ConstantTerm(1) ∈ FixedEffectModels.eachterm(formula.rhs)

    iterations = 0
    converged = false
    # get two dimensions

    id = FixedEffects.group(df[esample, m.id])
    time = FixedEffects.group(df[esample, m.time])

    ##############################################################################
    ##
    ## Construict vector y and matrix X
    ##
    ##############################################################################
    subdf = columntable(df[esample, unique(vcat(vars))])

    formula_schema = apply_schema(formula, schema(formula, subdf, contrasts), StatisticalModel)

    y = convert(Vector{Float64}, response(formula_schema, subdf))
    tss_total = FixedEffectModels.tss(y, has_intercept || has_fes_intercept, weights)

    X = convert(Matrix{Float64}, modelmatrix(formula_schema, subdf))

    # change default if has_regressors
    has_regressors = size(X, 2) > 0
    if save == nothing
        save = !has_regressors
    end

    # Compute demeaned X
    yname, coef_names = coefnames(formula_schema)
    if !isa(coef_names, Vector)
        coef_names = [coef_names]
    end
    yname = Symbol(yname)
    coef_names = Symbol.(coef_names)

    # demean variables
    if has_fes
        FixedEffectModels.solve_residuals!(y, feM)
        FixedEffectModels.solve_residuals!(X, feM)
    end

    ##############################################################################
    ##
    ## Estimate Model on Matrix
    ##
    ##############################################################################

    # initialize factor models at 0.1
    idpool = fill(0.1, length(id.pool), m.rank)
    timepool = fill(0.1, length(time.pool), m.rank)

    y .= y .* sqrtw
    X .= X .* sqrtw
    if !has_regressors
        # factor model
        fp = FactorModel(y, sqrtw, id.refs, time.refs, m.rank)
        fs = FactorSolution(idpool, timepool)
        (fs, iterations, converged) =
            fit!(Val{method}, fp, fs; maxiter = maxiter, tol = tol, lambda = lambda)
    else
        # interactive fixed effect
        coef = X \ y
        fp = FactorModel(y - X * coef, sqrtw, id.refs, time.refs, m.rank)
        fs = FactorSolution(idpool, timepool)
        fit!(Val{:levenberg_marquardt}, fp, fs; maxiter = 100, tol = 1e-3, lambda = lambda)

        fs = InteractiveFixedEffectsSolution(coef, fs.idpool, fs.timepool)
        fp = InteractiveFixedEffectsModel(y, sqrtw, X, id.refs, time.refs, m.rank)

        ym = copy(y)
        Xm = copy(X)

        while true
            # estimate the model
           (fs, iterations, converged) =
                fit!(Val{method}, fp, fs; maxiter = maxiter, tol = tol, lambda = lambda)
            # check that I obtain the same coefficient if I solve
            # y ~ x + γ1 x factors + γ2 x loadings
            # if not, this means fit! ended up on a a local minimum.
            # restart with randomized coefficients, factors, loadings
            newfeM = FixedEffectModels.AbstractFixedEffectSolver{Float64}(getfactors(fp, fs), weights, Val{:cpu})
            ym .= ym ./ sqrtw
            FixedEffectModels.solve_residuals!(ym, newfeM, tol = tol, maxiter = maxiter)
            ym .= ym .* sqrtw
            Xm .= Xm ./ sqrtw
            FixedEffectModels.solve_residuals!(Xm, newfeM, tol = tol, maxiter = maxiter)
            Xm .= Xm .* sqrtw
            ydiff = Xm  * (fs.b - Xm \ ym)
            if iterations >= maxiter || norm(ydiff)  <= 0.01 * norm(y)
                break
            end
            @info "Algorithm ended up on a local minimum. Restarting from a new, random, x0."
            map!(x -> randn() * x, fs, fs)
            copyto!(ym, y)
            copyto!(Xm, X)
        end
    end

    ##############################################################################
    ##
    ## Compute residuals
    ##
    ##############################################################################

    # compute residuals
    fp = FactorModel(copy(y), sqrtw, id.refs, time.refs, m.rank)
    if has_regressors
        LinearAlgebra.BLAS.gemm!('N', 'N', -1.0, X, fs.b, 1.0, fp.y)
    end
    subtract_factor!(fp, fs)
    fp.y .= fp.y ./ sqrtw
    residuals = fp.y
    ##############################################################################
    ##
    ## Compute errors
    ##
    ##############################################################################
    if !has_regressors
        rss = sum(abs2, residuals)
    else
        residualsm = ym .- Xm * fs.b
        crossxm = cholesky!(Symmetric(Xm' * Xm))
        ## compute the right degree of freedom
        df_absorb_fe = 0
        if has_fes
            for fe in fes
                df_absorb_fe += length(unique(fe.refs))
            end
        end
        dof_residual = max(size(X, 1) - size(X, 2) - df_absorb_fe, 1)

        ## estimate vcov matrix
        vcov_data = FixedEffectModels.VcovData(Xm, crossxm, residualsm, dof_residual)
        matrix_vcov = StatsBase.vcov(vcov_data, vcov_method_data)
        # compute various r2
        nobs = sum(esample)
        rss = sum(abs2, residualsm)
        _tss = FixedEffectModels.tss(ym ./ sqrtw, has_intercept || has_fes_intercept, weights)
        r2_within = 1 - rss / _tss

        rss = sum(abs2, residuals)
        r2 = 1 - rss / tss_total
        r2_a = 1 - rss / tss_total * (nobs - has_intercept) / dof_residual
    end

    ##############################################################################
    ##
    ## Save factors and loadings in a dataframe
    ##
    ##############################################################################

    if !save
        augmentdf = DataFrame()
    else
        augmentdf = DataFrame(fp, fs, esample)
        # save residuals in a dataframe
        if all(esample)
            augmentdf[!, :residuals] = residuals
        else
            augmentdf[!, :residuals] =  Vector{Union{Float64, Missing}}(missing, size(augmentdf, 1))
            augmentdf[esample, :residuals] = residuals
        end

        # save fixed effects in a dataframe
        if has_fes
            # residual before demeaning
             oldresiduals = convert(Vector{Float64}, response(formula_schema, subdf))
             oldresiduals .= oldresiduals .* sqrtw
             oldX = convert(Matrix{Float64}, modelmatrix(formula_schema, subdf))
             oldX .= oldX .* sqrtw
            if has_regressors
                LinearAlgebra.BLAS.gemm!('N', 'N', -1.0, oldX, coef, 1.0, oldresiduals)
            end
            fp = FactorModel(oldresiduals, sqrtw, id.refs, time.refs, m.rank)
            subtract_factor!(fp, fs)
            axpy!(-1.0, residuals, oldresiduals)
            # get fixed effect
            newfes, b, c = FixedEffectModels.solve_coefficients!(oldresiduals, feM; tol = tol, maxiter = maxiter)
            for j in 1:length(fes)
                augmentdf[!, ids[j]] = Vector{Union{Float64, Missing}}(missing, length(esample))
                augmentdf[esample, ids[j]] = newfes[j]
            end
        end
    end

    if !has_regressors
        return FactorResult(esample, augmentdf, rss, iterations, converged)
    else
        return InteractiveFixedEffectModel(fs.b, matrix_vcov, vcov, esample, augmentdf,
            coef_names, yname, f, nobs, dof_residual, r2, r2_a, r2_within,
            rss, sum(iterations), all(converged))
    end
end


text:
using Downhill

"""
    rosenbrock(x; b=2)

Return the value of Rosenbrock function
    ```math
    \\sum_{i=1}^n \\right[(1 - x_{i-1})^2 + b (x_i - x_{i-1})^2 \\left]
    ```
    in `n = length(x)` dimensions.
"""
function rosenbrock(x::AbstractVector{T}; b=2) where {T}
    f = zero(T) * zero(b) # type-generic zero
    for i in firstindex(x)+1:lastindex(x)
        f += (one(T) - x[i-1])^2 + b * (x[i] - x[i-1]^2)^2
    end
    return f
end

"""
    drosenbrock(x; b=2)

Return the gradient of Rosenbrock function
    ```math
    \\sum_{i=1}^n \\right[(1 - x_{i-1})^2 + b (x_i - x_{i-1})^2 \\left]
    ```
    in `n = length(x)` dimensions.
"""
function drosenbrock(x::AbstractVector{T}; b=2) where {T}
    g = fill!(similar(x, promote_type(T, typeof(b))), 0)
    for i in firstindex(x)+1:lastindex(x)
        g[i-1] += 2 * (x[i-1] - 1) + 4 * b * x[i-1] * (x[i-1]^2 - x[i])
        g[i] += 2 * b * (x[i] - x[i-1]^2)
    end
    return g
end

"""
    rosenbrock!(x, g; b=2)

Return the value of Rosenbrock function in `length(x)` dimensions.
"""
function rosenbrock!(x::AbstractVector, g::AbstractVector; b=2)
    f = zero(eltype(g))
    fill!(g, 0)
    inds = eachindex(x, g)
    for i in 2:last(inds)
        f += (1 - x[i-1])^2 + b * (x[i] - x[i-1]^2)^2
        g[i-1] += 2 * (x[i-1] - 1) + 4 * b * x[i-1] * (x[i-1]^2 - x[i])
        g[i] += 2 * b * (x[i] - x[i-1]^2)
    end
    return f, g
end

ros(x, g) = (rosenbrock(x), drosenbrock(x))

ros!(x, g) = rosenbrock!(x, g)

ans_nonmutating = let x0 = zeros(2)
    opt = BFGS(x0)
    optresult = optimize!(ros, opt, x0; maxiter=1000)
    println("""
        ==Nonmutating gradient evaluation==

        Optimization converged: $(optresult.converged)
        Final argument: $(optresult.argument)
        Final gradient: $(optresult.gradient)
        Iteration count: $(optresult.iterations)
        Call count: $(optresult.calls)
    """)
    optresult.argument
end

ans_mutating = let x0 = zeros(2)
    opt = BFGS(x0)
    optresult = optimize!(ros!, opt, x0; maxiter=1000)
    println("""
        ==Mutating gradient evaluation==

        Optimization converged: $(optresult.converged)
        Final argument: $(optresult.argument)
        Final gradient: $(optresult.gradient)
        Iteration count: $(optresult.iterations)
        Call count: $(optresult.calls)
    """)
    optresult.argument
end

@assert ans_nonmutating ≈ ans_mutating


text:
module MLJTime

# IMPORTS
using IndexedTables, CSVFiles, ZipFile
using  MLJBase: fit!, predict, machine, partition, l1, CV, range, cross_entropy,
       predict_mode, accuracy, predict_mean
import MLJModelInterface
import MLJModelInterface: @mlj_model, fit, predict, matrix
using StableRNGs

# EXPORTS
export RandomForestClassifierFit, InvFeatureGen, predict_new, predict_mean,
       fit!, predict, machine, fitted_params, predict_mode, accuracy

export univariate_datasets, load_dataset, load_from_tsfile_to_NDArray, dwt,
       load_ts_file, ts_dataset, partition, matrix, dtw_distance, load_NDdataset

export l1, StableRNG, CV, range, cross_entropy
#CONSTANTS
# the directory containing this file: (.../src/)
const MODULE_DIR = dirname(@__FILE__)
const MMI = MLJModelInterface

# Includes
include("interval_based_forest.jl")
include("datapath.jl")
include("interface.jl")
include("distances_based_classifier.jl")
include("rocket.jl")

end # module


text:
function DataFrame(x::TQueryResult)

    #get column name
    colnames = [desc.col_name for desc in x.row_set.row_desc]

    #get column types and nullable
    coltypes = [(getcolumntype(desc.col_type._type), desc.col_type.nullable) for desc in x.row_set.row_desc]

    #collapse vectors into a single Vector{T, Missing} vector
    mergecols = [squashbitmask(y, tup) for (y, tup) in zip(x.row_set.columns, coltypes)]

    #convert to DataFrame
    df = DataFrame(mergecols)
    rename!(df, Symbol.(colnames))

    return df

end

#iterate over propertynames and build dict, make parametric?
DataFrame(x::TDBInfo) = DataFrame(Dict(:db_name => x.db_name, :db_owner => x.db_owner))

DataFrame(x::TTableMeta) = DataFrame(Dict(:table_name => x.table_name,
                                          :num_cols => x.num_cols,
                                          :is_view => x.is_view,
                                          :is_replicated => x.is_replicated,
                                          :shard_count => x.shard_count,
                                          :max_rows => x.max_rows
                                          )
                                    )

DataFrame(x::TDashboard) = DataFrame(Dict(:dashboard_name => x.dashboard_name,
                                          :dashboard_state => x.dashboard_state,
                                          :image_hash => x.image_hash,
                                          :update_time => x.update_time,
                                          :dashboard_metadata => x.dashboard_metadata,
                                          :dashboard_id => x.dashboard_id,
                                          :dashboard_owner => x.dashboard_owner,
                                          :is_dash_shared => x.is_dash_shared
                                          )
                                    )

#This one digs into col_type, so slightly different than others that might be generated
DataFrame(x::TColumnType) = DataFrame(Dict(:col_name => x.col_name,
                                           :col_type => x.col_type._type,
                                           :encoding => x.col_type.encoding,
                                           :nullable => x.col_type.nullable,
                                           :is_array => x.col_type.is_array,
                                           :precision => x.col_type.precision,
                                           :scale => x.col_type.scale,
                                           :comp_param => x.col_type.comp_param,
                                           :size => x.col_type.size,
                                           :is_reserved_keyword => x.is_reserved_keyword,
                                           :src_name => x.src_name,
                                           :is_system => x.is_system,
                                           :is_physical => x.is_physical
                                          )
                                      )

#This is slightly awkward, in that the added columns make it seem like these are row-level info
#In reality, these fields are table-level
#TODO: think about how to better represent
function DataFrame(x::TTableDetails)

    tmp = DataFrame(x.row_desc)
    tmp[!, :fragment_size] .= x.fragment_size
    tmp[!, :page_size] .= x.page_size
    tmp[!, :max_rows] .= x.max_rows
    tmp[!, :view_sql] .= x.view_sql
    tmp[!, :shard_count] .= x.shard_count
    tmp[!, :key_metainfo] .= x.key_metainfo
    tmp[!, :is_temporary] .= x.is_temporary
    tmp[!, :partition_detail] .= x.partition_detail

    return tmp
end

#After individual types have DataFrame methods written, this method concatenates
DataFrame(x::Vector{<:Union{TDBInfo, TTableMeta, TDashboard, TColumnType}})  = vcat(DataFrame.(x)...)

#Need to validate this always returns a Vector length 1
#Figure out how to parameterize if possible
function Base.show(io::IO, ::MIME"text/plain", m::Vector{TServerStatus})

    props = propertynames(m[1])

    println(io, "$(eltype(m))\n")
    for sym in props
        println(io, """  $(sym): $(getfield(m[1], sym))""")
    end

end

#REPL display; show method for Juno uses inline tree display
Base.show(io::IO, ::MIME"text/plain", m::OmniSciConnection) = println(io, "Connected to $(m.c.p.t.host):$(m.c.p.t.port)")


text:
abstract type Neo4jBoltError <: Exception end
abstract type AbstractSessionError <: Neo4jBoltError end

struct SessionError <: AbstractSessionError
    session
    args
end

struct SessionExpired <: AbstractSessionError
    session
    args
end

struct TransactionError <: Neo4jBoltError
    transaction
    args
end

struct ProtocolError <: Neo4jBoltError
    driver
    args
end


text:
abstract type AbstractAggregator end

#=
INSTRUCTIONS FOR EXTENDING ABSTRACT CES AGGREGATOR

Subtypes should have:
- A type `shares` whose length determines the number of aggregated entities (e.g. goods)
- If you  want to use `aggr`, `gradient_aggr`, `partial_aggr` defined for AbstractCESAggregator:
+ `_aggr(x, t::SubType, NoSkip)`  should evaluate the aggregator at (iterable) point x.
                                  `aggr` will default to this specification.
+ `_aggr(x, t::SubType, leaveout::Int64)` should evaluate the aggregator leaving out element `j`

+ `_partial_aggr(x, j, Z, t::SubType)` should evaluate the aggregator's j-th partial derivative at,
                                       point x.
                                       The argument Z is the evaluated aggregator. Even if it's not
                                       not necessary, leave it there anyhow.
Otherwise, you have to define those methods yourself.

Moreover, users might expect the following:
+ `invpartial_aggr(x,j,k,t::SubType)` should solve the equation ∂ⱼF(ξ, x₋ⱼ) = k for ξ
=#
abstract type AbstractCESAggregator <: AbstractAggregator end
ngoods(aa::AbstractCESAggregator)=length(aa.shares)
share(j, aa::AbstractCESAggregator)=aa.shares[j]

aggr(x, aa::AbstractCESAggregator)=_aggr(x, aa, NoSkip)
# aggr(x, aa::AbstractCESAggregator, leaveout)=_aggr(x, aa, leaveout)

" Partial derivative of aggregator wrt coordinate `j`"
partial_aggr(x,j,aa::AbstractCESAggregator)=_partial_aggr(x,j,aggr(x, aa),aa)

function gradient_aggr!(gr, x, aa::AbstractCESAggregator)
    ev = aggr(x, aa)
    @inbounds for g in 1:length(aa)
        gr[g] = _partial_aggr(x, g, ev, aa)
    end
end

"Aggregator gradient (allocates) "
function gradient_aggr(x,aa::AbstractCESAggregator)
    gr = fill(0., ngoods(aa))
    gradient_aggr!(gr, x, aa)
    return gr
end

"Aggregator gradient (returns tuple)"
function gradient_aggr_t(x,aa::AbstractCESAggregator)
    ev = aggr(x, aa)
    map(g -> _partial_aggr(x, g, ev, aa), 1:ngoods(aa))
end

struct NoSkip end
# Inlining breaks the code. Seems to make execution (a bit) slower. Perhaps
# revert to regular checking?
@noinline cond_include(x, g, neutral, ::Type{NoSkip})=x
@noinline cond_include(x, g, neutral, n::Int64)=ifelse(g == n, neutral, x)
# @noinline cond_include(x, g, neutral, nn::NTuple{N, Int64}) where {N} = ifelse(g in nn, neutral, x)

# ARMINGTON AGGREGATOR (longest section) ---------------------------------------
struct Armington{N} <: AbstractCESAggregator
    es::Float64 # elasticity of substitution
    shares::NTuple{N, Float64}
    function Armington(es, shares)
        es == 1 && error("Please use Cobb-Douglas aggregator")
        any(shares .< 0.) && error("Negative shares not allowed")
        new{length(shares)}(es, Tuple(shares))
    end
end

function _aggr(x, aa::Armington, leaveout)
    @unpack shares, es = aa
    pow = 1 - 1/es

    r = 0.
    @inbounds for g in 1:ngoods(aa)
        # Include only if not left out
        r += cond_include(shares[g]*x[g]^pow, g, 0., leaveout)
    end
    return r^(1/pow)
end

doc"""
    invpartial_aggr(x,j,k,a::Armington)
Solve the equation
```math
\partial_j A(x_j, x_{-j}) = k
```
for xⱼ.

Note that entry x[j] is redundant.
"""
function invpartial_aggr(x,j,k,aa::Armington)
    pow = 1 - 1/aa.es
    vv = aggr(x, aa, j)

    t1 = (k/share(j, aa))^(aa.es-1)-share(j, aa)
    t1 ^= -1/pow

    return vv*t1
end

" Specialize inverse partial aggregate for 2 goods "
function invpartial_aggr(x, j, k, aa::Armington{2})
    pow = aa.es/(aa.es - 1)
    pow2 =  aa.es - 1

    w_mine = share(j, aa)
    w_other = share(3-j, aa)
    x_other = x[3-j]

    # Just check the formula... :(
    r = w_other/((k/w_mine)^pow2 - w_mine)
    r ^= pow
    r *= x_other
end

doc"""

Take advantage of the fact that each partial is

```math
\partial_j \mathcal A = \omega_j x_j^{-\frac{1}{\xi}}
\left[ \sum \omega_i x_i^{1-\frac{1}{\xi}}
      \right]^{\frac{1/\xi}{1 - 1/\xi}}
```

The term in square brackets corresponds to ``\mathcal A(x)^{\frac{1}{\xi}}``. In
`_partial_aggr`, we name this `aggr_val_e` (e stands for "along with exponent").

The extra argument for the function value is useful when computing the gradient
(need only compute it once).

"""
function _partial_aggr(x, j, aggr_val_e, aa::Armington)
    pow = 1/aa.es
    r = x[j]^(-pow) * share(j, aa)
    r *= aggr_val_e^pow
end

# COBB-DOUGLAS AGGREGATOR ------------------------------------------------------

# This is a special case of Armington when shares add up to one, and the
# elasticity of substitution parameter equals 1.

struct CobbDouglas{N} <: AbstractCESAggregator
    shares::NTuple{N, Float64}
    function CobbDouglas(shares)
        !all(shares .>= 0) && error("Please specify only positive shares")
        # !(sum(shares) == 1.) && error("Shares must add up to one")
        new{length(shares)}(shares)
    end
end

CobbDouglas(a::Float64)=CobbDouglas((a, 1-a))

function _aggr(x, cbd::CobbDouglas, leaveout)
    r = one(eltype(x))
    for i in 1:ngoods(cbd)
        r *= cond_include(x[i]^cbd.shares[i], i, 1., leaveout)
    end
    return r
end

function invaggr(x,j,y,cbd::CobbDouglas)
    r = y
    r /=  _aggr(x, cbd, j) # leave out entry j
    r ^= 1/share(j, cbd)
end

function invpartial_aggr(x,j,k,cbd::CobbDouglas)
    (k/shares(j, cbd)/_aggr(x, cbd, j))^(1/(shares(j, cbd)-1))
end

_partial_aggr(x, j, aggr_val, cbd::CobbDouglas)=shares(j, cbd)*aggr_val/x[j]

# Specialize 2 goods
aggr(x, cbd::CobbDouglas{2})=x[1]^share(1, cbd)*x[2]^share(2, cbd)

function invaggr(x, j, y, cbd::CobbDouglas{2})
    j_other = 3 - j
    r = y/x[j_other]^share(j_other, cbd)
    r ^= 1/share(j, cbd)
end

function partial_aggr(x, j, cbd::CobbDouglas{2})
    j_other = 3 - j
    share_mine = share(j, cbd)
    share_other = share(j_other, cbd)
    x_other = x[j_other]

    share_mine * x[j]^(share_mine-one(share_mine)) * x_other^share_other
end

function invpartial_aggr(x,j,k,cbd::CobbDouglas{2})
    aux = share(j, cbd) * x[3-j]^share(3-j, cbd)
    (k / aux)^(1/(share(j,cbd) - 1))
end

# LINEAR AGGREGATOR   ----------------------------------------------------------
# This when the elasticity of substitution is infinite.

struct AggLinear{N} <: AbstractCESAggregator
    shares::NTuple{N, Float64}
end

# aggr(x, la::AggLinear)=sum(la.shares.*x)
function aggr(x, la::AggLinear)
    for g in 1:ngoods(la)
        la.shares.*x
    end
end

invpartial_aggr(x, j, k, la::AggLinear)=error("Can't invert constant function")
partial_aggr(x,j,la::AggLinear)=la.shares[j]
gradient_aggr(x, la::AggLinear)=collect(la.shares)


text:
function Sonum8(x::Float64)
    if isfinite(x)
        if signbit(x)
            return -Sonum8(UInt8(findFirstSmaller(-x,bounds8)-2))
        else
            return Sonum8(UInt8(findFirstSmaller(x,bounds8)-2))
        end
    else    # Inf,-Inf, and NaN
        return notareal(Sonum8)
    end
end

function Sonum16(x::Float64)
    if isfinite(x)
        if signbit(x)
            return -Sonum16(UInt16(findFirstSmaller(-x,bounds16)-2))
        else
            return Sonum16(UInt16(findFirstSmaller(x,bounds16)-2))
        end
    else    # Inf, -Inf, and NaN
        return notareal(Sonum16)
    end
end

#TODO Float64(notareal(Sonum8)) is -Inf, problem?
Float64(x::Sonum8) = signbit(x) ? -sonum8[Int(-x)+1] : sonum8[Int(x)+1]
Float64(x::Sonum16) = signbit(x) ? -sonum16[Int(-x)+1] : sonum16[Int(x)+1]

Float32(x::Sonum8) = Float32(Float64(x))
Float32(x::Sonum16) = Float32(Float64(x))

Sonum8(x::Float32) = Sonum8(Float64(x))
Sonum16(x::Float32) = Sonum16(Float64(x))

Sonum8(x::Float16) = Sonum8(Float64(x))
Sonum16(x::Float16) = Sonum16(Float64(x))

Sonum8(x::Int64) = Sonum8(Float64(x))
Sonum16(x::Int64) = Sonum16(Float64(x))

promote_rule(::Type{Int64},::Type{Sonum16}) = Sonum16
promote_rule(::Type{Int64},::Type{Sonum8}) = Sonum8

promote_rule(::Type{Int32},::Type{Sonum16}) = Sonum16
promote_rule(::Type{Int32},::Type{Sonum8}) = Sonum8


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

relaxed_tols = (default_tol_relax = 100,)
insts = OrderedDict()
insts["minimal"] = [
    ((2, 2),),
    ((2, 2), :SOCExpPSD),
    ]
insts["fast"] = [
    ((:motzkin2,),),
    ((:motzkin2,), :SOCExpPSD),
    ((:motzkin3,),),
    ((:CS16ex8_13,),),
    ((:CS16ex8_14,),),
    ((:CS16ex18,),),
    ((:CS16ex12,),),
    ((:CS16ex13,),),
    ((:MCW19ex1_mod,),),
    ((:MCW19ex8,),),
    ((:MCW19ex8,), :SOCExpPSD),
    ((3, 2),),
    ((3, 2), :SOCExpPSD),
    ((6, 6),),
    ((20, 3),),
    ((20, 3), :SOCExpPSD),
    ]
insts["various"] = [
    ((:motzkin2,),),
    ((:motzkin3,),),
    ((:CS16ex8_13,),),
    ((:CS16ex8_14,),),
    ((:CS16ex18,),),
    ((:CS16ex12,),),
    ((:CS16ex13,),),
    ((:MCW19ex1_mod,),),
    ((:MCW19ex8,), nothing, relaxed_tols),
    ((6, 3),),
    ((6, 3), :SOCExpPSD),
    ((20, 3), nothing, relaxed_tols),
    ((20, 3), :SOCExpPSD, relaxed_tols),
    ]
return (SignomialMinJuMP, insts)


text:
using Printf
using Dates
using Base.Threads
#=
    记录、更新当前的最新实时状态，并记录当前状态的时间；
    如果是第一次调用此函数，则更新Lxtime、final_time到当前时间。
=#
function l0_handle(aobj::timeAnalsys ,curtime, reatime_eventOn::Bool,FrameCnt::UInt )
    # curtime = UInt( DT2UInt( dt ) )

    aobj.AvgValue = 0
    aobj.L0Status = reatime_eventOn

    # curtime = DT2UInt( dt )
    if aobj.L0time == UInt( 0 )
        # println( "-----",aobj.L0time ," DT:",dt )
        aobj.L0time = curtime
        aobj.L1time = curtime
        aobj.L2time = curtime
        aobj.final_time = curtime
        # aobj.FrameCnt

        aobj.L1Frame=FrameCnt
        aobj.L2Frame=FrameCnt
        aobj.L3Frame=FrameCnt
        aobj.L4Frame=FrameCnt
        # exit(-1)
    else
        aobj.L0time = curtime
        aobj.L0Frame=FrameCnt
        # println( "---Not Zero  --",aobj.L0time ," DT:",dt )
        # exit(1)
    end
# println( curtime ," ,,", aobj.final_time)
end
#=
    如果L1的状态与L0状态不同，则更新L1时间到L0时间，并同步更新L1状态
    L1 时间用于做后续的L1~ Lx、final处理
=#
function l1_handle( aobj::timeAnalsys )
    if aobj.L1Status  != aobj.L0Status
        # @debug  aobj.L1Status ,aobj.L0time  ,"l1_handle"
        aobj.L1time = aobj.L0time
        aobj.L1Status  = aobj.L0Status
        aobj.L1Frame=aobj.L0Frame
    end
end

#=
    处理Micro时间的过滤、合并。主要处理对象为 L1 => L2
    合并： 指在告警(L2 ==  true)为真时，在之后的某段时间：
       如果非告警状态(L1 == false )连续持续时间小于该合并时间,则认为告警状态一直持续，
       如果持续时间超过该时间，则认为告警停止。
    过滤：告警(L2 ==  true)为 否 时，在之后的某段时间：
       如果非告警状态(L1 == 真 )连续持续时间 大于 该过滤时间,则认为告警状态需要更新为 真，
        否则就不改变原L2状态
=#
function l2_handle( aobj::timeAnalsys )
    MergetimePeriod ,filterPeriod = getSegMicroTimeParaValue(glbSegAlarmParasDictObj,aobj.segID)
    L1_persistTime =  aobj.L0Frame - aobj.L1Frame
    #
    if aobj.L1Status == true  && L1_persistTime  >= filterPeriod
        if aobj.L2Status == true
            ;
        else  ## L2 状态转换
            aobj.L2Status = true
            aobj.L2time   = aobj.L1time
            aobj.L2Frame=aobj.L1Frame
            aobj.L3Status = false # 解除向下发展态势  L3 状态和时间 代表了 L1 向下发展的状态和开始时间，
        end
    end

    if aobj.L1Status == false  && aobj.L2Status  == true # 向下发展态势 L3 状态和时间 代表了 L1 向下发展的状态和开始时间，
        if  aobj.L3Status == false
            aobj.L3Status = true
            aobj.L3time = aobj.L1time
            aobj.L3Frame=aobj.L1Frame
        end
    end

    willDownLastingTime  = aobj.L0time -  aobj.L3time   #L3 状态和时间 代表了 L1 向下发展的状态和开始时间，
    if willDownLastingTime >= MergetimePeriod ## 不再等待 Merge
        if aobj.L2Status == true  && aobj.L1Status == false  ## L2 状态转换
            aobj.L2Status = false
            aobj.L2time   = aobj.L3time
            aobj.L2Frame=aobj.L3Frame
        end
    end
end

#=  l3_handle:
 虚拟处理
=#
function l3_handle( aobj::timeAnalsys )
     # L3 状态和时间 代表了 L1 向下发展的状态和开始时间，
end

function finalStatu( aobj::timeAnalsys ) ##
    aobj.final_Status = aobj.L4Status
    aobj.final_time = aobj.L4time
    # $(Main.machineName),$(aobj.segID), $(aobj.L0time), $(aobj.L4time) ,$( aobj.L4time-aobj.L0time )
end

function Alarmhappened( aobj::timeAnalsys  )
    finalStatu( aobj )
    if aobj.segID > 1 && aobj.segID < 4096
        # @info "AlarmOn,$(UInt2DTStr(aobj.L0time)),$(Main.machineName),$(aobj.L0Frame),$(aobj.segID), $(aobj.L0time), $(aobj.L2time) ,$( aobj.L0time-aobj.L2time ) "
    end
end
function AlarmStopped( aobj::timeAnalsys  )
    finalStatu( aobj )
    # @debug  "AlarmOff,$(Main.machineName),$(aobj.segID), $(aobj.L0time)"
    # exit(-1)
end
#=  l4_handle:
    处理Macro时间的处理。主要处理对象为 L3 => L4
    对于L3 为告警状态的：
        如果L3告警状态持续超过 最小告警时间 ，而且小于最大告警时间，才更新L4为告警状态。
        否则L4为无告警状态
    对于L3为 无告警状态的，实时将L4同步为L3状态。
=#
function l4_handle(aobj::timeAnalsys )
    L2_persistTime =  aobj.L0Frame -aobj.L2Frame
    # aobj.L3Frame=aobj.L1Frame
    # @info L2_persistTime , aobj.L2Status, aobj.L4Status
    if aobj.L2Status == true
        maxTimePeriod ,minTimePeriod = getSegMacroTimeParaValue(glbSegAlarmParasDictObj,aobj.segID)
        if L2_persistTime <  minTimePeriod
            # aobj.L4Status = false
            # aobj.L4time =  aobj.L0time
           # EventL1.ReportEventL1(UInt(aobj.segID),UInt(aobj.segID), aobj.final_time ,aobj.L0Frame ,"MacroTimeShort")
        elseif L2_persistTime  >= minTimePeriod && L2_persistTime <= maxTimePeriod &&  aobj.L4Status == false ## report the Alarm
                aobj.L4Status = true
                aobj.L4time =  aobj.L0time
                aobj.L4Frame =  aobj.L0Frame
                # L1Frame
                Alarmhappened( aobj  )
        elseif L2_persistTime > maxTimePeriod  && aobj.L4Status == true   ## cancel the Alarm
            aobj.L4Status = false
            aobj.L4time   = aobj.L0time
            aobj.L4Frame =  aobj.L0Frame

            AlarmStopped( aobj)
        elseif L2_persistTime > maxTimePeriod
            EventL1.ReportEventL1(UInt(aobj.segID),UInt(aobj.segID), aobj.final_time,aobj.L0Frame ,"MacroTimeLong")
        end
    else #L2Status == false
        if aobj.L4Status == true  ## stop the Alarm
            aobj.L4Status = false
            aobj.L4time   = aobj.L0time
            aobj.L4Frame =  aobj.L0Frame
            AlarmStopped( aobj )
        end
    end
end

function updateStatus(aobj::timeAnalsys,curtime,reatime_eventOn::Bool,FrameCnt::UInt )
    # @info " reatime_eventOn ",reatime_eventOn ,curDTStr
    l0_handle( aobj ,curtime, reatime_eventOn ,FrameCnt )
    l1_handle( aobj )
    l2_handle( aobj )
    l3_handle( aobj )
    l4_handle( aobj )
end

function ShowtimeSchmidt( timeAnalsysObj,dt)
    if timeAnalsysObj.final_Status == true
        # @debug ( timeAnalsysObj.final_Status ,timeAnalsysObj.final_time , dt,timeAnalsysObj.segID)
    end
end

# function proloop1(AlarmDynStatus,curAlarmsinSEG,FrameCnt,curtime)
#     maxl=  length( AlarmDynStatus )
#      @threads for segID = 1 : maxl
#           reatime_eventOn  = AlarmDynStatus[ segID ]
#         timeAnalsysObj   = curAlarmsinSEG[ segID ]
#         updateStatus( timeAnalsysObj, curtime , reatime_eventOn,FrameCnt   )
#      end
#  end
 function subproloop2(AlarmDynStatus,curAlarmsinSEG,FrameCnt,curtime,rng)
        for segID = rng
            reatime_eventOn  = AlarmDynStatus[ segID ]
            timeAnalsysObj   = curAlarmsinSEG[ segID ]
            updateStatus( timeAnalsysObj, curtime , reatime_eventOn,FrameCnt   )
        end
 end
 function proloop2(AlarmDynStatus,curAlarmsinSEG,FrameCnt,curtime)
    maxl=  length( AlarmDynStatus )
    parCnt = 4-1;
    loopmax =floor( Int, maxl/parCnt)

#    println( " spawn loope in  proloop2 ")
    rng1=1:loopmax
     a1=Base.Threads.@spawn subproloop2(AlarmDynStatus,curAlarmsinSEG,FrameCnt,curtime,rng1)

     rng2=1+loopmax : 2*loopmax
     #rng2=1+loopmax : 1+loopmax
     a2=Base.Threads.@spawn subproloop2(AlarmDynStatus,curAlarmsinSEG,FrameCnt,curtime,rng2)

     rng3=1+loopmax*2:3*loopmax
     a3=Base.Threads.@spawn subproloop2(AlarmDynStatus,curAlarmsinSEG,FrameCnt,curtime,rng3)
    if maxl  >= (1+loopmax*3)
        rng4=1+loopmax*3:maxl
        a4=Base.Threads.@spawn subproloop2(AlarmDynStatus,curAlarmsinSEG,FrameCnt,curtime,rng4)
        wait(a4);
    end

    #println("-----------------",rng1 ," ---",rng2," ",rng3," ",rng4)
     wait(a1);
     wait(a2);
     wait(a3);

 end
function FeaturetimeSchmidt( dt , FrameCnt::UInt)
    # @info "Enter into FeaturetimeSchmidt"
    curtime = UInt( DT2UInt( dt ) )

    curAlarmsinSEG[ 1 ].L0Frame = FrameCnt
    curAlarmsinSEG[ 1 ].L0time  = curtime
    # println( " spawn loope in  proloop2 ")
    # @time
      proloop2(AlarmDynStatus,curAlarmsinSEG,FrameCnt,curtime)

end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
# -----------------------------------------------------
# This script gives the methods for maniupulating the
# Cartesian indices corresponding to multilevels that
# are either kept or cut in various schemes of
# "sparsification"
# -----------------------------------------------------

# Efficiency Criticality: LOW
# This is not evaluated very often

# Accuracy Criticality: N/A
# No float manipulation

# Gives the appropriate boolean cutoff corresponding
# to a given scheme, e.g. sparse basis, full basis,
# and in the future possibly the energy basis of
# Bungartz and Griebel

function cutoff(::Val{:sparse}, x::CartesianIndex{D}, n::Int) where {D}
    sum(x.I) > n+D
end

function cutoff(::Val{:full}, x::CartesianIndex{D}, n::Int) where {D}
    false
end