text:
module Analysis
using Sessions, TaskData, SpikeTools, InformationTools, JLD, Distributions,
      DataFrames, DataStructures, StatsBase, HypothesisTests

"""
    mapcells(cellfunc, sessfunc, name, [extradict, [reduce]])

The main analysis driver routine. Does the following things:

1. Calls `sessresult = sessfunc(session::Session)` on each session.
2. Calls
   `cellfunc(session::Session, sessresult, spikes::Vector{Float64}, unit::ClusterMetadata)`
   for each unit. If this method returns `nothing`, the unit is skipped. Otherwise it
   must return a Dict of data corresponding to the given unit. The keys of the Dict must
   be the same for all units.
3. For each key of the Dict returned by `cellfunc`, concatenates values across units into
   vectors.
4. Saves vectors in file `name.jld` named by key, along with `chs::Vector{ClusterMetadata}`
   (the metadata objects for analyzed units), and contents of `extradict` if provided.

If a `reduce` Dict is provided, for keys present in `reduce`, values are reduced across
units by calling the value (a function) for that key, rather than concatenated.
"""
function mapcells{F1,F2}(cellfunc::F1, sessfunc::F2, name::Symbol,
                         extradict::Union{Associative,Void}=nothing,
                         reduce::Dict{Symbol}=Dict{Symbol,Any}())
    sess = sessionnames()
    res = pmap(sess) do s
        sdata = readsession(s)
        sessres = sessfunc(sdata)
        spikes = sdata.spikes
        chs = ClusterMetadata[]
        dicts = Dict{Symbol}[]
        local reduceddict = Dict{Symbol,Any}()

        for ich = 1:length(spikes)
            isdefined(spikes, ich) || continue
            clusters = spikes[ich]
            for icluster = 1:length(clusters)
                sdata.unit_type[ich][icluster] == :A && continue
                ch = ClusterMetadata(sdata, ich, icluster)
                println(ch)
                local res = cellfunc(sdata, sessres, clusters[icluster], ch)
                isa(res, Void) && continue
                push!(chs, ch)
                push!(dicts, res)
                if !isempty(reduce)
                    for (k, v) in reduce
                        reduceddict[k] = haskey(reduceddict, k) ? v(reduceddict[k], res[k]) : res[k]
                        delete!(res, k)
                    end
                end
            end
        end
        (chs, dicts, reduceddict)
    end
    filter!(x->!isempty(x[1]), res)

    dictout::Dict{Symbol} = Dict(k=>typeof(v)[] for (k, v) in res[1][2][1])
    reduceddict = Dict{Symbol,Any}()
    for x in res
        if !isempty(reduce)
            for (k, v) in reduce
                reduceddict[k] = haskey(reduceddict, k) ? v(reduceddict[k], x[3][k]) : x[3][k]
            end
        end
        for d::Dict{Symbol} in x[2]
            for k::Symbol in keys(d)
                push!(dictout[k], d[k])
            end
        end
    end
    !isempty(reduce) && (dictout = merge(dictout, reduceddict))

    j = jldopen(joinpath(PROCESSED_DATA_DIR, "$(name).jld"), "w")
    try
        write(j, "chs", vcat([x[1] for x in res]...))
        for (k, v) = dictout
            write(j, string(k), v)
        end
        if !isa(extradict, Void)
            for (k, v) = extradict
                write(j, string(k), v)
            end
        end
    finally
        close(j)
    end
end

"""
    compute_permpval()

Computes p-values for selectivity at image presentation, using Poisson likelihood ratio
permutation test.
"""
function compute_permpval()
    NPERM = 10000
    BINS = 100:10:1000
    function sessfunc(session::Session)
        trials = stimulusframe(session.trials)
        onsets = trials[:onset]
        df = DataFrame(y=pool(trials[:stimulus]), onset=trials[:stimulus])
        dmat = ModelMatrix(ModelFrame(onset ~ y, df)).m
        perms = makeperms(size(df, 1), NPERM)

        stimulus = convert(Vector{Int}, trials[:stimulus])
        inds = vcat([find(stimulus .== i) for i = 1:maximum(stimulus)]...)
        ningroup = [sum(stimulus .== i) for i = 1:maximum(stimulus)]
        (onsets[inds], perms, ningroup)
    end

    mapcells(sessfunc, :permpval, Dict(:nperm => NPERM, :bins => BINS)) do session, x, spikes, ch
        (onsets, perms, ningroup) = x
        rast = raster(spikes, onsets, BINS/1000)'
        @time pval = poisson_deviance_grid_pval(rast, ningroup, perms)
        Dict(:pval => pval)
    end
end

