text:
#=
    Logistic regresiion
    Inspired by WebPPL's examples/logisticRegression.wppl

    Cf logistic_regression.jl

    In this model we actually use logicstic().

      parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
          Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

               m    1.0285    0.5627     0.0178    0.0254    687.0390    1.0014      166.1119
               b    0.1912    0.9294     0.0294    0.0342    705.2296    0.9995      170.5100
  labels_post[1]    0.0180    0.1330     0.0042    0.0044    821.3831    0.9990      198.5936
  labels_post[2]    0.0610    0.2395     0.0076    0.0087    873.0023    0.9990      211.0741
  labels_post[3]    0.8440    0.3630     0.0115    0.0075    858.1636    0.9993      207.4864
  labels_post[4]    0.9740    0.1592     0.0050    0.0048    847.3995    0.9992      204.8838
  labels_post[5]    0.9880    0.1089     0.0034    0.0029   1028.5041    0.9993      248.6712

=#

using Turing, StatsPlots, DataFrames
using StatsFuns

include("jl_utils.jl")

@model function logistic_regression2(x, labels)
    m ~ Normal(0, 1)
    b ~ Normal(0, 1)

    n = length(x)
    v = tzeros(n)
    for i in 1:n
        v[i] = logistic(b + m*x[i])
        labels[i] ~ Bernoulli(v[i])

    end

    # Post predictive of labels
    labels_post = tzeros(n)
    for i in 1:n
        labels_post[i] ~ Bernoulli(v[i])
    end
end

xs = [-10, -5, 2, 6, 10]
labels = [0, 0, 1, 1, 1]
model = logistic_regression2(xs,labels)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(15), 1_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
export SparseHDF5Matrix, extractdense, extractsparse, sparse, issparse
import HDF5
import SparseArrays

"""
This class implements a `AbstractArray` wrapper around a sparse matrix in a HDF5 file.
Instances of this class can be used to represent a HDF5 dataset inside collections that expect an `AbstractArray`.
Most importantly, this class does not load any data into memory on construction (and only minimal loading in its `show()` method).
This allows users to manipulate arbitrarily large datasets in low-memory environments.

HDF5 does not actually implement a native sparse format (at least, not of time of writing).
Instead, we expect the 10X Genomics format where the matrix contents are stored in a group containing:

- `shape`, an integer dataset of length 2 containing the number of rows and columns in the matrix.
- `data`, a numeric dataset containing all of the non-zero elements in compressed sparse column (CSC) order.
- `indices`, an integer dataset containing the 0-based row indices of all non-zero elements in CSC order.
- `indptr`, an integer dataset containing the 0-based pointers into `indices` for the start and end of each column.

As with the [`DenseHDF5Array`](@ref), we recommend using instances of this class as a placeholder for a real matrix, rather than computing directly on it.
Once values are explicitly needed, the entire array can be loaded into memory with [`extractsparse`](@ref).
"""
struct SparseHDF5Matrix{Tv,Ti} <: HDF5Array{Tv, 2}
    file::String
    name::String
    dims::Tuple{Integer,Integer}
    ptrs::Vector{Integer}

    # Mutable cache. Technically we need to protect this from multiple threads,
    # but whatever, it's not like we can do parallel HDF5 access anyway.
    cache::Dict{Integer,SparseArrays.SparseVector{Tv,Ti}}
end

"""
    SparseHDF5Matrix(file, name)

Create a `HDF5Array` for a sparse matrix inside a HDF5 file `file` at group `name`.
The group is expected to follow the 10X Genomics format, see [`SparseHDF5Matrix`](@ref) for more details.

# Examples
```jldoctest
julia> using HDF5Arrays

julia> tmp = tempname();

julia> examplesparse(tmp, "stuff", (20, 10), 0.2)

julia> x = SparseHDF5Matrix(tmp, "stuff");

julia> size(x)
(20, 10)

julia> using SparseArrays

julia> SparseArrays.issparse(x)
true
```
"""
function SparseHDF5Matrix(file::String, name::String)
    handle = HDF5.h5open(file)
    group = handle[name]

    if !isa(group, HDF5.Group)
        throw(ErrorException("'" * name * "' in '" * file * "' should be a group"))
    end

    dims = HDF5.read(group["shape"])
    dtype = HDF5.get_jl_type(group["data"])
    itype = HDF5.get_jl_type(group["indices"])
    ptrs = HDF5.read(group["indptr"])

    return SparseHDF5Matrix{dtype,itype}(file, name, (dims...,), ptrs, Dict{Integer,SparseArrays.SparseVector{dtype,itype}}())
end

"""
    size(x::SparseHDF5Matrix{Tv,Ti})

Get the size of a `DenseHDF5Array{T,N}` as a `N`-tuple.
"""
function Base.size(x::SparseHDF5Matrix{Tv,Ti}) where {Tv,Ti}
    return x.dims
end

"""
    issparse(x::SparseHDF5Matrix{Tv,Ti})

Returns true.
"""
function SparseArrays.issparse(x::SparseHDF5Matrix{Tv, Ti}) where {Tv, Ti}
    return true
end

"""
    getindex(x::SparseHDF5Matrix{Tv,Ti}, i, j)

Get the value of `x[i,j]`.
It would be unwise to use this function for anything other than `show()` -
we suggest using [`extractsparse`](@ref) to obtain larger blocks of data.

# Examples
```jldoctest
julia> using HDF5Arrays

julia> tmp = tempname();

julia> examplesparse(tmp, "stuff", (20, 10), 0.2)

julia> x = SparseHDF5Matrix(tmp, "stuff");

julia> typeof(getindex(x, 1, 1))
Float64
```
"""
function Base.getindex(x::SparseHDF5Matrix{Tv,Ti}, i::Integer, j::Integer) where {Tv,Ti}
    if !haskey(x.cache, j)
        handle = HDF5.h5open(x.file)
        iset = handle[x.name * "/indices"]
        dset = handle[x.name * "/data"]

        start = Int(x.ptrs[j]) + 1 # 1-based indexing
        finish = Int(x.ptrs[j + 1]) # technically needs a -1, but this is cancelled by the +1 above

        dvals = dset[start:finish]
        ivals = iset[start:finish]
        for i in 1:length(ivals)
            ivals[i] += 1 # 1-based indexing.
        end
        vec = SparseArrays.sparsevec(ivals, dvals, x.dims[1])
        x.cache[j] = vec
    end

    return x.cache[j][i]
end

# TODO: use blockdim to determine whether we can extract multiple columns at
# once to reduce the number of calls.
function general_sparse_extractor(x::SparseHDF5Matrix{Tv,Ti}, indices, store::Function) where {Tv,Ti}
    handle = HDF5.h5open(x.file)
    iset = handle[x.name * "/indices"]
    dset = handle[x.name * "/data"]

    # Indices should be sorted and unique on input.
    allowed_rows = Set(indices[1])

    for ci in 1:length(indices[2])
        c = indices[2][ci]
        start = Int(x.ptrs[c]) + 1  # 1-based indexing
        finish = Int(x.ptrs[c + 1]) # technically needs a -1, but this is cancelled by the +1 above

        colvalues = dset[start:finish]
        colindices = iset[start:finish]
        for i in 1:length(colvalues)
            r = colindices[i] + 1 # 1-based indexing
            if r in allowed_rows
                store(r, c, colvalues[i])
            end
        end
    end
end

"""
    extractsparse(x::SparseHDF5Matrix{Tv,Ti}, i, j; blockdim = nothing)

Extract an in-memory sparse matrix from a `SparseHDF5Matrix` `x`.
The returned matrix contains the same values as `x[i, j]`.

`blockdim` is currently ignored and is only provided for consistency with the method for `DenseHDF5Array`s.

# Examples
```jldoctest
julia> using HDF5Arrays

julia> tmp = tempname();

julia> examplesparse(tmp, "stuff", (20, 10), 0.2)

julia> x = SparseHDF5Matrix(tmp, "stuff");

julia> y = extractsparse(x, 1:5, [2,3,6,7]);

julia> size(y)
(5, 4)
```
"""
function extractsparse(x::SparseHDF5Matrix{Tv,Ti}, i, j; blockdim = nothing) where {Tv,Ti}
    indices = to_indices(x, (i, j))
    NR = length(indices[1])
    NC = length(indices[2])

    collated = create_back_mapping(indices)
    my_indices = sorted_unique_indices(collated)

    # Use of Int is a hack, we can't store the actual index here.
    store = Vector{Dict{Int,Tv}}(undef, x.dims[2])
    for d in my_indices[2]
        store[d] = Dict{Int,Tv}()
    end

    FUN = function (r, c, val)
        store[c][r] = val
    end

    general_sparse_extractor(x, (my_indices...,), FUN)

    return create_csc_matrix(NR, NC, collated[1], collated[2], store)
end

"""
    extractdense(x::SparseHDF5Matrix{Tv,Ti}, i, j; blockdim = nothing)

Extract an in-memory dense `Matrix` from a `SparseHDF5Matrix`.
The returned matrix contains the same values as `x[i, j]`.

`blockdim` is currently ignored and is only provided for consistency with the method for `DenseHDF5Array`s.

# Examples
```jldoctest
julia> using HDF5Arrays

julia> tmp = tempname();

julia> examplesparse(tmp, "stuff", (20, 10), 0.2)

julia> x = SparseHDF5Matrix(tmp, "stuff");

julia> y = extractdense(x, 1:5, [2,3,6,7]);

julia> size(y)
(5, 4)
```
"""
function extractdense(x::SparseHDF5Matrix{Tv,Ti}, i, j; blockdim = nothing) where {Tv,Ti}
    indices = to_indices(x, (i, j))
    NR = length(indices[1])
    NC = length(indices[2])

    collated = create_back_mapping(indices)
    my_indices = sorted_unique_indices(collated)

    # Use of Int is a hack, we can't store the actual index here.
    output = Array{Tv}(undef, (NR, NC))
    fill!(output, 0)

    FUN = function (r, c, val)
        output[collated[1][r], collated[2][c]] .= val
    end

    general_sparse_extractor(x, (my_indices...,), FUN)

    return output
end

ext: jl


text: include("dbtest.jl")


text:
#=##############################################################################
# DESCRIPTION
Test rotor model simulating an isolated 10in propeller in forward flight. The
rotor and configuration approximates the APC Thin Electric 10x7 propeller as
used in McCrink, M. H., & Gregory, J. W. (2017), *Blade Element Momentum
Modeling of Low-Reynolds Electric Propulsion Systems*.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Dec 2019
  * License   : MIT
=###############################################################################

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

function singleprop(; xfoil=true,
                        # OUTPUT OPTIONS
                        save_path=nothing,
                        run_name="singlerotor",
                        prompt=true,
                        verbose=true, v_lvl=0)

    # TODO: Wake removal ?

    # ------------ PARAMETERS --------------------------------------------------

    # Rotor geometry
    rotor_file = "apc10x7.csv"          # Rotor geometry
    data_path = uns.def_data_path       # Path to rotor database
    pitch = 0.0                         # (deg) collective pitch of blades
    n = 10                              # Number of blade elements
    CW = false                          # Clock-wise rotation
    # xfoil = false                     # Whether to run XFOIL

    # Read radius of this rotor and number of blades
    R, B = uns.read_rotor(rotor_file; data_path=data_path)[[1,3]]

    # Simulation parameters
    J = 0.6                             # Advance ratio Vinf/(nD)
    ReD07 = 1.5e6                       # Diameter-based Reynolds at 70% span
    ReD = ReD07/0.7                     # Diameter-based Reynolds
    rho = 1.225                         # (kg/m^3) air density
    mu = 1.81e-5                        # (kg/ms) air dynamic viscosity
    nu = mu/rho
    sound_spd = 343                     # (m/s) speed of sound
    RPM = uns.calc_RPM(ReD, J, R, 2*R, nu) # RPM
    magVinf = J*RPM/60*2*R              # (m/s) freestream velocity
    Minf = magVinf / sound_spd
    Mtip = 2*pi*RPM/60*R / sound_spd

    if verbose
        println("\t"^(v_lvl+1)*"OPERATION PARAMETERS")
        println("\t"^(v_lvl+2)*"J:\t\t$(J)")
        println("\t"^(v_lvl+2)*"ReD07:\t\t$(ReD07)")
        println("\t"^(v_lvl+2)*"RPM:\t\t$(ceil(Int, RPM))")
        println("\t"^(v_lvl+2)*"Mtip:\t\t$(round(Mtip, digits=3))")
        println("\t"^(v_lvl+2)*"Minf:\t\t$(round(Minf, digits=3))")
    end

    Vinf(X,t) = magVinf*[1.0,0,0]       # (m/s) freestream velocity

    # Solver parameters
    nrevs = 8                           # Number of revolutions in simulation
    # nsteps_per_rev = 72               # Time steps per revolution
    nsteps_per_rev = 36
    # p_per_step = 2                    # Sheds per time step
    p_per_step = 1
    ttot = nrevs/(RPM/60)               # (s) total simulation time
    nsteps = nrevs*nsteps_per_rev       # Number of time steps
    lambda = 2.125                      # Core overlap
    overwrite_sigma = lambda * 2*pi*R/(nsteps_per_rev*p_per_step) # Smoothing core size
    surf_sigma = R/10                   # Smoothing radius of lifting surface
    # vlm_sigma = surf_sigma            # Smoothing radius of VLM
    vlm_sigma = -1
    shed_unsteady = true                # Shed particles from unsteady loading

    max_particles = ((2*n+1)*B)*nrevs*nsteps_per_rev*p_per_step # Max particles for memory pre-allocation
    plot_disc = true                    # Plot blade discretization for debugging

    # ------------ SIMULATION SETUP --------------------------------------------
    # Generate rotor
    rotor = uns.generate_rotor(rotor_file; pitch=pitch,
                                            n=n, CW=CW, ReD=ReD,
                                            verbose=verbose, xfoil=xfoil,
                                            data_path=data_path,
                                            plot_disc=plot_disc)
    # ----- VEHICLE DEFINITION
    # System of all FLOWVLM objects
    system = vlm.WingSystem()
    vlm.addwing(system, run_name, rotor)

    # Systems of rotors
    rotors = vlm.Rotor[rotor]   # Defining this rotor as its own system
    rotor_systems = (rotors,)

    # Wake-shedding system
    wake_system = vlm.WingSystem()
    vlm.addwing(wake_system, run_name, rotor)

    # FVS's Vehicle object
    vehicle = uns.VLMVehicle(   system;
                                rotor_systems=rotor_systems,
                                wake_system=wake_system
                             )

    # ----- MANEUVER DEFINITION
    RPM_fun(t) = 1.0                # RPM (normalized by reference RPM) as a
                                    # function of normalized time

    angle = ()                      # Angle of each tilting system (none in this case)
    sysRPM = (RPM_fun, )              # RPM of each rotor system
    Vvehicle(t) = zeros(3)          # Translational velocity of vehicle over Vcruise
    anglevehicle(t) = zeros(3)      # (deg) angle of the vehicle

    # FVS's Maneuver object
    maneuver = uns.KinematicManeuver(angle, sysRPM, Vvehicle, anglevehicle)

    # Plot maneuver path and controls
    uns.plot_maneuver(maneuver; vis_nsteps=nsteps)

    # ----- SIMULATION DEFINITION
    RPMref = RPM
    Vref = 0.0
    simulation = uns.Simulation(vehicle, maneuver, Vref, RPMref, ttot)

    monitor = generate_monitor_prop(J, rho, RPM, nsteps; save_path=save_path,
                                                            run_name=run_name)

    # ------------ RUN SIMULATION ----------------------------------------------
    pfield = uns.run_simulation(simulation, nsteps;
                                      # SIMULATION OPTIONS
                                      Vinf=Vinf,
                                      # SOLVERS OPTIONS
                                      p_per_step=p_per_step,
                                      overwrite_sigma=overwrite_sigma,
                                      vlm_sigma=vlm_sigma,
                                      surf_sigma=surf_sigma,
                                      max_particles=max_particles,
                                      shed_unsteady=shed_unsteady,
                                      extra_runtime_function=monitor,
                                      # OUTPUT OPTIONS
                                      save_path=save_path,
                                      run_name=run_name,
                                      prompt=prompt,
                                      verbose=verbose, v_lvl=v_lvl,
                                      save_code=splitdir(@__FILE__)[1],
                                      )
    return pfield, rotor
end

"""
Generate monitor for rotor performance parameters
"""
function generate_monitor_prop(J, rho, RPM, nsteps; save_path=nothing,
                                run_name="singlerotor",
                                figname="monitor_rotor", disp_conv=true,
                                nsteps_savefig=10)

    fcalls = 0                  # Number of function calls

    colors="rgbcmy"^100
    stls = "o^*.px"^100

    # Name of convergence file
    if save_path!=nothing
        fname = joinpath(save_path, run_name*"_convergence.csv")
    end

    # Function for run_vpm! to call on each iteration
    function extra_runtime_function(sim::uns.Simulation{V, M, R},
                                    PFIELD::uns.vpm.ParticleField,
                                    T::Real, DT::Real
                                   ) where{V<:uns.AbstractVLMVehicle, M, R}

        rotors = vcat(sim.vehicle.rotor_systems...)
        angle = T*360*RPM/60

        # Call figure
        if disp_conv; fig = figure(figname, figsize=(7*3,5*2)); end;

        if fcalls==0
            # Format subplots
            if disp_conv
                subplot(231)
                title("Circulation Distribution")
                xlabel("Element index")
                ylabel(L"Circulation $\Gamma$ (m$^2$/s)")
                grid(true, color="0.8", linestyle="--")
                subplot(232)
                title("Plane-of-rotation Normal Force")
                xlabel("Element index")
                ylabel(L"Normal Force $N_p$ (N)")
                grid(true, color="0.8", linestyle="--")
                subplot(233)
                title("Plane-of-rotation Tangential Force")
                xlabel("Element index")
                ylabel(L"Tangential Force $T_p$ (N)")
                grid(true, color="0.8", linestyle="--")
                subplot(234)
                xlabel(L"Age $\psi$ ($^\circ$)")
                ylabel(L"Thrust $C_T$")
                grid(true, color="0.8", linestyle="--")
                subplot(235)
                xlabel(L"Age $\psi$ ($^\circ$)")
                ylabel(L"Torque $C_Q$")
                grid(true, color="0.8", linestyle="--")
                subplot(236)
                xlabel(L"Age $\psi$ ($^\circ$)")
                ylabel(L"Propulsive efficiency $\eta$")
                grid(true, color="0.8", linestyle="--")
            end

            # Convergence file header
            if save_path!=nothing
                f = open(fname, "w")
                print(f, "age (deg),T,DT")
                for (i, rotor) in enumerate(rotors)
                    print(f, ",RPM_$i,CT_$i,CQ_$i,eta_$i")
                end
                print(f, "\n")
                close(f)
            end
        end

        # Write rotor position and time on convergence file
        if save_path!=nothing
            f = open(fname, "a")
            print(f, angle, ",", T, ",", DT)
        end

        # Plot circulation and loads distributions
        if disp_conv

            cratio = PFIELD.nt/nsteps
            cratio = cratio > 1 ? 1 : cratio
            clr = fcalls==0 && false ? (0,0,0) : (1-cratio, 0, cratio)
            stl = fcalls==0 && false ? "o" : "-"
            alpha = fcalls==0 && false ? 1 : 0.5

            # Circulation distribution
            subplot(231)
            this_sol = []
            for rotor in rotors
                this_sol = vcat(this_sol, [vlm.get_blade(rotor, j).sol["Gamma"] for j in 1:rotor.B]...)
            end
            plot(1:size(this_sol,1), this_sol, stl, alpha=alpha, color=clr)

            # Np distribution
            subplot(232)
            this_sol = []
            for rotor in rotors
                this_sol = vcat(this_sol, rotor.sol["Np"]["field_data"]...)
            end
            plot(1:size(this_sol,1), this_sol, stl, alpha=alpha, color=clr)

            # Tp distribution
            subplot(233)
            this_sol = []
            for rotor in rotors
                this_sol = vcat(this_sol, rotor.sol["Tp"]["field_data"]...)
            end
            plot(1:size(this_sol,1), this_sol, stl, alpha=alpha, color=clr)
        end

        # Plot performance parameters
        for (i,rotor) in enumerate(rotors)
            CT, CQ = vlm.calc_thrust_torque_coeffs(rotor, rho)
            eta = J*CT/(2*pi*CQ)

            if disp_conv
                subplot(234)
                plot([angle], [CT], "$(stls[i])", alpha=alpha, color=clr)
                subplot(235)
                plot([angle], [CQ], "$(stls[i])", alpha=alpha, color=clr)
                subplot(236)
                plot([angle], [eta], "$(stls[i])", alpha=alpha, color=clr)
            end

            if save_path!=nothing
                print(f, ",", rotor.RPM, ",", CT, ",", CQ, ",", eta)
            end
        end

        if disp_conv
            # Save figure
            if fcalls%nsteps_savefig==0 && fcalls!=0 && save_path!=nothing
                savefig(joinpath(save_path, run_name*"_convergence.png"),
                                                            transparent=false)
            end
        end

        # Close convergence file
        if save_path!=nothing
            print(f, "\n")
            close(f)
        end

        fcalls += 1

        return false
    end

    return extra_runtime_function
end


text:
# load packages
using HDF5
using PyPlot
using DSP
using Interpolations
using DelimitedFiles
using FFTW

#H1_16_32 = h5open("/dataset/H-H1_LOSC_16_V1-1126259446-32.hdf5", "r");
H1_4_32 = h5open("dataset/H-H1_LOSC_4_V1-1126259446-32.hdf5", "r");
#L1_16_32 = h5open("/dataset/L-L1_LOSC_16_V1-1126259446-32.hdf5", "r");
L1_4_32 = h5open("dataset/L-L1_LOSC_4_V1-1126259446-32.hdf5", "r");
GW=readdlm("dataset/GW150914_4_NR_waveform.txt");

#s1=read(H1_16_32["strain"])["Strain"]
s2=read(H1_4_32["strain"])["Strain"]
#s3=read(L1_16_32["strain"])["Strain"]
s4=read(L1_4_32["strain"])["Strain"]
m2=read(H1_4_32["meta"])["Duration"]
t2=read(H1_4_32["meta"])["GPSstart"]

#close(H1_16_32)
close(H1_4_32);
#close(L1_16_32)
close(L1_4_32);

#Plot the detector strains and the template for comparison
nt = length(s2)
interval = m2/nt
#xlabel("Time")
#ylabel("Strain")
#title("Raw")
#x=collect(t2:interval:t2+m2-interval)
#plot(s2)
#plot(s4)
#plot(GW[:,1],GW[:,2])

#Removing the noise

nft = round(Int,1/interval)
h_psd,freq1 = psd(s2,nft,nft)         #psd1
l_psd,freq2 = psd(s4,nft,nft)         #psd2
temp_psd,freq3 = psd(GW[:,2],nft,nft) #temp_psd
#loglog(freq1,sqrt(h_psd))  #asd1
#loglog(freq2,sqrt(l_psd))  #asd2
#loglog(freq3,sqrt(temp_psd))  #temp_asd

h_itp = interpolate((freq1,), h_psd, Gridded(Linear()))
l_itp = interpolate((freq2,), l_psd, Gridded(Linear()))

#Whitening
function whiten(strain,psd_,i)
  strain_len = length(strain)
  freq = rfftfreq(strain_len,i)
  hf = rfft(strain)
  white_ft = Complex{Float64}[]
  asd = sqrt.(psd_[freq] /i/2.)
  for i=1:length(hf)
    push!(white_ft,hf[i]/asd[i])
  end
  white_strain = irfft(white_ft,strain_len)
  return white_strain
end
hstrain_whiten = whiten(s2, h_itp, interval)
lstrain_whiten = whiten(s4, l_itp, interval)
temp_whiten = whiten(GW[:,2], h_itp, interval)

#Plot Whitened data
xlabel("Time")
ylabel("Strain")
title("Whitened")
responsetype = Bandpass(20/(nft/2.0), 300/(nft/2.0))
designmethod = Butterworth(4)
h_red = filt(digitalfilter(responsetype, designmethod), hstrain_whiten)
#plot(h_red)
l_red = filt(digitalfilter(responsetype, designmethod), lstrain_whiten)
#plot(l_red)
temp_red = filt(digitalfilter(responsetype, designmethod), temp_whiten)
#plot(temp_red)

#Finding the signal

#Plot the correlation between each detection strain and the reference template
hcorr = xcorr(h_red,temp_red)
lcorr = xcorr(l_red,temp_red)
plot(hcorr)
plot(lcorr)

#Plot the whitened H1 strain-template strain correlation between 15 and 17 seconds
startind = round(Int, 127000*(131072/262143))
endind = round(Int, 135000*(131072/262143))
hcorr = xcorr(h_red[startind:endind], temp_red)
plot(hcorr)
lcorr = xcorr(l_red[startind:endind], temp_red)
plot(lcorr)

#Zero in on the event
startind = 1400
endind = 6600
plot(hcorr[startind:endind])
plot(lcorr[startind:endind])

#Using the time of the event provided by LIGO, see if you found it
xlabel("Time")
ylabel("Strain")
title("Whitened")
plot(h_red[startind:endind])
plot(l_red[startind:endind])
plot(temp_red)

max_stars_repo_name: bingining/noisygravsignals

max_issues_repo_name: bingining/noisygravsignals

max_forks_repo_name: bingining/noisygravsignals


text:
#####
##### Move Action
#####
struct MoveAction
    consumers::Vector{XNode}
    location::TensorLocation
    replace_incumbent::Bool
    # Additional data needed for "asynchronous" moves
    concurrent::Union{Nothing, XNode}