"""
    compute_permpval_eachisi()

Computes p-values for selectivity at image presentation separately for each ISI, using
Poisson likelihood ratio permutation test.
"""
function compute_permpval_eachisi()
    NPERM = 10000
    BINS = 100:10:1000
    function sessfunc(session::Session)
        trials = stimulusframe(session.trials)
        onsets = trials[:onset].data
        isi = convert(Vector{Float64}, trials[:isi])
        perms = makeperms(sum(isi .== ISI[1]), NPERM)
        stimulus = convert(Vector{Int}, trials[:stimulus])
        inds = [vcat([find((stimulus .== i) & (isi .== x)) for i = 1:maximum(stimulus)]...) for x in ISI]
        ningroup = [[sum((stimulus .== i) & (isi .== x)) for i = 1:maximum(stimulus)] for x in ISI]
        (onsets, perms, inds, ningroup)
    end

    mapcells(sessfunc, :permpval_eachisi, Dict(:nperm => NPERM, :bins => BINS)) do session, x, spikes, ch
        ch.unit_type != :A || return
        (onsets, perms, inds, ningroup) = x
        rast = raster(spikes, onsets, BINS/1000)
        pval = zeros(length(ISI))
        for iisi = 1:length(ISI)
            @time pval[iisi] = poisson_deviance_grid_pval(rast[:, inds[iisi]]', ningroup[iisi], perms)
        end
        Dict(:pval => pval)
    end
end

"""
    compute_latency()

Computes latencies and durations of spiking responses by finding the window yielding the
maximal Poisson likelihood ratio.
"""
function compute_latency()
    LASTPT = 1500
    LASTPT_OPTWIN = 1000
    BASELINE_EPOCH = -1000:1000:0
    function sessfunc(session::Session)
        trials = stimulusframe(session.trials)
        @assert sort(unique(trials[:isi])) == ISI
        isitrials = [find(trials[:isi] .== isi) for isi in ISI]
        onset = trials[:onset].data
        offset = trials[:offset].data
        pos = convert(Vector{Int}, trials[:pos])
        df = DataFrame(y=pool(trials[:stimulus]), onset=onset)
        dmat = ModelMatrix(ModelFrame(onset ~ y, df)).m
        stimulus = convert(Vector{Int}, trials[:stimulus])
        trial_start_time = [session.trials[:onsets][i][1] for i = 1:size(session.trials, 1)]
        inds = vcat([find(stimulus .== i) for i = 1:maximum(stimulus)]...)
        ningroup = [sum(stimulus .== i) for i = 1:maximum(stimulus)]
        (isitrials, onset, offset, stimulus, inds, ningroup, trial_start_time, pos)
    end

    mapcells(sessfunc, :latency, Dict(:lastpt_optwin => LASTPT_OPTWIN, :lastpt => LASTPT)) do session, x, spikes, ch
        (isitrials, onset, offset, stimulus, inds, ningroup, trial_start_time, pos) = x
        nstim = maximum(stimulus)

        rast = raster(spikes, onset, (1:LASTPT)/1000)
        eta2_grid = poisson_deviance_grid(rast[1:LASTPT_OPTWIN, inds]', ningroup)
        ends, starts = ind2sub(size(eta2_grid), find(eta2_grid .== maximum(eta2_grid)))
        fr_optimal_window = zeros(nstim)
        if isempty(ends)
            stim_window = (-1, -1)
        else
            best = indmin(ends - starts)
            stim_window = (starts[best], ends[best]+1)
        end

        isi_windows = fill((-1, -1), length(isitrials))
        isi_windows_notfirst = fill((-1, -1), length(isitrials))

        for iisi = 1:length(isitrials)
            sortedisitrials = Int[]
            sortedisitrials_notfirst = Int[]
            ningroupisi = Int[]
            ningroupisi_notfirst = Int[]
            for istim = 1:nstim
                stimsel = isitrials[iisi][stimulus[isitrials[iisi]] .== istim]
                append!(sortedisitrials, stimsel)
                push!(ningroupisi, length(stimsel))
                stimsel_notfirst = isitrials[iisi][(stimulus[isitrials[iisi]] .== istim) & (pos[isitrials[iisi]] .!= 1)]
                append!(sortedisitrials_notfirst, stimsel_notfirst)
                push!(ningroupisi_notfirst, length(stimsel_notfirst))
            end

            eta2_grid = poisson_deviance_grid(rast[:, sortedisitrials]', ningroupisi)
            ends, starts = ind2sub(size(eta2_grid), find(eta2_grid .== maximum(eta2_grid)))
            if !isempty(ends)
                best = indmin(ends - starts)
                isi_windows[iisi] = (starts[best], ends[best]+1)
            end

            eta2_grid = poisson_deviance_grid(rast[:, sortedisitrials_notfirst]', ningroupisi_notfirst)
            ends, starts = ind2sub(size(eta2_grid), find(eta2_grid .== maximum(eta2_grid)))
            if !isempty(ends)
                best = indmin(ends - starts)
                isi_windows_notfirst[iisi] = (starts[best], ends[best]+1)
            end
        end

        bsdata = raster(spikes, trial_start_time, BASELINE_EPOCH/1000)
        bsscale = 1000/(last(BASELINE_EPOCH)-first(BASELINE_EPOCH))

        Dict(:window => stim_window,
             :isi_windows => isi_windows,
             :isi_windows_notfirst => isi_windows_notfirst,
             :baseline_mean => mean(bsdata)*bsscale,
             :baseline_std => std(bsdata)*bsscale
            )
    end
end

"""
    compute_alignedinfoisi()

Computes ω^2 for stimulus and interaction in different windows, along with the sum of
permutations across units.
"""
function compute_alignedinfoisi(; pos1::Bool=false)
    @load joinpath(PROCESSED_DATA_DIR, "latency.jld") window chs
    SMOOTHBY = 200
    NPERM = 1000
    PLOT_EPOCH = (-200, 1001)
    selch = allisisigunits()
    function sessfunc(session::Session)
        trials = stimulusframe(session.trials)
        nstim = maximum(trials[:stimulus])::Int

        trial_number = convert(Vector{Int}, trials[:trial_number])
        stimulus = convert(Vector{Int}, trials[:stimulus])
        dmat = zeros(length(stimulus), nstim)
        prevstim_dmat = zeros(length(stimulus), nstim)
        nextstim_dmat = zeros(length(stimulus), nstim)
        prevstim = zeros(length(stimulus))
        nextstim = zeros(length(stimulus))
        for i = 1:length(stimulus)
            if i+1 <= length(stimulus) && trial_number[i+1] == trial_number[i]
                nextstim_dmat[i, stimulus[i+1]] = 1
                nextstim[i] = stimulus[i+1]
            end
            if i-1 > 0 && trial_number[i-1] == trial_number[i]
                prevstim[i] = stimulus[i-1]
                prevstim_dmat[i, stimulus[i-1]] = 1
            end
            dmat[i, stimulus[i]] = 1
        end
        if pos1
            keep = trials[:pos] .== 1
            prevstim_dmat = zeros(sum(keep), 0)
            stimulus = stimulus[keep]
            dmat = dmat[keep, :]
            nextstim_dmat = nextstim_dmat[keep, 2:end]
            trials = trials[keep, :]
        end

        isitrials = Vector{Vector{Int}}(length(ISI))
        isiningroup = Vector{Vector{Int}}(length(ISI))
        for (iisi, isi) = enumerate(ISI)
            isitrials[iisi] = curisitrials = Int[]
            isiningroup[iisi] = curningroup = zeros(Int, nstim)
            for igrp = 1:nstim
                idx = find((trials[:isi] .== isi) & (trials[:stimulus] .== igrp))
                append!(curisitrials, idx)
                curningroup[igrp] = length(idx)
            end
        end

        isis = Vector{Float64}(trials[:isi])
        isiperms = fill(makeperms(length(isitrials[1]), NPERM), length(ISI))

        isi2perms = [begin
            perms = zeros(Int, length(isitrials[iisi])*2, NPERM)
            isistim = stimulus[(isis .== ISI[iisi]) | (isis .== ISI[iisi+1])]
            for stim = 1:nstim
                idx = find(isistim .== stim)
                shuffidx = copy(idx)
                for ibs = 1:NPERM
                    perms[idx, ibs] = shuffle!(shuffidx)
                end
            end
            perms
        end for iisi in 1:length(ISI)-1]

        (trials[:onset].data, isitrials, isiningroup, isis, dmat, prevstim_dmat, nextstim_dmat, isiperms, isi2perms)
    end

    mapcells(sessfunc, Symbol("alignedinfoisi$(pos1 ? "_pos1" : "")"),
             Dict(:fwhm => SMOOTHBY, :plot_epoch => PLOT_EPOCH, :nperm => NPERM),
             Dict(:alignedinfoperms => (a, b)->broadcast!(+, a, a, b),
                  :alignedinfoisidiffperms => (a, b)->broadcast!(+, a, a, b))) do session, x, spikes, ch
        ch in selch || return nothing
        (onset, isitrials, isiningroup, isi, dmat, prevstim_dmat, nextstim_dmat, isiperms, isi2perms) = x
        win = window[findfirst(chs, ch)]

        extract_times, extract_range, times = paddedbins(PLOT_EPOCH, 1, 5*SMOOTHBY)
        rast = boxcarsmooth(raster(spikes, onset, extract_times+win[1]/1000), SMOOTHBY, Int)[extract_range, :]

        all_dmat = [dmat nextstim_dmat prevstim_dmat]
        null_dmat = [ones(size(dmat, 1)) nextstim_dmat prevstim_dmat]
        alignedinfo = zeros(size(rast, 1), length(isitrials))
        alignedinfoperms = zeros(size(rast, 1), NPERM, length(isitrials))
        alignedinfoisidiff = zeros(size(rast, 1), length(isitrials)-1)
        alignedinfoisidiffperms = zeros(size(rast, 1), NPERM, length(isitrials)-1)
        @time for iisi = 1:length(isitrials)
            sel = isi .== ISI[iisi]
            isinulldmat = null_dmat[sel, :]
            isinulldmatqr = qrfactchk(isinulldmat)
            isidmat = all_dmat[sel, :]
            isidmatqr = qrfactchk(isidmat)
            isirast = rast[:, sel]'
            fullrss = computessr(isidmatqr, isirast)
            reducedrss = computessr(isinulldmatqr, isirast)
            fulldfe = size(isirast, 1) - size(all_dmat, 2)
            reduceddfe = size(isirast, 1) - size(isinulldmat, 2)
            omega2!(view(alignedinfo, :, iisi), fullrss, fulldfe, reducedrss, reduceddfe, size(isirast, 1))

            permisirast = zeros(Float64, size(isirast))
            isirasttmp = zeros(Float64, size(isirast))
            @assert size(isiperms[iisi], 1) == size(isirast, 1)
            for iperm = 1:NPERM
                copy!(permisirast, view(isirast, view(isiperms[iisi], :, iperm), :))
                computessr!(fullrss, isidmatqr, copy!(isirasttmp, permisirast))
                computessr!(reducedrss, isinulldmatqr, permisirast)
                omega2!(view(alignedinfoperms, :, iperm, iisi), fullrss, fulldfe, reducedrss, reduceddfe, size(isirast, 1))
            end

            if iisi <= length(isitrials) - 1
                sel = find((isi .== ISI[iisi]) | (isi .== ISI[iisi+1]))
                isirast = rast[:, sel]'
                nulldmat = zeros(length(sel), (size(null_dmat, 2)-1)*2)
                dmat2 = zeros(length(sel), size(dmat, 2)*2)
                for ipt = 1:length(sel)
                    if isi[sel[ipt]] == ISI[iisi]
                        nulldmat[ipt, 1:size(null_dmat, 2)-1] = view(null_dmat, sel[ipt], 2:size(null_dmat, 2))
                        dmat2[ipt, 1:size(dmat, 2)] = view(dmat, sel[ipt], :)
                    else
                        nulldmat[ipt, size(null_dmat, 2):end] = view(null_dmat, sel[ipt], 2:size(null_dmat, 2))
                        dmat2[ipt, size(dmat, 2)+1:end] = view(dmat, sel[ipt], :)
                    end
                end
                nulldmatqr = qrfactchk([isi[sel] .== ISI[iisi] nulldmat dmat[sel, :]])
                dmatqr = qrfactchk([nulldmat dmat2])
                permisirast = zeros(Float64, size(isirast))
                isirasttmp = zeros(Float64, size(isirast))
                reduceddfe2 = size(isirast, 1) - size(nulldmatqr, 2)
                fulldfe2 = size(isirast, 1) - size(dmatqr, 2)
                computessr!(fullrss, dmatqr, copy!(isirasttmp, isirast))
                computessr!(reducedrss, nulldmatqr, copy!(isirasttmp, isirast))
                omega2!(view(alignedinfoisidiff, :, iisi), fullrss, fulldfe2, reducedrss, reduceddfe2, size(isirast, 1))
                @assert size(isi2perms[iisi], 1) == size(isirast, 1)
                for iperm = 1:NPERM
                    copy!(permisirast, view(isirast, view(isi2perms[iisi], :, iperm), :))
                    computessr!(fullrss, dmatqr, copy!(isirasttmp, permisirast))
                    computessr!(reducedrss, nulldmatqr, permisirast)
                    omega2!(view(alignedinfoisidiffperms, :, iperm, iisi), fullrss, fulldfe2, reducedrss, reduceddfe2, size(isirast, 1))
                end
            end
        end
        alignedinfo[isnan(alignedinfo)] = 0
        alignedinfoperms[isnan(alignedinfoperms)] = 0
        alignedinfoisidiff[isnan(alignedinfoisidiff)] = 0
        alignedinfoisidiffperms[isnan(alignedinfoisidiffperms)] = 0

        Dict(:window => win,
             :alignedinfo => alignedinfo,
             :alignedinfoperms => alignedinfoperms,
             :alignedinfoisidiff => alignedinfoisidiff,
             :alignedinfoisidiffperms => alignedinfoisidiffperms)
    end
end

"""
    shuffleby(x::AbstractVector, NPERM)::Matrix{Int}

Create a `length(x)` × `NPERM` matrix of indices for permutations of `x`, permuting within
levels of `x`.
"""
function shuffleby(x::AbstractVector, NPERM)
    perms = zeros(Int, length(x), NPERM)
    for val in unique(x)
        inds = find(x .== val)
        shuffled = copy(inds)
        for i = 1:NPERM
            perms[inds, i] = shuffle!(shuffled)
        end
    end
    perms
end

"""
    rearrange(x::AbstractVector)

Returns indices that rearrange `x` so that levels are adjacent, as well as counts for each
group.
"""
function rearrange(x::AbstractVector)
    uniques = sort!(unique(x))
    inds = vcat([find(x .== i) for i = uniques]...)
    ningroup = [sum(x .== i) for i = uniques]
    (inds, ningroup)
end

"""
    compute_permpval_isiwin()

Performs permutation test for image selectivity in different bins, for different groups of
ISIs.
"""
function compute_permpval_isiwin()
    selch = permselunits()
    @load joinpath(PROCESSED_DATA_DIR, "latency.jld") chs window
    NPERM = 10000
    WIN = [(0, 200), (200, 400), (400, 600), (700, 900)]
    function sessfunc(session::Session)
        trials = stimulusframe(session.trials)
        isi = convert(Vector{Float64}, trials[:isi])
        stimulus = convert(Vector{Float64}, trials[:stimulus])
        tgrp = Vector{Vector{Int}}(length(ISI))
        ningroup = Vector{Vector{Int}}(length(ISI))
        for i = 1:length(ISI)
            idx = find(isi .>= ISI[i])
            inds, ningroup[i] = rearrange(stimulus[idx])
            tgrp[i] = idx[inds]
        end
        next_stimulus = convert(Vector{Int}, trials[:next_stimulus])
        perms = [shuffleby(next_stimulus[x], NPERM) for x in tgrp]
        (trials[:onset].data, tgrp, ningroup, perms)
    end

    mapcells(sessfunc, :permpval_isiwin, Dict(:nperm => NPERM, :win => WIN, :minisi => ISI)) do session, x, spikes, ch
        ch in selch || return nothing
        stim_window = window[findfirst(chs, ch)]

        (onset, tgrp, ningroup, perms) = x
        pval_cumulative = ones(length(ISI), length(WIN))
        @time for (i, win) = enumerate(WIN)
            rast = vec(raster(spikes, onset+stim_window[1]/1000, (win[1]:win[2]-win[1]:win[2])/1000))
            for j = 1:length(ningroup)
                y = rast[tgrp[j]]
                pval_cumulative[j, i] = poisson_deviance_pval(y, ningroup[j], perms[j])
            end
        end
        Dict(
            :pval_cumulative => pval_cumulative
        )
    end
end

"""
    compute_permpval_maintenance()

Performs permutation test for maintenance period stimulus selectivity.
"""
function compute_permpval_maintenance()
    NPERM = 100000
    WIN = (300, 2400)
    function sessfunc(session::Session)
        trials = session.trials
        nstim = maximum(map(maximum, trials[:stimuli]))
        X = zeros(size(trials, 1), nstim)
        for i = 1:size(trials, 1)
            X[i, trials[i, :stimuli]::Vector{Int}] = 1
        end
        perms = makeperms(size(X, 1), NPERM)
        (trials[:maint].data, qrfactchk(X), perms)
    end

    mapcells(sessfunc, :permpval_maintenance, Dict(:nperm => NPERM, :win => WIN)) do session, x, spikes, ch
        (maint, qrf, perms) = x
        null = qrfact(ones(size(qrf, 1), 1))
        rast = map(Float64, raster(spikes, maint, (WIN[1]:WIN[2]-WIN[1]:WIN[2])/1000)')
        trueval = calculate_f(rast, qrf, null)[1]
        if isnan(trueval)
            pval = NaN
        else
            permvals = calculate_f(rast[perms], qrf, null)
            pval = mean(x->x >= trueval, permvals)
        end
        Dict(
            :pval => pval
        )
    end
end

"""
    compute_modulation_maintenance()

Computes ω^2 at maintenance, as well as mean responses by stimulus.
"""
function compute_modulation_maintenance()
    WIN = (300, 2400)
    function sessfunc(session::Session)
        trials = session.trials
        nstim = maximum(map(maximum, trials[:stimuli]))
        X = zeros(size(trials, 1), nstim)
        for i = 1:size(trials, 1)
            X[i, trials[i, :stimuli]::Vector{Int}] = 1
        end
        stimuli = trials[:stimuli]
        (trials[:maint].data, stimuli.data, nstim, qrfactchk(X), qrfactchk(ones(size(X, 1), 1)))
    end

    mapcells(sessfunc, :modulation_maintenance, Dict(:win => WIN)) do session, x, spikes, ch
        (maint, stimuli, nstim, X, null) = x
        rast = map(Float64, raster(spikes, maint, (WIN[1]:WIN[2]-WIN[1]:WIN[2])/1000)')
        means = [mean(rast[[i for i = 1:length(rast) if stim in stimuli[i]]]) for stim in 1:nstim]
        Dict(:pve => calculate_omega2(rast, X, null)[1],
             :means => means)
    end
end

"""
    compute_modulation_presentation()

Computes mean responses by stimulus at presentation, as well as baseline firing rate and
standard deviation.
"""
function compute_modulation_presentation()
    @load joinpath(PROCESSED_DATA_DIR, "latency.jld") chs window
    chs = chs::Vector{ClusterMetadata}
    window = window::Vector{Tuple{Int,Int}}

    function sessfunc(session::Session)
        trials = stimulusframe(session.trials)
        @assert sort(unique(trials[:isi])) == ISI
        isi = convert(Vector{Float64}, trials[:isi])
        onset = trials[:onset].data
        stimulus = convert(Vector{Int}, trials[:stimulus])::Vector{Int}
        trial_start_time = [session.trials[:onsets][i][1] for i = 1:size(session.trials, 1)]
        (isi, onset, stimulus, trial_start_time)
    end

    mapcells(sessfunc, :modulation_presentation) do session, x, spikes, ch
        (isi, onset, stimulus, trial_start_time) = x
        chind = findfirst(chs, ch)
        win = window[chind]
        win_rg = (win[1]:win[2]-win[1]:win[2])/1000
        @assert length(win_rg) == 2
        nstim = maximum(stimulus)
        stim_fr_win = zeros(nstim)
        for stim = 1:nstim
            stim_onset = onset[stimulus .== stim]
            stim_fr_win[stim] = mean(raster(spikes, stim_onset, win_rg))
        end
        baseline_fr_all, baseline_std_all = mean_and_std(raster(spikes, trial_start_time, -1:0))
        baseline_fr_win, baseline_std_win = mean_and_std(raster(spikes, trial_start_time, win_rg-1))
        Dict(:stim_fr_win => stim_fr_win,
             :baseline_fr_all => baseline_fr_all,
             :baseline_std_all => baseline_std_all,
             :baseline_fr_win => baseline_fr_win,
             :baseline_std_win => baseline_std_win,
             :win => win)
    end
end

"""
    compute_maintenance_topstim_holdout()

Computes normalized modulation for best stimulus at presentation during maintenance,
determining best stimulus and maintenance period modulation on different subsets of trials.
"""
function compute_maintenance_topstim_holdout()
    WIN = (300, 2400)
    BSWIN = (300, 1000)
    LASTPT = 1000
    @load joinpath(PROCESSED_DATA_DIR, "modulation_presentation.jld") chs baseline_std_all

    function sessfunc(session::Session)
        s1 = 1:4:size(session.trials, 1)
        s2 = [!(x in s1) for x in 1:size(session.trials, 1)]
        strials = stimulusframe(session.trials[s1, :])
        onset = strials[:onset].data
        stimulus = convert(Vector{Int}, strials[:stimulus])::Vector{Int}
        trial_start_time = [x[1] for x in session.trials[s1, :onsets]]
        inds = vcat([find(stimulus .== i) for i = 1:maximum(stimulus)]...)
        ningroup = [sum(stimulus .== i) for i = 1:maximum(stimulus)]

        trials = session.trials[s2, :]
        (onset, stimulus, inds, ningroup, trial_start_time, getindex.(trials[:onsets].data, 1), trials[:maint].data, trials[:stimuli].data)
    end

    mapcells(sessfunc, :maintenance_topstim_holdout, Dict(:win => WIN)) do session, x, spikes, ch
        chind = findfirst(chs, ch)
        chind == 0 && return
        (ponsets, pstimulus, inds, ningroup, pbaseline, mbaseline, maint, stimuli) = x

        rast = raster(spikes, ponsets, (1:LASTPT)/1000)
        eta2_grid = poisson_deviance_grid(rast[:, inds]', ningroup)
        ends, starts = ind2sub(size(eta2_grid), find(eta2_grid .== maximum(eta2_grid)))
        isempty(ends) && return
        best = indmin(ends - starts)
        pwin = (starts[best], ends[best]+1)

        rast = raster(spikes, ponsets, (pwin[1]:pwin[2]-pwin[1]:pwin[2])/1000)
        stimresp = [mean(rast[pstimulus .== i]) for i = 1:maximum(pstimulus)]
        baseline = mean(raster(spikes, pbaseline, (pwin[1]:pwin[2]-pwin[1]:pwin[2])/1000-1))
        beststim = indmax(abs(stimresp - baseline))
        ispositive = indmax(stimresp) == beststim

        rast = raster(spikes, maint, (WIN[1]:WIN[2]-WIN[1]:WIN[2])/1000)'
        g1 = find(x->beststim in x, stimuli)
        g2 = find(x->!(beststim in x), stimuli)
        pval_mwu = pvalue(MannWhitneyUTest(rast[g1], rast[g2]))
        mbsrast = raster(spikes, maint, (BSWIN[1]:BSWIN[2]-BSWIN[1]:BSWIN[2])/1000)'
        bsrast = raster(spikes, [pbaseline; mbaseline], (BSWIN[1]:BSWIN[2]-BSWIN[1]:BSWIN[2])/1000-1)'
        sdest = baseline_std_all[chind]*(WIN[2]-WIN[1])/1000

        mg2 = mean(rast[g2])
        sep = (mean(rast[g1]) - mg2)/sdest*ifelse(ispositive, 1, -1)
        sep_by_pos = [mean(rast[find(x->x[pos] == beststim, stimuli)]) - mg2 for pos = 1:length(stimuli[1])]/sdest*ifelse(ispositive, 1, -1)
        Dict(
            :pwin => pwin,
            :beststim => beststim,
            :ng1 => length(g1),
            :ng2 => length(g2),
            :ispositive => ispositive,
            :pval_mwu => pval_mwu,
            :sep => sep,
            :sep_by_pos => sep_by_pos
        )
    end
end

"""
    compute_maintenance_topstim_correct()

Computes normalized modulation for best stimulus at presentation during maintenance on
correct trials relative to incorrect trials, when the stimulus was probed.
"""
function compute_maintenance_topstim_correct()
    individual_unit_test(rast, g1, g2) =
        isempty(g1) || isempty(g2) ? NaN : pvalue(MannWhitneyUTest(rast[g1], rast[g2]))

    modulation(rast, g1, g2, sdest, ispositive) =
        ifelse(ispositive, 1, -1)*(mean(rast[g1]) - mean(rast[g2]))/sdest

    WIN = (300, 2400)
    EARLY_WIN = (0, 200)
    PERSIST_WIN = (300, 700)

    sel = permselunits()
    @load joinpath(PROCESSED_DATA_DIR, "modulation_presentation.jld") chs stim_fr_win baseline_fr_all baseline_fr_win win baseline_std_all
    stim_chs = chs
    ispositive = [maxabs(fr-bs) == maximum(fr-bs) for (fr, bs) in zip(stim_fr_win, baseline_fr_win)]
    beststim = [indmax(ifelse(ip, 1, -1)*x) for (x, ip) in zip(stim_fr_win, ispositive)]

    function sessfunc(session::Session)
        trials = session.trials
        (trials[:onsets].data, trials[:maint].data, trials[:probes].data, trials[:stimuli].data, trials[:probe1].data, trials[:probe2].data,
         convert(Vector{Float64}, trials[:isi]), convert(Vector{Bool}, trials[:correct]))
    end

    mapcells(sessfunc, :maintenance_topstim_correct, Dict(:win => WIN)) do session, x, spikes, ch
        ch in sel || return
        chind = findfirst(chs, ch)
        chind == 0 && return
        (onsets, maint, probes, stimuli, probe1, probe2, isi, correct) = x
        bs = beststim[chind]
        intrial = [bs in x for x in stimuli]
        baselines = [x[1]-1 for x in onsets]

        rast = vec(raster(spikes, maint, (WIN[1]:WIN[2]-WIN[1]:WIN[2])/1000))
        probed = (probe1 .== bs) | (probe2 .== bs)
        g1 = find(intrial & probed & correct)
        ng1 = length(g1)
        g2 = find(intrial & probed & !correct)
        ng2 = length(g2)
        pval_mwu_maint = individual_unit_test(rast, g1, g2)
        sep_maint = modulation(rast, g1, g2, baseline_std_all[chind]*(WIN[2]-WIN[1])/1000, ispositive[chind])

        pres = [y[findfirst(x, bs)] for (y, x) in zip(onsets[intrial], stimuli[intrial])]
        cwin = win[chind]
        g1 = find(probed[intrial] & correct[intrial])
        g2 = find(probed[intrial] & !correct[intrial])

        rast = vec(raster(spikes, pres, (cwin[1]+EARLY_WIN[1]:EARLY_WIN[2]-EARLY_WIN[1]:cwin[1]+EARLY_WIN[2])/1000))
        pval_mwu_pres = individual_unit_test(rast, g1, g2)
        sep_pres = modulation(rast, g1, g2, baseline_std_all[chind]*(EARLY_WIN[2]-EARLY_WIN[1])/1000, ispositive[chind])

        use = isi[intrial] .>= .5
        rast = vec(raster(spikes, pres, (cwin[1]+PERSIST_WIN[1]:PERSIST_WIN[2]-PERSIST_WIN[1]:cwin[1]+PERSIST_WIN[2])/1000))
        sdest = baseline_std_all[chind]*(PERSIST_WIN[2]-PERSIST_WIN[1])/1000
        g1 = find(probed[intrial] & correct[intrial] & use)
        g2 = find(probed[intrial] & !correct[intrial] & use)
        ng1_persist = length(g1)
        ng2_persist = length(g2)
        pval_mwu_persist = individual_unit_test(rast, g1, g2)
        sep_persist = modulation(rast, g1, g2, baseline_std_all[chind]*(PERSIST_WIN[2]-PERSIST_WIN[1])/1000, ispositive[chind])

        Dict(
            :ng1 => ng1,
            :ng2 => ng2,
            :ng1_persist => ng1_persist,
            :ng2_persist => ng2_persist,
            :ispositive => ispositive[chind],
            :pval_mwu_maint => pval_mwu_maint,
            :pval_mwu_pres => pval_mwu_pres,
            :pval_mwu_persist => pval_mwu_persist,
            :sep_maint => sep_maint,
            :sep_persist => sep_persist,
            :sep_pres => sep_pres
        )
    end
end

"""
    bs95ci(X)

Bootstrap 95% confidence interval; trials along first dimension of X
"""
function bs95ci(X::AbstractArray, bs=rand(1:size(X, 1), size(X, 1), 1000))
    @assert size(X, 1) == size(bs, 1)
    tmp = zeros(size(bs, 2))
    ci = zeros(2, size(X, 2))
    for ip = 1:size(X, 2)
        for i = 1:size(bs, 2)
            mu = 0.0
            @inbounds @simd for j = 1:size(bs, 1)
                mu += X[bs[j, i], ip]
            end
            mu /= size(X, 1)
            tmp[i] = mu
        end
        ci[:, ip] = quantile!(tmp, [0.025, 0.975])
    end
    ci
end

"""
    bsvs0(X)

Single-sample bootstrap test that X is > 0
"""
function bsvs0(X::AbstractVector, nbs=10000)
    n = 0
    gen = Base.Random.RangeGenerator(1:size(X, 1))
    generator = Base.Random.GLOBAL_RNG
    for i = 1:nbs
        mu = 0.0
        @inbounds for j = 1:size(X, 1)
            mu += X[rand(generator, gen)]
        end
        n += mu <= 0
    end
    n/nbs
end

"""
    permselunits(X)::Vector{ClusterMetadata}

ClusterMetadata for units significant at presentation.
"""
function permselunits(; alpha=0.001)
    @load joinpath(PROCESSED_DATA_DIR, "permpval.jld") pval chs
    chsel = find(x->(x.unit_type == :MU || x.unit_type == :SU), chs)
    selchs = chs[chsel][pval[chsel] .< alpha]
end

"""
    allisisigunits(X)::Vector{ClusterMetadata}

ClusterMetadata for units significant at presentation separately at each ISI.
"""
function allisisigunits(; alpha=0.05)
    @load joinpath(PROCESSED_DATA_DIR, "permpval_eachisi.jld") chs pval
    chsel = find(x->(x.unit_type == :MU || x.unit_type == :SU), chs)
    chs[chsel[[all(x->x < alpha, x) for x in pval[chsel]]]]
end

"""
    maintselunits(X)::Vector{ClusterMetadata}

ClusterMetadata for units significant at maintenance.
"""
function maintselunits(; alpha=0.01)
    @load joinpath(PROCESSED_DATA_DIR, "permpval_maintenance.jld") pval chs
    chsel = find(x->(x.unit_type == :MU || x.unit_type == :SU), chs)
    selchs = chs[chsel][pval[chsel] .< alpha]
end
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
using Combinatorics
using DataFrames
export Walsh, k_bit_epistasis, total_epistasis

# See the functions k_bit_epistasis and total_epistasis for Walsh transform computation

# count_ones(x) is a built-in Julia function with the same purpose
function count1bits( x::MyInt, numinputs::Int64 )
  shift = (2^numinputs-1)
  mask = MyInt(1) << shift
  ssum = 0
  while mask != MyInt(0)
    ssum += (mask & x) >> shift
    #println("mask: ",mask,"  shift: ",shift,"  ssum: ",ssum)
    shift -= 1
    mask >>= 1
  end
  ssum
end

function my_ones( output::MyInt, numinputs::Int64 )
  mask = One()
  result = [ count1bits( output & MyInt(m), numinputs )/count1bits( MyInt(m), numinputs ) for m = mask:-1:1 ]
  push!(result,0.0)
  result
end

function my_ones( output::MyInt )
  mask = One()
  result = [ count_ones( output & MyInt(m) )/count_ones( MyInt(m) ) for m = mask:-1:1 ]
  push!(result,0.0)
  result
end

function my_mask( position::Int64 )
  one(MyInt) << (position-1)
end

function my_mask( positions::Vector{Int64} )
  result = zero(MyInt)
  for p in positions
    result |= my_mask( p )
  end
  result
end

function my_size()
  i = 1
  while my_mask(i) != zero(MyInt)
    i+=1
  end
  i-1
end

function my_size( numinputs::Int64 )
  m = One()
  i = 1
  while (my_mask(i) & m) != zero(MyInt)
    i+=1
  end
  i-1
end

function epistasis( x::Integer, p::Vector{Int64}, numinputs::Int64 )
    #println("p: ",p,"  ", [[s for s in combinations(p,k)]  for k = 1:order ] )
    x = MyInt(x)
    summand = 0
    count = 0
    for q in combinations(p,length(p)-1)
      #println("  q: ",q,"  perm_output( x, p, numinputs ): ",perm_output( x, q, numinputs ))
      summand += perm_output( x, q, numinputs )
      count += 1
    end
    #println("perm_output(x,p,numinputs): ",perm_output(x,p,numinputs),"  summand/count: ",summand/count )
    abs(perm_output(x,p,numinputs) - summand/count)
end

function epistasis( x::Integer, order::Int64, numinputs::Int64 )
  x = MyInt(x)
  ssum = 0
  for p in combinations( collect(1:numinputs), order )
    summand = epistasis( x, p, numinputs )
    ssum += summand
    #println("p: ",p,"  epistasis(p): ", summand )
  end
  ssum
end

#[ (p, my_mask(Int64(my_mask( p ))), output & my_mask(Int64(my_mask( p )))) for p in combinations(collect(1:8),2) ]

function perm_output( output::MyInt, p::Vector{Int64}, numinputs::Int64 )
  context = construct_context(numinputs)
  (output & my_mask(get_bits(context,numinputs)[my_mask(p)]+1)) != zero(MyInt) ? 1 : -1
end

function Walsh( L::Int64 )
  denom = 1.0/2.0^L
  W = zeros(Int64,L,L)
  for i = 0:(L-1)
    for j = 0:(L-1)
      W[i+1,j+1] = iseven(count_ones( i & j )) ? 1 : -1
    end
  end
  W
end

# Not correct
function Wtrans( L::Int64, x::Integer )
  result = zeros(Int64,L)
  xb = to_binary(x,L)
  for i = 0:(L-1)
    result[i+1] += xb[i+1]*(iseven( count_ones( i & xb[i+1] )) ? 1 : -1 )
  end
  result
end

# returns a list of the 1-based indices i such that count_ones(i-1) == k
function k_bit_indices( L::Int64, k::Int64 )
  findall( x->count_ones(x)==k, 0:(L-1))
end

function k_bit_epistasis( W::DenseMatrix{Int64}, k::Int64, x::MyInt )
  L = size(W)[1]
  k_inds = k_bit_indices( L, k )
  wx = W*to_binary(x,L)
  #println("k_inds: ",k_inds,"  wx: ",wx)
  sum( abs(wx[i]) for i in k_inds )
end

function total_epistasis( W::DenseMatrix{Int64}, x::MyInt )
  L = size(W)[1]
  wx = W*to_binary(x,L)
  sum( abs(wx[i]) for i = 2:L )
end

function reciprocal_sign_epistasis()
end

function add_epistasis_to_dataframe( df::DataFrame, W::DenseMatrix{Int64} )
  numinputs = Int(df[1,:numinputs])
  all_epistasis = zeros(Int64,size(df)[1],numinputs+1)
  for i = 1:size(df)[1]
    goal = parse(MyInt,df[i,:goal])
    for k = 1:numinputs
      all_epistasis[i,k] = k_bit_epistasis( W, k, goal )
    end
    all_epistasis[i,numinputs+1] = total_epistasis( W, goal )
  end
  #all_epistasis
  field_names = vcat([ Symbol("epistasis$k") for k = 1:numinputs ], [Symbol("total_epistasis")] )
  for k = 1:numinputs
    df[!,Symbol("epistasis$k")] = all_epistasis[:,k]
  end
  df[!,Symbol("total_epistasis")] = all_epistasis[:,numinputs+1]
  df
end

function testw( W, x )
  L=size(W)[1]
  #println( [ k_bit_epistasis(W,k,x) for k = 1:Int(log2(L)) ])
  @assert sum( k_bit_epistasis(W,k,x) for k = 1:Int(log2(L))) == total_epistasis( W,x)
end

function testw( W )
  L = size(W)[1]
  for x = zero(MyInt):MyInt(Int(log2(L)))-MyInt(1)
    testw(W,x)
  end
end

ext: jl


text:
struct Solver{S,T,P,PS}
    strength::S
    CF::T
    presmoother::P
    postsmoother::PS
    max_levels::Int
    max_coarse::Int
end

function ruge_stuben(_A::Union{TA, Symmetric{Ti, TA}, Hermitian{Ti, TA}},
                ::Type{Val{bs}}=Val{1};
                strength = Classical(0.25),
                CF = RS(),
                presmoother = GaussSeidel(),
                postsmoother = GaussSeidel(),
                max_levels = 10,
                max_coarse = 10,
                coarse_solver = Pinv) where {Ti,Tv,bs,TA<:SparseMatrixCSC{Ti,Tv}}

    s = Solver(strength, CF, presmoother,
                postsmoother, max_levels, max_levels)

    if _A isa Symmetric && Ti <: Real || _A isa Hermitian
        A = _A.data
        At = A
        symmetric = true
        @static if VERSION < v"0.7-"
            levels = Vector{Level{TA, TA}}()
        else
            levels = Vector{Level{TA, Adjoint{Ti, TA}, TA}}()
        end
    else
        symmetric = false
        A = _A
        At = adjoint(A)
        @static if VERSION < v"0.7-"
            levels = Vector{Level{TA, TA, TA}}()
        else
            levels = Vector{Level{TA, Adjoint{Ti, TA}, TA}}()
        end
    end
    w = MultiLevelWorkspace(Val{bs}, eltype(A))
    residual!(w, size(A, 1))

    while length(levels) + 1 < max_levels && size(A, 1) > max_coarse
        A = extend_heirarchy!(levels, strength, CF, A, symmetric)
        coarse_x!(w, size(A, 1))
        coarse_b!(w, size(A, 1))
        residual!(w, size(A, 1))
    end

    MultiLevel(levels, A, coarse_solver(A), presmoother, postsmoother, w)
end

function extend_heirarchy!(levels, strength, CF, A::SparseMatrixCSC{Ti,Tv}, symmetric) where {Ti,Tv}
    if symmetric
        At = A
    else
        At = adjoint(A)
    end
    S, T = strength(At)
    splitting = CF(S)
    P, R = direct_interpolation(At, T, splitting)
    push!(levels, Level(A, P, R))
    A = R * A * P
end

function direct_interpolation(At, T, splitting)
    fill!(T.nzval, eltype(At)(1))
    T .= At .* T

    Pp = rs_direct_interpolation_pass1(T, splitting)
    Px, Pj, Pp = rs_direct_interpolation_pass2(At, T, splitting, Pp)
    R = SparseMatrixCSC(maximum(Pj), size(At, 1), Pp, Pj, Px)
    P = R'

    P, R
end

# calculates the number of nonzeros in each column of the interpolation matrix
function rs_direct_interpolation_pass1(T, splitting)
    n = size(T, 2)
    Bp = ones(Int, n+1)
    nnzplus1 = 1
    for i = 1:n
        if splitting[i] == C_NODE
            nnzplus1 += 1
        else
            for j in nzrange(T, i)
                row = T.rowval[j]
                if splitting[row] == C_NODE
                    nnzplus1 += 1
                end
            end
        end
        Bp[i+1] = nnzplus1
    end
    Bp
 end

function rs_direct_interpolation_pass2(At::SparseMatrixCSC{Tv,Ti},
                                                T::SparseMatrixCSC{Tv, Ti},
                                                splitting::Vector{Ti},
                                                Bp::Vector{Ti}) where {Tv,Ti}

    Bx = zeros(Tv, Bp[end] - 1)
    Bj = zeros(Ti, Bp[end] - 1)

    n = size(At, 1)

    for i = 1:n
        if splitting[i] == C_NODE
            Bj[Bp[i]] = i
            Bx[Bp[i]] = 1
        else
            sum_strong_pos = zero(Tv)
            sum_strong_neg = zero(Tv)
            for j in nzrange(T, i)
                row = T.rowval[j]
                sval = T.nzval[j]
                if splitting[row] == C_NODE
                    if sval < 0
                        sum_strong_neg += sval
                    else
                        sum_strong_pos += sval
                    end
                end
            end
            sum_all_pos = zero(Tv)
            sum_all_neg = zero(Tv)
            diag = zero(Tv)
            for j in nzrange(At, i)
                row = At.rowval[j]
                aval = At.nzval[j]
                if row == i
                    diag += aval
                else
                    if aval < 0
                        sum_all_neg += aval
                    else
                        sum_all_pos += aval
                    end
                end
            end

            if sum_strong_pos == 0
                beta = zero(diag)
                if diag >= 0
                    diag += sum_all_pos
                end
            else
                beta = sum_all_pos / sum_strong_pos
            end

            if sum_strong_neg == 0
                alpha = zero(diag)
                if diag < 0
                    diag += sum_all_neg
                end
            else
                alpha = sum_all_neg / sum_strong_neg
            end

            if isapprox(diag, 0, atol=eps(Tv))
                neg_coeff = Tv(0)
                pos_coeff = Tv(0)
            else
                neg_coeff = alpha / diag
                pos_coeff = beta / diag
            end

            nnz = Bp[i]
            for j in nzrange(T, i)
                row = T.rowval[j]
                sval = T.nzval[j]
                if splitting[row] == C_NODE
                    Bj[nnz] = row
                    if sval < 0
                        Bx[nnz] = abs(neg_coeff * sval)
                    else
                        Bx[nnz] = abs(pos_coeff * sval)
                    end
                    nnz += 1
                end
            end
        end
    end

    m = zeros(Ti, n)
    sum = zero(Ti)
    for i = 1:n
        m[i] = sum
        sum += splitting[i]
    end
    Bj .= m[Bj] .+ 1

    Bx, Bj, Bp
end


text:
using Test
using EMIRT
using EMIRT.IOs
@testset "test IO" begin

fileName = "$(tempname()).h5"
# test image IO
img = EMImage(rand(UInt8, 300,300,10))
saveimg(fileName, img)
@test img == readimg(fileName)
rm(fileName)

# test segmentation IO
seg = Segmentation(rand(UInt32, 300,300,10))
saveseg(fileName, seg)
@test seg == readseg(fileName)
rm(fileName)

# test affinity IO
aff = AffinityMap(rand(Float32, 300,300,30,3))
saveaff(fileName, aff)
@test aff == readaff(fileName)
rm(fileName)

# test sgm IO
sgm = EMIRT.Segmentations.seg2sgm(seg)
savesgm(fileName, sgm)
sgm2 = readsgm(fileName)
# Note! we can only compare by internal field
# julia == function only compare the memory address for mutable objects.
# see https://github.com/JuliaLang/julia/issues/5340
@test sgm2.segmentation == sgm.segmentation
@test sgm2.segmentPairs == sgm.segmentPairs
@test sgm2.segmentPairAffinities == sgm.segmentPairAffinities
rm(fileName)

end # end of IO test set


text:
"""
    FiniteDifferenceGrids

A simple 1-dimensional uniform grid
for finite difference method.
"""
module FiniteDifferenceGrids

export Grid, UniformGrid, over_elems, over_elems_real, over_elems_ghost
export ZBoundary, Zmin, Zmax, n_hat, binary, ghost_vec, ghost_dual
export first_interior, boundary, second_interior, first_ghost, boundary_points

abstract type Grid{FT <: AbstractFloat, IT <: Integer} end

"""
    UniformGrid{T, I}

A simple 1-dimensional uniform grid of
type `T` for finite difference method.
"""
struct UniformGrid{FT, IT} <: Grid{FT, IT}
    """ Number of physical elements """
    n_elem::IT

    """ Number of real elements """
    n_elem_real::IT

    """ Number of ghost points per side """
    n_ghost::IT

    """ Element size """
    Δz::FT

    """ Inverse of element size squared """
    Δzi::FT

    """ Inverse of element size squared """
    Δzi2::FT

    """ z-coordinate at cell edges """
    ze::Vector{FT}

    """ z-coordinate at cell centers """
    zc::Vector{FT}

    """ z-coordinate at bottom of boundary node """
    zn_min::FT

    """ z-coordinate at top of boundary node """
    zn_max::FT
end

"""
    UniformGrid(z_min::T, z_max::T, n_elem_real::Int, n_ghost::Int = 1)

A simple grid implementation that accepts the domain interval
`z_min` and `z_max`, the number of elements `n_elem_real` and
the number of ghost points per side `n_ghost`.
"""
function UniformGrid(
    z_min::T,
    z_max::T,
    n_elem_real::I,
    n_ghost::I = 1,
) where {T, I}
    n_elem = n_elem_real + 2 * n_ghost
    Δz = (z_max - z_min) / n_elem_real
    Δzi = inv(Δz)
    Δzi2 = Δzi^2
    zc = [(i - n_ghost) * Δz - Δz / 2 for i in 1:n_elem]
    ze = [(i - n_ghost) * Δz - Δz for i in 1:(n_elem + 1)]
    zn_min = (zc[n_ghost] + zc[1 + n_ghost]) / 2
    zn_max = (zc[end - n_ghost] + zc[end - n_ghost + 1]) / 2
    return UniformGrid{T, I}(
        n_elem,
        n_elem_real,
        n_ghost,
        Δz,
        Δzi,
        Δzi2,
        ze,
        zc,
        zn_min,
        zn_max,
    )
end

function Base.show(io::IO, grid::Grid)
    println(io, "----------------- Grid")
    println(io, "zc      = ", grid.zc)
    println(io, "ze      = ", grid.ze)
    println(io, "n_elem  = ", grid.n_elem)
    println(io, "n_ghost = ", grid.n_ghost)
    println(io, "zn_min  = ", grid.zn_min)
    println(io, "zn_max  = ", grid.zn_max)
    println(io, "-----------------")
end

abstract type ZBoundary end

"""
    Zmax
Type to dispatch function to particular location
"""
struct Zmax <: ZBoundary end

"""
    Zmin
Type to dispatch function to particular location
"""
struct Zmin <: ZBoundary end

"""
    n_hat(::ZBoundary)

The outward normal vector to the boundary
"""
n_hat(::Zmin) = -1
n_hat(::Zmax) = 1

"""
    binary(::ZBoundary)

Returns 0 for Zmin and 1 for Zmax
"""
binary(::Zmin) = 0
binary(::Zmax) = 1

"""
    ghost_vec(::ZBoundary)

A 3-element vector near the boundary with 1's on ghost cells and 0's on interior cells
"""
ghost_vec(::Zmin) = [1, 0, 0]
ghost_vec(::Zmax) = [0, 0, 1]

"""
    ghost_dual(::ZBoundary)

A 2-element vector near the boundary with 1's on ghost cells and 0's on interior cells
"""
ghost_dual(::Zmin) = [1, 0]
ghost_dual(::Zmax) = [0, 1]

"""
    over_elems(grid::Grid)

Get the range of indexes to traverse real and ghost grid elements
"""
@inline over_elems(grid::Grid) = 1:(grid.n_elem)

"""
    over_elems_real(grid::Grid)

Get the range of indexes to traverse only real grid elements
"""
@inline over_elems_real(grid::Grid) =
    (1 + grid.n_ghost):(grid.n_elem - grid.n_ghost)

"""
    first_interior(grid::Grid, ::ZBoundary)

Get the first element index in the physical domain
"""
@inline first_interior(grid::Grid, ::Zmax) = grid.n_elem - grid.n_ghost
@inline first_interior(grid::Grid, ::Zmin) = 1 + grid.n_ghost

"""
    boundary(grid::Grid, ::ZBoundary)

Get the element index for edge data at the boundary
"""
@inline boundary(grid::Grid, ::Zmax) = grid.n_elem - grid.n_ghost
@inline boundary(grid::Grid, ::Zmin) = 1 + grid.n_ghost

"""
    second_interior(grid::Grid, ::ZBoundary)

Get the second element index in the physical domain
"""
@inline second_interior(grid::Grid, ::Zmax) = first_interior(grid, Zmax()) - 1
@inline second_interior(grid::Grid, ::Zmin) = first_interior(grid, Zmin()) + 1

"""
    first_ghost(grid::Grid, ::ZBoundary)

Get the first element index in the ghost domain
"""
@inline first_ghost(grid::Grid, ::Zmax) = first_interior(grid, Zmax()) + 1
@inline first_ghost(grid::Grid, ::Zmin) = first_interior(grid, Zmin()) - 1

"""
    boundary_points(grid::Grid, ::ZBoundary)

Get the first element index in the ghost domain
"""
@inline boundary_points(grid::Grid, b::ZBoundary) =
    first_ghost(grid, b), first_interior(grid, b), second_interior(grid, b)

"""
    over_elems_ghost(grid::Grid)

Get the range of indexes to traverse only ghost grid elements
"""
@inline over_elems_ghost(grid::Grid) =
    setdiff(over_elems(grid), over_elems_real(grid))

"""
    over_elems_ghost(grid::Grid, ::ZBoundary)

Get the range of indexes to traverse only ghost grid elements below ``z_{min}`` and above ``z_{max}``.
"""
@inline over_elems_ghost(grid::Grid, ::Zmin) = 1:(grid.n_ghost)
@inline over_elems_ghost(grid::Grid, ::Zmax) =
    (grid.n_elem - grid.n_ghost + 1):(grid.n_elem)

Base.eltype(::Grid{T}) where {T} = T

include("discrete_operators.jl")

end


text:
normal_kurtosis(x) = kurtosis(x, false)

function ks_test(a,b)
    ecdf_a = ecdf(a[:, 1])
    ecdf_b = ecdf(b[:, 1])
    all_x = vcat(a[:,1],b[:,1])
    return maximum(abs.(ecdf_a(all_x) .- ecdf_b(all_x)))
end

acf_lag_first(x) = autocor(x, [1])[1]

function generalized_hurts_exp(obs)
    q = 1
    max_τ = 19

    k = 0
    H = zeros(15)
    L = size(obs, 1)

    for iter_τ = 5:max_τ
        k += 1

        x = collect(1:iter_τ)
        k_q_t = zeros(iter_τ)
        for τ = 1:iter_τ
            numer = obs[(τ+1):τ:L] - obs[1:τ:(L-τ)]
            denom = obs[1:τ:L]

            # Determine Drift
            N = size(numer, 1) + 1
            X = collect(1:N)
            Y = denom
            mx = sum(X)/N
            SSxx = sum(X.^2) - N * mx^2
            my = sum(Y) / N
            SSxy = sum(X .* Y) - N * mx * my
            cc = [SSxy/SSxx, my - SSxy/SSxx*mx]

            # Subtract Drift
            numer = numer .- cc[1]
            denom = denom .- cc[1] * collect(1:N) .- cc[2]

            k_q_t[τ] = mean(abs.(numer).^q) / mean(abs.(denom).^q)
        end

        # Calculate Hurst Exponent for current iteration
        log_10_x = log10.(x)
        mx = mean(log_10_x)
        SSxx = sum(log_10_x.^2) - iter_τ * mx^2
        log_10_k_q_t = log10.(k_q_t)
        my = mean(log_10_k_q_t)
        SSxy = sum(log_10_x .* log_10_k_q_t) - iter_τ * mx * my
        H[k] = SSxy/SSxx
    end

    H = mean(H) / q
    return H
end

function get_summary_stats(simulated_obs, obs)

    num_replications = size(simulated_obs,2)
    stats_mat = zeros(num_replications,6)

    stats_mat[:,1] = mean(simulated_obs, dims=1)
    stats_mat[:,2] = std(simulated_obs, dims=1)
    stats_mat[:,3] = mapslices(normal_kurtosis, simulated_obs, dims=1)

    ks_test_stat(x) = ks_test(obs, x)
    stats_mat[:,4] = mapslices(ks_test_stat, simulated_obs, dims=1)
    stats_mat[:,5] = mapslices(generalized_hurts_exp, simulated_obs,
        dims=1)
    stats_mat[:, 6] = mapslices(acf_lag_first, simulated_obs, dims=1)
    return stats_mat
end


text:
import SmoothingSplines: fit,predict

function idx_cycle(voltage)

    _voltage = voltage .+1000
    if _voltage[1] > _voltage[10]
        _idx=findnext(x->x>= _voltage[1],_voltage,2)
        idx=findnext(x->x<= _voltage[1],_voltage,_idx+1)

    elseif _voltage[1] < _voltage[10]
        _idx=findnext(x->x<= _voltage[1],_voltage,2)
        idx=findnext(x->x>= _voltage[1],_voltage,_idx+1)
    end
    return idx
end

function cycle(voltage,cycle::Int64)

    idx = idx_cycle(voltage)
    total_cycle = round(length(voltage)/idx) |> Int

    if cycle > total_cycle
        error("There are only $(total_cycle) full cycles in your Array.")
    else
        _start = 1+(cycle-1)*idx
        _end = cycle*idx
        return voltage[_start:_end]
    end
end

"""
cycle(voltage::Array{Any,1},current::Array{Any,1},cycle::Int64)

Return the Input voltage and current Arrays for the given cycle argument.
"""
function cycle(voltage,current,cycle::Int64)

    idx = idx_cycle(voltage)
    total_cycle = round(length(voltage)/idx) |> Int

    if cycle > total_cycle
        error("There are only $(total_cycle) full cycles in your Array.")
    else
        _start = 1+(cycle-1)*idx
        _end   = cycle*idx

        return voltage[_start:_end],current[_start:_end]

    end
end

"""
total_cycles(voltage::Array{Number,1})

Returns the total number of full cycles in voltage.
"""
function total_cycles(voltage)

    idx = idx_cycle(voltage)
    total_cycle = round(length(voltage)/idx) |> Int

    return total_cycle
end

"""Returns the anodic and cathodic current at the given potential value (val).Works only if c and v consist of only one cycle. For Value you can put any potential in the scan range.
get_eq_anodickathodiccurrents(v,c,val)

julia> get_eq_anodickathodiccurrents(v_cycle1[1],c_cycle1[1],-0.03)
    (-1.3481340532268844, 0.9398593965871291)
"""
function get_eq_anodickathodiccurrents(v,c,val)
    idx_1,idx_2 = get_ind_eq_anodickathodicpotentials(v,c,val)

    val_p1 = v[idx_1]
    val_p2 = v[idx_2]
    val_p3 = v[idx_2+1]

    inc_idx = abs(val_p1-val_p2 )/abs(val_p3-val_p2)

    itp = interpolate(c, BSpline(Quadratic(Reflect(OnCell()))))

    if val_p1 < val_p2

        return c[idx_1],itp(idx_2+inc_idx)

    elseif val_p1 > val_p2

        return c[idx_1],itp(idx_2+(1-inc_idx))
    end
end

"""Returns the anodic and cathodic current indices at the given potential value (val).Works only if c and v consist of only one cycle. For Value you can put any potential in the scan range.
get_ind_eq_anodickathodicpotentials(v,c,val)

julia> get_eq_anodickathodicpotentials(v_cycle1[1],c_cycle1[1],-0.03)
    (42,97)
"""
function get_ind_eq_anodickathodicpotentials(v,c,val)

    p=sortperm(v)
    v_sorted = v[p]

    _v = v_sorted .- val

    idx_1=findnext(x->x>= 0,_v,1)

    idx_2=findpotential(v,p,idx_1,c)

    idxs= [p[idx_1],p[idx_2]]

    return idxs

end

function findpotential(v,p,idx,c;sample_rate=10)
    val = v[p[idx]]
    pot_halfcycle= abs(maximum(v) - minimum(v))

    time4halfcycle=  abs((argmin(v) - argmax(v))/sample_rate)

    scan_rate = pot_halfcycle/time4halfcycle

    if c[p[idx]] > 0
        current= "anodic"
    else
        current= "cathodic"
    end

    idx_2= 0

    for i in idx-10:idx+10
        if current == "anodic"
            if abs(v[p][i] - val) < (scan_rate/sample_rate) && c[p[i]] < 0 && c[p[i]] < 0
                idx_2= i
                return idx_2
            end
        elseif current == "cathodic"
            if abs(v[p][i] - val) < (scan_rate/sample_rate) && c[p[i]] > 0
                idx_2= i
                return idx_2
            end
        end
    end

end

max_stars_repo_name: lambcount/EC

max_issues_repo_name: lambcount/EC

max_forks_repo_name: lambcount/EC


text:
"""
$(SIGNATURES)

BitBuffer to buffer bits
"""
struct BitBuffer
    buffer::UInt64
    length::Int
end

function BitBuffer()
    BitBuffer(0, 0)
end

@inline get_bits(bit_buffer::BitBuffer) = bit_buffer.buffer
@inline length(bit_buffer::BitBuffer) = bit_buffer.length

"""
$(SIGNATURES)

Buffer data bits based on the prompt accumulation and the current prompt value.
"""
function buffer(
    system::AbstractGNSS,
    bit_buffer,
    prompt_accumulator,
    secondary_code_or_bit_found,
    prev_code_phase,
    code_phase,
    integration_time,
    prompt_correlator
)
    prompt_accumulator = prompt_accumulator + secondary_code_or_bit_found *
        prompt_correlator

    if secondary_code_or_bit_found &&
        (code_phase - prev_code_phase < 0 || integration_time == 1 / get_data_frequency(system))
        bit = real(prompt_accumulator) > 0
        bit_buffer = BitBuffer(
            get_bits(bit_buffer) << 1 + UInt64(bit),
            length(bit_buffer) + 1
        )
        prompt_accumulator = zero(prompt_accumulator)
    end
    bit_buffer, prompt_accumulator
end


text:
"Simple start-end gradient structure."
struct Gradient{T<:Colorant}
    start::T
    stop::T
end

"Lineary interpolates between two colors of a `Gradient`, `x` between `xmin`
and `xmax`."
function grad(gradient::Gradient{T}, x::Real, xmin::Real=0, xmax::Real=1) where T
    mapc(
        (start, stop) -> clamp(interp(x, xmin, xmax, start, stop), 0, 1),
        gradient.start, gradient.stop
    )
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

# Julia does not allow identifiers starting with a digit, so we can't
# call this just "4ti2"
name = "lib4ti2"
version = v"1.6.9"

# Collection of sources required to build 4ti2
sources = [
    ArchiveSource("https://github.com/4ti2/4ti2/releases/download/Release_$(version.major)_$(version.minor)_$(version.patch)/4ti2-$(version).tar.gz",
                  "3053e7467b5585ad852f6a56e78e28352653943e7249ad5e5174d4744d174966"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/4ti2-*

# Remove misleading libtool files
rm -f ${prefix}/lib/*.la
rm -f /opt/${target}/${target}/lib*/*.la
rm -f /opt/${MACHTYPE}/${MACHTYPE}/lib*/*.la

# Patch #1 for fixing cross-compilation: The correctness of the patch
# relies on us using clang or GCC (in a new enough version) as compiler;
# otherwise we ought to turn off ftrapv support to be on the safe side.
atomic_patch -p1 ../patches/ftrapv.patch

# Patch #2 for fixing cross-compilation: The code detecting presence and
# usability of GMP is broken for cross compilation; fix that
atomic_patch -p1 ../patches/gmp.patch

# Patch to fix compilation on mingw32: add missing #include <time.h>
atomic_patch -p1 ../patches/time.patch

./configure \
    --prefix=$prefix \
    --build=${MACHTYPE} \
    --host=$target \
    --with-gcc-arch=pentium4 \
    --with-gmp=${prefix} \
    --with-glpk=${prefix} \
    --enable-shared \
    --disable-static
make -j${nproc}
make install

rm -rf ${WORKSPACE}/destdir/${target}

# On Windows, make sure non-versioned filename exists...
if [[ ${target} == *mingw* ]]; then
    cp -v ${prefix}/bin/lib4ti2common-*.dll ${prefix}/bin/lib4ti2common.dll
    cp -v ${prefix}/bin/lib4ti2gmp-*.dll ${prefix}/bin/lib4ti2gmp.dll
    cp -v ${prefix}/bin/lib4ti2int32-*.dll ${prefix}/bin/lib4ti2int32.dll
    cp -v ${prefix}/bin/lib4ti2int64-*.dll ${prefix}/bin/lib4ti2int64.dll
    cp -v ${prefix}/bin/lib4ti2util-*.dll ${prefix}/bin/lib4ti2util.dll
    cp -v ${prefix}/bin/libzsolve-*.dll ${prefix}/bin/libzsolve.dll
fi
"""

# Build for all platforms
platforms = supported_platforms()

# 4ti2 contains std::string values; to avoid incompatibilities across
# the GCC 4/5 version boundary, we need the following:
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
   LibraryProduct("lib4ti2gmp", :lib4ti2),
   LibraryProduct("lib4ti2int32", :lib4ti2int32),
   LibraryProduct("lib4ti2int64", :lib4ti2int64),
   LibraryProduct("libzsolve", :libzsolve),

   # Julia identifiers can't start with a digit, so we have to add
   # a prefix to the symbol used to refer to some of the executables.
   ExecutableProduct("4ti2gmp", :exe4ti2gmp),
   ExecutableProduct("4ti2int32", :exe4ti2int32),
   ExecutableProduct("4ti2int64", :exe4ti2int64),
   ExecutableProduct("genmodel", :genmodel),
   ExecutableProduct("gensymm", :gensymm),
   ExecutableProduct("output", :output),
   ExecutableProduct("ppi", :ppi),
   ExecutableProduct("zsolve", :zsolve),

   # The following are shell scripts; we can't declare them as
   # ExecutableProduct because that would break on Windows (where only .exe
   # files are accepted as executables). But people still are interested in
   # having these scripts, it seems, at least on POSIX systems.
   FileProduct("bin/circuits", :circuits),
   FileProduct("bin/graver", :graver),
   FileProduct("bin/groebner", :groebner),
   FileProduct("bin/hilbert", :hilbert),
   FileProduct("bin/markov", :markov),
   FileProduct("bin/minimize", :minimize),
   FileProduct("bin/normalform", :normalform),
   FileProduct("bin/qsolve", :qsolve),
   FileProduct("bin/rays", :rays),
   FileProduct("bin/walk", :walk),
   FileProduct("bin/zbasis", :zbasis),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("GMP_jll", v"6.1.2"),
    Dependency("GLPK_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using ExprTools

@testset "parsing" begin

    # kwargs
    let kwargs = splitdef(:( foo(a; b=2, c::Int=1, d::Int, e) = a + b + c))[:kwargs]
        @test Pretend.arg_names(kwargs) == [:b, :c, :d, :e]
    end

end


text:
include("middleware/Static.jl")
include("middleware/Rendering.jl")
include("middleware/Storage.jl")
include("middleware/Filter.jl")
include("middleware/Logging.jl")
include("middleware/Security.jl")


text:
using Test
using JuMP
using Polyhedra

@testset "Model to Polyhedra with default library" begin
    @testset "1 variable" begin
        model = Model()
        @variable(model, x >= -1, Int)
        @constraint(model, x ≤ 1)
        p = polyhedron(model)
        @test p isa Interval{Float64}
        @test nhalfspaces(p) == 2
    end
    @testset "2 variables" begin
        model = Model()
        @variable(model, x[1:2] >= -1, integer = true)
        @constraint(model, sum(x) ≤ 1)
        p = polyhedron(model)
        @test p isa DefaultPolyhedron{Float64}
        @test nhalfspaces(p) == 3
    end
    # TODO add test with binary variables once https://github.com/JuliaOpt/MathOptInterface.jl/issues/704 is done.
end


text:
lines = readlines(ARGS[1])

function mode(mat)
    mode_array = []
    for col in mat
        val_dict = Dict()
        for i in col
            if i in keys(val_dict)
                val_dict[i] += 1
            else
                val_dict[i] = 1
            end
        end
        max_count = nothing
        max_val = nothing
        for i in keys(val_dict)
            if isnothing(max_count)
                max_count = val_dict[i]
                max_val = i
            elseif val_dict[i] > max_count
                max_count = val_dict[i]
                max_val = i
            end
        end
        append!(mode_array, parse(Int, max_val))
    end
    return mode_array
end

function flip_bit(val_array)
    eps_array = []
    for i in val_array
        if i == 0
            append!(eps_array, 1)
        elseif i == 1
            append!(eps_array, 0)
        end
    end
    return eps_array
end

function bin_array_to_int(bin_array)
    output = 0
    for i = length(bin_array):-1:1
        output += bin_array[i] * 2^(length(bin_array) - i)
    end
    return output
end

function string_array_to_int(bin_array)
    output = 0
    for i = length(bin_array):-1:1
        output += parse(Int, bin_array[i]) * 2^(length(bin_array) - i)
    end
    return output
end

function calculate_product(val_array)
    eps_array = flip_bit(val_array)
    gamma = bin_array_to_int(val_array)
    epsilon = bin_array_to_int(eps_array)
    return epsilon * gamma
end

function valmatrix(lines)
    line_length = nothing
    val_matrix = []
    for line in lines
        if isnothing(line_length)
            line_length = length(line)
            for i = 1:line_length
                append!(val_matrix, [[]])
            end
        end
        for i = 1:line_length
            append!(val_matrix[i], line[i])
        end

    end
    return val_matrix
end

function solution1(lines)
    val_matrix = valmatrix(lines)
    return calculate_product(mode(val_matrix))
end

function solution2(lines)
    val_matrix = valmatrix(lines)
    co2_rating = nothing
    o2_rating = nothing
    co2_rows = range(1, length(val_matrix[1]))
    o2_rows = range(1, length(val_matrix[1]))
    for i = 1:12

        count_0 = length([val_matrix[i][j] for j in co2_rows if val_matrix[i][j] == '0'])
        count_1 = length([val_matrix[i][j] for j in co2_rows if val_matrix[i][j] == '1'])
        if count_0 > count_1
            co2_rows = [j for j in co2_rows if val_matrix[i][j] == '1']
        elseif count_1 > count_0
            co2_rows = [j for j in co2_rows if val_matrix[i][j] == '0']
        else
            co2_rows = [j for j in co2_rows if val_matrix[i][j] == '0']
        end
        if length(co2_rows) == 1
            co2_rating = string_array_to_int([col[co2_rows[1]] for col in val_matrix])
        end
    end

    for i = 1:12

        count_0 = length([val_matrix[i][j] for j in o2_rows if val_matrix[i][j] == '0'])
        count_1 = length([val_matrix[i][j] for j in o2_rows if val_matrix[i][j] == '1'])
        if count_0 > count_1
            o2_rows = [j for j in o2_rows if val_matrix[i][j] == '0']
        elseif count_1 > count_0
            o2_rows = [j for j in o2_rows if val_matrix[i][j] == '1']
        else
            o2_rows = [j for j in o2_rows if val_matrix[i][j] == '1']
        end
        if length(o2_rows) == 1
            o2_rating = string_array_to_int([col[o2_rows[1]] for col in val_matrix])
        end
    end
    return o2_rating * co2_rating

end

println(solution1(lines))
println(solution2(lines))


text:
module AnalogDiscoveryNeuromore

using AnalogDiscovery

include("RingBuffer.jl")

export go

function float32msg!(buffer::Vector{UInt8},ch::Int64,val::Vector{Float32})
  buffer[1:5]=(ch==1) ? "/in/0".data : "/in/1".data
  @inbounds for i in 6:8
    buffer[i]=UInt8(0)
  end
  buffer[9:10]=",f".data
  @inbounds for i in 11:12
    buffer[i]=UInt8(0)
  end
  @inbounds for i in 1:4
    buffer[12+i]=unsafe_load(Ptr{UInt8}(pointer(val)),4-i+1)
  end
  #@show buffer
  return buffer
end

function go(fs,runForMin,twoChannels,lowVolt)
  cDev = enumDevices()
  @show ("Available",cDev)
  myDev = Int32(0)
  if cDev > 0
    if enumDeviceIsOpened(myDev)
      error("Error opening device, already in use.")
    end
    hdwf = deviceOpen(myDev)
  else
    error("No devices available.")
  end

  bmin,bmax = analogInBufferSizeInfo(hdwf)
  adbufsize = Int32(min(fs,bmax))
  adbufsize64 = Int64(adbufsize)
  analogInFrequencySet(hdwf,fs)
  analogInBufferSizeSet(hdwf,Int32(bmax)) # must be set at bmax for acqmodeRecord!
  analogInAcquisitionModeSet(hdwf,ACQMODERECORD)
  analogInRecordLengthSet(hdwf,0.0) # run indefinitely

  chVoltMin,chVoltMax,chSteps = analogInChannelRangeInfo(hdwf)
  chCount = analogInChannelCount(hdwf)
  for ch in Int32(0):(twoChannels?Int32(chCount-1):Int32(0))
    analogInChannelEnableSet(hdwf,ch,true)
    analogInChannelFilterSet(hdwf,ch,FILTERAVERAGE)
    analogInChannelRangeSet(hdwf,ch,lowVolt?chVoltMin:chVoltMax)
  end

  sock = UDPSocket(); const localh = ip"127.0.0.1"
  const fs_Int32 = round(Int32,fs)
  # Collect s-sec worth of data
  const s = 60*runForMin
  tBuf = Vector{Float64}(ceil(Int32,s*fs))

  # RingBuffers
  ch1buf = RingBuffer(Vector{Float64}(adbufsize),1,1)
  ch2buf = RingBuffer(Vector{Float64}(adbufsize),1,1)
  # Single-element buffer to reinterpret Float32 value later using unsafe_load
  valBuf::Vector{Float32} = Vector{Float32}(1)
  # Buffer for OSC packets sent over UDP
  udpBuf = Array{UInt8,1}(16)

  sleep(1.0) # gets rid of offset estimation transient
  analogInConfigure(hdwf,false,true) # Start acquistion

  # Prefill
  status = analogInStatus(hdwf,true)
  while status != DWFSTATETRIGGEREDRUNNING
    status = analogInStatus(hdwf,true)
  end
  cSamples::Int32 = 0
  chDataBuf = Vector{Float64}(adbufsize)
  sleep(0.25) # Deliberately introduce latency to improve throughput
  cAvailable,cLost,cCorrupted = analogInStatusRecord(hdwf)
  cSamples += cLost
  @inbounds for ch in Int32(0):(twoChannels?Int32(chCount-1):Int32(0))
    analogInStatusData!(chDataBuf,hdwf,ch,Int32(cAvailable))
    for n in 1:Int64(cAvailable)
      if ch==Int32(0)
        enque!(ch1buf,chDataBuf[n])
      else
        enque!(ch2buf,chDataBuf[n])
      end
    end
  end
  cSamples += cAvailable

  voltVal::Float64 = 0.0
  nSamples::Int32 = ceil(Int32,s*fs)
  startT = Libc.time(); t = -1
  display("Streaming data...")
  while t < nSamples-1
    if used(ch1buf) == 0
      status = analogInStatus(hdwf,true)
      if (cSamples == 0) && ((status == DWFSTATECONFIG) || (status == DWFSTATEPREFILL) || (status == DWFSTATEARMED))
        # Acquisition not yet started
        continue
      end
      cAvailable,cLost,cCorrupted = analogInStatusRecord(hdwf)
      # @show (cSamples,cAvailable,cLost,cCorrupted)
      cSamples += cLost
      if cAvailable==0
        continue
      end
      if cSamples+cAvailable > nSamples
        cAvailable = nSamples-cSamples
      end
      # get samples
      @inbounds for ch in Int32(0):(twoChannels?Int32(chCount-1):Int32(0))
        analogInStatusData!(chDataBuf,hdwf,ch,Int32(cAvailable))
        for n in 1:cAvailable
          if ch==0
            enque!(ch1buf,chDataBuf[n])
          else
            enque!(ch2buf,chDataBuf[n])
          end
        end
      end
      cSamples += cAvailable
    end

    # Now send UDP samples
    currT = Libc.time()
    theorT = t/fs+startT
    if (currT>=theorT)
      if used(ch1buf) > 0
        # @show used(ch1buf)
        t += 1
        tBuf[t+1] = theorT-currT
        for ch in Int32(0):(twoChannels?Int32(chCount-1):Int32(0))
          if ch==0
            voltVal = deque!(ch1buf)
          else
            voltVal = deque!(ch2buf)
          end
          valBuf[1] = Float32(voltVal)
          float32msg!(udpBuf,ch+1,valBuf)
          Base._send(sock,localh,UInt16(4545),udpBuf) # internal _send is faster and we don't care about failures
        end
      end
    end
  end

  deviceClose(hdwf)
  return tBuf
end

end # module


text:
#2021/7/18
#Edited by shwh
using Plots
using Fmt
function parameter()
    xmax = Float64(pi)
    xmin = 0.0
    xsp = 100
    dx = (xmax-xmin)/xsp

    ymax = Float64(pi)
    ymin = 0.0
    ysp = 100
    dy = (ymax-ymin)/ysp

    x = collect(xmin:dx:xmax)
    y = collect(ymin:dy:ymax)

    tstep = 600
    dt = 0.01
    c = 1.0
    alpha = (c*dt/dx)
    u = zeros(tstep, length(x), length(y))
    return x,y,dx,dy,dt,alpha,u,tstep
end

#set up initial condition
function IC(x,y,u,alpha)
    for i = 1:length(x)
        for j = 1:length(y)
            #=
            u[1,i,j] = exp(-10*(i*dx)^2)*exp(-10*(j*dy)^2)*2
            =#
            u[1,i,j] = exp(-10*(i*dx-2)^2)*exp(-10*(j*dy-2)^2)*2
        end
    end
    ubn = []
    push!(ubn,u[1,1,:])
    push!(ubn,u[1,end,:])
    push!(ubn,u[1,:,1])
    push!(ubn,u[1,:,end])

    p1 = surface(x,y,u[1,:,:],color=cgrad(:bwr),title = ("0.0s"),zlims=(-1,1),bg = :lightgray, camera = (30,80),colorbar = false)
    plot(p1)
    savefig("wave_"*string(Int(0))*".png")
    for i = 2:length(x)-1
        for j = 2:length(y)-1
            u[2,i,j] = u[1,i,j] + 0.5*alpha^2*(u[1,i+1,j] + u[1,i-1,j] + u[1,i,j+1] + u[1,i,j-1] -4*u[1,i,j])
        end
    end
    bound_d(x,y,u,2)
    #bound_n(x,y,u,2,ubn)
    return u,ubn
end

#bound
function bound_d(x,y,u,t)
    u[t,1,:] = u[t,2,:]
    u[t,length(x),:] = u[t,length(x)-1,:]
    u[t,:,1] = u[t,:,2]
    u[t,:,length(y)] = u[t,:,length(y)-1]
end

function bound_n(x,y,u,t,ubn)
    u[t,1,:] = ubn[1]
    u[t,length(x),:] = ubn[2]
    u[t,:,1] = ubn[3]
    u[t,:,length(y)] = ubn[4]
end

function main()
    x,y,dx,dy,dt,alpha,u,tstep = parameter()
    u,ubn = IC(x,y,u,alpha)
    #println(ubn)
    for t = 3:tstep
        tt = dt*(t-1)
        for i = 2:length(x)-1
            for j = 2:length(y)-1
                u[t,i,j] = 2*u[t-1,i,j] - u[t-2,i,j] + alpha^2*(u[t-1,i+1,j] + u[t-1,i-1,j] + u[t-1,i,j+1] + u[t-1,i,j-1] -4*u[t-1,i,j])
            end
        end
        bound_d(x,y,u,t)
        #bound_n(x,y,u,t,ubn)
        if (t-1)%10 == 0
            p2 = surface(x,y,u[t,:,:],color=cgrad(:bwr),title = (((Fmt.format(f"{:.2f}",tt))*"s")),zlims=(-1,1),bg = :lightgray, camera = (30,80),colorbar = false)
            plot(p2)
            savefig("wave_"*string((Int(t-1)))*".png")
        end
    end
end

main()


text:
using LinearTransforms, Arpack
using SparseArrays, LinearAlgebra
using Test

@testset "opMask" begin

    m, n = 123, 321
    M = sprand(Bool,m,n,0.5)
    I, J = findnz(M)[1:2]
    u, v = randn(m,2), randn(n,2)
    X = u*v'

    @testset "full matrices" begin
        op = opMask(M)

        # forward operator
        w = op*X
        wt = [X[i,j] for (i,j) in zip(I,J)]
        @test norm(w-wt, Inf) ≤ 1e-10*(1+norm(wt,Inf))

        # adjoint operator
        V = op'*w
        Vt = sparse(I, J, w, m, n)
        @test norm(V-Vt) ≤ 1e-10

        # dot test
        @test dot(op*X,w) ≈ dot(X,op'*w)

        # test eltypes
        op = opMask(M)
        @test eltype(op) == Float64
        op = opMask{Complex{Float64}}(M)
        @test eltype(op) == Complex{Float64}

        # test size
        @test size(op) == (nnz(M), size(M))

    end

    @testset "factored matrices" begin
        op = opMaskFactored(M)

        # forward operator
        w = op*(u,v)
        wT = [X[i,j] for (i,j) in zip(I,J)]
        @test norm(w-wT, Inf) ≤ 1e-10*(1+norm(wT,Inf))

        # adjoint operator
        V = op'*w
        Vt = sparse(I, J, w, m, n)
        @test norm(V-Vt) ≤ 1e-10

        # dot test
        @test dot(op*(u,v),w) ≈ dot(u*v',op'*w)

        # eltypes
        op = opMaskFactored(M)
        @test eltype(op) == Float64
        op = opMaskFactored{Int64}(M)
        @test eltype(op) == Int64

        # test size
        @test size(op) == (nnz(M), size(M))

    end
end


text:
"""
al_solver.jl
"""

import Base: show

@doc raw""" ```julia
struct AugmentedLagrangianSolver <: ConstrainedSolver{T}
```
Augmented Lagrangian (AL) is a standard tool for constrained optimization.
For a trajectory optimization problem of the form:
```math
\begin{aligned}
  \min_{x_{0:N},u_{0:N-1}} \quad & \ell_f(x_N) + \sum_{k=0}^{N-1} \ell_k(x_k, u_k, dt) \\
  \textrm{s.t.}            \quad & x_{k+1} = f(x_k, u_k), \\
                                 & g_k(x_k,u_k) \leq 0, \\
                                 & h_k(x_k,u_k) = 0.
\end{aligned}
```
AL methods form the following augmented Lagrangian function:
```math
\begin{aligned}
    \ell_f(x_N) + &λ_N^T c_N(x_N) + c_N(x_N)^T I_{\mu_N} c_N(x_N) \\
           & + \sum_{k=0}^{N-1} \ell_k(x_k,u_k,dt) + λ_k^T c_k(x_k,u_k)
+ c_k(x_k,u_k)^T I_{\mu_k} c_k(x_k,u_k)
\end{aligned}
```
This function is then minimized with respect to the primal variables using any
unconstrained minimization solver (e.g. iLQR).
After a local minima is found, the AL method updates the Lagrange multipliers λ
and the penalty terms μ and repeats the unconstrained minimization.
AL methods have superlinear convergence as long as the penalty term μ
is updated each iteration.
"""
struct AugmentedLagrangianSolver{T,S} <: ConstrainedSolver{T}
    opts::SolverOptions{T}
    stats::SolverStats{T}
    solver_uncon::S
    solver_name::Symbol
end

"""$(TYPEDSIGNATURES)
Form an augmented Lagrangian cost function from a Problem and AugmentedLagrangianSolver.
Does not allocate new memory for the internal arrays, but points to the arrays in the solver.
"""
function AugmentedLagrangianSolver(prob::Problem{IR,T}, opts::SolverOptions,
                                   stats::SolverStats; solver_uncon=iLQRSolver) where {IR,T}
    # set up unconstrained solver
    solver_uncon = solver_uncon(prob, opts, stats)
    # put it all together
    S = typeof(solver_uncon)
    solver = AugmentedLagrangianSolver{T,S}(opts, stats, solver_uncon, :AugmentedLagrangian)
    reset!(solver)
    return solver
end

# options methods
Base.show(io::IO, solver::AugmentedLagrangianSolver) = (
    print(io, "ALSolver")
)

Base.show(io::IO, m::MIME"text/plain", solver::AugmentedLagrangianSolver) = (
    print(io, "ALSolver")
)

function set_verbosity!(solver::AugmentedLagrangianSolver)
    llevel = log_level(solver)
    if is_verbose(solver)
        set_logger()
        Logging.disable_logging(LogLevel(llevel.level-1))
        logger = global_logger()
        if is_verbose(solver.solver_uncon)
            freq = 1
        else
            freq = 5
        end
        logger.leveldata[llevel].freq = freq
    else
        Logging.disable_logging(llevel)
    end
end

function reset!(solver::AugmentedLagrangianSolver)
    # reset duals and penalties
    for convals in solver.solver_uncon.obj.convals
        for conval in convals
            conval.λ .= 0
            conval.μ .= conval.con.params.μ0
            conval.a .= false
        end
    end
    # reset unconstrained solver
    reset!(solver.solver_uncon)
    return nothing
end


text:
module Levels

import ..ED
using .. SimLib
using ..SimLib: FArray
using SharedArrays: sdata

export Energies, LevelDataDescriptor, LevelData, load_levels

### Descriptor

struct LevelDataDescriptor <: ED.EDDerivedDataDescriptor
    derivedfrom::ED.EDDataDescriptor
end

LevelDataDescriptor(args...; kwargs...) = LevelDataDescriptor(EDDataDescriptor(args...; kwargs...))

### Data obj

struct LevelData{N} <: SimLib.AbstractSimpleData
    descriptor::LevelDataDescriptor
    data::FArray{N}
end

ED._default_folder(::LevelDataDescriptor) = "levels"
#_filename_addition(::LevelData) = "" # default

"""
    load_levels(edd)
    load_levels(model[, diagtype][, location])
"""
load_levels(args...; kwargs...) = load(LevelDataDescriptor(args...; kwargs...))

### Task

mutable struct LevelTask <: ED.EDTask
    data
end

Energies() = LevelTask(nothing)

function ED.initialize!(task::LevelTask, arrayconstructor, spectral_size)
    task.data = arrayconstructor(Float64, spectral_size)
end

function ED.compute_task!(task::LevelTask, evals, evecs, inds...)
    n = min(size(task.data,1), length(evals))
    task.data[1:n, inds...] .= view(evals, 1:n)
end

function ED.failed_task!(task::LevelTask, inds...)
    task.data[:, inds...] .= NaN64
end

function ED.assemble(task::LevelTask, edd)
    LevelData(LevelDataDescriptor(edd), sdata(task.data))
end

end # module


text:
using ServerMetrics
using HTTP

@testcase "Simple counter" begin
    clear_registry!(get_default_registry())
    c = register!(Counter(), "my_counter")
    inc!(c, 2.0)

    resp = handle_metrics(HTTP.Request("GET", "/metrics"), nothing)
    @test String(resp.body) == """
    # TYPE my_counter counter
    my_counter 2.0

    """

    inc!(c, 3.0)
    resp = handle_metrics(HTTP.Request("GET", "/metrics"), nothing)
    @test String(resp.body) == """
    # TYPE my_counter counter
    my_counter 5.0

    """
end

@testcase "Simple gauge" begin
    clear_registry!(get_default_registry())
    g = register!(Gauge(1.5), "my_gauge")
    resp = handle_metrics(HTTP.Request("GET", "/metrics"), nothing)

    @test String(resp.body) == """
    # TYPE my_gauge gauge
    my_gauge 1.5

    """
end

@testcase "Counter with labels" begin
    clear_registry!(get_default_registry())
    c = register!(Counter(;action=String,response_code=Int64), "requests")
    inc!(c; action="get", response_code=404)
    inc!(c; action="put", response_code=500)  # What a bad day for prod :-(
    resp = handle_metrics(HTTP.Request("GET", "/metrics"), nothing)

    @test String(resp.body) == """
    # TYPE requests counter
    requests{action="get",response_code="404"} 1.0
    requests{action="put",response_code="500"} 1.0

    """
end

@testcase "Gauge with labels" begin
    clear_registry!(get_default_registry())
    temp = register!(Gauge(;location=String, hour=Int64), "temperature")
    set!(temp, 36.0; location="outside", hour=6)
    set!(temp, 40.0; location="outside", hour=8)
    set!(temp, 60.0; location="inside", hour=8)
    resp = handle_metrics(HTTP.Request("GET", "/metrics"), nothing)

    @test String(resp.body) == """
    # TYPE temperature gauge
    temperature{hour="6",location="outside"} 36.0
    temperature{hour="8",location="inside"} 60.0
    temperature{hour="8",location="outside"} 40.0

    """
end

@testcase "Exported metrics are sorted" begin
    clear_registry!(get_default_registry())
    bbb = register!(Counter(), "bbb")
    zzz = register!(Gauge(0.2), "zzz")
    aaa = register!(Counter(), "aaa")
    resp = handle_metrics(HTTP.Request("GET", "/metrics"), nothing)

    @test String(resp.body) == """
    # TYPE aaa counter
    aaa 0.0

    # TYPE bbb counter
    bbb 0.0

    # TYPE zzz gauge
    zzz 0.2

    """
end


text:
# Singing voice separation demo using Robust PCA

using WAV
using RobustPCA

include(joinpath(dirname(@__FILE__), "stft.jl"))

filename = "titon_2_07_SNR5.wav"
x, fs = wavread(joinpath(dirname(@__FILE__), filename))
x = vec(x)

framelen = 1024
hopsize = div(framelen, 4)
win = hanning(framelen)
Xᶜ = stft(x, framelen, hopsize, win)

X, P = abs(Xᶜ), angle(Xᶜ)
@show size(X)

# performs RPCA
# A: row-rank matrix
# S: sparse matrix
elapsed_rpca = @elapsed begin
    s = 1.0/sqrt(maximum(size(X)))
    A, E = inexact_alm_rpca(X; verbose=true, error_tol=1.0e-7,
                            sparseness=s)
end
println("Elapsed time in Robust PCA: $(elapsed_rpca) sec.")

# Frequency masking
const binary_mask = true
if binary_mask
    m = abs(E) .> abs(A)
    Emask = zeros(size(E))
    Emask[m] = X[m]
    Amask = X - Emask
else
    Emask = E
    Amask = A
end

# back to time domain
a = istft(Amask .* exp(im * P), framelen, hopsize, win)
e = istft(Emask .* exp(im * P), framelen, hopsize, win)

name, suffix = splitext(filename)
wavwrite(float(a), string(name, "_A", suffix), Fs=fs)
wavwrite(float(e), string(name, "_E", suffix), Fs=fs)


text:
struct NotImplementedOperation <: AbstractGroupOperation end

struct NotImplementedManifold <: AbstractManifold{ℝ} end

struct NotImplementedGroupDecorator{M} <:
       AbstractDecoratorManifold{ℝ,TransparentGroupDecoratorType}
    manifold::M
end

struct DefaultTransparencyGroup{M,A<:AbstractGroupOperation} <:
       AbstractGroupManifold{ℝ,A,DefaultGroupDecoratorType}
    manifold::M
    op::A
end


text:
# ---
# title: 387. First Unique Character in a String
# id: problem387
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Hash Table, String
# link: <https://leetcode.com/problems/first-unique-character-in-a-string/description/>
# hidden: true
# ---
#
# Given a string, find the first non-repeating character in it and return its
# index. If it doesn't exist, return -1.
#
# **Examples:**
#
#
#
#     s = "leetcode"
#     return 0.
#
#     s = "loveleetcode"
#     return 2.
#
#
#
#
# **Note:** You may assume the string contains only lowercase English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Test
using NumericallyIntegrateArrays

@testset "simps" begin
    x = range(0,stop=π,length=10); dx = step(x)

    @testset "dx" begin
	    @testset "1D" begin
		    f = sin.(x)
		    @test isapprox(simps(f,step(x)),2,rtol=1e-3)
	    end
	    @testset "2D" begin
	        f = sin.(x) * cos.(x)'
	       	@test isapprox(simps(f,step(x)),2cos.(x),rtol=1e-3)
	    end
    end
    @testset "x array" begin
	    @testset "1D" begin
		    f = sin.(x)
		    @test isapprox(simps(f,x),2,rtol=1e-3)
	    end
	    @testset "2D" begin
	        f = sin.(x) * cos.(x)'
	       	@test isapprox(simps(f,x),2cos.(x),rtol=1e-3)
	    end
    end
end

@testset "trapz" begin
    x = range(0,stop=π,length=10); dx = step(x)

    @testset "dx" begin
	    @testset "1D" begin
		    f = sin.(x)
		    @test isapprox(trapz(f,step(x)),2,rtol=1e-1)
	    end
	    @testset "2D" begin
	        f = sin.(x) * cos.(x)'
	       	@test isapprox(trapz(f,step(x)),2cos.(x),rtol=1e-1)
	    end
    end
    @testset "x array" begin
	    @testset "1D" begin
		    f = sin.(x)
		    @test isapprox(trapz(f,x),2,rtol=1e-1)
	    end
	    @testset "2D" begin
	        f = sin.(x) * cos.(x)'
	       	@test isapprox(trapz(f,x),2cos.(x),rtol=1e-1)
	    end
    end
end


text:
######################################################################
# ekf.jl
# handles kalman filter stuff
######################################################################
#include("kalman.jl")

# length is size of one side of search domain
type EKF <: GaussianFilter
	mu::Vector{Float64}			# belief
	Sigma::Matrix{Float64}		# belief
	length::Float64

	initialized::Bool
	initializer::Initializer

	function EKF(m::SearchDomain)
		mu = m.length/2 * ones(2)
		S = 1e9 * eye(2)
		initer = NaiveInitializer(m)
		return new(mu, S, m.length,false,initer)
	end
end

# have to include search domain because jammer location factors in
# Really, I should just fold that into the state
# TODO: just subtracting is not ok, need circle distance (can be negative)
function update!(ekf::EKF, x::Vehicle, o::Float64)

	if !ekf.initialized
		check_initialization(ekf, x, o)
		return nothing
	end

	xr = ekf.mu[1] - x.x
	yr = ekf.mu[2] - x.y
	if xr == 0.0 && yr === 0.0
		xr = 1e-6
		yr = 1e-6
	end
	Ht = [yr, -xr]' * (180.0 / pi) / (xr^2 + yr^2)
	Kt = ekf.Sigma * Ht' * inv(Ht * ekf.Sigma * Ht' + x.sensor.noise_sigma^2)

	o_predict = true_bearing(x, ekf.mu)
	# Computing o - o_predict... both will be in range 0 to 360
	# Make sure this value is within -180,180
	o_diff = fit_180(o - o_predict)
	#print("o = ", round(o,1), ", op = ", round(o_predict,1), ", od = ", round(o_diff,1))
	#println("od = ", o_diff)
	#mu_t = ekf.mu + Kt * (o - o_predict)
	mu_t = ekf.mu + Kt * o_diff
	Sigma_t = (eye(2)  - Kt * Ht) * ekf.Sigma
	#println(": mu = ($(round(mu_t[1],1)), $(round(mu_t[2],1)))")

	# TODO: use copy here
	ekf.mu = vec(mu_t)
	ekf.Sigma = Sigma_t

	return nothing
end

centroid(ekf::EKF) = (ekf.mu[1], ekf.mu[2])
function reset!(f::EKF)
	f.mu = f.length/2 * ones(2)
	f.Sigma = 1e9 * eye(2)
end


text:
using JLD2, PyPlot, JUDI.TimeModeling

PyPlot.rc("font", family="serif")
PyPlot.rc("xtick", labelsize=20)
PyPlot.rc("ytick", labelsize=20)

@load "./data/gl_tti_g.jld"

n = model_true.n
d = model_true.d
extent = [0, (n[1] -1) * d[1],  (n[2] -2) * d[2], 0]

### Plot model
figure(figsize=(30, 30))
subplot(221)
imshow(model_true.m', vmin=1/4.5^2, vmax=.44, cmap="jet", extent=extent, aspect="auto")
ax1 = gca()
ax1.set_xticks([])
ylabel("Depth (m)", fontsize=20)
title("m", fontsize=20)
colorbar(ax=ax1)

subplot(222)
imshow(model_true.epsilon', vmin=0.0, vmax=.3, cmap="jet", extent=extent, aspect="auto")
ax1 = gca()
ax1.set_xticks([])
ax1.set_yticks([])
title("ε", fontsize=20)
colorbar(ax=ax1)

subplot(223)
imshow(model_true.delta', vmin=0.0, vmax=.3, cmap="jet", extent=extent, aspect="auto")
ax1 = gca()
xlabel("X (m)", fontsize=20)
ylabel("Depth (m)", fontsize=20)
title("δ", fontsize=20)
colorbar(ax=ax1)

subplot(224)
imshow(model_true.theta', vmin=-.25, vmax=.25, cmap="jet", extent=extent, aspect="auto")
xlabel("X (m)", fontsize=20)
ax1 = gca()
ax1.set_yticks([])
title("θ", fontsize=20)
colorbar(ax=ax1)

savefig("./data/gl_model.png", bbox_inches="tight")

### Plot inital model and perturbationgradients
plt.figure(figsize=(30, 30))
subplot(121)
imshow(model0_tti.m', vmin=1/4.5^2, vmax=.44, cmap="jet", extent=extent, aspect="auto")
xlabel("X (m)", fontsize=20)
ylabel("Depth (m)", fontsize=20)
ax1 = gca()
title("m0", fontsize=20)
ax1 = gca()
colorbar(ax=ax1)

subplot(122)
imshow(model0_tti.m' - model_true.m', vmin=-.05, vmax=.05, cmap="jet", extent=extent, aspect="auto")
xlabel("X (m)", fontsize=20)
ax1 = gca()
ax1.set_yticks([])
title("δ m", fontsize=20)
colorbar(ax=ax1)

savefig("./data/gl_dm.png", bbox_inches="tight")

### Plot gradients

# True tti params
@load "./data/gl_tti_g.jld"
figure(figsize=(30, 30))
subplot(321)
imshow(g_wri_tti', vmin=-1, vmax=1, cmap="jet", extent=extent, aspect="auto")
ylabel("Depth (m)", fontsize=20)
ax1 = gca()
ax1.set_xticks([])
# ax1.set_yticks([])
title("WRI true anisotropy", fontsize=20)
colorbar(ax=ax1)
subplot(322)
imshow(g_fwi_tti', vmin=-.5, vmax=.5, cmap="jet", extent=extent, aspect="auto")
ax1 = gca()
ax1.set_xticks([])
ax1.set_yticks([])
title("FWI true anisotropy", fontsize=20)
colorbar(ax=ax1)

# Error in tti params
@load "./data/gl_tti_g_w.jld"
subplot(323)
imshow(g_wri_tti_w', vmin=-1, vmax=1, cmap="jet", extent=extent, aspect="auto")
ylabel("Depth (m)", fontsize=20)
ax1 = gca()
ax1.set_xticks([])
title("WRI inaccurate anisotropy", fontsize=20)
colorbar(ax=ax1)
subplot(324)
imshow(g_fwi_tti_w', vmin=-.5, vmax=.5, cmap="jet", extent=extent, aspect="auto")
ax1 = gca()
ax1.set_xticks([])
ax1.set_yticks([])
title("FWI inaccurate anisotropy", fontsize=20)
colorbar(ax=ax1)

# Acosutic WE gradient
@load "./data/gl_acou_g.jld"
subplot(325)
imshow(g_wri_a', vmin=-.1, vmax=.1, cmap="jet", extent=extent, aspect="auto")
xlabel("X (m)", fontsize=20)
ylabel("Depth (m)", fontsize=20)
ax1 = gca()
title("WRI acoustic", fontsize=20)
colorbar(ax=ax1)
subplot(326)
imshow(g_fwi_a', vmin=-.25, vmax=.25, cmap="jet", extent=extent, aspect="auto")
xlabel("X (m)", fontsize=20)
ax1 = gca()
ax1.set_yticks([])
title("FWI acoustic", fontsize=20)
colorbar(ax=ax1)

savefig("./data/gl_wrong.png", bbox_inches="tight")


text:
using VideoCaptureWrap
using Test

@testset "VideoCaptureWrap.jl" begin
    # Write your tests here.
end


text:
# Autogenerated wrapper script for CALCEPH_jll for x86_64-apple-darwin
export libcalceph

JLLWrappers.@generate_wrapper_header("CALCEPH")
JLLWrappers.@declare_library_product(libcalceph, "@rpath/libcalceph.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcalceph,
        "lib/libcalceph.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module DrawFunc

using TargetFunction
using PyCall
@pyimport matplotlib.pyplot as plt
@pyimport matplotlib.cm as cm

function draw2D(fun::TargetFunction.ObjFun)

    dx = [-3.0:.02:2.5]
    dy = [-4.5:.02:4.5]'

    m, n = length(dx), length(dy)
    dx = repmat(dx, 1, n)[:]
    dy = repmat(dy, m, 1)[:]
    z = [fun.f_x([i;j]) for (i, j) in zip(dx, dy)]

    dx, dy, z = reshape(dx, m, n), reshape(dy, m, n), reshape(z, m, n)
    cs = plt.contourf(dx, dy, z, 20, cmap=cm.get_cmap("YlGnBu"))
end

end

max_stars_repo_name: wyli/gradient

max_issues_repo_name: wyli/gradient

max_forks_repo_name: wyli/gradient


text:
using Test, Sandbox, Base.BinaryPlatforms, LazyArtifacts

@testset "multiarch" begin
    # This is a set of multiarch platforms that is _not_ our current platform
    native_arch = arch(HostPlatform())
    sub_arch = native_arch
    if native_arch == "x86_64"
        sub_arch = "i686"
    elseif native_arch == "aarch64"
        sub_arch = "armv7l"
    end
    alien_arch = native_arch ∈ ("x86_64", "i686") ? "aarch64" : "x86_64"
    @testset "argument parsing" begin
        # Test that our `multiarch` kwarg is correctly parsed
        config = SandboxConfig(
            Dict("/" => Sandbox.multiarch_rootfs());
            multiarch = [
                Platform(alien_arch, "linux"; libc="glibc"),
                Platform(alien_arch, "linux"; libc="musl"),
                Platform(native_arch, "linux"; libgfortran_version=v"4"),
                Platform(native_arch, "linux"; libgfortran_version=v"5"),
                Platform(sub_arch, "linux"),
            ],
        )

        # qemu doesn't care about `libc` or `libgfortran_version` or anything like that.
        # Also, native architectures (and sub-architectures such as `i686` for `x86_64`,
        # or `armv7l` for `aarch64`) get ignored, so we end up with only one multiarch
        # format from all that above, which is just `alien_arch`.
        @test length(config.multiarch_formats) == 1
        @test occursin(alien_arch, config.multiarch_formats[1].name)
    end

    # Of our available executors, let's check to see if each can be used to run multiarch workloads
    multiarch_executors = filter(executor_available, Sandbox.all_executors)
    old_binfmt_misc_regs = nothing

    if get(ENV, "SANDBOX_TEST_MULTIARCH", "true") != "true"
        @warn("Refusing to test multiarch because SANDBOX_TEST_MULTIARCH set to $(ENV["SANDBOX_TEST_MULTIARCH"])")
        multiarch_executors = Sandbox.SandboxExecutor[]
    elseif Sys.islinux()
        # On Linux, we need passwordless sudo to be able to register things
        if !success(`sudo -k -n true`)
            @warn("Refusing to test multiarch on a system without passwordless sudo!")
            multiarch_executors = Sandbox.SandboxExecutor[]
        end

        # Otherwise, let's save the current set of binfmt_misc registrations
        old_binfmt_misc_regs = Sandbox.read_binfmt_misc_registrations()
    end

    for executor in multiarch_executors
        if Sys.islinux()
            # Start by clearing out the binfmt_misc registrations, so that each executor
            # has to set things up from scratch.
            Sandbox.clear_binfmt_misc_registrations!()
        end

        @testset "HelloWorldC_jll" begin
            multiarch = [
                Platform("x86_64", "linux"; libc="glibc"),
                Platform("x86_64", "linux"; libc="musl"),
                Platform("i686", "linux"; libc="glibc"),
                Platform("i686", "linux"; libc="musl"),
                Platform("aarch64", "linux"; libc="glibc"),
                Platform("aarch64", "linux"; libc="musl"),
                Platform("armv7l", "linux"; libc="glibc"),
                Platform("armv7l", "linux"; libc="musl"),
                Platform("powerpc64le", "linux"; libc="glibc"),
                # We don't have this one yet
                #Platform("powerpc64le", "linux"; libc="musl"),
            ]
            stdout = IOBuffer()
            stderr = IOBuffer()
            config = SandboxConfig(
                Dict(
                    "/" => Sandbox.multiarch_rootfs(),
                    "/apps" => LazyArtifacts.ensure_artifact_installed("multiarch-testing", joinpath(dirname(@__DIR__), "Artifacts.toml")),
                );
                multiarch,
                stdout,
                stderr,
            )
            with_executor(executor) do exe
                for platform in multiarch
                    @test success(exe, config, `/apps/hello_world.$(triplet(platform))`)
                    @test String(take!(stdout)) == "Hello, World!\n";
                    @test isempty(take!(stderr))
                end
            end
        end
    end

    if old_binfmt_misc_regs !== nothing && !isempty(old_binfmt_misc_regs)
        # Restore old binfmt_misc registrations so that our test suite isn't clobbering things for others
        Sandbox.clear_binfmt_misc_registrations!()
        Sandbox.write_binfmt_misc_registration!.(old_binfmt_misc_regs)
    end
end


text:
mutable struct EnhancedDirectMethod{S<:RatesCacheTrait} <: DecoupledSearchAlgorithm
  rates::Vector{Float64}
  total_rate::Float64
end

RatesCache(::EnhancedDirectMethod{S}) where S = S()

@inline search_jumps(algorithm::EnhancedDirectMethod) = search_jump_rates(RatesCache(algorithm), algorithm.rates, algorithm.total_rate)

@inline update!(algorithm::EnhancedDirectMethod, state, model, j) = update_jump_rates!(algorithm, state, model, j)


text:
include("utils/utils.jl")
include("preprocessors.jl")
include("approximators/approximators.jl")
include("environment_models/environment_models.jl")
include("learners/learners.jl")
include("policies/policies.jl")
include("iteration_methods.jl")