end
isasync(M::MoveAction) = !isnothing(M.concurrent)

# Define this at the XNode level because we have to correctly maintain our "shadow" copy
# of the nGraph graph
function insert_move_node!(
        this_producer::XNode,
        producer_output::Integer,
        consumers::Vector{XNode},
        consumer_inputs::Vector{T},
        action,
    ) where {T <: Integer}

    # First, mutate the nGraph graph underneath all the Julia types
    #
    # Dispatch to the correct implementation depending on whether or not this is an
    # asynchronous move of just a normal synchronous move.
    if !isasync(action)
        move_node = Utils.insert_move_node!(
            unx(this_producer),
            producer_output,
            unx.(consumers),
            consumer_inputs,
        )
    else
        move_node = Utils.insert_move_node!(
            unx(this_producer),
            producer_output,
            unx.(consumers),
            consumer_inputs,
            unx(action.concurrent),
        )

        # If this is not a move node, make sure that the concurrent node is scheduled after
        # the producer node
        if !ismove(unx(this_producer))
            @assert this_producer.index < action.concurrent.index
        end
    end

    # Now, create an appropriate xnode
    #
    # Set its index to 0 temporarily
    xnode = XNode(move_node, 0)

    # Remove all consumer tensors from the producer's output tensor
    producer_xtensor = outputs(this_producer)[producer_output]
    push!(xnode.inputs, producer_xtensor)
    deleteat!(producer_xtensor.users, findall(in(consumers), producer_xtensor.users))

    # Create an output xtensor for the newly created node
    xtensor = XTensor(first(outputs(move_node)), nGraph.CPU)
    push!(xtensor.users, xnode)
    push!(xnode.outputs, xtensor)

    for (input, consumer) in zip(consumer_inputs, consumers)
        consumer.inputs[input] = xtensor
    end

    return xnode
end

#####
##### inner configure!
#####

function configure!(fn::nGraph.NFunction, schedule, data::Profiler.FunctionData)
    # Get the locations of the tensors currently in the graph
    #
    # NOTE: Leave this as TensorDescriptor since it will ultimately be used for ngraphj
    # configuration.
    config = Dict{TensorDescriptor, TensorLocation}()

    # We find all nodes that are targets of a move to DRAM and insert a synchronization
    # barrier.
    #
    # We only do this once for each target because a synchronization barrier will
    # synchronize ALL asynchronous moves
    #
    # KEY: The node that has at least one input coming from an async move
    # VALUE: An async move to this node.
    synced_nodes = Set{XNode}()

    for (tensor, (initial_location, actions)) in schedule
        config[unx(tensor)] = initial_location

        this_producer = producer(tensor)
        producer_output = findonly(isequal(tensor), outputs(this_producer))
        incumbent = tensor

        for action in actions
            consumers = action.consumers
            consumer_inputs = [findonly(isequal(incumbent), inputs(n)) for n in consumers]

            move_node = insert_move_node!(
                this_producer,
                producer_output,
                consumers,
                consumer_inputs,
                action
            )

            # Quick debug
            if action.location == PMEM && !isasync(action)
                @assert initial_location == DRAM
            end

            # Add this move node to `node_dict` and assign its output tensor to the config.
            output_tensor = first(outputs(move_node))
            config[unx(output_tensor)] = action.location

            if action.replace_incumbent
                this_producer = move_node
                # Since we're just inserting move nodes, the output index will now always
                # be 1
                producer_output = 1
                incumbent = output_tensor
            end
        end
    end

    #####
    ##### Apply the config
    #####

    # Make a map from tensor descriptors to XTensors.
    # This will be returned to notify upstream processes which inputs should be in PMM
    arg_descriptors = Set(t.tensor for t in tensors(data) if isarg(t))
    remote_args = Set{TensorDescriptor}()

    # Iterate over each node and each output tensor for each node. Each output tensor should
    # have an assigned location
    for node in fn, output in outputs(NodeDescriptor(node))
        if config[output] == PMEM
            make_persistent(output)

            # Add the xtensor to the list of persistent args if applicable
            if in(output, arg_descriptors)
                push!(remote_args, output)
            end
        end
    end

    # Run priority pass after configuration
    priority_pass!(fn)

    # Set algorithms and workspaces
    return remote_args
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
"""
Internal data for the worker process of the parallel evaluator.
"""
struct ParallelEvaluatorWorker{P<:OptimizationProblem}
    problem::P

    ParallelEvaluatorWorker(problem::P) where {P<:OptimizationProblem} =
        new{P}(problem)
end

fitness(params::Individual, worker::ParallelEvaluatorWorker) =
    fitness(params, worker.problem)

const ChannelRef{T} = RemoteChannel{Channel{T}}
const ParallelEvaluatorWorkerRef{P} = ChannelRef{ParallelEvaluatorWorker{P}}

fitness(params::Individual, worker_ref::ParallelEvaluatorWorkerRef{P}) where P =
    fitness(params, fetch(fetch(worker_ref))::ParallelEvaluatorWorker{P})

"""
Current state of fitness function evaluation for the vector of candidates.
"""
mutable struct ParallelEvaluationState{F, FS}
    fitness_scheme::FS
    candidates::Vector{Candidate{F}}  # candidates to calculate fitness for

    worker_busy::Vector{Bool}         # if the worker is busy calculating
    retry_queue::Vector{Int}          # queue to candidate indices to retry

    worker_finished::Condition        # gets notified each time worker is done
                                      # fitness calculation
    next_index::Int                   # index of the next candidate to calculate fitness

    function ParallelEvaluationState(fitness_scheme::FS,
                                     nworkers::Int) where {FS<:FitnessScheme}
        F = fitness_type(fitness_scheme)
        new{F,FS}(fitness_scheme, Vector{Candidate{F}}(),
                  fill(false, nworkers), Vector{Int}(), Condition(), 0)
    end
end

is_stopped(estate::ParallelEvaluationState) = estate.next_index == 0

abort!(estate::ParallelEvaluationState) = (estate.next_index = 0)

"""
Reset the current `ParallelEvaluationState` and the vector
of candidates that need fitness evaluation.
"""
function reset!(estate::ParallelEvaluationState{F}, candidates::Vector{Candidate{F}}) where F
    estate.candidates = candidates
    fill!(estate.worker_busy, false)
    empty!(estate.retry_queue)
    estate.next_index = 1
    return estate
end

"""
Get the index of the next candidate for evaluation
based on the Base.pmap() code.
"""
function next_candidate!(estate::ParallelEvaluationState, worker_ix::Int)
    @assert !estate.worker_busy[worker_ix]

    task_ix = 0
    if !is_stopped(estate) && estate.next_index <= length(estate.candidates)
        # find the next candidate with unevaluated fitness
        while !is_stopped(estate) && estate.next_index <= length(estate.candidates)
            if isnafitness(estate.candidates[estate.next_index].fitness, estate.fitness_scheme)
                task_ix = estate.next_index
                estate.next_index += 1
                break;
            else
                estate.next_index += 1
            end
        end
    end
    if task_ix == 0
        if !isempty(estate.retry_queue)
            task_ix = popfirst!(estate.retry_queue)
        else
            # Handles the condition where we have finished processing the requested
            # lsts as well as any retryqueue entries, but there are still some jobs
            # active that may result in an error and have to be retried.
            while any(estate.worker_busy)
                wait(estate.worker_finished)
                if !isempty(estate.retry_queue)
                    task_ix = popfirst!(estate.retry_queue)
                    break
                end
            end
        end
    end
    estate.worker_busy[worker_ix] = task_ix != 0
    return task_ix
end

"""
Notify that the worker process is finished and reset its busy flag.
"""
function worker_finished!(estate::ParallelEvaluationState, worker_ix::Int)
    @assert estate.worker_busy[worker_ix]
    estate.worker_busy[worker_ix] = false
    notify(estate.worker_finished; all=true)
    return estate
end

"""
Fitness evaluator that distributes candidates fitness calculation
among several worker processes.
"""
mutable struct ParallelEvaluator{F, FS, P<:OptimizationProblem} <: Evaluator{P}
    problem::P
    archive::Archive
    num_evals::Int
    last_fitness::F

    worker_refs::Vector{ParallelEvaluatorWorkerRef{P}}
    eval_state::ParallelEvaluationState{F, FS}
end

function ParallelEvaluator(
        problem::P, archive::Archive;
        pids::Vector{Int} = workers()) where {P<:OptimizationProblem}
    fs = fitness_scheme(problem)
    ParallelEvaluator{fitness_type(fs), typeof(fs), P}(problem,
        archive,
        0, nafitness(fs),
        [RemoteChannel(function ()
                     # create fake channel and put problem there
                     ch = Channel{ParallelEvaluatorWorker{P}}(1)
                     put!(ch, ParallelEvaluatorWorker(copy(problem)))
                     ch
                   end, pid) for pid in pids],
        ParallelEvaluationState(fs, length(pids))
    )
end

ParallelEvaluator(problem::OptimizationProblem;
        pids::Vector{Int} = workers(),
        archiveCapacity::Integer = 10) =
    ParallelEvaluator(problem, TopListArchive(fitness_scheme(problem), numdims(problem), archiveCapacity),
                      pids=pids)

num_evals(e::ParallelEvaluator) = e.num_evals

function update_fitness!(e::ParallelEvaluator{F}, candidates::Vector{Candidate{F}}) where F
    reset!(e.eval_state, candidates)

    # based on pmap() code
    @sync begin
        for (widx, wref) in enumerate(e.worker_refs)
            @async begin
                while (candi_ix = next_candidate!(e.eval_state, widx)) != 0
                    try
                        candi = candidates[candi_ix]
                        candi.fitness = remotecall_fetch(fitness, wref.where, candi.params, wref)
                        e.last_fitness = candi.fitness
                        e.num_evals += 1
                        add_candidate!(e.archive, candi.fitness, candi.params, e.num_evals)
                    catch ex
                        abort!(e.eval_state) # when one worker fails, the whole process is aborted
                        rethrow(ex)
                    finally
                        worker_finished!(e.eval_state, widx)
                    end
                end
            end
        end
    end

    return candidates
end

# FIXME it's not efficient to calculate fitness like that with `ParallelEvaluator`
function fitness(params::Individual, e::ParallelEvaluator{F}) where F
    candi = Candidate{F}(params)
    update_fitness!(e, [candi])
    return candi.fitness
end


text:
@testset "Bimodal Cauchy Distributions" begin

    @testset "sample bimodal cauchy with truncation above 1" begin
        bc = BlackBoxOptim.BimodalCauchy(0.65, 0.1, 1.0, 0.1, clampBelow0=false)

        n_gt0 = 0
        n_le1 = 0
        for _ in 1:10000
            v = rand(bc)
            (v > 0.0) && (n_gt0 += 1) # Very unlikely to be 0.0 so should be ok...
            (v <= 1.0) && (n_le1 += 1)
        end
        @test n_gt0 == 10000
        @test n_le1 == 10000
    end

    @testset "sample bimodal cauchy with truncation below 0 and above 1" begin
        bc = BlackBoxOptim.BimodalCauchy(0.1, 0.1, 0.95, 0.1)

        n_ge0 = 0
        n_le1 = 0
        for _ in 1:10000
            v = rand(bc)
            (v >= 0.0) && (n_ge0 += 1)
            (v <= 1.0) && (n_le1 += 1)
        end
        @test n_ge0 == 10000
        @test n_le1 == 10000
    end

end


text:
# 3D field tracing on a regular grid.

"""
    trilin_reg(x, y, z, Q)

Trilinear interpolation for x1,y1,z1=(0,0,0) and x2,y2,z2=(1,1,1)
Q's are surrounding points such that Q000 = F[0,0,0], Q100 = F[1,0,0], etc.
"""
function trilin_reg(x, y, z, Q)
   fout =
      Q[1]*(1.0-x)*(1.0-y)*(1.0-z) +
      Q[2]* x *    (1.0-y)*(1.0-z) +
      Q[3]* y *    (1.0-x)*(1.0-z) +
      Q[4]* x * y * (1.0-z) +
      Q[5]*(1.0-x)*(1.0-y)*z +
      Q[6]* x *    (1.0-y)*z +
      Q[7]* y *    (1.0-x)*z +
      Q[8]* x * y * z
end

# Extension from 2d case
function DoBreak(iloc, jloc, kloc, iSize, jSize, kSize)
   ibreak = false
   if iloc ≥ iSize-1 || jloc ≥ jSize-1 || kloc ≥ kSize-1; ibreak = true end
   if iloc < 0 || jloc < 0 || kloc < 0; ibreak = true end
   ibreak
end

"Create unit vectors of field in normalized coordinates."
function normalize_field(ux, uy, uz, dx, dy, dz)
   fx, fy, fz = similar(ux), similar(uy), similar(uz)
   dxInv, dyInv, dzInv = 1/dx, 1/dy, 1/dz
   @inbounds for i = 1:length(ux)
      uInv = 1.0 / hypot(ux[i]*dxInv, uy[i]*dyInv, uz[i]*dzInv)
      fx[i] = ux[i] * dxInv * uInv
      fy[i] = uy[i] * dyInv * uInv
      fz[i] = uz[i] * dzInv * uInv
   end
   fx, fy, fz
end

"""
    grid_interp(x, y, z, field, ix, iy, iz, xsize, ysize)

Interpolate a value at (x,y,z) in a field. `ix`,`iy` and `iz` are indexes for x, y and z
locations (0-based). `xsize` and `ysize` are the sizes of field in X and Y.
"""
grid_interp(x, y, z, field, ix, iy, iz) =
   trilin_reg(x-ix, y-iy, z-iz,
   [
   field[ix+1, iy+1, iz+1],
   field[ix+2, iy+1, iz+1],
   field[ix+1, iy+2, iz+1],
   field[ix+2, iy+2, iz+1],
   field[ix+1, iy+1, iz+2],
   field[ix+2, iy+1, iz+2],
   field[ix+1, iy+2, iz+2],
   field[ix+2, iy+2, iz+2]
   ])

"""
    euler(maxstep, ds, startx, starty, startz, xGrid, yGrid, zGrid, ux, uy, uz)

Fast 3D tracing using Euler's method. It takes at most `maxstep` with step size `ds` tracing
the vector field given by `ux,uy,uz` starting from  `(startx,starty,startz)` in the
Cartesian grid specified by ranges `xGrid`, `yGrid` and `zGrid`.
Return footprints' coordinates in (`x`,`y`,`z`).
"""
function euler(maxstep, ds, startx, starty, startz, xGrid, yGrid, zGrid, ux, uy, uz)

   @assert size(ux) == size(uy) == size(uz) "field array sizes must be equal!"

   x = Vector{eltype(startx)}(undef, maxstep)
   y = Vector{eltype(starty)}(undef, maxstep)
   z = Vector{eltype(startz)}(undef, maxstep)

   iSize, jSize, kSize = size(xGrid,1), size(yGrid,1), size(zGrid,1)

   # Get starting points in normalized/array coordinates
   dx = xGrid[2] - xGrid[1]
   dy = yGrid[2] - yGrid[1]
   dz = zGrid[2] - zGrid[1]
   x[1] = (startx - xGrid[1]) / dx
   y[1] = (starty - yGrid[1]) / dy
   z[1] = (startz - zGrid[1]) / dz

   # Create unit vectors from full vector field
   f1, f2, f3 = normalize_field(ux, uy, uz, dx, dy, dz)

   nstep = 0
   # Perform tracing using Euler's method
   for n = 1:maxstep-1
      # Find surrounding points
      ix = floor(Int, x[n])
      iy = floor(Int, y[n])
      iz = floor(Int, z[n])

      # Break if we leave the domain
      if DoBreak(ix, iy, iz, iSize, jSize, kSize)
         nstep = n; break
      end

      # Interpolate unit vectors to current location
      fx = grid_interp(x[n], y[n], z[n], f1, ix,iy,iz)
      fy = grid_interp(x[n], y[n], z[n], f2, ix,iy,iz)
      fz = grid_interp(x[n], y[n], z[n], f3, ix,iy,iz)

      if any(isnan,[fx, fy, fz]) || any(isinf, [fx, fy, fz])
         nstep = n
         break
      end

      @muladd x[n+1] = x[n] + ds * fx
      @muladd y[n+1] = y[n] + ds * fy
      @muladd z[n+1] = z[n] + ds * fz

      nstep = n
   end

   # Convert traced points to original coordinate system.
   @inbounds for i = 1:nstep
      @muladd x[i] = x[i]*dx + xGrid[1]
      @muladd y[i] = y[i]*dy + yGrid[1]
      @muladd z[i] = z[i]*dz + zGrid[1]
   end
   x[1:nstep], y[1:nstep], z[1:nstep]
end

"""
    rk4(maxstep, ds, startx, starty, startz, xGrid, yGrid, zGrid, ux, uy, uz)

Fast and reasonably accurate 3D tracing with 4th order Runge-Kutta method and constant step
size `ds`. See also [`euler`](@ref).
"""
function rk4(maxstep, ds, startx, starty, startz, xGrid, yGrid, zGrid,
   ux, uy, uz)

   @assert size(ux) == size(uy) == size(uz) "field array sizes must be equal!"

   x = Vector{eltype(startx)}(undef, maxstep)
   y = Vector{eltype(starty)}(undef, maxstep)
   z = Vector{eltype(startz)}(undef, maxstep)

   iSize, jSize, kSize = size(xGrid,1), size(yGrid,1), size(zGrid,1)

   # Get starting points in normalized/array coordinates
   dx = xGrid[2] - xGrid[1]
   dy = yGrid[2] - yGrid[1]
   dz = zGrid[2] - zGrid[1]
   x[1] = (startx - xGrid[1]) / dx
   y[1] = (starty - yGrid[1]) / dy
   z[1] = (startz - zGrid[1]) / dz

   # Create unit vectors from full vector field
   fx, fy, fz = normalize_field(ux, uy, uz, dx, dy, dz)

   nstep = 0
   # Perform tracing using RK4
   for n = 1:maxstep-1
      # SUBSTEP #1
      ix = floor(Int, x[n])
      iy = floor(Int, y[n])
      iz = floor(Int, z[n])
      if DoBreak(ix, iy, iz, iSize, jSize, kSize); nstep = n; break end

      f1x = grid_interp(x[n],y[n],z[n], fx, ix,iy,iz)
      f1y = grid_interp(x[n],y[n],z[n], fy, ix,iy,iz)
      f1z = grid_interp(x[n],y[n],z[n], fz, ix,iy,iz)
      if any(isnan,[f1x, f1y, f1z]) || any(isinf, [f1x, f1y, f1z])
         nstep = n; break
      end
      # SUBSTEP #2
      xpos = x[n] + f1x*ds/2.0
      ypos = y[n] + f1y*ds/2.0
      zpos = z[n] + f1z*ds/2.0
      ix = floor(Int, xpos)
      iy = floor(Int, ypos)
      iz = floor(Int, zpos)
      if DoBreak(ix, iy, iz, iSize, jSize, kSize); nstep = n; break end

      f2x = grid_interp(xpos,ypos,zpos, fx, ix,iy,iz)
      f2y = grid_interp(xpos,ypos,zpos, fy, ix,iy,iz)
      f2z = grid_interp(xpos,ypos,zpos, fz, ix,iy,iz)
      if any(isnan,[f2x, f2y, f2z]) || any(isinf, [f2x, f2y, f2z])
         nstep = n; break
      end
      # SUBSTEP #3
      xpos = x[n] + f2x*ds/2.0
      ypos = y[n] + f2y*ds/2.0
      zpos = z[n] + f2z*ds/2.0
      ix = floor(Int, xpos)
      iy = floor(Int, ypos)
      iz = floor(Int, zpos)
      if DoBreak(ix, iy, iz, iSize, jSize, kSize); nstep = n; break end

      f3x = grid_interp(xpos,ypos,zpos, fx, ix,iy,iz)
      f3y = grid_interp(xpos,ypos,zpos, fy, ix,iy,iz)
      f3z = grid_interp(xpos,ypos,zpos, fz, ix,iy,iz)
      if any(isnan,[f3x, f3y, f3z]) || any(isinf, [f3x, f3y, f3z])
         nstep = n; break
      end

      # SUBSTEP #4
      xpos = x[n] + f3x*ds
      ypos = y[n] + f3y*ds
      zpos = z[n] + f3z*ds
      ix = floor(Int, xpos)
      iy = floor(Int, ypos)
      iz = floor(Int, zpos)
      if DoBreak(ix, iy, iz, iSize, jSize, kSize); nstep = n; break end

      f4x = grid_interp(xpos,ypos,zpos, fx, ix,iy,iz)
      f4y = grid_interp(xpos,ypos,zpos, fy, ix,iy,iz)
      f4z = grid_interp(xpos,ypos,zpos, fz, ix,iy,iz)
      if any(isnan,[f4x, f4y, f4z]) || any(isinf, [f4x, f4y, f4z])
         nstep = n; break
      end

      # Perform the full step using all substeps
      @muladd x[n+1] = x[n] + ds/6.0 * (f1x + f2x*2.0 + f3x*2.0 + f4x)
      @muladd y[n+1] = y[n] + ds/6.0 * (f1y + f2y*2.0 + f3y*2.0 + f4y)
      @muladd z[n+1] = z[n] + ds/6.0 * (f1z + f2z*2.0 + f3z*2.0 + f4z)

      nstep = n
   end

   # Convert traced points to original coordinate system.
   @inbounds for i = 1:nstep
      @muladd x[i] = x[i]*dx + xGrid[1]
      @muladd y[i] = y[i]*dy + yGrid[1]
      @muladd z[i] = z[i]*dz + zGrid[1]
   end
   x[1:nstep], y[1:nstep], z[1:nstep]
end

"""
	 trace3d_euler(fieldx, fieldy, fieldz, startx, starty, startz, gridx, gridy, gridz;
       maxstep=20000, ds=0.01)

Given a 3D vector field, trace a streamline from a given point to the edge of the vector
field. The field is integrated using Euler's method. Only valid for regular grid with
coordinates `gridx`, `gridy`, `gridz`.
The field can be in both `meshgrid` or `ndgrid` (default) format.
Supporting `direction` of {"both","forward","backward"}.
"""
function trace3d_euler(fieldx, fieldy, fieldz, startx, starty, startz, gridx, gridy, gridz;
   maxstep=20000, ds=0.01, gridtype="ndgrid", direction="both")

   if gridtype == "ndgrid"
      fx, fy, fz = fieldx, fieldy, fieldz
   else # meshgrid
      fx, fy, fz = permutedims(fieldx), permutedims(fieldy), permutedims(fieldz)
   end

   if direction == "forward"
      xt, yt, zt = euler(maxstep, ds, startx,starty,startz, gridx,gridy,gridz, fx, fy, fz)
   elseif direction == "backward"
      xt, yt, zt = euler(maxstep, ds, startx,starty,startz, gridx,gridy,gridz, -fx,-fy,-fz)
   else
      x1, y1, z1 = euler(floor(Int, maxstep/2), ds, startx, starty, startz,
         gridx, gridy, gridz, -fx, -fy, -fz)
      blen = length(x1)
      x2, y2, z2 = euler(maxstep-blen, ds, startx, starty, startz,
         gridx, gridy, gridz, fx, fy, fz)
      # concatenate with duplicates removed
      xt = vcat(reverse!(x1), x2[2:end])
      yt = vcat(reverse!(y1), y2[2:end])
      zt = vcat(reverse!(z1), z2[2:end])
   end

   xt, yt, zt
end

"""
	 trace3d_rk4(fieldx, fieldy, fieldz, startx, starty, startz, gridx, gridy, gridz;
       maxstep=20000, ds=0.01)

Given a 3D vector field, trace a streamline from a given point to the edge of the vector
field. The field is integrated using Euler's method. Only valid for regular grid with
coordinates `gridx`, `gridy`, `gridz`.
The field can be in both `meshgrid` or `ndgrid` (default) format.
See also [`trace3d_euler`](@ref).
"""
function trace3d_rk4(fieldx, fieldy, fieldz, startx, starty, startz, gridx, gridy, gridz;
   maxstep=20000, ds=0.01, gridtype="ndgrid", direction="both")

   if gridtype == "ndgrid"
      fx, fy, fz = fieldx, fieldy, fieldz
   else # meshgrid
      fx, fy, fz = permutedims(fieldx), permutedims(fieldy), permutedims(fieldz)
   end

   if direction == "forward"
      xt, yt, zt = rk4(maxstep, ds, startx,starty,startz, gridx,gridy,gridz, fx, fy, fz)
   elseif direction == "backward"
      xt, yt, zt = rk4(maxstep, ds, startx,starty,startz, gridx,gridy,gridz, -fx,-fy,-fz)
   else
      x1, y1, z1 = rk4(floor(Int,maxstep/2), ds, startx, starty, startz,
         gridx, gridy, gridz, -fx, -fy, -fz)
      blen = length(x1)
      x2, y2, z2 = rk4(maxstep-blen, ds, startx, starty, startz,
         gridx, gridy, gridz, fx, fy, fz)
      # concatenate with duplicates removed
      xt = vcat(reverse!(x1), x2[2:end])
      yt = vcat(reverse!(y1), y2[2:end])
      zt = vcat(reverse!(z1), z2[2:end])
   end

   xt, yt, zt
end

"""
    trace3d_euler(fieldx, fieldy, fieldz, startx, starty, startz, grid::CartesianGrid;
		 maxstep=20000, ds=0.01, gridtype="ndgrid", direction="both")

See also [`trace3d_rk4`](@ref).
"""
function trace3d_euler(fieldx, fieldy, fieldz, startx, starty, startz, grid::CartesianGrid;
   kwargs...)

   gridmin = coordinates(minimum(grid))
   gridmax = coordinates(maximum(grid))
   Δx = spacing(grid)

   gridx = range(gridmin[1], gridmax[1], step=Δx[1])
   gridy = range(gridmin[2], gridmax[2], step=Δx[2])
   gridz = range(gridmin[3], gridmax[3], step=Δx[3])

   trace3d_euler(fieldx, fieldy, fieldz, startx, starty, startz, gridx, gridy, gridz;
      kwargs...)
end


text:
using CoordinateTransformations, Rotations, StaticArrays
using LinearAlgebra
using TransformUtils
using Plots
using PyPlot;
const plt = PyPlot;
pygui(true)

include("lie_groups.jl")

D = [
    1 0 0
    0 1 0
    0 0 1
    0 0 0
]

# Define batch estimator structure
mutable struct BatchEstimator

    # Vectors
    time_steps::Vector{Float32}

    # Rotational velocity vk_vk
    ω_vₖ_vₖ_i::Matrix{Float32}

    # Linear velocity vk_vk
    v_vₖ_vₖ_i::Matrix{Float32}

    # GT rotation (axis-angle)
    θ_vₖ_i::Matrix{Float32}

    # GT position
    rᵢ_vₖ_i::Matrix{Float32}

    # landmark position
    rho_i_pj_i::Matrix{Float32}

    # observations
    yₖ_j::Array{Float32,3}

    # Transformation from vehicle to camera frame
    T_c_v::Matrix{Float32}

    # Estimated states for each timestep
    states::Array{Float32,3}

    # Ground truth
    ground_truth::Vector{Matrix{Float32}}

    # Variances
    v_var::Matrix{Float32}
    w_var::Matrix{Float32}
    y_var::Matrix{Float32}

    # Scalars
    num_landmarks::Int
    fu::Float32
    fv::Float32
    cu::Float32
    cv::Float32
    b::Float32
    k1::Int
    k2::Int
    batch_size::Int
    max_iterations::Int

    BatchEstimator(data::Dict) = create_batch_estimator!(new(), data)
end

# Define batch estimator functions
function create_batch_estimator!(estimator::BatchEstimator, data::Dict)

    # Extract raw data
    estimator.ω_vₖ_vₖ_i = data["w_vk_vk_i"]
    estimator.v_vₖ_vₖ_i = data["v_vk_vk_i"]
    estimator.θ_vₖ_i = data["theta_vk_i"]
    estimator.rᵢ_vₖ_i = data["r_i_vk_i"]
    estimator.rho_i_pj_i = data["rho_i_pj_i"]
    estimator.yₖ_j = data["y_k_j"]

    # Variances
    estimator.v_var = data["v_var"] * 3
    estimator.w_var = data["w_var"]
    estimator.y_var = data["y_var"]

    # Constants
    estimator.num_landmarks = 20
    estimator.fu = data["fu"]
    estimator.fv = data["fv"]
    estimator.cu = data["cu"]
    estimator.cv = data["cv"]
    estimator.b = data["b"]

    # Time steps
    times = data["t"]
    estimator.time_steps = [
        times[1, k] - times[1, k-1] for k in 2:size(times, 2)
    ]

    # Set vehicle to camera transform
    estimator.T_c_v = [data["C_c_v"] -data["C_c_v"]*data["rho_v_c_v"]; 0 0 0 1]

    # Set the period we're looking at
    estimator.k1 = 1218
    estimator.k2 = 1426
    estimator.batch_size = estimator.k2 - estimator.k1

    # Set ground truth
    estimator.ground_truth = ground_truth(estimator)

    # Set gauss newton parameters
    estimator.max_iterations = 3

    return estimator
end

function ground_truth_k(estimator::BatchEstimator, k)
    C = so3_exp(estimator.θ_vₖ_i[:, k])
    d = estimator.rᵢ_vₖ_i[:, k]
    return [C -C*d; 0 0 0 1]
end

function ground_truth(estimator::BatchEstimator)
    return [
        inv(ground_truth_k(estimator, k))
        for k ∈ estimator.k1:estimator.k2
    ]
end

function forward(estimator::BatchEstimator, k::Int64, T::Matrix{<:AbstractFloat})

    # Rotation = Angular velocity * time period
    psi = estimator.ω_vₖ_vₖ_i[:, k] * estimator.time_steps[k]

    # Distance = Linear velocity * time period
    d = estimator.v_vₖ_vₖ_i[:, k] * estimator.time_steps[k]

    # Change in pose
    dT = inv(se3_exp(d, -psi))

    return dT * T
end

function dead_reconing(estimator::BatchEstimator)

    # Get initial state from ground truth
    curr_state = ground_truth_k(estimator, estimator.k1)

    # Set initial state
    states = Array{Float32,3}(undef, estimator.batch_size, 4, 4)
    states[1, :, :] = curr_state

    for k ∈ estimator.k1:estimator.k2-2

        # Compute next state
        curr_state = forward(estimator, k, curr_state)

        # update current state
        states[k-estimator.k1+2, :, :] = curr_state
    end

    return states
end

# Given a point p in the camera frame, project it
# on the images such that pixels are returned
function g(estimator::BatchEstimator, p::Vector{<:AbstractFloat})
    """Given a point and pose project into camera."""
    return [
        estimator.fu * p[1] / p[3] + estimator.cu
        estimator.fv * p[2] / p[3] + estimator.cv
        estimator.fu * (p[1] - estimator.b) / p[3] + estimator.cu
        estimator.fv * p[2] / p[3] + estimator.cv
    ]
end

# Given the state {r,C} of the vehicle, project the
# i-th point expressed wrt the inertial frame
# into the camera frame and return its [x,y,z] coords.
function p_i(estimator::BatchEstimator, T::Matrix{<:AbstractFloat}, i::Int64)
    """Project point into camera frame."""
    return D' * estimator.T_c_v * T * [estimator.rho_i_pj_i[:, i]; 1]
end

# Jacobian of the motion model wrt the perturbation
# If you have the previous pose and the current pose,
# this can be computed as Ad( T_curr inv(T_prev) ).
# The provided poses needs to be transformation matrices.
function motion_jacobian(T_prev::Matrix{<:AbstractFloat}, T_curr::Matrix{<:AbstractFloat})
    return se3_adjoint(T_curr * se3_inv(T_prev))
end

# Return the error between two poses expressed as
# transformation matrices. The error is a 6x1 vector.
function pose_error(T1::Matrix{<:AbstractFloat}, T2::Matrix{<:AbstractFloat})
    error = se3_log(T1 * se3_inv(T2))
    error[error.<1e-4] .= 0.0
    return error
end

# Return the error, in pixels between an observation
# and the predicted pixel position (4x1) of a landmark
function measurement_error(obs_1::Vector{<:AbstractFloat}, obs_2::Vector{<:AbstractFloat})
    return obs_1 - obs_2
end

function measurement_to_point(self, p)
    ul = p[0]
    vl = p[1]
    ur = p[2]
    vr = p[3]
    z = self.fu * self.b / (ul - ur)
    x = (ul - self.cu) * z / self.fu
    y = (vl - self.cv) * z / self.fv
    return np.array([x, y, z, 1])
end

# Jacobian of g() wrt to point p
function dg_dp(estimator::BatchEstimator, p::Vector{<:AbstractFloat})
    """Camera jacobian."""
    x = p[1]
    y = p[2]
    z = p[3]
    z_squared = z^2
    return [
        estimator.fu/z 0 -estimator.fu*x/z_squared
        0 estimator.fv/z -estimator.fv*y/z_squared
        estimator.fu/z 0 -estimator.fu*(x-estimator.b)/z_squared
        0 estimator.fv/z -estimator.fv*y/z_squared
    ]
end

function dp_dx(estimator::BatchEstimator, T::Matrix{<:AbstractFloat}, p::Vector{<:AbstractFloat})
    """
    C: Orientation of the vehicle wrt inertial frame
    r: Position of the vehicle wrt inertial frame, expressed in inertial frame
    p: Position of the observed point wrt inertial frame
    """
    return D' * estimator.T_c_v * se3_odot(T * [p; 1])
end

function gauss_newton(estimator::BatchEstimator, initial_estimate::Array{<:AbstractFloat,3})

    # Initalize guass newton updated states
    T_op = copy(initial_estimate)
    T0 = T_op[1, :, :]

    # Holds the perturbations used on each iteration
    delta_x = zeros(estimator.batch_size, 6)

    # Width of H jacobian matrix
    width = 6 * estimator.batch_size

    # Iterate refinement until convergence
    for i ∈ 1:estimator.max_iterations

        # Initialize lists of arrays
        ev = Vector{Float32}()  # Motion errors throughout trajectory
        ey = Vector{Float32}()          # Measurement errors for observations
        Gs = Vector{Matrix{Float32}}()  # Lower part of the H jacobian matrix

        # Variance P0 of the initial pose
        sigma_init = 1

        # Wv_inv is the inverse of the first part of the W matrix
        # with diagonal entries related to process noise.
        Wv_inv = ones(6, 1) * sigma_init

        # Wy_inv is the inverse of the second part of the W matrix
        # with diagonal entries related to measurement noise.
        Wy_inv = Vector{Matrix{Float32}}()

        # push!(ev, pose_error(T0, T_op[1, :, :]))
        ev = [ev; pose_error(T0, T_op[1, :, :])]

        # Each "row" of Fs consists in a 6 x 6K array
        # The first one has only np.eye(6) at the front of the row
        F = zeros(6, width)
        F[:, 1:6] = I(6)
        Fs = F

        # This loop generates the upper part of the H jacobian matrix
        # and the upper part of the error matrix.
        for k ∈ estimator.k1:estimator.k2-2

            # Relative k index starting at 0 when k = k1
            k_rel = k - estimator.k1 + 1

            # Current pose T_curr and previous pose T_prev
            T_curr = T_op[k_rel+1, :, :]
            T_prev = T_op[k_rel, :, :]

            # Given the previous pose, use the motion model to predict
            # the current pose.
            T_pred = forward(estimator, k, T_prev)

            # Compare predicted pose with current pose
            error = pose_error(T_pred, T_curr)
            # Append the error to the list
            # push!(ev, error)
            ev = [ev; error]

            # Since the Q matrix is diagonal, the inverse of this matrix
            # is obtained by taking the reciprocal of the diagonal elements.
            # The variance of the linear and angular process noise must be multiplied
            # by the squared period of the timestep since the variance is expressed
            # in [m/s] and [rad/s] but we want [m] and [rad].
            qr_inv = 1.0 ./ (estimator.v_var * (estimator.time_steps[k]^2))
            qc_inv = 1.0 ./ (estimator.w_var * (estimator.time_steps[k]^2))
            q_inv = [qr_inv; qc_inv]

            # Append the inverse of the Q matrix to the list
            Wv_inv = [Wv_inv; q_inv]

            # Motion jacobian of the previous timestep F_{k-1}
            F = motion_jacobian(T_prev, T_curr)

            # Most of the elements of the row are zeros
            F_row = zeros(6, width)
            # Horizontal offset where to place the -1*F matrix
            offset = 6 * (k_rel - 1) + 1
            F_row[:, offset:offset+5] = -1 * F
            # Horizontal offset where to place the identity matrix
            offset = offset + 6
            F_row[:, offset:offset+5] = I(6)

            # Append the produced row to the list
            Fs = vcat(Fs, F_row)
        end

        # This loop generates the lower part of the H jacobian matrix
        # and the lower part of the error matrix.
        for k ∈ estimator.k1:estimator.k2-1

            # Produced matrices for the k-th timestep
            Gk = Vector{Matrix{Float32}}()      # Measurement jacobian
            eyk = Vector{Float32}()     # Observation errors
            Rk_inv = Vector{Matrix{Float32}}()  # Inverted measurement covariance matrix

            # Relative k index starting at 0 when k = k1
            k_rel = k - estimator.k1 + 1

            # For convenience
            T = T_op[k_rel, :, :]

            # Iterate through all of the possible landmarks
            for j ∈ 1:estimator.num_landmarks

                # Pixels values (4x1) of the j-th landmark
                # observed at the k-th timestep.
                landmark = estimator.yₖ_j[:, k, j]

                # If the landmark is invalid, each pixel will be set to -1
                if landmark[1] > 0

                    # Project the j-th point in camera frame
                    P_cam = p_i(estimator, T, j)
                    # Project a point in space to a point in images
                    P_image = g(estimator, P_cam)

                    # Compute the measurement error and append to list
                    meas_error = measurement_error(landmark, P_image)
                    eyk = [eyk; meas_error]

                    # Append the inverse of R
                    # push!(Rk_inv, 1 ./ estimator.y_var)
                    if size(Rk_inv, 1) > 0
                        Rk_inv = vcat([Rk_inv; 1 ./ estimator.y_var])
                    else
                        Rk_inv = 1 ./ estimator.y_var
                    end

                    # Jacobian of g() wrt perturbations given pose and landmark
                    G = dg_dp(estimator, P_cam) * dp_dx(estimator, T, estimator.rho_i_pj_i[:, j])
                    push!(Gk, G)
                end
            end

            # Add the optimization variables to the problem
            # if you were able to observe more than a certain
            # number of landmarks.
            if size(Gk, 1) > 0

                ey = [ey; eyk]
                if size(Wy_inv, 1) > 0
                    Wy_inv = [Wy_inv; Rk_inv]
                else
                    Wy_inv = Rk_inv
                end

                # Each G jacobian is a 4x6 matrix
                G_row = zeros(4 * size(Gk, 1), width)

                # Offset each jacobian below and right of the previous one
                for (i, Gki) ∈ enumerate(Gk)
                    vert_offset = 4 * (i - 1) + 1       # Vertical offset
                    horz_offset = 6 * (k_rel - 1) + 1   # Horizontal offset
                    G_row[
                        vert_offset:vert_offset+3, horz_offset:horz_offset+5
                    ] = Gki
                end

                if size(Gs, 1) > 1
                    Gs = vcat(Gs, G_row)
                else
                    Gs = G_row
                end
            end
        end

        # Finally, the H matrix can be produced from its upper and lower parts
        H = [Fs; Gs]
        # Build the inverse covariance matrix
        W_inv = Diagonal([Wv_inv; Wy_inv][:, 1])
        # Build the error matrix
        e = [ev; ey]

        # For Ax=b and solving for x
        A = H' * W_inv * H
        b = H' * W_inv * e

        # We solve for the optimal perturbations \delta_x
        # for the whole trajectory. Each perturbation is a
        # 6x1 vector.
        x_opt = A \ b

        # Reshape the very 6*N long vector in a Nx6 matrix where
        # N is the number of timesteps in the trajectory.

        # TODO: FIGURE THIS OUT
        delta_x = reshape(x_opt, 6, estimator.batch_size)'

        # To update the operating point of the optimization we
        # need to unstack all optimal perturbations and then
        # use T = exp(x^)T to update the operating point T
        # which represents our best guess in terms of pose
        for k ∈ 1:estimator.batch_size
            d = delta_x[k, 1:3]
            psi = delta_x[k, 4:6]
            # d = clamp!(d, -0.2, 0.2)
            # psi = clamp!(psi, -0.2, 0.2)
            T_op[k, :, :] = se3_exp(d, -psi) * T_op[k, :, :]
        end

    end

    return T_op
end

function plot_trajectory(estimator::BatchEstimator, dead_reconing::Array{<:AbstractFloat,3}, gauss_newton::Array{<:AbstractFloat,3})
    num_steps = estimator.batch_size

    # Create plot
    ax = plt.axes(projection="3d")
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')

    states = [inv(dead_reconing[k, :, :]) for k ∈ 1:num_steps]

    # Get state estimates
    x = [states[k][1, 4] for k ∈ 1:num_steps]
    y = [states[k][2, 4] for k ∈ 1:num_steps]
    z = [states[k][3, 4] for k ∈ 1:num_steps]
    ax.plot3D(x, y, z, c="red", label="Dead Reconing")

    states = [inv(gauss_newton[k, :, :]) for k ∈ 1:num_steps]

    # Get state estimates
    x = [states[k][1, 4] for k ∈ 1:num_steps]
    y = [states[k][2, 4] for k ∈ 1:num_steps]
    z = [states[k][3, 4] for k ∈ 1:num_steps]
    ax.plot3D(x, y, z, c="purple", label="Gauss Newton")

    # Get ground truth
    gt_x = [estimator.ground_truth[k][1, 4] for k ∈ 1:num_steps]
    gt_y = [estimator.ground_truth[k][2, 4] for k ∈ 1:num_steps]
    gt_z = [estimator.ground_truth[k][3, 4] for k ∈ 1:num_steps]
    ax.plot3D(gt_x, gt_y, gt_z, c="blue", label="GT")

    # Landmark positions
    l_x = [estimator.rho_i_pj_i[1, i] for i ∈ 1:estimator.num_landmarks]
    l_y = [estimator.rho_i_pj_i[2, i] for i ∈ 1:estimator.num_landmarks]
    l_z = [estimator.rho_i_pj_i[3, i] for i ∈ 1:estimator.num_landmarks]
    ax.scatter(l_x, l_y, l_z, c="r", label="landmarks")

    title("Estimate vs Ground Truth")
    plt.legend()
    plt.show()

end

max_stars_repo_name: jandaa/state-estimation

max_issues_repo_name: jandaa/state-estimation

max_forks_repo_name: jandaa/state-estimation


text:
println("Workout")

for pk in keys(Pkg.installed())
       (pk != "PackageCompiler") && (eval(parse("import $pk")))
end
import CommonTools
using CommonTools: parseini,cached_load
using Base.Test
using PyPlot
import Reactive
# imshow(rand(100,100))
# write your own tests here
# @test 1 == 1
using StaticArrays
using GLVisualize
using Rotations,Parameters
window = GLVisualize.glscreen()

using BenchmarkTools
import Optim

using FileIO
using MeshIO
using StaticArrays
using GeometryTypes
ply_1 = string(@__DIR__ ,"/sawbone.ply")
ply = load(ply_1)
ply

ini_1 = string(@__DIR__ ,"/HMD.ini")
CommonTools.parseini(ini_1)

_view(visualize(ply),window)
# GLVisualize.renderloop(window)
# @async GLWindow.renderloop(window)
# sleep(10)
# addprocs(1)
# @spawn rand(100,100)
println("Workout")

for i=1:10
            GLVisualize.render_frame(window)
            GLVisualize.swapbuffers(window)
            GLVisualize.poll_glfw()
            yield()
end
GLVisualize.destroy!(window)
ply = nothing

println(nworkers())

import LsqFit
println("Workout")

# a two-parameter exponential model
# x: array of independent variables
# p: array of model parameters
model(x, p) = p[1]*exp.(-x.*p[2])

# some example data
# xdata: independent variables
# ydata: dependent variable
xdata = linspace(0,10,20)
ydata = model(xdata, [1.0 2.0]) + 0.01*randn(length(xdata))
p0 = [0.5, 0.5]

fit = LsqFit.curve_fit(model, xdata, ydata, p0)
# fit is a composite type (LsqFitResult), with some interesting values:
#	fit.dof: degrees of freedom
#	fit.param: best fit parameters
#	fit.resid: residuals = vector of residuals
#	fit.jacobian: estimated Jacobian at solution

# We can estimate errors on the fit parameters,
# to get 95% confidence error bars:
errors = LsqFit.estimate_errors(fit, 0.95)

# The finite difference method is used above to approximate the Jacobian.
# Alternatively, a function which calculates it exactly can be supplied instead.
function jacobian_model(x,p)
    J = Array{Float64}(length(x),length(p))
    J[:,1] = exp.(-x.*p[2])    #dmodel/dp[1]
    J[:,2] = -x.*p[1].*J[:,1]  #dmodel/dp[2]
    J
end
fit = LsqFit.curve_fit(model, jacobian_model, xdata, ydata, p0)
println(fit)


text:
# wrappers of low-level functionality

using CUDAapi: cudaDataType
using CUDAdrv: CuDefaultStream, CuStream

function version()
    ver = cutensorGetVersion()
    major, ver = divrem(ver, 10000)
    minor, patch = divrem(ver, 100)

    VersionNumber(major, minor, patch)
end

function cuda_version()
  ver = cutensorGetCudartVersion()
  major, ver = divrem(ver, 1000)
  minor, patch = divrem(ver, 10)

  VersionNumber(major, minor, patch)
end

const ModeType = AbstractVector{<:Union{Char, Integer}}

is_unary(op::cutensorOperator_t) =
    (op ∈ (CUTENSOR_OP_IDENTITY, CUTENSOR_OP_SQRT, CUTENSOR_OP_RELU, CUTENSOR_OP_CONJ,
            CUTENSOR_OP_RCP))
is_binary(op::cutensorOperator_t) =
    (op ∈ (CUTENSOR_OP_ADD, CUTENSOR_OP_MUL, CUTENSOR_OP_MAX, CUTENSOR_OP_MIN))

mutable struct CuTensorDescriptor
    desc::Ref{cutensorTensorDescriptor_t}

    function CuTensorDescriptor(a; size = size(a), strides = strides(a), eltype = eltype(a),
                                   op = CUTENSOR_OP_IDENTITY)
        sz = collect(Int64, size)
        st = collect(Int64, strides)
        desc = Ref{cutensorTensorDescriptor_t}()
        cutensorInitTensorDescriptor(handle(), desc, length(sz), sz, st,
                                     cudaDataType(eltype), op)
        obj = new(desc)
        return obj
    end
end

Base.cconvert(::Type{Ptr{cutensorTensorDescriptor_t}}, obj::CuTensorDescriptor) = obj.desc

function elementwiseTrinary!(
    alpha::Number, A::CuArray, Ainds::ModeType, opA::cutensorOperator_t,
    beta::Number,  B::CuArray, Binds::ModeType, opB::cutensorOperator_t,
    gamma::Number, C::CuArray{T}, Cinds::ModeType, opC::cutensorOperator_t,
    D::CuArray{T}, Dinds::ModeType, opAB::cutensorOperator_t,
    opABC::cutensorOperator_t; stream::CuStream=CuDefaultStream()) where {T}

    !is_unary(opA)    && throw(ArgumentError("opA must be a unary op!"))
    !is_unary(opB)    && throw(ArgumentError("opB must be a unary op!"))
    !is_unary(opC)    && throw(ArgumentError("opC must be a unary op!"))
    !is_binary(opAB)  && throw(ArgumentError("opAB must be a binary op!"))
    !is_binary(opABC) && throw(ArgumentError("opABC must be a binary op!"))
    descA = CuTensorDescriptor(A; op = opA)
    descB = CuTensorDescriptor(B; op = opB)
    descC = CuTensorDescriptor(C; op = opC)
    @assert size(C) == size(D) && strides(C) == strides(D)
    descD = descC # must currently be identical
    #typeCompute = cudaDataType(T)
    typeCompute = cudaDataType(T)
    modeA = collect(Cint, Ainds)
    modeB = collect(Cint, Binds)
    modeC = collect(Cint, Cinds)
    modeD = modeC
    cutensorElementwiseTrinary(handle(),
                                T[alpha], A, descA, modeA,
                                T[beta],  B, descB, modeB,
                                T[gamma], C, descC, modeC,
                                          D, descD, modeD,
                                opAB, opABC, typeCompute, stream)
    return D
end

function elementwiseTrinary!(
    alpha::Number, A::Array, Ainds::ModeType, opA::cutensorOperator_t,
    beta::Number, B::Array, Binds::ModeType, opB::cutensorOperator_t,
    gamma::Number, C::Array{T}, Cinds::ModeType, opC::cutensorOperator_t,
    D::Array{T}, Dinds::ModeType, opAB::cutensorOperator_t,
    opABC::cutensorOperator_t; stream::CuStream=CuDefaultStream()) where {T}

    !is_unary(opA)    && throw(ArgumentError("opA must be a unary op!"))
    !is_unary(opB)    && throw(ArgumentError("opB must be a unary op!"))
    !is_unary(opC)    && throw(ArgumentError("opC must be a unary op!"))
    !is_binary(opAB)  && throw(ArgumentError("opAB must be a binary op!"))
    !is_binary(opABC) && throw(ArgumentError("opABC must be a binary op!"))
    descA = CuTensorDescriptor(A; op = opA)
    descB = CuTensorDescriptor(B; op = opB)
    descC = CuTensorDescriptor(C; op = opC)
    @assert size(C) == size(D) && strides(C) == strides(D)
    descD = descC # must currently be identical
    typeCompute = cudaDataType(T)
    modeA = collect(Cint, Ainds)
    modeB = collect(Cint, Binds)
    modeC = collect(Cint, Cinds)
    modeD = modeC
    cutensorElementwiseTrinary(handle(),
                               T[alpha], A, descA, modeA,
                               T[beta],  B, descB, modeB,
                               T[gamma], C, descC, modeC,
                                         D, descD, modeD,
                               opAB, opABC, typeCompute, stream)
    return D
end

function elementwiseBinary!(
    alpha::Number, A::CuArray, Ainds::ModeType, opA::cutensorOperator_t,
    gamma::Number, C::CuArray{T}, Cinds::ModeType, opC::cutensorOperator_t,
    D::CuArray{T}, Dinds::ModeType, opAC::cutensorOperator_t;
    stream::CuStream=CuDefaultStream()) where {T}

    !is_unary(opA)    && throw(ArgumentError("opA must be a unary op!"))
    !is_unary(opC)    && throw(ArgumentError("opC must be a unary op!"))
    !is_binary(opAC)  && throw(ArgumentError("opAC must be a binary op!"))
    descA = CuTensorDescriptor(A; op = opA)
    descC = CuTensorDescriptor(C; op = opC)
    @assert size(C) == size(D) && strides(C) == strides(D)
    descD = descC # must currently be identical
    typeCompute = cudaDataType(T)
    modeA = collect(Cint, Ainds)
    modeC = collect(Cint, Cinds)
    modeD = modeC
    cutensorElementwiseBinary(handle(),
                              T[alpha], A, descA, modeA,
                              T[gamma], C, descC, modeC,
                                        D, descD, modeD,
                              opAC, typeCompute, stream)
    return D
end

function elementwiseBinary!(
    alpha::Number, A::Array, Ainds::ModeType, opA::cutensorOperator_t,
    gamma::Number, C::Array{T}, Cinds::ModeType, opC::cutensorOperator_t,
    D::Array{T}, Dinds::ModeType, opAC::cutensorOperator_t;
    stream::CuStream=CuDefaultStream()) where {T}

    !is_unary(opA)    && throw(ArgumentError("opA must be a unary op!"))
    !is_unary(opC)    && throw(ArgumentError("opC must be a unary op!"))
    !is_binary(opAC)  && throw(ArgumentError("opAC must be a binary op!"))
    descA = CuTensorDescriptor(A; op = opA)
    descC = CuTensorDescriptor(C; op = opC)
    @assert size(C) == size(D) && strides(C) == strides(D)
    descD = descC # must currently be identical
    typeCompute = cudaDataType(T)
    modeA = collect(Cint, Ainds)
    modeC = collect(Cint, Cinds)
    modeD = modeC
    cutensorElementwiseBinary(handle(),
                              T[alpha], A, descA, modeA,
                              T[gamma], C, descC, modeC,
                                        D, descD, modeD,
                              opAC, typeCompute, stream)
    return D
end

function elementwiseBinary!(
    alpha::Number, A::CuTensor, opA::cutensorOperator_t, gamma::Number, C::CuTensor{T}, opC::cutensorOperator_t, D::CuTensor{T}, opAC::cutensorOperator_t; stream::CuStream=CuDefaultStream()) where {T}

    !is_unary(opA)    && throw(ArgumentError("opA must be a unary op!"))
    !is_unary(opC)    && throw(ArgumentError("opC must be a unary op!"))
    !is_binary(opAC)  && throw(ArgumentError("opAC must be a binary op!"))
    descA = CuTensorDescriptor(A; op = opA)
    descC = CuTensorDescriptor(C; op = opC)
    @assert size(C) == size(D) && strides(C) == strides(D)
    descD = descC # must currently be identical
    typeCompute = cudaDataType(T)
    cutensorElementwiseBinary(handle(),
                              T[alpha], A.data, descA, A.inds,
                              T[gamma], C.data, descC, C.inds,
                                        D.data, descD, C.inds,
                              opAC, typeCompute, stream)
    return D
end

function permutation!(alpha::Number, A::CuArray, Ainds::ModeType,
                      B::CuArray, Binds::ModeType; stream::CuStream=CuDefaultStream())
    #!is_unary(opPsi)    && throw(ArgumentError("opPsi must be a unary op!"))
    descA = CuTensorDescriptor(A)
    descB = CuTensorDescriptor(B)
    T = eltype(B)
    typeCompute = cudaDataType(T)
    modeA = collect(Cint, Ainds)
    modeB = collect(Cint, Binds)
    cutensorPermutation(handle(), T[alpha], A, descA, modeA, B, descB, modeB, typeCompute,
                        stream)
    return B
end
function permutation!(alpha::Number, A::Array, Ainds::ModeType,
                      B::Array, Binds::ModeType; stream::CuStream=CuDefaultStream())
    #!is_unary(opPsi)    && throw(ArgumentError("opPsi must be a unary op!"))
    descA = CuTensorDescriptor(A)
    descB = CuTensorDescriptor(B)
    T = eltype(B)
    typeCompute = cudaDataType(T)
    modeA = collect(Cint, Ainds)
    modeB = collect(Cint, Binds)
    cutensorPermutation(handle(), T[alpha], A, descA, modeA, B, descB, modeB, typeCompute,
                        stream)
    return B
end

function contraction!(
    alpha::Number, A::CuArray, Ainds::ModeType, opA::cutensorOperator_t,
                   B::CuArray, Binds::ModeType, opB::cutensorOperator_t,
    beta::Number,  C::CuArray, Cinds::ModeType, opC::cutensorOperator_t,
                                                opOut::cutensorOperator_t;
    pref::cutensorWorksizePreference_t=CUTENSOR_WORKSPACE_RECOMMENDED,
    algo::cutensorAlgo_t=CUTENSOR_ALGO_DEFAULT, stream::CuStream=CuDefaultStream(),
    compute_type::Type=eltype(C), plan::Union{cutensorContractionPlan_t, Nothing}=nothing)

    !is_unary(opA)    && throw(ArgumentError("opA must be a unary op!"))
    !is_unary(opB)    && throw(ArgumentError("opB must be a unary op!"))
    !is_unary(opC)    && throw(ArgumentError("opC must be a unary op!"))
    !is_unary(opOut)  && throw(ArgumentError("opOut must be a unary op!"))
    descA = CuTensorDescriptor(A; op = opA)
    descB = CuTensorDescriptor(B; op = opB)
    descC = CuTensorDescriptor(C; op = opC)
    # for now, D must be identical to C (and thus, descD must be identical to descC)
    T = eltype(C)
    computeType = cutensorComputeType(T) #CUTENSOR_R_MIN_64F #TODO cudaDataType(T)
    modeA = collect(Cint, Ainds)
    modeB = collect(Cint, Binds)
    modeC = collect(Cint, Cinds)

    alignmentRequirementA = Ref{UInt32}(C_NULL)
    cutensorGetAlignmentRequirement(handle(), A, descA, alignmentRequirementA)
    alignmentRequirementB = Ref{UInt32}(C_NULL)
    cutensorGetAlignmentRequirement(handle(), B, descB, alignmentRequirementB)
    alignmentRequirementC = Ref{UInt32}(C_NULL)
    cutensorGetAlignmentRequirement(handle(), C, descC, alignmentRequirementC)
    desc = Ref(cutensorContractionDescriptor_t(ntuple(i->0, Val(256))))
    cutensorInitContractionDescriptor(handle(),
                                      desc,
                   descA, modeA, alignmentRequirementA[],
                   descB, modeB, alignmentRequirementB[],
                   descC, modeC, alignmentRequirementC[],
                   descC, modeC, alignmentRequirementC[],
                   computeType)
    find = Ref(cutensorContractionFind_t(ntuple(i->0, Val(64))))
    cutensorInitContractionFind(handle(), find, algo)

    @workspace fallback=1<<27 size=@argout(
            cutensorContractionGetWorkspace(handle(), desc, find, pref,
                                            out(Ref{UInt64}(C_NULL)))
        )[] workspace->begin
            plan_ref = Ref(cutensorContractionPlan_t(ntuple(i->0, Val(640))))
            if isnothing(plan)
                cutensorInitContractionPlan(handle(), plan_ref, desc, find, sizeof(workspace))
            else
                plan_ref = Ref(plan)
            end
            cutensorContraction(handle(), plan_ref,
                                T[convert(T, alpha)], A, B,
                                T[convert(T, beta)],  C, C,
                                workspace, sizeof(workspace), stream)
        end
    return C
end

function plan_contraction(
    A::CuArray, Ainds::ModeType, opA::cutensorOperator_t,
    B::CuArray, Binds::ModeType, opB::cutensorOperator_t,
    C::CuArray, Cinds::ModeType, opC::cutensorOperator_t,
                                                opOut::cutensorOperator_t;
    pref::cutensorWorksizePreference_t=CUTENSOR_WORKSPACE_RECOMMENDED,
    algo::cutensorAlgo_t=CUTENSOR_ALGO_DEFAULT, compute_type::Type=eltype(C))

    !is_unary(opA)    && throw(ArgumentError("opA must be a unary op!"))
    !is_unary(opB)    && throw(ArgumentError("opB must be a unary op!"))
    !is_unary(opC)    && throw(ArgumentError("opC must be a unary op!"))
    !is_unary(opOut)  && throw(ArgumentError("opOut must be a unary op!"))
    descA = CuTensorDescriptor(A; op = opA)
    descB = CuTensorDescriptor(B; op = opB)
    descC = CuTensorDescriptor(C; op = opC)
    # for now, D must be identical to C (and thus, descD must be identical to descC)
    computeType = cutensorComputeType(compute_type)
    T     = sizeof(compute_type) < sizeof(eltype(C)) ? eltype(C) : compute_type
    modeA = collect(Cint, Ainds)
    modeB = collect(Cint, Binds)
    modeC = collect(Cint, Cinds)

    alignmentRequirementA = Ref{UInt32}(C_NULL)
    cutensorGetAlignmentRequirement(handle(), A, descA, alignmentRequirementA)
    alignmentRequirementB = Ref{UInt32}(C_NULL)
    cutensorGetAlignmentRequirement(handle(), B, descB, alignmentRequirementB)
    alignmentRequirementC = Ref{UInt32}(C_NULL)
    cutensorGetAlignmentRequirement(handle(), C, descC, alignmentRequirementC)
    desc = Ref(cutensorContractionDescriptor_t(ntuple(i->0, Val(256))))
    cutensorInitContractionDescriptor(handle(),
                                      desc,
                   descA, modeA, alignmentRequirementA[],
                   descB, modeB, alignmentRequirementB[],
                   descC, modeC, alignmentRequirementC[],
                   descC, modeC, alignmentRequirementC[],
                   computeType)

    find = Ref(cutensorContractionFind_t(ntuple(i->0, Val(64))))
    cutensorInitContractionFind(handle(), find, algo)
    plan = Ref(cutensorContractionPlan_t(ntuple(i->0, Val(640))))
    workspace_size = Ref{UInt64}(C_NULL)
    cutensorContractionGetWorkspace(handle(), desc, find, pref, workspace_size)
    cutensorInitContractionPlan(handle(), plan, desc, find, workspace_size[])
    return plan[]
end

function reduction!(
    alpha::Number, A::CuArray, Ainds::ModeType, opA::cutensorOperator_t,
    beta::Number,  C::CuArray, Cinds::ModeType, opC::cutensorOperator_t,
    opReduce::cutensorOperator_t; stream::CuStream=CuDefaultStream())

    !is_unary(opA)    && throw(ArgumentError("opA must be a unary op!"))
    !is_unary(opC)    && throw(ArgumentError("opC must be a unary op!"))
    !is_binary(opReduce)  && throw(ArgumentError("opReduce must be a binary op!"))
    descA = CuTensorDescriptor(A; op = opA)
    descC = CuTensorDescriptor(C; op = opC)
    # for now, D must be identical to C (and thus, descD must be identical to descC)
    T = eltype(C)
    typeCompute = cutensorComputeType(T)
    modeA = collect(Cint, Ainds)
    modeC = collect(Cint, Cinds)

    @workspace fallback=1<<13 size=@argout(
            cutensorReductionGetWorkspace(handle(),
                A, descA, modeA,
                C, descC, modeC,
                C, descC, modeC,
                opReduce, typeCompute,
                out(Ref{UInt64}(C_NULL)))
        )[] workspace->begin
            cutensorReduction(handle(),
                T[alpha], A, descA, modeA,
                T[beta],  C, descC, modeC,
                        C, descC, modeC,
                opReduce, typeCompute,
                workspace, sizeof(workspace), stream)
        end

    return C
end

function cutensorComputeType(T::DataType)
    if T == Float32
        return CUTENSOR_R_MIN_32F
    elseif T == ComplexF32
        return CUTENSOR_C_MIN_32F
    elseif T == Float16
        return CUTENSOR_R_MIN_16F
    elseif T == ComplexF16
        return CUTENSOR_C_MIN_16F
    elseif T == Float64
        return CUTENSOR_R_MIN_64F
    elseif T == ComplexF64
        return CUTENSOR_C_MIN_64F
    elseif T == Int8
        return CUTENSOR_R_MIN_8I
    elseif T == Int32
        return CUTENSOR_R_MIN_32I
    elseif T == UInt8
        return CUTENSOR_R_MIN_8U
    elseif T == UInt32
        return CUTENSOR_R_MIN_32U
    else
        throw(ArgumentError("cutensorComputeType equivalent for input type $T does not exist!"))
    end
end


text:
#using Base.Cartesian

"""
    auc(scores, events, increasing = true)

compute the Area Under the receiver operator Curve (AUC), given some output `scores` array and some ground truth (`events`).
By default, it is assumed, that the `scores` are ordered increasingly (`increasing = true`), i.e. high scores represent events.
AUC is computed according to Fawcett, T. (2006). An introduction to ROC analysis. Pattern Recognition Letters, 27(8), 861–874. http://doi.org/10.1016/j.patrec.2005.10.010

# Examples
```
julia> scores = rand(10, 2)
julia> events = rand(0:1, 10, 2)
julia> auc(scores, events)
julia> auc(scores, boolevents(events))
```
"""
function auc(scores,events; increasing::Bool = true)
    s = sortperm(reshape(scores,length(scores)),rev=increasing);
    length(scores) == length(events) || error("Scores and events must have same number of elements")
    f=scores[s]
    L=events[s]
    fp=0
    tp=0
    fpprev=0
    tpprev=0
    A=0.0
    fprev=-Inf
    P=sum(L)
    N=length(L)-P
    for i=1:length(L)
        if f[i]!=fprev
            A+=trap_area(fp,fpprev,tp,tpprev)
            @inbounds fprev=f[i]
            fpprev=fp
            tpprev=tp
        end
        if isextreme(L[i])
            tp+=1
        else
            fp+=1
        end
    end
    A+=trap_area(N,fpprev,P,tpprev)
    A=A/(P*N)
end

function trap_area(x1,x2,y1,y2)
    b=abs(x1-x2)
    h=0.5*(y1+y2)
    return b*h
end

isextreme(l::Bool)=l
isextreme(l::Integer)=l>0

###################################
# end


text:
module JuliaDBMeta

using IndexedTables, MacroTools, Reexport
using IndexedTables: transform
import IndexedTables: select
import IterTools

import JuliaDB: Dataset, DDataset, fromchunks
import Dagger: delayedmap

@reexport using JuliaDB
export @with, @map, @byrow!, @filter
export @transform, @transform_vec, transformcol
export @where, @where_vec
export @groupby
export @apply, @applychunked
export @=>, @select

include("utils.jl")
include("byrow.jl")
include("with.jl")
include("map.jl")
include("filter.jl")
include("where.jl")
include("transform.jl")
include("groupby.jl")
include("apply.jl")
include("pair.jl")

end


text:
## Show the effect of seeded PageRank on MNIST

## Load the MNIST Graph and draw it
##
using DelimitedFiles
using SparseArrays
using LinearAlgebra
ei,ej = copy.(eachcol(Int.(readdlm("mnist-train-4.edges"))))
xy = readdlm("mnist-train-4.xy")
n = maximum(ei)
A = sparse(ei,ej,1,n,n)
@assert issymmetric(A)

##
using Random
function simple_spectral_eigenvector(A,nsteps=500)
  @assert issymmetric(A) # expensive, but useful...
  Random.seed!(1)
  n = size(A,1)
  d = vec(sum(A,dims=1))
  nd2 = d'*d
  x = randn(n)
  # project x orthogonal to d
  x .-= ((d'*x)/(d'*d)).*d
  x ./= x'*(d.*x) # normalize
  for i=1:nsteps
    x = (A*x)./d          # matrix-vector
    x .-= ((d'*x)/nd2).*d
    x ./= sqrt(x'*(d.*x)) # normalize
    # make sure the first component is positive
    x .*= sign(x[1])
  end

  return x
end
f = simple_spectral_eigenvector(A)
##
using MatrixNetworks
f2 = fiedler_vector(A)[1]
# make sure the first component is positive
f2 .*= sign(f2[1])
##
norm(f-f2)
## why different? convergence is slow :(
f = simple_spectral_eigenvector(A,2500)
norm(f-f2)
## why different? convergence is slow :(
f = simple_spectral_eigenvector(A,4000)
norm(f-f2)
##
using Plots
x = f
p = sortperm(x) # this gets things to draw in the right z-order
scatter(xy[p,1],xy[p,2], marker_z=x[p],colorbar=false,alpha=0.5,
  color=:magma,
  size=(1200,1200),
  markerstrokecolor=:white,markerstrokewidth=0.0,label="",framestyle=:none,
  background=nothing,dpi=220)
savefig("spectral-mnist-500.png")
##
using Plots
x = f2
p = sortperm(x) # this gets things to draw in the right z-order
scatter(xy[p,1],xy[p,2], marker_z=x[p],colorbar=false,alpha=0.5,
  color=:magma,
  size=(1200,1200),
  markerstrokecolor=:white,markerstrokewidth=0.0,label="",framestyle=:none,
  background=nothing,dpi=220)
savefig("spectral-mnist-converged.png")
## Let's watch spectra converge
using Plots
using GraphRecipes
Random.seed!(1)
function spectral_convergence_animation(A,xy,nsteps_per_frame,nframes)
  @assert issymmetric(A) # expensive, but useful...
  Random.seed!(1)
  n = size(A,1)
  d = vec(sum(A,dims=1))
  sumd = sum(d)
  x = randn(n)
  # project x orthogonal to d
  x .-= (d'*x)/sum(d)
  x ./= x'*(d.*x) # normalize

  anim = @animate for i=1:nframes

    f = zeros(n)
    for j=1:nsteps_per_frame
      x = (A*x)./d          # matrix-vector
      x .-= (d'*x)/sum(d)   # project out
      x ./= sqrt(x'*(d.*x)) # normalize
      # make sure the first component is positive
      #x .*= sign(x[1])
      f .+= x./nsteps_per_frame
    end

    p = sortperm(f) # this gets things to draw in the right z-order
    scatter(xy[p,1],xy[p,2], marker_z=f[p],colorbar=false,alpha=0.5,
      color=:magma,
      size=(1200,1200),dpi=100,
      markerstrokecolor=nothing,markerstrokewidth=0,label="",framestyle=:none,background=:black)
  end
  return anim
end
anim = spectral_convergence_animation(A,xy,10,500)
gif(anim,"spectral-convergence-mnist.gif",fps=30)
##
anim = spectral_convergence_animation(A,xy,1,500)
gif(anim,"spectral-convergence-mnist-1-500.gif",fps=30)

## Now show the entire graph
using GraphRecipes
using LinearAlgebra
using Plots

ei, ej = findnz(A)
graphplot(ei, ej, x =xy[:,1], y=xy[:,2],
  markercolor=:black, markerstrokecolor=:white,
  size=(1200,1200),dpi=300,
  markersize=0, linecolor=1, linealpha=0.2, linewidth=0.5,
  markeralpha=0.2,colorbar=false,
  axis_buffer=0.02, background=nothing)
p=sortperm(pr)
scatter!(xy[p,1],xy[p,2],marker_z=-abs.(log10.(pr[p])).^0.5,markerstrokewidth=0,alpha=0.5)
savefig("pagerank-mnist.png")


text:
function _count_by_state(events::Events{T},
                         state::DiseaseState,
                         time::Float64) where T <: EpidemicModel
  local n_ids
  if state ∉ _state_progressions[T]
    @error "Invalid state specified"
  end
  if state == _state_progressions[T][1] # S
    nextstate = advance(state, T) # Either E or I
    n_ids = sum(events[nextstate] .> Ref(time)) # E/I after `time`
    n_ids += sum(isnan.(events[nextstate])) # Never E/I
  elseif state in _state_progressions[T][2:end-1]
    nextstate = advance(state, T) # Either I or R
    n_ids = sum((events[state] .<= Ref(time)) .& (events[nextstate] .> Ref(time))) # E/I at or before time and I/R after time
    n_ids += sum((events[state] .<= Ref(time)) .& isnan.(events[nextstate])) # E/I at or before time and never I/R
  elseif state == _state_progressions[T][end] # I or R
    n_ids = sum(events[state] .<= Ref(time)) # I/R at or before time
  end
  @debug "$n_ids individual(s) in state $state at t = $time"
  return n_ids
end

function _epidemic_curve(events::Events{T},
                         state::DiseaseState,
                         min::Float64,
                         max::Float64) where T <: EpidemicModel
  if min >= max
    @error "Minimum time must be less than maximum time"
  end
  local times
  if state == _state_progressions[T][1]
    nextstate = _state_progressions[T][2]
    times = events[nextstate]
  elseif state in _state_progressions[T][2:end-1]
    nextstate = advance(state, T)
    times = events[[state; nextstate]][:]
  elseif state == _state_progressions[T][end]
    times = events[state]
  else
    @error "Invalid state specified"
  end
  times = times[Ref(min) .< times .< Ref(max)]
  sort!(times)
  insert!(times, 1, min)
  push!(times, max)
  counts = [_count_by_state(events, state, t) for t in times]
  return times, counts
end

@recipe function f(events::Events{T},
                   state::DiseaseState,
                   min::Float64,
                   max::Float64) where T<: EpidemicModel
  xguide --> "Time"
  yguide --> "N"
  xlims --> (min - 1.0, max + 1.0)
  linewidth --> 2.0
  linecolor --> :cornflowerblue
  label --> ""
  seriestype --> :steppost
  _epidemic_curve(events, state, min, max)
end

@recipe function f(events::Events{T},
                   state::DiseaseState) where T <: EpidemicModel
  events, state, 0.0, maximum(events)
end

@recipe function f(events::Events{T},
                   min::Float64,
                   max::Float64) where T<: EpidemicModel
  @series begin
    linecolor --> :purple
    label --> "S"
    events, State_S, min, max
  end
  if T in [SEIR; SEI]
    @series begin
      linecolor --> :lightblue4
      label --> "E"
      events, State_E, min, max
    end
  end
  @series begin
    linecolor --> :lightgreen
    label --> "I"
    events, State_I, min, max
  end
  if T in [SEIR; SIR]
    @series begin
      linecolor --> :yellow
      label --> "R"
      events, State_R, min, max
    end
  end
end

@recipe function f(events::Events{T}) where T <: EpidemicModel
  events, 0.0, maximum(events)
end


text:
macro def_fit_sepp(params_to_fit)
    l = (qn -> qn.value).(params_to_fit.args)
    n = length(l)
    x0_expr = Expr[]
    fitted_params = :(setproperty!(sepp, $(l[1]), Optim.minimizer(res)[1]))

    h_it = Expr[]

    kw_expr_it = Expr[]

    for i in 1:n
        push!(x0_expr, :(θ[$(l[i])]))
        push!(kw_expr_it, :($(l[i]) => x[$i]))
        push!(h_it, :(getproperty(sepp, $(l[i]))))
    end

    for i in 2:n
        fitted_params = :($fitted_params ; $(:(setproperty!(sepp, $(l[i]), Optim.minimizer(res)[$i]))))
    end

    kw_expr = :(kw = Dict($(kw_expr_it...)))

    to_min_expr = Expr(:function, Expr(:call, :to_min, :x), quote
        default_kw = Dict{Symbol, Real}(:μ => θ[:μ], :ϕ => θ[:ϕ], :γ => θ[:γ])
        $kw_expr
        final_kw = merge(default_kw, kw)
        return negloglik(ts; final_kw...)
    end)

    h = :([$(h_it...)])

    x0_expr = :([$(x0_expr...)])

    m = esc.(l)

    return Expr(:function, Expr(:call, :($(esc(:fit!))), Expr(Symbol("::"), :sepp, :SEPPExpKern), m...),quote

        ts = sepp.data
        isnothing(ts) && error("No data in model, can't fit")

        θ = params(sepp)

        $to_min_expr

        lower = zeros($n)
        upper = fill(Inf, $n)

        res = optimize(to_min, lower, upper, $x0_expr)

        $fitted_params

        hess = ForwardDiff.hessian(to_min, $h)

        sepp.cov_mat = inv(hess)

        return res

    end)
end

macro def_fit_sempp(params_to_fit)
    l = (qn -> qn.value).(params_to_fit.args)
    n = length(l)
    x0_expr = Expr[]
    fitted_params = :(setproperty!(sempp, $(l[1]), Optim.minimizer(res)[1]))

    h_it = Expr[]

    kw_expr_it = Expr[]

    for i in 1:n
        push!(x0_expr, :(θ[$(l[i])]))
        push!(kw_expr_it, :($(l[i]) => x[$i]))
        push!(h_it, :(getproperty(sempp, $(l[i]))))
    end

    for i in 2:n
        fitted_params = :($fitted_params ; $(:(setproperty!(sempp, $(l[i]), Optim.minimizer(res)[$i]))))
    end

    kw_expr = :(kw = Dict($(kw_expr_it...)))

    to_min_expr = Expr(:function, Expr(:call, :to_min, :x), quote
        default_kw = Dict{Symbol, Real}(:μ => θ[:μ], :ϕ => θ[:ϕ], :γ => θ[:γ])
        $kw_expr
        final_kw = merge(default_kw, kw)
        return negloglik(mts, markdens, impact_func; final_kw...)
    end)

    h = :([$(h_it...)])

    x0_expr = :([$(x0_expr...)])

    m = esc.(l)

    return Expr(:function, Expr(:call, :($(esc(:fit!))), Expr(:parameters, Expr(:kw, :lower, zeros(n)), Expr(:kw, :upper, fill(Inf, n))), Expr(Symbol("::"), :sempp, :SEMPPExpKern), m...),quote

        mts = sempp.data
        isnothing(mts) && error("No data in model, can't fit")

        θ = params(sempp)
        markdens = θ[:markdens]
        impact_func = θ[:impact_function]

        $to_min_expr

        res = optimize(to_min, lower, upper, $x0_expr)

        $fitted_params

        hess = ForwardDiff.hessian(to_min, $h)

        sempp.cov_mat = inv(hess)

        return res

    end)
end