text:
module Summary
using Sessions, TaskData, Analysis, SpikeTools, PlotTools, PyPlot, PyCall, JLD,
      HypothesisTests, Distributions, DataFrames, DataStructures, StatsBase, RCall
@pyimport matplotlib.gridspec as gridspec
@pyimport matplotlib.patches as patches
const slice = pybuiltin(:slice)

## Results

unit_mask(chs, region=:MTL, su_only::Bool=false) =
    [(ch.unit_type == :SU || (!su_only && ch.unit_type == :MU)) && ismatch(REGIONS[region], ch.site)
     for ch in chs]

function to_pct(x::AbstractArray{Bool})
    @sprintf "%.2g%% (%d/%d)" mean(x)*100 sum(x) length(x)
end

function region_breakdown(chs, sig::AbstractArray{Bool})
    println("    $(to_pct(sig[[ismatch(REGIONS[:PHG], ch.site) for ch in chs]])) Parahippocampal Units")
    println("    $(to_pct(sig[[ismatch(REGIONS[:EC], ch.site) for ch in chs]])) Entorhinal Units")
    println("    $(to_pct(sig[[ismatch(REGIONS[:H], ch.site) for ch in chs]])) Hippocampal Units")
    println("    $(to_pct(sig[[ismatch(REGIONS[:A], ch.site) for ch in chs]])) Amygdala Units")
end

function stimsel_info()
    regions = [:PHG, :EC, :H, :A]
    @load joinpath(PROCESSED_DATA_DIR, "permpval.jld") chs pval
    selch = unit_mask(chs, :MTL)
    chs = chs[selch]
    pval = pval[selch]
    sig = pval .< 0.001

    println("$(to_pct(sig)) units responsive")
    region_breakdown(chs, sig)
    proportions = [begin
        maskregion = unit_mask(chs, region)
        (sum(sig[maskregion]), sum(maskregion))
    end for region in regions]
    for ir1 = 1:length(regions), ir2 = ir1+1:length(regions)
        p1 = proportions[ir1]
        p2 = proportions[ir2]
        @printf "%3s vs. %3s p = %.3g\n" regions[ir1] regions[ir2] pvalue(FisherExactTest(p1[1], p1[2]-p1[1], p2[1], p2[2]-p2[1]))
    end
end

function showstimulus(f)
    text(0.5, 0.5, f, horizontalalignment="center", verticalalignment="center", fontsize=9)
    axis("off")
end

function plot_stimulus_all(epoch, session, channel_number, cluster_number; header::Bool=true, highlight=nothing, out="pdf")
    const LIMS = epoch == :maintenance ? (-1000, 2400) : (-200, 800)
    spikes = session.spikes[channel_number][cluster_number]
    sesstrials = session.trials
    trials = stimulusframe(sesstrials)
    nstim = length(session.stimuli)
    ticks = epoch == :maintenance ? (0:1000:2000) : 0:600:1000
    gcf()[:set_size_inches](4.488, 2.1)
    if header
        suptitle("$(session.folder) channel $(channel_number) cluster $(cluster_number) $(session.sites[channel_number]) $epoch")
        top = 0.93
    else
        top = 0.995
    end
    gs = gridspec.GridSpec(1, nstim, left=0.06, right=0.98, bottom=0.135, top=top, wspace=0.2, hspace=0.0)
    local ax
    nthimg = 4
    nthraster = 3
    rasterbin = epoch == :maintenance ? 200 : 50
    if isa(highlight, Void)
        highlighttrials = Int[]
    else
        stims = sesstrials[:stimuli].data::Vector{Vector{Int}}
        highlighttrials = find(Bool[highlight in x for x in stims])
    end

    for stim = 1:nstim
        gs_item = gs[:__getitem__]((0, stim-1))
        inner_gs = gridspec.GridSpecFromSubplotSpec(nthimg, 1, subplot_spec=gs_item, hspace=0.02)
        subplot(inner_gs[:__getitem__](0))
        showstimulus(stim)

        gs_item = inner_gs[:__getitem__](slice(1, nthimg))
        inner_gs = gridspec.GridSpecFromSubplotSpec(nthraster, 1, subplot_spec=gs_item, hspace=0.0)
        axa = subplot(inner_gs[:__getitem__]((slice(0, nthraster-1))), zorder=1000)
        local usetrials::Vector{Int}
        if epoch == :sample
            stims = sesstrials[:stimuli].data::Vector{Vector{Int}}
            usetrials = find([stim in x for x in stims])
            onsets = trials[:onset][trials[:stimulus] .== stim].data
        elseif epoch == :maintenance
            stims = sesstrials[:stimuli].data::Vector{Vector{Int}}
            usetrials = find([stim in x for x in stims])
            onsets = sesstrials[:maint][usetrials]
        end
        t1 = findin(usetrials, setdiff(usetrials, highlighttrials))
        t2 = findin(usetrials, highlighttrials)

        plotraster(spikes, onsets[[t1; t2]], LIMS)
        if !isa(highlight, Void)
            axhline(length(t1)+0.5, color=(1, 0, 0, 0.5), linewidth=0.5)
        end
        axis("off")
        axvline(0, ymin=-1/(nthraster-1), ymax=1, linestyle=":", color="darkcyan", dashes=(1, 1),
                   linewidth=0.5, clip_on=false)

        if stim == 1
            ylabel("Presentation Number")
            bax = ax = subplot(inner_gs[:__getitem__](nthraster-1))
        else
            gca()[:set_yticklabels](())
            bax = subplot(inner_gs[:__getitem__](nthraster-1), sharey=ax)
        end
        rast = raster(spikes, onsets, (LIMS[1]:rasterbin:LIMS[2])/1000)
        bar(LIMS[1]:rasterbin:LIMS[2]-rasterbin, mean(rast, 2)*1000/rasterbin, rasterbin, color="k", edgecolor="none", clip_on=false)
        if !isa(highlight, Void)
            bar(LIMS[1]:rasterbin:LIMS[2]-rasterbin, mean(rast[:, t1], 2)*1000/rasterbin,
                rasterbin, color=(1, 0, 0, 0.5), edgecolor="none", clip_on=false)
        end
        if epoch == :maintenance
            xticks(ticks, ["0"; [string(Int(x/1000)) for x in ticks[2:end]]], fontsize=7)
            stim == 1 && xlabel("Time After Mask Offset (Seconds)", fontsize=7, labelpad=2, horizontalalignment="left", x=-0.6)
        else
            xticks(ticks, fontsize=7)
            stim == 1 && xlabel("Time After Stimulus Onset (ms)", fontsize=7, labelpad=2, horizontalalignment="left", x=-0.6)
        end
        if stim == 1
            ylabel("Firing Rate (Hz)", fontsize=7, labelpad=3)
            removespines(["top", "right"])
            gca()[:yaxis][:set_ticks_position]("left")
            gca()[:xaxis][:set_ticks_position]("bottom")
        else
            bax[:axes][:get_yaxis]()[:set_visible](false)
            removespines(["left", "top", "right"])
            gca()[:xaxis][:set_ticks_position]("bottom")
        end
        tick_params(axis="y", pad=2)
        tick_params(axis="x", direction="out", length=3, pad=1)
        xlim(LIMS[1], LIMS[2])
    end
    sca(ax)
    yl = ax[:get_ylim]()
    yticks((0:max(div(yl[2], 3), 1):yl[2]), fontsize=7)
    savefig("$(session.folder)_ch$(channel_number)c$(cluster_number)_$(session.sites[channel_number])_$(session.unit_type[channel_number][cluster_number])_all_$(epoch).$(out)", dpi=300)
end

## Responses Persist Until Subsequent Image Presentations

function read_resptime(about; notfirst::Bool=false)
    @load joinpath(PROCESSED_DATA_DIR, "latency.jld") chs isi_windows isi_windows_notfirst
    if notfirst
        isi_windows = isi_windows_notfirst
    end
    sel = findin(chs, intersect(Analysis.allisisigunits()))
    if about == :latency
        x = [isi_windows[isel][iisi][1] for iisi = 1:4, isel = sel]
    elseif about == :duration
        x = [isi_windows[isel][iisi][2]-isi_windows[isel][iisi][1] for iisi = 1:4, isel = sel]
    else
        error("invalid about $about")
    end
    x, chs[sel]
end

function resptime_anova(about::Symbol=:duration; notfirst::Bool=false)
    x, chs = read_resptime(about, notfirst=notfirst)
    isi = vec(repmat([1, 2, 3, 4], 1, size(x, 2)))
    region = vec(repmat(reshape([string(siteregion(x)) for x in chs], 1, length(chs)), size(x, 1), 1))
    cell = vec(repmat(reshape(1:length(chs), 1, length(chs)), size(x, 1), 1))
    [@assert sum((isi .== y) & (cell .== x)) == 1 for x in cell, y in isi]
    println()
    println(R"""
    library(ez)
    df = data.frame(x=$(vec(x)), isi=factor($isi), region=factor($region), cell=factor($cell))
    ezANOVA(df, x, cell, within=isi, between=region, type=2, detailed=TRUE)
    """)
    for i = 1:length(ISI), j = i+1:length(ISI)
        test = OneSampleTTest(x[i, :], x[j, :])
        @printf("ISI %4d vs. %4d  t(%d) = %.1f p = %.3g\n", ISI[i]*1000, ISI[j]*1000, test.df, test.t, pvalue(test))
    end
end

function resptime_table(about::Symbol=:duration)
    x, chs = read_resptime(about)
    regions = [:PHG, :EC, :H, :A, :All]
    txt = Array(String, length(regions), length(ISI)+1)
    for i = 1:length(regions)
        masked = x[:, unit_mask(chs, regions[i])]
        for j = 1:length(ISI)
            txt[i, j] = @sprintf "%d (%d)" mean(masked[j, :]) std(masked[j, :])
        end
        txt[i, end] = @sprintf "%d (%d)" mean(mean(masked, 1)) std(mean(masked, 1))
    end
    [["" reshape(["ISI $(x)" for x in round(Int, ISI*1000)], 1, 4) "Average"]; regions txt]
end

function plot_isi_single(session, channel_number, cluster_number, stim)
    @load joinpath(PROCESSED_DATA_DIR, "latency.jld") chs isi_windows
    isiwins = isi_windows[findfirst(chs, ClusterMetadata(session, channel_number, cluster_number))]

    const LIMS = (-200, 1500)
    spikes = session.spikes[channel_number][cluster_number]
    trials = stimulusframe(session.trials)
    nstim = maximum(trials[:stimulus])::Int
    byval::Vector{Float64} = sort(unique(trials[:isi]))
    gcf()[:set_size_inches](3.34, 1.5)
    # ioff()
    gs = gridspec.GridSpec(2, 4, left=0.1, right=0.99, bottom=0.2, top=0.9, wspace=0.15, hspace=0.0)
    rasterbin = 100
    local ax
    for ibv = 1:length(byval)
        subplot(gs[:__getitem__]((0, ibv-1)), zorder=1000)
        onsets = trials[:onset][(trials[:stimulus] .== stim) & (trials[:isi] .== byval[ibv])].data
        plotraster(spikes, onsets, LIMS)
        xlim(LIMS[1], LIMS[2])
        axis("off")
        isi = Int(byval[ibv]*1000)
        title("ISI = $(isi) ms", y=0.98, fontsize=8)
        axvline(isiwins[ibv][1], ymin=-1, ymax=1, linestyle=":", color="darkcyan", dashes=(1, 1),
                linewidth=0.5, clip_on=false)
        axvline(isiwins[ibv][2], ymin=-1, ymax=1, linestyle=":", color="darkcyan", dashes=(1, 1),
                linewidth=0.5, clip_on=false)
        println("ISI = $(isi) ms onset = $(isiwins[ibv][1]) offset = $(isiwins[ibv][2]) duration = $(isiwins[ibv][2]-isiwins[ibv][1])")

        if ibv == 1
            setp(gca()[:get_yticklabels](), fontsize=9)
        else
            gca()[:set_yticklabels](())
        end

        if ibv == 1
            bax = ax = subplot(gs[:__getitem__]((1, ibv-1)))
        else
            bax = subplot(gs[:__getitem__]((1, ibv-1)), sharey=ax)
        end
        rast = raster(spikes, onsets, (LIMS[1]:rasterbin:LIMS[2])/1000)
        bar(LIMS[1]:rasterbin:LIMS[2]-rasterbin, mean(rast, 2)*1000/rasterbin, rasterbin, color="k", edgecolor="none", clip_on=false)
        xticks(0:500:1500, ("0", "", "1000", ""), fontsize=7)
        removespines(["top", "right"])
        bax[:yaxis][:set_ticks_position]("left")
        bax[:xaxis][:set_ticks_position]("bottom")
        if ibv == 1
            setp(bax[:get_yticklabels](), fontsize=7)
            setp(bax[:get_xticklabels](), fontsize=7)
            ylabel("Firing Rate (Hz)", fontsize=7, labelpad=3)
            removespines(["top", "right"])
            gca()[:yaxis][:set_ticks_position]("left")
            gca()[:xaxis][:set_ticks_position]("bottom")
            xlabel("Time After Stimulus Onset (ms)", fontsize=7, labelpad=2, horizontalalignment="left", x=-0.4)
        else
            bax[:set_xticklabels](())
            bax[:axes][:get_yaxis]()[:set_visible](false)
            removespines(["left", "top", "right"])
            gca()[:xaxis][:set_ticks_position]("bottom")
        end
        tick_params(axis="y", pad=3)
        tick_params(axis="x", direction="out", length=3, pad=1)
        xlim(LIMS[1], LIMS[2])
    end
    yl = ax[:get_ylim]()
    yticks((0:max(div(yl[2], 4), 1):yl[2]), fontsize=7)
    savefig("$(session.folder)_ch$(channel_number)c$(cluster_number)_isi_stim$stim.pdf", dpi=300)
end

function plot_resptime(about::Symbol=:duration)
    gcf()[:set_size_inches](3.34, 1.5)
    x, chs = read_resptime(about)
    scatter(repeat(ISI*1000, outer=size(x, 2))+(rand(length(x))-0.5)*100, vec(x), 20, marker=".", alpha=0.4, color="k", edgecolor="none")
    means = zeros(4)
    stds = zeros(4)
    for i = 1:4
        xi = x[i, :]
        means[i] = mean(xi)
        stds[i] = std(xi)/sqrt(length(xi))
    end
    errorbar(ISI*1000, means, stds, color="k")
    removespines()
    gca()[:yaxis][:set_ticks_position]("left")
    gca()[:xaxis][:set_ticks_position]("bottom")
    xticks(ISI*1000, fontsize=7)
    yticks(fontsize=7)
    xlim(-100, 900)
    ylim(0, 1085)
    xlabel("Inter-Stimulus Interval (ms)", fontsize=7, labelpad=3)
    ylabel("Response Duration (ms)", fontsize=7, labelpad=2, y=0.45)
    gap = 10
    ht1 = 1025
    len1 = 25
    makeline(ISI[1]*1000, ISI[2]*1000-gap, ht1, ht1, len1, pvalue(OneSampleTTest(x[1, :], x[2, :])))
    makeline(ISI[2]*1000+gap, ISI[3]*1000-gap, ht1, ht1, len1, pvalue(OneSampleTTest(x[2, :], x[3, :])))
    makeline(ISI[3]*1000+gap, ISI[4]*1000, ht1, ht1, len1, pvalue(OneSampleTTest(x[3, :], x[4, :])))
    tight_layout(rect=[-.03, -.07, 1.04, 1.0])
    savefig("$(about)_by_isi.pdf")
end

function resptime_region_differences(about::Symbol=:duration)
    x, chs = read_resptime(about)
    regions = [:PHG, :EC, :H, :A]
    region = [siteregion(ch) for ch in chs]
    regionnum = [findfirst(regions, x) for x in region]
    meandur = vec(mean(x, 1))
    means = zeros(4)
    stds = zeros(4)
    for i = 1:4
        xi = meandur[regionnum .== i]
        means[i] = mean(xi)
        stds[i] = std(xi)/sqrt(length(xi))
    end
    PHG = meandur[regionnum .== 1]
    EC = meandur[regionnum .== 2]
    H = meandur[regionnum .== 3]
    A = meandur[regionnum .== 4]
    @printf("PHG: %d ms\n", mean(PHG))
    @printf("EC: %d ms\n", mean(EC))
    @printf("H: %d ms\n", mean(H))
    @printf("A: %d ms\n", mean(A))
    for ir1 = 1:length(regions), ir2 = ir1+1:length(regions)
        md1 = meandur[regionnum .== ir1]
        md2 = meandur[regionnum .== ir2]
        uet = UnequalVarianceTTest(md1, md2)
        @printf("%3s vs. %3s t(%d) = %.2g p = %.2g\n",
            regions[ir1],
            regions[ir2],
            uet.df,
            uet.t,
            pvalue(uet))
    end
end

function plot_alignedisi(; ch=nothing, compare::Bool=false, pos1::Bool=false)
    @load(joinpath(PROCESSED_DATA_DIR, "alignedinfoisi$(pos1 ? "_pos1" : "").jld"), chs, plot_epoch, alignedinfo,
          alignedinfoperms, alignedinfoisidiff, alignedinfoisidiffperms)
    if compare
        alignedinfo = alignedinfoisidiff
        alignedinfoperms = alignedinfoisidiffperms
    end
    _, sel = read_resptime(:duration)
    @assert sel == chs
    info = cat(3, alignedinfo...)
    alignedinfoperms ./= length(chs)

    clf()
    gcf()[:set_size_inches](3.34, 2)
    win = plot_epoch
    xv = win[1]:win[2]-1
    xvsel = 1:length(xv)
    xv = xv[xvsel]
    info = info[xvsel, :, :]
    mu = mean(info, 3)
    sd = std(info, 3)
    lim = 0.0
    limmin = 0.0
    colors = ["b", "g", "r", "c"]
    compare && shift!(colors)
    for i = 1:size(info, 2)
        y = vec(mu[:, i])
        yerr = vec(sd[:, i])/sqrt(size(info, 3))
        lim = max(lim, maximum(y+yerr))
        limmin = min(limmin, minimum(y-yerr))
        errbarpatch(xv, y-yerr, y+yerr, color=colors[i])
        plot(xv, y, linewidth=0.5, color=colors[i])
    end
    xlabel("Time after Response Onset (ms)", fontsize=7, labelpad=3)
    if compare
        lstrings = ["ISI $(round(Int, ISI[i+1]*1000-200)) vs. $(round(Int, ISI[i]*1000-200)) ms" for i = 1:length(ISI)-1]
    else
        lstrings = ["ISI $(round(Int, x)) ms" for x in ISI*1000-200]
    end
    if pos1 && compare
        ylabel("Proportion Variance Explained (\$\\omega^2\$)", fontsize=7, labelpad=1)
    else
        ylabel("Proportion Variance Explained (\$\\omega^2\$)", fontsize=7)
    end
    removespines()
    gca()[:yaxis][:set_ticks_position]("left")
    gca()[:xaxis][:set_ticks_position]("bottom")
    tick_params(axis="x", length=2)
    xticks(fontsize=7)
    yticks(fontsize=7)
    if pos1 && compare
        tight_layout(rect=[-.06, -.06, 1.04, 1.047])
    else
        tight_layout(rect=[-.07, -.06, 1.04, 1.047])
    end
    xlim((-200, 1000))
    limmin = pos1 || compare ? floor(limmin, 4) : 0.0
    plot([0, 0], [limmin, lim*.99], color="k", linestyle=":")
    plot([200, 200], [limmin, lim*.99], color="b", linestyle=":")
    plot([400, 400], [limmin, lim*.99], color="g", linestyle=":")
    plot([700, 700], [limmin, lim*.99], color="r", linestyle=":")
    plot([1000, 1000], [limmin, lim*.99], color="c", linestyle=":", clip_on=false)
    legend(lstrings, fontsize=7, loc="upper right", bbox_to_anchor=(1, 0.9), borderaxespad=0)

    lim = lim+(lim - limmin)*0.09
    ylim(limmin, lim)
    spacing = lim/50

    maxes = maximum(alignedinfoperms, 1)
    thr = mapslices(x->quantile(x, 0.95), maxes, 2)
    for iisi = 1:size(info, 2)
        pv = vec(mean(maxes[1:1, :, iisi] .>= mu[:, iisi], 2))
        sig = pv .< 0.05
        sigs = find(sig)::Vector{Int}
        biggestsig = (-1, -1)
        for i = 1:length(sigs)
            j = i + 1
            while j <= length(sigs) && sigs[j] == sigs[j-1]+1
                sigs[j] - sigs[i] > biggestsig[2] - biggestsig[1] && (biggestsig = (sigs[i], sigs[j]))
                j += 1
            end
        end
        if compare
            println("ISI $(ISI[iisi+1]*1000) - $(ISI[iisi]*1000) min p = $(minimum(pv)) @ $(xv[indmin(pv)])")
            sum(sig) > 0 && println("ISI $(ISI[iisi+1]*1000) - $(ISI[iisi]*1000) significant in $((xv[biggestsig[1]], xv[biggestsig[2]]))")
        else
            println("ISI $(ISI[iisi]*1000) min p = $(minimum(pv)) @ $(xv[indmin(pv)])")
            sum(sig) > 0 && println("ISI $(ISI[iisi]*1000) significant in $((xv[biggestsig[1]], xv[biggestsig[2]]))")
        end
        scatter(xv[find(sig)], fill(lim-(4.5-iisi)*spacing, sum(sig)), 1, edgecolor="none", color=colors[iisi], zorder=10000, marker="s")
    end
    (pos1 || compare) && axhline(0, color="k", zorder=-1)
    savefig("alignedisi$(compare ? "_compare" : "")$(pos1 ? "_pos1" : "").pdf")
end

## Maintenance Period Activity

function wm_units(epoch=:all)
    @load joinpath(PROCESSED_DATA_DIR, "permpval.jld") chs pval
    chs_pres = chs
    pval_pres = pval
    if epoch == :all
        @load joinpath(PROCESSED_DATA_DIR, "permpval_maintenance.jld") pval chs
    elseif epoch == :last
        @load joinpath(PROCESSED_DATA_DIR, "permpval_maintenance_pos.jld") pval_last chs
        pval = pval_last
    elseif epoch == :notlast
        @load joinpath(PROCESSED_DATA_DIR, "permpval_maintenance_pos.jld") pval_notlast chs
        pval = pval_notlast
    end
    mask = unit_mask(chs)
    pval = pval[mask]
    chs = chs[mask]
    pval_pres = pval_pres[indexin(chs, chs_pres)]
    presalpha = 0.001

    for alpha2 in (0.01, 0.05, 0.001)
        selsel2 = pval[pval_pres .< presalpha] .< alpha2
        @printf("%s of stimulus-selective units significant at maintenance at α = %s (p = %.2g)\n",
        to_pct(selsel2), alpha2, pvalue(BinomialTest(selsel2, alpha2), tail=:right))
    end

    for alpha2 in (0.01, 0.05, 0.001)
        selsel3 = pval[pval_pres .>= presalpha] .< alpha2
        @printf("%s of non-stimulus-selective units significant at maintenance at α = %s (p = %.2g)\n",
        to_pct(selsel3), alpha2, pvalue(BinomialTest(selsel3, alpha2), tail=:right))
    end

    alpha = 0.01
    selmask = pval .< alpha
    rgd = Dict(region => selmask[(pval_pres .< presalpha) & [ismatch(Sessions.REGIONS[region], ch.site) for ch in chs]] for region in (:PHG, :EC, :H, :A))
    for region in (:PHG, :EC, :H, :A)
        selmask = pval .< alpha
        selrg = rgd[region]
        @printf("%s of units in %s show significant maintenance period activity at α = %s (p = %.2g)\n",
                to_pct(selrg), region, alpha, pvalue(BinomialTest(selrg, alpha), tail=:right))
    end
end

const FIGS2_UNITS = [
    "399e11sb-2007-7-21_16-34-42 LPHG ch46 clu1 MU",
    "405e23sb-2008-3-2_10-52-20 LPHG ch54 clu1 MU",
    "409e4sb-2008-7-23_16-8-10 LPHG ch1 clu1 SU",
    "409e19sb-2008-7-29_3-9-13 LPHG ch8 clu1 SU",
    "409e32sb-2008-8-1_3-2-50 LPHG ch6 clu1 MU",
    "409e32sb-2008-8-1_3-2-50 LPHG ch8 clu1 MU",
    "398e12sb-2007-5-16_15-55-53 REC ch2 clu2 SU",
    "396e33sb-2006-12-11_15-56-47 LMH ch43 clu2 SU",
    "399e11sb-2007-7-21_16-34-42 LAH ch36 clu2 MU",
    "399e11sb-2007-7-21_16-34-42 RA ch5 clu1 MU",
    "399e11sb-2007-7-21_16-34-42 RA ch7 clu2 SU",
    "399e11sb-2007-7-21_16-34-42 RA ch7 clu5 SU",
    "399e11sb-2007-7-21_16-34-42 RA ch8 clu3 SU",
    "407e13sb-2008-4-11_14-36-46 RA ch1 clu2 MU",
    "407e13sb-2008-4-11_14-36-46 RA ch5 clu2 SU"
]

function plot_stimulus_both(session, channel_number, cluster_number; highlight=nothing, out="pdf", small::Bool=false)
    small ? gcf()[:set_size_inches](2, 2.5) : gcf()[:set_size_inches](4.488, 3.5)
    fs = small ? 5 : 7
    nstim = length(session.stimuli)
    gs = gridspec.GridSpec(1, nstim, left=0.093, right=0.98, bottom=0.09, top=1.006, wspace=0.2, hspace=0.0)
    nthimg = 10
    nthraster = 3
    gss = Array{Any}(nstim)
    for stim = 1:nstim
        gs_item = gs[:__getitem__]((0, stim-1))
        inner_gs = gridspec.GridSpecFromSubplotSpec(nthimg, 1, subplot_spec=gs_item, hspace=0.02)
        subplot(inner_gs[:__getitem__](0))
        showstimulus(stim)
        gs_item = inner_gs[:__getitem__](slice(1, nthimg))
        gss[stim] = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=gs_item, hspace=0.3)
    end
    for (iepoch, epoch) in enumerate((:sample, :maintenance))
        const LIMS = epoch == :maintenance ? (0, 2400) : (-200, 800)
        spikes = session.spikes[channel_number][cluster_number]
        trials = stimulusframe(session.trials)
        ticks = epoch == :maintenance ? (0:1000:2000) : 0:600:1000
        local ax
        rasterbin = epoch == :maintenance ? 200 : 50
        if isa(highlight, Void)
            highlighttrials = Int[]
        else
            stims = session.trials[:stimuli].data::Vector{Vector{Int}}
            highlighttrials = find(Bool[highlight in x for x in stims])
        end

        for stim = 1:nstim
            inner_gs = gridspec.GridSpecFromSubplotSpec(nthraster, 1, subplot_spec=gss[stim][:__getitem__]((iepoch-1, 0)), hspace=0.0)
            axa = subplot(inner_gs[:__getitem__]((slice(0, nthraster-1))), zorder=1000)
            local usetrials::Vector{Int}
            if epoch == :sample
                stims = session.trials[:stimuli].data::Vector{Vector{Int}}
                usetrials = find([stim in x for x in stims])
                onsets = trials[:onset][trials[:stimulus] .== stim].data
            elseif epoch == :maintenance
                stims = session.trials[:stimuli].data::Vector{Vector{Int}}
                usetrials = find([stim in x for x in stims])
                onsets = session.trials[:maint][usetrials]
            end
            t1 = findin(usetrials, setdiff(usetrials, highlighttrials))
            t2 = findin(usetrials, highlighttrials)

            plotraster(spikes, onsets[[t1; t2]], LIMS)
            if !isa(highlight, Void)
                axhline(length(t1)+0.5, color=(1, 0, 0, 0.5), linewidth=0.5)
            end
            axis("off")
            axvline(0, ymin=-1/(nthraster-1), ymax=1, linestyle=":", color="darkcyan", dashes=(1, 1),
                       linewidth=0.5, clip_on=false)

            if stim == 1
                ylabel("Presentation Number")
                bax = ax = subplot(inner_gs[:__getitem__](nthraster-1))
            else
                gca()[:set_yticklabels](())
                bax = subplot(inner_gs[:__getitem__](nthraster-1), sharey=ax)
            end
            rast = raster(spikes, onsets, (LIMS[1]:rasterbin:LIMS[2])/1000)
            bar(LIMS[1]:rasterbin:LIMS[2]-rasterbin, mean(rast, 2)*1000/rasterbin, rasterbin, color="k", edgecolor="none", clip_on=false)
            if !isa(highlight, Void)
                bar(LIMS[1]:rasterbin:LIMS[2]-rasterbin, mean(rast[:, t1], 2)*1000/rasterbin,
                    rasterbin, color=(1, 0, 0, 0.5), edgecolor="none", clip_on=false)
            end
            if epoch == :maintenance
                xticks(ticks, ["0"; [string(Int(x/1000)) for x in ticks[2:end]]], fontsize=fs)
                stim == 1 && xlabel("Time After Mask Offset (Seconds)", fontsize=fs, labelpad=small ? 1 : 2, horizontalalignment="left", x=-0.6)
            else
                xticks(ticks, fontsize=fs)
                stim == 1 && xlabel("Time After Stimulus Onset (ms)", fontsize=fs, labelpad=small ? 1 : 2, horizontalalignment="left", x=-0.6)
            end
            if stim == 1
                ylabel("Firing Rate (Hz)", fontsize=fs, labelpad=small ? 1 : 3)
                removespines(["top", "right"])
                gca()[:yaxis][:set_ticks_position]("left")
                gca()[:xaxis][:set_ticks_position]("bottom")
            else
                bax[:axes][:get_yaxis]()[:set_visible](false)
                removespines(["left", "top", "right"])
                gca()[:xaxis][:set_ticks_position]("bottom")
            end
            tick_params(axis="y", pad=2)
            tick_params(axis="x", direction="out", length=3, pad=1)
            xlim(LIMS[1], LIMS[2])
        end
        sca(ax)
        yl = ax[:get_ylim]()
        yticks((0:max(div(yl[2], 3), 1):yl[2]), fontsize=fs)
    end
    savefig("$(session.folder)_ch$(channel_number)c$(cluster_number)_$(session.sites[channel_number])_$(session.unit_type[channel_number][cluster_number])_samplemaintenance$(small ? "_small" : "").$(out)", dpi=300)
end

function selunits_table()
    patientid(ch) = parse(Int, ch.session[1:3])
    function nperpatient(patients, chs)
        n = zeros(Int, length(patients))
        for ch in chs
            n[searchsortedfirst(patients, patientid(ch))] += 1
        end
        n
    end

    @load joinpath("processed", "permpval_maintenance.jld") chs pval
    allunit = chs[Summary.unit_mask(chs)]
    chpatient = map(patientid, allunit)
    vissel = Analysis.permselunits()
    maintsel = intersect(vissel, chs[pval .< 0.01])
    patients = sort(unique(chpatient))
    nallunit = nperpatient(patients, allunit)
    nvissel = nperpatient(patients, vissel)
    nmaintsel = nperpatient(patients, maintsel)
    push!(nallunit, sum(nallunit))
    push!(nvissel, sum(nvissel))
    push!(nmaintsel, sum(nmaintsel))
    println("Patient\t# Units\t# Vis\t# Maint\t% Vis\t\t% Maint")
    for i = 1:length(patients)+1
        visci = rcopy(R"binom.test($(nvissel[i]), $(nallunit[i]))$conf.int")
        if nvissel[i] == 0
            maintci = [NaN, NaN]
        else
            maintci = rcopy(R"binom.test($(nmaintsel[i]), $(nvissel[i]))$conf.int")
        end
        patient = i > length(patients) ? "TOTAL" : patients[i]
        @printf("%s\t%d\t%d\t%d\t%.0f%% (%.0f%%-%.0f%%)\t%.0f%% (%.0f%%-%.0f%%)\n",
        patient, nallunit[i], nvissel[i], nmaintsel[i],
        nvissel[i]./nallunit[i]*100, visci[1]*100, visci[2]*100,
        nmaintsel[i]./nvissel[i]*100, maintci[1]*100, maintci[2]*100)
    end
end

function meancortest(x, y, nperm)
    cors = map(cor, x, y)
    meancor = tanh(mean(atanh(cors)))
    permcor = zeros(nperm)
    pvalunit = zeros(length(x))
    for ix = 1:length(x)
        permx = copy(x[ix])
        cury = y[ix]
        unitgt = 0
        for iperm = 1:nperm
            shuffle!(permx)
            v = cor(permx, cury)
            permcor[iperm] += atanh(v)
            unitgt += v >= cors[ix]
        end
        pvalunit[ix] = unitgt/nperm
    end
    map!(x->tanh(x/length(cors)), permcor)
    (meancor, cors, mean(permcor .>= meancor) + 1/nperm, pvalunit)
end

function corscatter()
    wd = 6.5
    ht = 4.55
    regions = [:PHG, :EC, :H, :A]
    region_names = ["PHC", "EC", "Hippocampus", "Amygdala"]
    gcf()[:set_size_inches](wd, ht)
    subplots_adjust(wspace=0.35, hspace=0.4, left=0.072, right=0.995, bottom=0.08, top=0.965)
    selselch = intersect(Analysis.maintselunits(), Analysis.permselunits())
    @load joinpath("processed", "modulation_maintenance.jld") means chs
    sel = findin(chs, selselch)
    chsel = chs
    n = map(length, means[sel])
    @load joinpath("processed", "modulation_presentation.jld") stim_fr_win win chs baseline_fr_all baseline_fr_win
    sel2 = findin(chs, selselch)
    sel2region = [findfirst(regions, siteregion(x)) for x in chs[sel2]]
    sp = sortperm(selselch, lt=(x,y)-> siteregion(x) == siteregion(y) ? string(x.site)[1] < string(y.site)[1] :
                                                                        findfirst(regions, siteregion(x)) < findfirst(regions, siteregion(y)))
    sel = sel[sp]
    sel2 = sel2[sp]
    sel2region = sel2region[sp]
    @assert chsel[sel] == chs[sel2] == selselch[sp]
    fig3ch = findfirst(x->x.session == "399e11sb-2007-7-21_16-34-42" && x.channel == 8 && x.cluster == 2, selselch[sp])
    for (i, isel) = enumerate([1:fig3ch-1; fig3ch+1:length(sel); fig3ch])
        subplot(4, 6, i)
        x = stim_fr_win[sel2[isel]]./(win[sel2[isel]][2] - win[sel2[isel]][1])*1000
        y = means[sel[isel]]./(2400 - 300)*1000
        plot(x, y, color="k", marker="o", markerfacecolor="none", markeredgewidth=0.5, linestyle="none", markersize=2)
        rgx = maximum(x) - minimum(x)
        rgy = maximum(y) - minimum(y)
        stepx = rgx > 20 ? 10 : rgx > 10 ? 5 : rgx > 6 ? 3 : rgx > 6 ? 2 : rgx > 2 ? 1 : rgx > 1 ? 0.5 : rgx > 0.6 ? 0.3 : 0.1
        stepy = rgy > 10 ? 5 : rgy > 6 ? 2 : rgy > 3 ? 1 : rgy > 1.5 ? 0.5 : rgy > 0.6 ? 0.2 : rgy > 0.3 ? 0.1 : 0.05
        xticks(-100*stepx:stepx:100*stepx, fontsize=7)
        yticks(-100*stepy:stepy:100*stepy, fontsize=7)
        minx = minimum(x)-rgx/10
        maxx = maximum(x)+rgx/10
        miny = minimum(y)-rgy/10
        maxy = maximum(y)+rgy/10
        xlim(minx, maxx)
        ylim(miny, maxy)
        tick_params(axis="y", direction="out", pad=2, length=2)
        tick_params(axis="x", direction="out", pad=2, length=2)
        gca()[:yaxis][:set_ticks_position]("left")
        gca()[:xaxis][:set_ticks_position]("bottom")
        if (i - 1) % 6 == 0
            ylabel("Maintenance\nFiring Rate (Hz)", fontsize=7, labelpad=2)
        end
        if i > 18
            xlabel("Sample\nFiring Rate (Hz)", fontsize=7, labelpad=1)
        end
        left, bottom, right, top = gca()[:get_position]()[:extents]
        cortxt = "\$r = $(round(cor(x, y), 2))\$"
        idx = findfirst(FIGS2_UNITS, string(chsel[sel[isel]]))
        if idx != 0
            cortxt = "Fig. S2$('A'+idx-1)\n$cortxt"
        elseif isel == fig3ch
            cortxt = "Fig. 3\n$cortxt"
        end
        text(maxx-rgx/30, miny+rgy/50, cortxt, horizontalalignment="right", verticalalignment="baseline", fontsize=8)
        side = string(chs[sel2[isel]].site)[1] == 'L' ? "Left" : "Right"
        title(side*" "*region_names[sel2region[isel]], y=0.97, fontsize=8)
    end
    savefig("corscatter.pdf")
end

function responsestrength()
    function stimhist(stim_fr, baseline_fr)
        resps = [begin
            stimbs = stim .- bs
            bin = maximum(abs(stimbs)) == minimum(abs(stimbs)) ? indmax(stimbs) : indmax(abs(stimbs))
            stimbs[bin]
            end for (stim, bs) in zip(stim_fr, baseline_fr)]
        resps[!isfinite(resps)] = 0
        edges = -4:0.2:9
        respsc = copy(resps)
        respsc[resps .<= edges[1]] = edges[1]+0.01
        respsc[resps .>= edges[end]] = edges[end]-0.01
        h = fit(Histogram, respsc, edges, closed=:right)
        w = h.weights
        (resps, h.edges[1], w)
    end

    gcf()[:set_size_inches](3.5, 2)
    chs_maintenance = Analysis.maintselunits()
    chs_resp = Analysis.permselunits()
    @load joinpath(PROCESSED_DATA_DIR, "modulation_presentation.jld") chs stim_fr_win baseline_fr_win
    stim_fr = stim_fr_win
    baseline_fr = baseline_fr_win
    chsel = indexin(chs_resp, chs)
    vs, edge, w = stimhist(stim_fr[chsel], baseline_fr[chsel])
    bar(edge[1:end-1], w, step(edge), color=(0.5, 0.5, 0.5), linewidth=0.1)
    maintenancesel = findin(chs[chsel], chs_maintenance)
    selmod = vs[maintenancesel]
    nselmod = vs[setdiff(1:length(vs), maintenancesel)]
    selmod = sort(selmod)
    test = MannWhitneyUTest(abs(selmod), abs(nselmod))
    @printf "Median maintenance: %.1f; median no maintenance: %.1f; U = %.3g, p = %.2g\n" median(abs(selmod)) median(abs(nselmod)) test.U pvalue(test)
    showch = chs[chsel][sortperm(vs)]
    ylim(0.0, maximum(w)+1)
    selmod[selmod .>= edge[end]] = edge[end] - step(edge)/2
    selmod = sort!(selmod)
    y = zeros(length(selmod))
    for i = 2:length(selmod)
        i1 = i
        i2 = i-1
        while i2 > 0 && selmod[i1] - selmod[i2] < 0.1
            y[i] += 1
            i2 -= 1
        end
    end
    @assert length(y) == length(selmod)
    scatter(selmod, y/(maximum(w)+1)*25+1.5, edgecolor="k", facecolor="none", linewidth=0.5, clip_on=false, zorder=1000)
    removespines()
    gca()[:yaxis][:set_ticks_position]("left")
    gca()[:xaxis][:set_ticks_position]("bottom")

    xticks(Int(edge[1]):Int(edge[end]))
    xlim(first(edge), last(edge))
    ylabel("Number of Units")
    xlabel("Change in Firing Rate from Baseline (Number of Spikes)")
    tight_layout(pad=0.1)
    savefig("responsestrength.pdf")
end

function wm_correct(epoch=:maint)
    selch = intersect(Analysis.permselunits(), Analysis.maintselunits())
    jld = jldopen(joinpath(PROCESSED_DATA_DIR, "maintenance_topstim_correct.jld"))
    chs = read(jld, "chs")
    pval_mwu = read(jld, "pval_mwu_$epoch")
    sep = read(jld, "sep_$epoch")
    ng1 = read(jld, epoch == :persist ? "ng1_persist" : "ng1")
    ng2 = read(jld, epoch == :persist ? "ng2_persist" : "ng2")
    close(jld)
    alpha = 0.01

    for (title, mych) = [
            ("all visually selective", Analysis.permselunits()),
            ("visually selective, but not maintenance-selective",
            setdiff(Analysis.permselunits(), selch)),
            ("maintenance-selective", selch)
        ]
        println(title)
        sel = findin(chs, mych)
        sel = sel[(ng1[sel] .> 0) & (ng2[sel] .> 0)]
        @printf "  median %s correct trials, %s incorrect trials\n" median(ng1[sel]) median(ng2[sel])
        @printf("  at α = %s, %s significant (p = %.2g)\n", alpha,
                to_pct(pval_mwu[sel] .< alpha),
                pvalue(BinomialTest(pval_mwu[sel] .< alpha, alpha), tail=:right))
        test = OneSampleTTest(sep[sel])
        @printf("  population effect %.2f, t(%d) = %.3f, p = %.4g\n",
                test.xbar, test.df, test.t, pvalue(test))
        @printf "  bootstrap vs. 0 p = %.2g\n" 2*(min(Analysis.bsvs0(sep[sel], 1000000), Analysis.bsvs0(sep[sel], 1000000)))
    end
end
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
using Statistics

"""
    ReadData(f)

Reads a .mtx file that contains a
sparse undirected graph.

Returns the number of non zero elements per row
and the column index of the non zero elements.
"""
function ReadData(f)
    println("Reading: ",f)
    dfile = open(f)
    header = readline(dfile)
    # Skipping comments
    while true
        if header[1] == '%'
            header = readline(dfile)
        else
            break
        end
    end
    header =  map(x -> parse(Int32,x),split(header))
    col_indx = Array{Int32,1}(undef,header[3])
    csr_rows = zeros(Int32, header[1])
    for i  = 1:header[3]
        line = (readline(dfile))
        mid = findfirst(" ", line)[1]
        col_indx[i] = parse(Int32, line[1:(mid-1)])
        csr_rows[parse(Int32, line[(mid+1):end])] += 1
    end
    return csr_rows, col_indx
end

function SparseReadData(f)
    println("Reading: ",f)
    dfile = open(f)
    header = readline(dfile)
    # Skipping comments
    while true
        if header[1] == '%'
            header = readline(dfile)
        else
            break
        end
    end
    header =  map(x -> parse(Int32,x),split(header))
    col_indx = Array{Int32,1}(undef,header[3])
    row_indx = Array{Int32,1}(undef,header[3])
    for i  = 1:header[3]
        line = (readline(dfile))
        mid = findfirst(" ", line)[1]
        col_indx[i] = parse(Int32, line[1:(mid-1)])
        row_indx[i] = parse(Int32, line[(mid+1):end])
    end
    return row_indx, col_indx
end
function GetIntegralCsrRows(csr_rows)
    int_csr_rows = similar(csr_rows)
    int_csr_rows[1] = 1

    for i = 2:length(csr_rows)
        @inbounds int_csr_rows[i] = csr_rows[i-1] + int_csr_rows[i-1]
    end
    return int_csr_rows
end

#
# function SlowReadData(f="auto.mtx")
#     println("Slow Reading: ",f)
#     dfile = open(f)
#     header = map(x -> parse(Int32,x),split(readline(dfile)))
#     col_indx = Array{Int32,1}(undef,header[3])
#     csr_rows = zeros(Int32, header[1])
#     for i  = 1:header[3]
#         t = map(x -> parse(Int32,x),split(readline(dfile)))
#         col_indx[i] = t[1]
#         csr_rows[t[2]] += 1
#     end
#
#     return csr_rows, col_indx
# end

ext: jl


text:
using Dates
using Test

using MarketData

using TimeSeries

@testset "timearray" begin

@testset "field extraction methods work" begin
    @testset "timestamp, values, colnames and meta" begin
        for ta ∈ [cl, op, ohlc]
            @test timestamp(ta) isa Vector{Date}
            @test values(ta)    isa Array{Float64}
            @test colnames(ta)  isa Vector{Symbol}
            @test meta(mdata) == "Apple"
        end
    end
end

@testset "type constructors allow views" begin
    source_rows = 101:121
    source_cols = 1:size(values(AAPL), 2)
    tstamps = view(timestamp(AAPL), source_rows)
    tvalues = view(values(AAPL), source_rows, source_cols)

    AAPL1 = TimeArray(timestamp(AAPL)[source_rows],
                      values(AAPL)[source_rows, source_cols],
                      colnames(AAPL), meta(AAPL))

    AAPL2 = TimeArray(tstamps, tvalues, colnames(AAPL), meta(AAPL))

    @testset "match first date" begin
        @test timestamp(AAPL1)[1] == timestamp(AAPL2)[1]
    end

    @testset "match first values" begin
        @test values(AAPL1)[1] == values(AAPL2)[1]
    end

    @testset "match all values" begin
        @test values(AAPL1) == values(AAPL2)
    end
end

@testset "type constructors enforce invariants" begin
    mangled_stamp = vcat(timestamp(cl)[200:end], timestamp(cl)[1:199])
    dupe_stamp    = vcat(timestamp(cl)[1:499], timestamp(cl)[499])

    @testset "unequal length between values and timestamp fails" begin
        @test_throws(
            DimensionMismatch,
            TimeArray(timestamp(cl), values(cl)[2:end], [:Close]))
    end

    @testset "unequal length between colnames and array width fails" begin
        @test_throws(
            DimensionMismatch,
            TimeArray(timestamp(cl), values(cl), [:Close, :Open]))
    end

    @testset "duplicate timestamp values fails" begin
        @test_throws(
            ArgumentError,
            TimeArray(dupe_stamp, values(cl), [:Close]))
    end

    @testset "mangled order of timestamp values fails" begin
        @test_throws(
            ArgumentError,
            TimeArray(mangled_stamp, values(cl), [:Close]))
    end

    @testset "reverse occurs when needed" begin
        rev_timestamp = reverse(timestamp(cl), dims = 1)
        rev_values = reverse(values(cl), dims = 1)
        ta = TimeArray(rev_timestamp, rev_values, [:Close])
        @test timestamp(ta)[1] == Date(2000,1,3)
        @test values(ta)[1]    == 111.94
    end

    @testset "duplicate column names are enumerated by inner constructor" begin
        cols = [:a, :b, :c, :a, :a, :b, :d, :e, :e, :e, :e, :f]
        ta = TimeArray(timestamp(AAPL), values(AAPL), cols)
        @test colnames(ta)[1]  == :a
        @test colnames(ta)[2]  == :b
        @test colnames(ta)[3]  == :c
        @test colnames(ta)[4]  == :a_1
        @test colnames(ta)[5]  == :a_2
        @test colnames(ta)[6]  == :b_1
        @test colnames(ta)[7]  == :d
        @test colnames(ta)[8]  == :e
        @test colnames(ta)[9]  == :e_1
        @test colnames(ta)[10] == :e_2
        @test colnames(ta)[11] == :e_3
        @test colnames(ta)[12] == :f
    end

    @testset "and doesn't when unchecked" begin
        let
            ta = TimeArray(mangled_stamp, values(cl); unchecked = true)
            @test values(ta)    === values(cl)
            @test timestamp(ta) === mangled_stamp
        end

        let
            ta = TimeArray(dupe_stamp, values(cl); unchecked = true)
            @test timestamp(ta) === dupe_stamp
        end
    end
end

@testset "construction without colnames" begin
    one   = TimeArray(timestamp(cl), values(cl))
    multi = TimeArray(timestamp(AAPL), values(AAPL))
    more  = TimeArray(timestamp(cl)[1], collect(1:50)')

    @testset "default colnames" begin
        @test colnames(one)   == [:A]
        @test colnames(multi) == [:A, :B, :C, :D, :E, :F, :G, :H, :I, :J, :K, :L]
        @test colnames(more)  == [:A, :B, :C, :D, :E, :F, :G, :H, :I, :J, :K, :L, :M,
                                  :N, :O, :P, :Q, :R, :S, :T, :U, :V, :W, :X, :Y, :Z,
                                  :AA, :AB, :AC, :AD, :AE, :AF, :AG, :AH, :AI, :AJ, :AK,
                                  :AL, :AM, :AN, :AO, :AP, :AQ, :AR, :AS, :AT, :AU, :AV,
                                  :AW, :AX]
    end

    @testset "empty colnames forces meta to nothing" begin
        @test meta(one)   == nothing
        @test meta(multi) == nothing
        @test meta(more)  == nothing
    end
end

@testset "construct with StepRange{Date,Day}" begin
    drng = Date(2000,1,1):Day(1):Date(2000,1,5)
    ta = TimeArray(drng, 1:5)

    @test timestamp(ta)[1]   == first(drng)
    @test timestamp(ta)[end] == last(drng)

    @test values(ta)[1]   == 1
    @test values(ta)[end] == 5
end

@testset "construction from existing TimeArray" begin
    ts = Date(2018, 1, 1):Day(1):Date(2018, 1, 31)
    ta = TimeArray(ts, 1:31, [:x], :Meta)

    let
        ts′ = Date(2018, 3, 1):Day(1):Date(2018, 3, 31)
        ta′ = TimeArray(ta; timestamp = ts′)
        @test timestamp(ta′) == ts′
        @test values(ta′)    == values(ta)
        @test colnames(ta′)  == colnames(ta)
        @test meta(ta′)      == meta(ta)
    end

    let
        ta′ = TimeArray(ta; values = 2:32)
        @test timestamp(ta′) == timestamp(ta)
        @test values(ta′)    == 2:32
        @test colnames(ta′)  == colnames(ta)
        @test meta(ta′)      == meta(ta)
    end

    let
        ta′ = TimeArray(ta; colnames = [:y])
        @test timestamp(ta′) == timestamp(ta)
        @test values(ta′)    == values(ta)
        @test colnames(ta′)  == [:y]
        @test meta(ta′)      == meta(ta)
    end

    let
        ta′ = TimeArray(ta; meta = :Meta42)
        @test timestamp(ta′) == timestamp(ta)
        @test values(ta′)    == values(ta)
        @test colnames(ta′)  == colnames(ta)
        @test meta(ta′)      == :Meta42
    end
end

@testset "conversion methods" begin
    @testset "convert works " begin
        @test convert(TimeArray{Float64,1}, (cl.>op))                isa TimeArray{Float64,1}
        @test convert(TimeArray{Float64,2}, (merge(cl.<op, cl.>op))) isa TimeArray{Float64,2}
        @test convert(cl.>op)                                        isa TimeArray{Float64,1}
        @test convert(merge(cl.<op, cl.>op))                         isa TimeArray{Float64,2}
    end
end

@testset "copy methods" begin
    cop = copy(op)
    cohlc = copy(ohlc)

    @testset "copy works" begin
        @test timestamp(cop) == timestamp(op)
        @test values(cop)    == values(op)
        @test colnames(cop)  == colnames(op)
        @test meta(cop)      == meta(op)

        @test timestamp(cohlc) == timestamp(ohlc)
        @test values(cohlc)    == values(ohlc)
        @test colnames(cohlc)  == colnames(ohlc)
        @test meta(cohlc)      == meta(ohlc)
    end
end

@testset "index by integer works with both 1d and 2d time array" begin
    @testset "1d time array" begin
        ta = cl[1]
        @test timestamp(ta) == [Date(2000,1,3)]
        @test values(ta)    == [111.94]
        @test colnames(ta)  == [:Close]
        @test meta(ta)      == "AAPL"
    end

    @testset "2d time array" begin
        ta = ohlc[1]
        @test timestamp(ta) == [Date(2000,1,3)]
        @test values(ta)    == [104.88 112.5 101.69 111.94]
        @test colnames(ta)  == [:Open, :High, :Low, :Close]
        @test meta(ta)      == "AAPL"
    end
end

@testset "ordered collection methods" begin
    @testset "iteration protocol is valid" begin
        let  # single column
            i = 1
            for (t, val) ∈ cl[1:3]
                @test t   == timestamp(cl)[i]
                @test val == values(cl)[i]
                i += 1
            end
        end

        let  # multiple column
            i = 1
            for (t, val) ∈ ohlc[1:3]
                @test t   == timestamp(ohlc)[i]
                @test val == values(ohlc)[i, :]
                i += 1
            end
        end
    end

    @testset "end keyword returns correct index" begin
        @test timestamp(ohlc[end])[1] == timestamp(ohlc)[end]
    end

    @testset "getindex on single Int and Date" begin
        @test timestamp(ohlc[1])              == [Date(2000,1,3)]
        @test timestamp(ohlc[Date(2000,1,3)]) == [Date(2000,1,3)]
    end

    @testset "getindex on array of Int and Date" begin
        @test timestamp(ohlc[[1,10]])                           == [Date(2000,1,3), Date(2000,1,14)]
        @test timestamp(ohlc[[Date(2000,1,3),Date(2000,1,14)]]) == [Date(2000,1,3), Date(2000,1,14)]
    end

    @testset "getindex on range of Int and Date" begin
        irng = Int8(1):Int8(2):Int8(4)
        drng = Date(2000,1,3):Day(1):Date(2000,1,4)
        @test timestamp(ohlc[1:2])   == [Date(2000,1,3), Date(2000,1,4)]
        @test timestamp(ohlc[1:2:4]) == [Date(2000,1,3), Date(2000,1,5)]
        @test timestamp(ohlc[irng])  == [Date(2000,1,3), Date(2000,1,5)]
        @test timestamp(ohlc[drng])  == [Date(2000,1,3), Date(2000,1,4)]
    end

    @testset "getindex on range of DateTime when only Date is in timestamp" begin
        @test_throws(
            MethodError,
            ohlc[DateTime(2000,1,3,0,0,0)])
        @test_throws(
            MethodError,
            ohlc[[DateTime(2000,1,3,0,0,0),DateTime(2000,1,14,0,0,0)]])
        @test_throws(
            MethodError,
            ohlc[DateTime(2000,1,3,0,0,0):Day(1):DateTime(2000,1,4,0,0,0)])
    end

    @testset "getindex on range of Date" begin
        @test length(cl[Date(2000,1,1):Day(1):Date(2001,12,31)]) == 500
    end

    @testset "getindex on single column name" begin
        idx = Date(2000,1,3):Day(1):Date(2000,1,14)
        @test size(values(ohlc[:Open]), 2)      == 1
        @test size(values(ohlc[:Open][idx]), 1) == 10
    end

    @testset "getindex on multiple column name" begin
        @test values(ohlc[:Open, :Close])[1]   == 104.88
        @test values(ohlc[:Open, :Close])[2]   == 108.25
        @test values(ohlc[:Open, :Close])[501] == 111.94
    end

    @testset "getindex on 1d returns 1d object" begin
        @test cl[1]   isa TimeArray{Float64,1}
        @test cl[1:2] isa TimeArray{Float64,1}
    end

    @testset "getindex on a 1d Boolean TimeArray returns appropriate rows" begin
        @test values(ohlc[op .> cl][2])             == values(ohlc[4])
        @test timestamp(ohlc[op[300:end] .> cl][2]) == timestamp(ohlc[303])
        # MethodError, Bool must be 1D-TimeArray
        @test_throws MethodError ohlc[merge(op .> cl, op .< cl)]
    end

    @testset "Base.eachindex" begin
        @test eachindex(cl)   == 1:length(cl)
        @test eachindex(ohlc) == 1:length(ohlc)
    end
end

@testset "Base.ndims" begin
    @test ndims(cl)   == 1
    @test ndims(ohlc) == 2
end

@testset "Base.eltype" begin
    @test eltype(cl)   == Tuple{Date,Float64}
    @test eltype(ohlc) == Tuple{Date,Vector{Float64}}

    @test eltype(cl .> 1)   == Tuple{Date,Bool}
    @test eltype(ohlc .> 1) == Tuple{Date,Vector{Bool}}
end

@testset "Base.collect" begin
    @testset "cl" begin
        A = collect(cl)
        @test A[1] == (Date(2000, 01, 03), 111.94)
        @test A[2] == (Date(2000, 01, 04), 102.5)
        @test A[3] == (Date(2000, 01, 05), 104.0)
    end

    @testset "ohlc" begin
        A = collect(ohlc)
        @test A[1] == (Date(2000, 01, 03), [104.88, 112.5, 101.69, 111.94])
        @test A[2] == (Date(2000, 01, 04), [108.25, 110.62, 101.19, 102.5])
        @test A[3] == (Date(2000, 01, 05), [103.75, 110.56, 103.0, 104.0])
    end
end

@testset "Base.size" begin
    @test size(ohlc)    == (500, 4)
    @test size(ohlc, 1) == 500
    @test size(ohlc, 2) == 4

    @test size(cl)    == (500,)
    @test size(cl, 1) == 500
    @test size(cl, 2) == 1
end

@testset "equal" begin
    @test cl == copy(cl)
    @test cl ≠ ohlc  # rely on fallback definition
    @test cl ≠ lag(cl)

    @test  isequal(cl, copy(cl))
    @test !isequal(cl, ohlc)
    @test !isequal(cl, lag(cl))

    ds = DateTime(2017, 12, 25):Day(1):DateTime(2017, 12, 31) |> collect

    let  # diff colnames
        x = TimeArray(ds, 1:7, [:foo])
        y = TimeArray(ds, 1:7, [:bar])
        @test x != y
    end

    let  # Float vs Int
        x = TimeArray(ds, 1:7)
        y = TimeArray(ds, 1.0:7)
        @test x == y
    end

    let  # Date vs DateTime
        ds2 = Date(2017, 12, 25):Day(1):Date(2017, 12, 31) |> collect
        x = TimeArray(ds,  1:7, [:foo], :bar)
        y = TimeArray(ds2, 1:7, [:foo], :bar)
        @test x == y
    end

    let  # diff meta
        x = TimeArray(ds, 1:7, [:foo], :bar)
        y = TimeArray(ds, 1:7, [:foo], :baz)
        @test x != y
    end

    @testset "hash" begin
        @test hash(cl)   == hash(copy(cl))
        @test hash(ohlc) == hash(copy(ohlc))

        d = Dict(cl => 42)
        @test d[cl] == 42
        @test d[copy(cl)] == 42

        d = Dict(ohlc => 24)
        @test d[ohlc] == 24
        @test d[copy(ohlc)] == 24
    end
end

@testset "show methods don't throw errors" begin
    let str = sprint(show, cl)
        out = """500×1 TimeArray{Float64,1,Date,Array{Float64,1}} 2000-01-03 to 2001-12-31
│            │ Close  │
├────────────┼────────┤
│ 2000-01-03 │ 111.94 │
│ 2000-01-04 │ 102.5  │
│ 2000-01-05 │ 104.0  │
│ 2000-01-06 │ 95.0   │
│ 2000-01-07 │ 99.5   │
│ 2000-01-10 │ 97.75  │
│ 2000-01-11 │ 92.75  │
│ 2000-01-12 │ 87.19  │
│ 2000-01-13 │ 96.75  │
   ⋮
│ 2001-12-19 │ 21.62  │
│ 2001-12-20 │ 20.67  │
│ 2001-12-21 │ 21.0   │
│ 2001-12-24 │ 21.36  │
│ 2001-12-26 │ 21.49  │
│ 2001-12-27 │ 22.07  │
│ 2001-12-28 │ 22.43  │
│ 2001-12-31 │ 21.9   │"""
        @test str == out
    end

    let str = sprint(show, ohlc)
        out = """500×4 TimeArray{Float64,2,Date,Array{Float64,2}} 2000-01-03 to 2001-12-31
│            │ Open   │ High   │ Low    │ Close  │
├────────────┼────────┼────────┼────────┼────────┤
│ 2000-01-03 │ 104.88 │ 112.5  │ 101.69 │ 111.94 │
│ 2000-01-04 │ 108.25 │ 110.62 │ 101.19 │ 102.5  │
│ 2000-01-05 │ 103.75 │ 110.56 │ 103.0  │ 104.0  │
│ 2000-01-06 │ 106.12 │ 107.0  │ 95.0   │ 95.0   │
│ 2000-01-07 │ 96.5   │ 101.0  │ 95.5   │ 99.5   │
│ 2000-01-10 │ 102.0  │ 102.25 │ 94.75  │ 97.75  │
│ 2000-01-11 │ 95.94  │ 99.38  │ 90.5   │ 92.75  │
│ 2000-01-12 │ 95.0   │ 95.5   │ 86.5   │ 87.19  │
│ 2000-01-13 │ 94.48  │ 98.75  │ 92.5   │ 96.75  │
   ⋮
│ 2001-12-19 │ 20.58  │ 21.68  │ 20.47  │ 21.62  │
│ 2001-12-20 │ 21.4   │ 21.47  │ 20.62  │ 20.67  │
│ 2001-12-21 │ 21.01  │ 21.54  │ 20.8   │ 21.0   │
│ 2001-12-24 │ 20.9   │ 21.45  │ 20.9   │ 21.36  │
│ 2001-12-26 │ 21.35  │ 22.3   │ 21.14  │ 21.49  │
│ 2001-12-27 │ 21.58  │ 22.25  │ 21.58  │ 22.07  │
│ 2001-12-28 │ 21.97  │ 23.0   │ 21.96  │ 22.43  │
│ 2001-12-31 │ 22.51  │ 22.66  │ 21.83  │ 21.9   │"""
        @test str == out
    end

    let str = sprint(show, AAPL)
        out = """8336×12 TimeArray{Float64,2,Date,Array{Float64,2}} 1980-12-12 to 2013-12-31
│            │ Open   │ High   │ Low    │ Close  │ Volume    │ Ex-Dividend │
├────────────┼────────┼────────┼────────┼────────┼───────────┼─────────────┤
│ 1980-12-12 │ 28.75  │ 28.88  │ 28.75  │ 28.75  │ 2.0939e6  │ 0.0         │
│ 1980-12-15 │ 27.38  │ 27.38  │ 27.25  │ 27.25  │ 785200.0  │ 0.0         │
│ 1980-12-16 │ 25.38  │ 25.38  │ 25.25  │ 25.25  │ 472000.0  │ 0.0         │
│ 1980-12-17 │ 25.88  │ 26.0   │ 25.88  │ 25.88  │ 385900.0  │ 0.0         │
│ 1980-12-18 │ 26.62  │ 26.75  │ 26.62  │ 26.62  │ 327900.0  │ 0.0         │
│ 1980-12-19 │ 28.25  │ 28.38  │ 28.25  │ 28.25  │ 217100.0  │ 0.0         │
│ 1980-12-22 │ 29.62  │ 29.75  │ 29.62  │ 29.62  │ 166800.0  │ 0.0         │
│ 1980-12-23 │ 30.88  │ 31.0   │ 30.88  │ 30.88  │ 209600.0  │ 0.0         │
│ 1980-12-24 │ 32.5   │ 32.62  │ 32.5   │ 32.5   │ 214300.0  │ 0.0         │
   ⋮
│ 2013-12-19 │ 549.5  │ 550.0  │ 543.73 │ 544.46 │ 1.14396e7 │ 0.0         │
│ 2013-12-20 │ 545.43 │ 551.61 │ 544.82 │ 549.02 │ 1.55862e7 │ 0.0         │
│ 2013-12-23 │ 568.0  │ 570.72 │ 562.76 │ 570.09 │ 1.79038e7 │ 0.0         │
│ 2013-12-24 │ 569.89 │ 571.88 │ 566.03 │ 567.67 │ 5.9841e6  │ 0.0         │
│ 2013-12-26 │ 568.1  │ 569.5  │ 563.38 │ 563.9  │ 7.286e6   │ 0.0         │
│ 2013-12-27 │ 563.82 │ 564.41 │ 559.5  │ 560.09 │ 8.0673e6  │ 0.0         │
│ 2013-12-30 │ 557.46 │ 560.09 │ 552.32 │ 554.52 │ 9.0582e6  │ 0.0         │
│ 2013-12-31 │ 554.17 │ 561.28 │ 554.0  │ 561.02 │ 7.9673e6  │ 0.0         │

│            │ Split Ratio │ Adj. Open │ Adj. High │ Adj. Low │ Adj. Close │
├────────────┼─────────────┼───────────┼───────────┼──────────┼────────────┤
│ 1980-12-12 │ 1.0         │ 3.3766    │ 3.3919    │ 3.3766   │ 3.3766     │
│ 1980-12-15 │ 1.0         │ 3.2157    │ 3.2157    │ 3.2004   │ 3.2004     │
│ 1980-12-16 │ 1.0         │ 2.9808    │ 2.9808    │ 2.9655   │ 2.9655     │
│ 1980-12-17 │ 1.0         │ 3.0395    │ 3.0536    │ 3.0395   │ 3.0395     │
│ 1980-12-18 │ 1.0         │ 3.1264    │ 3.1417    │ 3.1264   │ 3.1264     │
│ 1980-12-19 │ 1.0         │ 3.3179    │ 3.3331    │ 3.3179   │ 3.3179     │
│ 1980-12-22 │ 1.0         │ 3.4788    │ 3.494     │ 3.4788   │ 3.4788     │
│ 1980-12-23 │ 1.0         │ 3.6267    │ 3.6408    │ 3.6267   │ 3.6267     │
│ 1980-12-24 │ 1.0         │ 3.817     │ 3.8311    │ 3.817    │ 3.817      │
   ⋮
│ 2013-12-19 │ 1.0         │ 546.2492  │ 546.7463  │ 540.5133 │ 541.239    │
│ 2013-12-20 │ 1.0         │ 542.2033  │ 548.3467  │ 541.5969 │ 545.7721   │
│ 2013-12-23 │ 1.0         │ 564.6398  │ 567.3437  │ 559.4308 │ 566.7174   │
│ 2013-12-24 │ 1.0         │ 566.5186  │ 568.4968  │ 562.6814 │ 564.3117   │
│ 2013-12-26 │ 1.0         │ 564.7392  │ 566.1309  │ 560.0471 │ 560.564    │
│ 2013-12-27 │ 1.0         │ 560.4845  │ 561.071   │ 556.1901 │ 556.7766   │
│ 2013-12-30 │ 1.0         │ 554.1621  │ 556.7766  │ 549.0525 │ 551.2395   │
│ 2013-12-31 │ 1.0         │ 550.8916  │ 557.9595  │ 550.7226 │ 557.7011   │"""
        @test str == out
    end

    let str = sprint(show, ohlc[1:4])
        out = """4×4 TimeArray{Float64,2,Date,Array{Float64,2}} 2000-01-03 to 2000-01-06
│            │ Open   │ High   │ Low    │ Close  │
├────────────┼────────┼────────┼────────┼────────┤
│ 2000-01-03 │ 104.88 │ 112.5  │ 101.69 │ 111.94 │
│ 2000-01-04 │ 108.25 │ 110.62 │ 101.19 │ 102.5  │
│ 2000-01-05 │ 103.75 │ 110.56 │ 103.0  │ 104.0  │
│ 2000-01-06 │ 106.12 │ 107.0  │ 95.0   │ 95.0   │"""
        @test str == out
    end

    let str = sprint(show, ohlc[1:0])
        @test str == "0×4 TimeArray{Float64,2,Date,Array{Float64,2}}"
    end

    let str = sprint(show, TimeArray(Date[], []))
        @test str == "0×1 TimeArray{Any,1,Date,Array{Any,1}}"
    end

    let str = sprint(show, lag(cl[1:2], padding=true))
        out = """2×1 TimeArray{Float64,1,Date,Array{Float64,1}} 2000-01-03 to 2000-01-04
│            │ Close  │
├────────────┼────────┤
│ 2000-01-03 │ NaN    │
│ 2000-01-04 │ 111.94 │"""
        @test str == out
    end
end  # @testset "show methods don't throw errors"

@testset "getproperty" begin
    let
        @test cl.Close  == cl[:Close]
        @test_throws KeyError cl.NotFound
    end

    let
        @test ohlc.Open  == ohlc[:Open]
        @test ohlc.High  == ohlc[:High]
        @test ohlc.Low   == ohlc[:Low]
        @test ohlc.Close == ohlc[:Close]
        @test_throws KeyError ohlc.NotFound
    end
end

end  # @testset "timearray"


text:
using BenchmarkTools, Plots, Plots.PlotMeasures, LaTeXStrings
using BenchmarkTools: minimum, median

SUITE = BenchmarkGroup()
model = "GEARBOX"
cases = ["GRBX01-MES01", "GRBX02-MES01", "GRBX01-MES01-discrete", "GRBX02-MES01-discrete"]

SUITE[model] = BenchmarkGroup()

include("gearbox.jl")
validation = []

# Template directions
const boxdirs = BoxDirections{Float64, Vector{Float64}}(6);
const octdirs = CustomDirections([Vector(vi) for vi in OctDirections(6)]);

vx, vy, px, py, n, θ = 1, 2, 3, 4, 6, 0.628318530717959
vec1 = sparsevec([px, py], [tan(θ), 1.], n)
vec2 = sparsevec([px, py], [tan(θ), -1.], n)
vec1perp = sparsevec([px, py], [1, -tan(θ)], n)
vec2perp = sparsevec([px, py], [1., tan(θ)], n)
dirs_plane_34 = [vec1, -vec1, vec2, -vec2, vec1perp, -vec1perp, vec2perp, -vec2perp]

vec3 = sparsevec([vx, vy], [-sin(θ), -cos(θ)], n)
vec4 = sparsevec([vx, vy], [-sin(θ), cos(θ)], n)
vec3perp = sparsevec([vx, vy], [cos(θ), -sin(θ)], n)
vec4perp = sparsevec([vx, vy], [cos(θ), sin(θ)], n)

dirs_plane_12 = [vec3, -vec3, vec4, -vec4, vec3perp, -vec3perp, vec4perp, -vec4perp]
const extdirs = vcat(collect(boxdirs), Vector.(dirs_plane_34), Vector.(dirs_plane_12)) |> CustomDirections

# Initial states (non-homogeneous system)
X0_GRBX01 = Hyperrectangle(low=[0, 0, -0.0168, 0.0029, 0, 0], high=[0, 0, -0.0166, 0.0031, 0, 0])
X0_GRBX02 = Hyperrectangle(low=[0, 0, -0.01675, 0.00285, 0, 0], high=[0, 0, -0.01665, 0.00315, 0, 0]);

# Initial states (homogeneized system)
X0_GRBX01h = Hyperrectangle(low=[0, 0, -0.0168, 0.0029, 0, 1.], high=[0, 0, -0.0166, 0.0031, 0, 1.])
X0_GRBX02h = Hyperrectangle(low=[0, 0, -0.01675, 0.00285, 0, 1.], high=[0, 0, -0.01665, 0.00315, 0, 1.]);

thull = TemplateHullIntersection(extdirs);

# ----------------------------------------
#  GRBX01 (dense time)
# ----------------------------------------
prob_GRBX01 = gearbox_homog(X0=X0_GRBX01h)

# GRBX01-MES01
sol_GRBX01 = solve(prob_GRBX01, max_jumps=100,
                   clustering_method=LazyClustering(1, convex=false),
                   intersect_source_invariant=true,
                   intersection_source_invariant_method=thull,
                   intersection_method=thull,
                   tspan = 0 .. 0.21,
                   fixpoint_check=true,
                   alg=LGG09(δ=0.0005, template=extdirs, approx_model=Forward()))
property = req1(sol_GRBX01) && req2(sol_GRBX01)
push!(validation, Int(property))
SUITE[model][cases[1]] = @benchmarkable solve($prob_GRBX01, max_jumps=100,
                   clustering_method=LazyClustering(1, convex=false),
                   intersect_source_invariant=true,
                   intersection_source_invariant_method=$thull,
                   intersection_method=$thull,
                   tspan = 0 .. 0.21,
                   fixpoint_check=true,
                   alg=LGG09(δ=0.0005, template=$extdirs, approx_model=Forward()))

# ----------------------------------------
#  GRBX02 (dense time)
# ----------------------------------------
prob_GRBX02 = gearbox_homog(X0=X0_GRBX02h)

# GRBX02-MES01
sol_GRBX02 = solve(prob_GRBX02, max_jumps=100,
                   clustering_method=LazyClustering(1, convex=false),
                   intersect_source_invariant=true,
                   intersection_source_invariant_method=thull,
                   intersection_method=thull,
                   tspan = 0 .. 0.21,
                   fixpoint_check=true,
                   alg=LGG09(δ=0.0008, template=extdirs, approx_model=Forward()))
property = req1(sol_GRBX02) && req2(sol_GRBX02)
push!(validation, Int(property))
SUITE[model][cases[2]] = @benchmarkable solve($prob_GRBX02, max_jumps=100,
                   clustering_method=LazyClustering(1, convex=false),
                   intersect_source_invariant=true,
                   intersection_source_invariant_method=$thull,
                   intersection_method=$thull,
                   tspan = 0 .. 0.21,
                   fixpoint_check=true,
                   alg=LGG09(δ=0.0008, template=$extdirs, approx_model=Forward()))

# ----------------------------------------
#  GRBX01 (discrete time)
# ----------------------------------------
prob_GRBX01 = gearbox_homog(X0=X0_GRBX01h)

# GRBX01-MES01-discrete
sol_GRBX01d = solve(prob_GRBX01, max_jumps=100,
                   clustering_method=LazyClustering(1, convex=false),
                   intersect_source_invariant=true,
                   intersection_source_invariant_method=thull,
                   intersection_method=thull,
                   tspan = 0 .. 0.21,
                   fixpoint_check=true,
                   alg=LGG09(δ=0.0001, template=extdirs, cache=true, approx_model=NoBloating()))
property = req1(sol_GRBX01d) && req2(sol_GRBX01d)
push!(validation, Int(property))
SUITE[model][cases[3]] = @benchmarkable solve($prob_GRBX01, max_jumps=100,
                   clustering_method=LazyClustering(1, convex=false),
                   intersect_source_invariant=true,
                   intersection_source_invariant_method=$thull,
                   intersection_method=$thull,
                   tspan = 0 .. 0.21,
                   fixpoint_check=true,
                   alg=LGG09(δ=0.0001, template=extdirs, cache=true, approx_model=NoBloating()))

# ----------------------------------------
#  GRBX02 (discrete time)
# ----------------------------------------
prob_GRBX02 = gearbox_homog(X0=X0_GRBX02h)

# GRBX02-MES01-discrete
sol_GRBX02d = solve(prob_GRBX02, max_jumps=100,
                   clustering_method=LazyClustering(1, convex=false),
                   intersect_source_invariant=true,
                   intersection_source_invariant_method=thull,
                   intersection_method=thull,
                   tspan = 0 .. 0.21,
                   fixpoint_check=true,
                   alg=LGG09(δ=0.0001, template=extdirs, cache=true, approx_model=NoBloating()))
property = req1(sol_GRBX02d) && req2(sol_GRBX02d)
push!(validation, Int(property))
SUITE[model][cases[4]] = @benchmarkable solve($prob_GRBX02, max_jumps=100,
                   clustering_method=LazyClustering(1, convex=false),
                   intersect_source_invariant=true,
                   intersection_source_invariant_method=$thull,
                   intersection_method=$thull,
                   tspan = 0 .. 0.21,
                   fixpoint_check=true,
                   alg=LGG09(δ=0.0001, template=extdirs, cache=true, approx_model=NoBloating()))

# ==============================================================================
# Execute benchmarks and save benchmark results
# ==============================================================================

# tune parameters
tune!(SUITE)

# run the benchmarks
results = run(SUITE, verbose=true)

# return the sample with the smallest time value in each test
println("minimum time for each benchmark:\n", minimum(results))

# return the median for each test
println("median time for each benchmark:\n", median(results))

# export runtimes
runtimes = Dict()
for (i, c) in enumerate(cases)
   t = median(results[model][c]).time * 1e-9
   runtimes[c] = round(t, digits=4)
end

for (i, c) in enumerate(cases)
   print(io, "JuliaReach, $model, $c, $(validation[i]), $(runtimes[c])\n")
end

# ==============================================================================
# Plot
# ==============================================================================

fig = Plots.plot()
Plots.plot!(fig, sol_GRBX01, vars=(3, 4), ε=1e-5,
           color=:blue, alpha=0.5, lw=1.0, linecolor=:black,
           tickfont=font(30, "Times"), guidefontsize=45,
           xlab=L"x_3",
           ylab=L"x_4",
           xtick=([-0.016, -0.009, -0.002],
                  [L"-0.016", L"-0.009", L"-0.002"]),
           ytick=([-0.008, -0.004, 0.0, 0.004],
                  [L"-0.008", L"-0.004", L"0", L"0.004"]),
           xlims=(-0.017, -0.0015), ylims=(-0.008, 0.004),
           bottom_margin=6mm, left_margin=5mm, right_margin=12mm, top_margin=3mm,
           size=(1000, 1000))
savefig(fig, joinpath(@__DIR__, "ARCH-COMP21-JuliaReach-Gearbox_GRBX01.pdf"))

sol_GRBX01 = nothing
sol_GRBX02 = nothing
GC.gc()


text:
module runtests

tests = [## "ttconstraints.jl",
         ## "moments.jl",
         ## "ttinvestments.jl",
         ## "portfolio.jl",
         "ttutils.jl",
         ## "ttstrategies.jl",
         "ttmuSigmaFunctions.jl",
         "ttuniverseModels.jl"
         ]

for t in tests
    include(string(Pkg.dir("AssetMgmt"), "/test/", t))
end

end


text:
"""
    Discontinuity(t, order::Int)

Object of discontinuity of order `order` at time `t`, i.e. discontinuity of
`order`th derivative at time `t`.
"""
struct Discontinuity{tType}
    t::tType
    order::Int
end

# ordering of discontinuities
Base.:<(a::Discontinuity, b::Discontinuity) = a.t < b.t || (a.t == b.t && a.order < b.order)
Base.isless(a::Discontinuity, b::Discontinuity) = isless(a.t, b.t) || (isequal(a.t, b.t) && a.order < b.order)
Base.isequal(a::Discontinuity, b::Discontinuity) = isequal(a.t, b.t) && a.order == b.order
Base.:(==)(a::Discontinuity, b::Discontinuity) = a.t == b.t && a.order == b.order

# ordering with numbers
Base.:<(a::Discontinuity, b::Number) = a.t < b
Base.:<(a::Number, b::Discontinuity) = a < b.t
Base.isless(a::Discontinuity, b::Number) = isless(a.t, b)
Base.isless(a::Number, b::Discontinuity) = isless(a, b.t)
Base.:(==)(a::Discontinuity, b::Number) = a.t == b
Base.:(==)(a::Number, b::Discontinuity) = a == b.t
Base.isequal(a::Discontinuity, b::Number) = isequal(a.t, b)
Base.isequal(a::Number, b::Discontinuity) = isequal(a, b.t)

# multiplication with numbers
Base.:*(a::Discontinuity, b::Number) = a.t * b
Base.:*(a::Number, b::Discontinuity) = a * b.t

# conversion to numbers
Base.convert(::Type{T}, d::Discontinuity) where {T<:Number} = convert(T, d.t)
Base.convert(::Type{T}, d::Discontinuity{T}) where {T<:Number} = d.t

# display
Base.show(io::IO, d::Discontinuity) = print(io, d.t, " (order ", d.order, ")")
Base.show(io::IO, ::MIME"text/plain", d::Discontinuity{T}) where {T} =
    print(io, "Discontinuity{", T, "}:\n   ", d)


text:
include("textureatlas.jl")

isnewline(x) = x == '\n'

type Text
    data
    text
    atlas
    cursors
    # default values
    font
    scale
    offset
    color
    startposition
    lineheight
end

function Sprite{N, T}(char, position::Point{N, T}, text)
    Sprite(
        char, position, text.scale, text.offset,
        text.color, text.font, text.atlas
    )
end
function Sprite{N, T}(
        char, position::Point{N, T}, scale, offset, color,
        font = defaultfont(),  atlas = get_texture_atlas()
    )
    Sprite{N, T}(
        position,
        glyph_bearing!(atlas, char, font, scale) + offset,
        glyph_scale!(atlas, char, font, scale),
        glyph_uv_width!(atlas, char, font),
        color
    )
end

function nextposition(sprite::Sprite, char, text)
    advance_x, advance_y = glyph_advance!(text.atlas, char, text.font, text.scale)
    position = sprite.position
    if isnewline(char)
        return Point2f0(text.startposition[1], position[2] - advance_y * text.lineheight) #reset to startx
    else
        return position + Point2f0(advance_x, 0)
    end
end

function printat(text::Text, idx::Integer, char::Char)
    position = if checkbounds(Bool, text.data, idx)
        sprite = text.data[idx]
        nextposition(sprite, text.text[idx], text)
    else
        text.startposition
    end
    nextsprite = Sprite(char, position, text)
    idx += 1
    insert!(text.data, idx, nextsprite)
    insert!(text.text, idx, char)
    idx
end

function printat(text::Text, idx::Int, str::String)
    sprite = text.data[idx]
    position = sprite.position
    for char in str
        char == '\r' && continue # stupid windows!
        idx = printat(text, idx, char)
    end
    idx
end

function Base.print(text::Text, char::Union{Char, String})
    map!(text.cursors, text.cursors) do idx
        idx = printat(text, idx, char)
        return idx
    end
    nothing
end
Base.String(text::Text) = join(text.text)


text:
using Unitful

import PipeLayout: weymouth

@testset "unit consistency in pressure loss" begin
    L = 100.0u"km"
    D = 1.0u"m"
    q = 100.0u"kg/s"
    p_i = 80.0u"bar"
    C = L/D^5*weymouth

    # pi^2 - po^2 = C*q|q|,  pressure should drop
    @test p_i^2 - C*q*abs.(q) <= p_i^2
end


text:
module ScalarTest
using Base.Test
using Quantum

xb = basis(:X,[1:5])
yb = basis(:Y,[1:5])

@test (1+(xb[1]'*yb[1])).ex.args==[:+, 1, (xb[1]'*yb[1])]
for op=(:+, :-, :*, :/)
	@test @eval qeval((b,k)->b==yb[1]' ? 2.5 : 3.4, ($op)(xb[1]'*yb[1], yb[1]'*xb[1])) == ($op)(3.4,2.5)
end

end


text:
using Pkg; Pkg.add("PlutoSliderServer")

using PlutoSliderServer, ClimateModels
pth=dirname(pathof(ClimateModels))

## Hector

if false
fil=joinpath(pth,"..","examples","Hector.jl")
PlutoSliderServer.export_notebook(fil)

tmp1=readdir(tempdir())
tst1=[isfile(joinpath(tempdir(),i,"hector","src","hector")) for i in tmp1]
tst2=[!islink(joinpath(tempdir(),i,"hector","src","hector")) for i in tmp1]
ii=findall(tst1.*tst2)[1]
symlink(joinpath(tempdir(),tmp1[ii],"hector","src","hector"),"hector")
end

## Speedy

if false
fil=joinpath(pth,"..","examples","Speedy.jl")
PlutoSliderServer.export_notebook(fil)

tmp1=readdir(tempdir())
tst1=[isfile(joinpath(tempdir(),i,"bin","speedy")) for i in tmp1]
ii=findall(tst1)[1]
symlink(joinpath(tempdir(),tmp1[ii],"bin","speedy"),"speedy")
end

## MITgcm

fil=joinpath(pth,"..","examples","MITgcm.jl")
PlutoSliderServer.export_notebook(fil)

max_stars_repo_name: crbinz/GlobalOceanNotebooks

max_issues_repo_name: crbinz/GlobalOceanNotebooks

max_forks_repo_name: crbinz/GlobalOceanNotebooks


text:
module ACEatoms

using StaticArrays, Reexport

include("ext_imports.jl")

include("configs.jl")

include("species_1pbasis.jl")

include("pop_1pbasis.jl")

include("pairpots/pair.jl")
@reexport using ACEatoms.PairPotentials

include("electrostatics.jl")

include("utils.jl")

include("siteenergy.jl")

# DONT UNCOMMENT ANYTHING BELOW HERE FOR NOW!!!
# include("ad.jl")
# include("models/models.jl")

end


text:
struct Chunks
    ctx::Context
    Chunks(ctx::Context) = new(ctx)
end

"""
    CSV.Chunks(source; ntasks::Integer=Threads.nthreads(), kwargs...) => CSV.Chunks

Returns a file "chunk" iterator. Accepts all the same inputs and keyword arguments as [`CSV.File`](@ref),
see those docs for explanations of each keyword argument.

The `ntasks` keyword argument specifies how many chunks a file should be split up into, defaulting to
the # of threads available to Julia (i.e. `JULIA_NUM_THREADS` environment variable) or 8 if Julia is
run single-threaded.

Each iteration of `CSV.Chunks` produces the next chunk of a file as a `CSV.File`. While initial file
metadata detection is done only once (to determine # of columns, column names, etc), each iteration
does independent type inference on columns. This is significant as different chunks may end up with
different column types than previous chunks as new values are encountered in the file. Note that, as
with `CSV.File`, types may be passed manually via the `type` or `types` keyword arguments.

This functionality is new and thus considered experimental; please
[open an issue](https://github.com/JuliaData/CSV.jl/issues/new) if you run into any problems/bugs.

$KEYWORD_DOCS
"""
function Chunks(source;
    # file options
    # header can be a row number, range of rows, or actual string vector
    header::Union{Integer, Vector{Symbol}, Vector{String}, AbstractVector{<:Integer}}=1,
    normalizenames::Bool=false,
    # by default, data starts immediately after header or start of file
    datarow::Integer=-1,
    skipto::Integer=-1,
    footerskip::Integer=0,
    transpose::Bool=false,
    comment::Union{String, Nothing}=nothing,
    ignoreemptyrows::Bool=true,
    ignoreemptylines=nothing,
    select=nothing,
    drop=nothing,
    limit::Union{Integer, Nothing}=nothing,
    ntasks::Union{Integer, Nothing}=nothing,
    tasks::Union{Nothing, Integer}=nothing,
    rows_to_check::Integer=DEFAULT_ROWS_TO_CHECK,
    lines_to_check=nothing,
    # parsing options
    missingstrings=String[],
    missingstring="",
    delim::Union{Nothing, Char, String}=nothing,
    ignorerepeated::Bool=false,
    quoted::Bool=true,
    quotechar::Union{UInt8, Char}='"',
    openquotechar::Union{UInt8, Char, Nothing}=nothing,
    closequotechar::Union{UInt8, Char, Nothing}=nothing,
    escapechar::Union{UInt8, Char}='"',
    dateformat::Union{String, Dates.DateFormat, Nothing, AbstractDict}=nothing,
    dateformats=nothing,
    decimal::Union{UInt8, Char}=UInt8('.'),
    truestrings::Union{Vector{String}, Nothing}=TRUE_STRINGS,
    falsestrings::Union{Vector{String}, Nothing}=FALSE_STRINGS,
    # type options
    type=nothing,
    types=nothing,
    typemap::Dict=Dict{Type, Type}(),
    pool::Union{Bool, Real, AbstractVector, AbstractDict}=DEFAULT_POOL,
    downcast::Bool=false,
    lazystrings::Bool=false,
    stringtype::StringTypes=DEFAULT_STRINGTYPE,
    strict::Bool=false,
    silencewarnings::Bool=false,
    maxwarnings::Int=DEFAULT_MAX_WARNINGS,
    debug::Bool=false,
    parsingdebug::Bool=false)

    ctx = Context(source, header, normalizenames, datarow, skipto, footerskip, transpose, comment, ignoreemptyrows, ignoreemptylines, select, drop, limit, nothing, ntasks, tasks, rows_to_check, lines_to_check, missingstrings, missingstring, delim, ignorerepeated, quoted, quotechar, openquotechar, closequotechar, escapechar, dateformat, dateformats, decimal, truestrings, falsestrings, type, types, typemap, pool, downcast, lazystrings, stringtype, strict, silencewarnings, maxwarnings, debug, parsingdebug, false)
    !ctx.threaded && throw(ArgumentError("unable to iterate chunks from input file source"))
    foreach(col -> col.lock = ReentrantLock(), ctx.columns)
    return Chunks(ctx)
end

Base.IteratorSize(::Type{Chunks}) = Base.HasLength()
Base.length(x::Chunks) = x.ctx.ntasks
Base.eltype(x::Chunks) = File{false}

function Base.iterate(x::Chunks, i=1)
    i > x.ctx.ntasks && return nothing
    names = copy(x.ctx.names)
    columns = [Column(col) for col in x.ctx.columns]
    datapos = x.ctx.chunkpositions[i]
    len = x.ctx.chunkpositions[i + 1] - 1
    rowsguess = cld(x.ctx.rowsguess, x.ctx.ntasks)
    threaded = false
    ntasks = 1
    limit = typemax(Int64)
    ctx = Context(x.ctx.transpose, x.ctx.name, names, rowsguess, x.ctx.cols, x.ctx.buf, datapos, len, 1, x.ctx.options, columns, x.ctx.pool, x.ctx.downcast, x.ctx.customtypes, x.ctx.typemap, x.ctx.stringtype, limit, threaded, ntasks, x.ctx.chunkpositions, x.ctx.strict, x.ctx.silencewarnings, x.ctx.maxwarnings, x.ctx.debug, x.ctx.streaming)
    f = File(ctx, true)
    return f, i + 1
end

Tables.partitions(x::Chunks) = x


text:
using OmegaCore
using Test
using Distributions
using OmegaTest

function samplemodel()
  x = 1 ~ Normal(0, 1)
  y = 2 ~ Normal(0, 1)
  z(ω) = (x(ω), y(ω), x(ω))
  randsample(z)
end

function test_ciid()
  @test isinferred(samplemodel)
  z_ = samplemodel()
  @test z_[1] != z_[2]
  @test z_[1] == z_[3]
end

function test_shared_parent()
  parent(ω) = Bool(Bernoulli(0.5)((1,), ω)) ? -100 : 100
  x(id, ω) = parent(ω) + Uniform(0, 1)(id, ω)
  x1 = 2 ~ x
  x2 = 3 ~ x
  z = ω -> (x1(ω), x2(ω))
  samples = randsample(z, 10)
  t1 = [(x1_ != x2_) for (x1_, x2_) in samples]
  t2 = [(abs(x1_ - x2_) <= 2) for (x1_, x2_) in samples]
  @test all(t1)
  @test all(t2)
end

@testset "ciid" begin
  test_ciid()
  test_shared_parent()
end


text:
# Low-level, isolated, functions tests.
include("unit.jl")

# High-level API tests.
# Ensure that everything works well together, for different kinds of HMMs.
include("integration.jl")

max_stars_repo_name: Shu-Lea-Lai/module-MultiHMMBase


text:
cd(@__DIR__)

using Pkg

using Test, TestSetExtensions, SafeTestsets
using Genie

isdir(joinpath(Genie.config.path_build, Genie.Renderer.BUILD_NAME)) && rm(joinpath(Genie.config.path_build, Genie.Renderer.BUILD_NAME), force = true, recursive = true)

@testset ExtendedTestSet "Genie tests" begin
  @includetests ARGS
end


text:
using Unwrap
using Test

@time @testset "Unwrap 1D" begin include("test_unwrap_1D.jl") end
@time @testset "Unwrap 2D" begin include("test_unwrap_2D.jl") end
@time @testset "Unwrap 3D" begin include("test_unwrap_3D.jl") end

max_forks_repo_name: platawiec/Unwrap


text:
module SuperModule
import AbstractModule.foo
include("SubModule1.jl")
include("SubModule2.jl")
export foo
end


text:
const SLICELENGTH = 100000

const MISSINGLEVEL = "."

const HEADLENGTH = 20

const JUMLDIR = "JuML"

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
function issquarefree{I<:Integer}(d::I)
  # factor(d) is a dictionary whose keys are primes
  # and whose values are corresponding prime powers.
  # An integer is square free iff none of its prime
  # factors have power 2 or greater. This is equivalent
  # to its Mobius function being nonzero.
  mobius(d) != 0
end

"""

[Link](https://proofwiki.org/wiki/Definition:M%C3%B6bius_Function)
"""
function mobius{I<:Integer}(n::I)
  if n==1 return I(1) end
  powers = collect(values(factor(n)))
  return all(powers.== 1) ? I(-1)^length(powers) : I(0)
end

#  Division by monic polynomials

function onestep!{I<:Integer}(num::Array{I,1}, den::Array{I,1}, quo::Array{I,1})
  n = findlast(num .!= 0)
  d = findlast(den .!= 0)
  den[d] == 1 || error("denominator is not monic")
  k = n-d
  k >= 0 || error("numerator not divisible by denominator")
  quo[1+k] += num[n]
  for i=1:d
    num[i+k] -= num[n]*den[i]
  end
end

function divide_by_monic{I<:Integer}(num::Array{I,1}, den::Array{I,1})
  cnum = deepcopy(num)
  quo = zeros(I,length(num))
  while any(cnum .!= 0)
    onestep!(cnum, den, quo)
  end
  quo
end

function cyclotomic{I<:Integer}(n::I)
  n >= 1 || error("There is no cyclotomic of degree < 1")
  if n==1
    return [I(-1), I(1)]
  end
  num = zeros(I,n+1)
  num[n+1]=I(1)
  num[1] = I(-1)
  for d in 1:(n-1)
    if n%d == 0
      num = divide_by_monic(num, cyclotomic(d))
    end
  end
  d = findlast(num .!= 0)
  num[1:d]
end

max_stars_repo_name: WilCrofter/ZAdjoin

max_issues_repo_name: WilCrofter/ZAdjoin

max_forks_repo_name: WilCrofter/ZAdjoin


text:
mutable struct BayesianRegressor
    _alpha::Float64
    _beta::Float64
    # prior/posterior distribution N(w | w_mean, w_precision^(-1))
    _w_mean::AbstractArray{Float64,1}
    _w_precision::AbstractArray{Float64,2}
    function BayesianRegressor(alpha::Float64, beta::Float64, n_features::Int64)
        w_mean = zeros(n_features)
        w_precision = alpha * Matrix(I, n_features, n_features)
        new(alpha, beta, w_mean, w_precision)
    end
end

function fitting(
    regressor::BayesianRegressor,
    Phi::AbstractArray{Float64,2},
    t::AbstractArray{Float64,1},
)
    # Phi is the design matrix [phi(x_1), phi(x_2), ,,, phi(x_N)]^T of size (n_samples, n_features)
    n_features = size(Phi)[2]
    n_samples = size(Phi)[1]
    @assert n_samples == size(t)[1]

    w_0 = regressor._w_mean
    S_0_inv = regressor._w_precision
    beta = regressor._beta

    tmp = S_0_inv + beta * transpose(Phi) * Phi
    # make S_N_inv Cholesky
    S_N_inv = (tmp + transpose(tmp)) / 2.0
    w_N = inv(S_N_inv) * (S_0_inv * w_0 + beta * transpose(Phi) * t)

    # update
    regressor._w_mean = w_N
    regressor._w_precision = S_N_inv
end

function predict(
    regressor::BayesianRegressor,
    Phi::AbstractArray{Float64,2},
    return_std::Bool = false,
)
    # Phi is the design matrix [phi(x_1), phi(x_2), ,,, phi(x_N)]^T of size (n_samples, n_features)
    N = size(Phi)[1]
    w_mean = regressor._w_mean
    w_cov = inv(regressor._w_precision)

    y = Phi * w_mean

    if return_std == true
        y_var =
            1.0 / regressor._beta .+ [transpose(Phi[i, :]) * w_cov * Phi[i, :] for i = 1:N]
        y_std = sqrt.(y_var)
        return y, y_std
    end

    return y
end

function predictSampling(
    regressor::BayesianRegressor,
    Phi::AbstractArray{Float64,2},
    n_sampling::Int64,
)
    # Phi is the design matrix [phi(x_1), phi(x_2), ,,, phi(x_N)]^T of size (n_samples, n_features)
    # returns [y(x_1), y(x_2),,,, y(x_N)], each y(x_i) is a AbstractArray of size n_sampling
    N = size(Phi)[1]
    ret = zeros(n_sampling, N)

    w_mean = regressor._w_mean
    w_cov = inv(regressor._w_precision)

    w_samples = rand(MvNormal(w_mean, w_cov), n_sampling)
    for i = 1:n_sampling
        ret[i, :] = transpose(w_samples[:, i]) * transpose(Phi)
    end

    return ret
end

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
using IndirectTrajOpt
using StaticArrays
using LinearAlgebra
using BenchmarkTools
using DifferentialEquations
using ModelingToolkit
using Plots

# Initialize parameters
ps = initCR3BPIndirectParams("Low Thrust 10 CR3BP")
ps.ϵ = 0.0

# Initial condition vector
y0 = @SVector [-0.019488511458668, -0.016033479812051, 0.0,
                8.918881923678198, -4.081793688818725, 0.0, 1.0,
                -24.3123620400503, -60.46766521717989, -0.2883477429418446, 0.18005551834799977, -0.09283727657419928, -0.00021162595932502253, 0.5613066482340945]

# Integrate
tspan = (0.0, 8.6404*24*3600/ps.crp.TU)
#@btime cr3bpOptIntegrate($y0, $tspan, $ps)
sol = cr3bpOptIntegrate(y0, tspan, ps)

# Check that utype is set appropriately
cSc = ps.sp.isp*9.81*ps.crp.TU / (ps.crp.LU*1000.0)
λv = norm(view(y0,11:13))
S = IndirectTrajOpt.computeS(y0, λv, cSc)
if S > ps.ϵ
    ps.utype = 0
elseif S < -ps.ϵ
    ps.utype = 2
else
    ps.utype = 1
end

# Callbacks
cb = ContinuousCallback(
    IndirectTrajOpt.cr3bpEomsConditionNoTerm,
    IndirectTrajOpt.cr3bpEomsAffectNoTerm!,
    IndirectTrajOpt.cr3bpEomsAffectNoTerm!;
    idxs = nothing,
    rootfind = true,
    interp_points = 10,
    abstol = 1e-14,
    reltol = 0.0,
    save_positions = (true, true))

# ODE Problem
ff = ODEFunction{true}(IndirectTrajOpt.cr3bpEomIndirect!)
prob = ODEProblem(ff, Vector{Float64}(y0), tspan, ps; callback=cb)

# Solve ode
sol = solve(
    prob,
    Vern9(),
    reltol = 1e-14,
    abstol = 1e-14,
    save_everystep = true,
    save_start = true,
    initialize_save = true,
    maxiters = 1e6
    )

display(plot(sol, vars=(1,2)))

#@btime IndirectTrajOpt.cr3bpEomIndirect($y0, $ps, 0.0)
#@btime IndirectTrajOpt.cr3bpCostateEom($y0, $ps, 1.0)
#at = @SVector [1.0, 0.0, 0.0]
#@btime IndirectTrajOpt.cr3bpEomControl($y0, $ps.crp,0.0,$at)
#dy0 = similar(y0)
#@btime IndirectTrajOpt.cr3bpEomIndirect!($dy0,$y0,$ps,0.0)


text:
@testset "AbstractStepRangeLen" begin
    for frange in (mrange, srange, range)
        @testset "$frange" begin
            r = frange(1.0, step=1, stop=10.0)
            b = range(1.0, step=1, stop=10.0)
            @test stephi(r) == 1
            @test steplo(r) == 0
            @test refhi(r) == 1
            @test reflo(r) == 0
            @test eltype(StepRangeLen{Float64}(r)) == Float64

            @test intersect(r, r[2]) == intersect(b, b[2])
            @test intersect(r, r[2]) == intersect(b, b[2])

            if frange == mrange
                setproperty!(r, :step, 2)
                @test r == range(1.0, step=2, length=10)

                setproperty!(r, :ref, 2)
                @test r == range(2.0, step=2, length=10)

                setproperty!(r, :len, 5)
                @test r == range(2.0, step=2, length=5)

                setproperty!(r, :offset, 2)
                @test r == StepRangeLen(2.0, 2.0, 5, 2)

                @test_throws ErrorException setproperty!(r, :anything, 3)
            end
        end
    end
end


text:
using Documenter, EnglishText

makedocs(
    format = :html,
    sitename = "EnglishText.jl",
    authors = "Fengyang Wang",
    analytics = "UA-68884109-1",
    pages = [
        "index.md"
    ]
)

deploydocs(
    julia  = "nightly",
    repo   = "github.com/TotalVerb/EnglishText.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing,
)


text:
immutable DisplacementMesh <: ImageTransformation
    input_vertices::Matrix{Float64}
    output_vertices::Matrix{Float64}
    indices::Matrix{Int}

    function DisplacementMesh(input_vertices::Matrix{Float64}, output_vertices::Matrix{Float64}, indices::Matrix{Int})
        @assert size(input_vertices) == size(output_vertices)
        new(input_vertices, output_vertices, indices)
    end
end

function DisplacementMesh(field::DisplacementField, img_width::Int, img_height::Int)
    input_vertices, output_vertices = _compute_vertices(field, Float64(img_width), Float64(img_height))
    indices = _compute_indices(field)
    DisplacementMesh(input_vertices, output_vertices, indices)
end

DisplacementMesh(field::DisplacementField, img_size::Tuple{Int,Int}) = DisplacementMesh(field, img_size[1], img_size[2])
DisplacementMesh(field::DisplacementField, img::AbstractImage) = DisplacementMesh(field, size(img, "x"), size(img, "y"))

function Base.show(io::IO, dm::DisplacementMesh)
    print(io, "DisplacementMesh (vertices: $(size(dm.input_vertices,1)), triangles: $(size(dm.indices,1)))")
end

@recipe function plot(dm::DisplacementMesh, img::Image)
    legend --> false
    yflip := true
    layout := 2

    @series begin
        title --> "Input"
        seriestype := :image
        subplot := 1
        img
    end

    @series begin
        title --> "Output"
        seriestype := :image
        subplot := 2
        transform(dm,img)
    end

    dm
end

@recipe function plot(dm::DisplacementMesh)
    fillcolor --> :transparent
    legend --> false
    yflip := true
    layout := 2

    n_indices = size(dm.indices,1)
    x_in = zeros(n_indices*4)
    y_in = zeros(n_indices*4)
    x_out = zeros(n_indices*4)
    y_out = zeros(n_indices*4)
    ti = 1
    for i = 1:n_indices
        for j = 1:3
            x_in[ti] = dm.input_vertices[dm.indices[i, j], 2]
            y_in[ti] = dm.input_vertices[dm.indices[i, j], 1]
            x_out[ti] = dm.output_vertices[dm.indices[i, j], 2]
            y_out[ti] = dm.output_vertices[dm.indices[i, j], 1]
            ti += 1
        end
        x_in[ti] = NaN
        y_in[ti] = NaN
        x_out[ti] = NaN
        y_out[ti] = NaN
        ti += 1
    end

    @series begin
        title --> "Input"
        seriestype := :shape
        subplot := 1
        x_in, y_in
    end
    @series begin
        title --> "Output"
        seriestype := :shape
        subplot := 2
        x_out, y_out
    end
end

function _transform{T}(dm::DisplacementMesh, img::Image{T})
    img_sep = separate(img)
    rawdata = convert(Array{Float64, 3}, data(img_sep))
    wrp = PiecewiseAffineTransforms.pa_warp(rawdata, dm.input_vertices, dm.output_vertices, dm.indices)
    rawdata_new = convert(Array{eltype(T),3}, wrp)
    img_sep_new = Image(rawdata_new, properties(img_sep))
    if isxfirst(img)
        copyproperties(img, permutedims(convert(Image{T}, img_sep_new), [2, 1]))
    else
        copyproperties(img, convert(Image{T}, img_sep_new))
    end
end

function transform{T<:AbstractImage}(dm::DisplacementMesh, img::T)
    result = _transform(dm, img)::T
    _log_operation!(result, dm)::T
end

function _compute_vertices(field::DisplacementField, img_width::Float64, img_height::Float64)
    height, width = size(field.delta_X)
    input_vertices  = zeros(height*width, 2)
    output_vertices = zeros(height*width, 2)
    i = 1
    for x = 1:width, y = 1:height
        input_vertices[i,1] = clamp(clamp(field.y[y], 0., 1.) * img_height, 1., img_height)
        input_vertices[i,2] = clamp(clamp(field.x[x], 0., 1.) * img_width,  1., img_width)
        output_vertices[i,1] = clamp(input_vertices[i,1] + field.delta_Y[y,x] * img_height, 1., img_height)
        output_vertices[i,2] = clamp(input_vertices[i,2] + field.delta_X[y,x] * img_width,  1., img_width)
        i = i + 1
    end
    input_vertices, output_vertices
end

function _compute_indices(field::DisplacementField)
    grid_size = size(field.delta_X)
    height, width = grid_size
    w_half = floor(Int, width/2)
    h_half = floor(Int, height/2)
    indices = zeros(Int, (height-1)*(width-1)*2, 3)
    i = 1
    for x = 1:(width-1), y = 1:(height-1)
        if (x <= w_half && y <= h_half) || (x > w_half && y > h_half)
            # upper left or lower right
            # *--*
            # |\ |
            # | \|
            # *--*
            indices[i, 1] = sub2ind(grid_size, y,   x  )
            indices[i, 2] = sub2ind(grid_size, y+1, x+1)
            indices[i, 3] = sub2ind(grid_size, y+1, x  )
            i = i + 1

            indices[i, 1] = sub2ind(grid_size, y,   x  )
            indices[i, 2] = sub2ind(grid_size, y,   x+1)
            indices[i, 3] = sub2ind(grid_size, y+1, x+1)
            i = i + 1
        else
            # lower left or upper right
            # *--*
            # | /|
            # |/ |
            # *--*
            indices[i, 1] = sub2ind(grid_size, y,   x  )
            indices[i, 2] = sub2ind(grid_size, y,   x+1)
            indices[i, 3] = sub2ind(grid_size, y+1, x  )
            i = i + 1

            indices[i, 1] = sub2ind(grid_size, y+1,   x)
            indices[i, 2] = sub2ind(grid_size, y,   x+1)
            indices[i, 3] = sub2ind(grid_size, y+1, x+1)
            i = i + 1
        end
    end
    indices
end


text:
module Quadratic

export quadratic, +, -, *, /, ≤

using IntervalSets
using Base: +, ≤

struct quadratic{S,T,U}
    a::S
    b::T
    c::U
end

(f::quadratic)(x) = f.a * x^2 + f.b * x + f.c

Base.:+(f::quadratic, g::quadratic) = quadratic(f.a + g.a, f.b + g.b, f.c + g.c)
Base.:-(f::quadratic, g::quadratic) = quadratic(f.a - g.a, f.b - g.b, f.c - g.c)
Base.:*(f::quadratic, c) = quadratic(f.a * c, f.b * c, f.c * c)
Base.:/(f::quadratic, c) = quadratic(f.a / c, f.b / c, f.c / c)

# return: region that f is less than g
function Base.:≤(f::quadratic, g::quadratic)
    # coefs of h ≔ f - g
    a = f.a - g.a
    b = f.b - g.b
    c = f.c - g.c

    # if h is linear or constant
    if iszero(a)
        iszero(b) && return c ≤ 0 ? [-Inf..Inf] : ClosedInterval{Float64}[]
        return b > 0 ? [-Inf..(-c / b)] : [(-c / b)..Inf]
    end

    # if h is quadratic
    D = b^2 - 4 * a * c
    if D > 0
        x₁, x₂ = (-b + √D) / (2a), (-b - √D) / (2a)
        if a > 0
            @assert x₁ > x₂
            return [x₂..x₁]
        else
            @assert x₁ < x₂
            return [-Inf..x₁, x₂..Inf]
        end
    else
        return a < 0 ? [-Inf..Inf] : ClosedInterval{Float64}[]
    end
end

end # module Quadratic

# # test
# using Plots
# f = quadratic(2, -2, 1)
# g = quadratic(-4, 5, 100)
# xs = -10:0.1:10
# plot(xs, f.(xs), label="f", c="blue", legend=:inline)
# plot!(xs, g.(xs), label="g", c="green", legend=:inline)
# region = f ≤ g
# for interval in region
#     xs = range(max(interval.left, -10), min(interval.right, 10); length=100)
#     plot!(xs, f.(xs), label=nothing, color="red", lw=3) |> display
# end


text:
"""
    struct Parser <: Any
Hold a pointer to a `clang::Parser` object.
"""
struct Parser
    ptr::CXParser
end

Base.unsafe_convert(::Type{CXParser}, x::Parser) = x.ptr
Base.cconvert(::Type{CXParser}, x::Parser) = x


text:
using Iterators
using RLESUtils, RunUtils, IFTTTUtils

const COMP = gethostname()

V_seeds = collect(1:20)
V_efs = [0.1, 0.25, 0.5, 0.6]
V_wnews = [0.1, 0.2, 0.4, 0.6, 0.8, 1.0]
V_wprior = [0.01, 0.05, 0.1, 0.2, 0.4]
V_maxsteps = [20, 40, 60]

template_ce(seed, elite_frac, w_new, w_prior, maxsteps) =
"""
using ExprSearch,Symbolic_CE,RLESUtils,IFTTTUtils
config=Dict{Symbol,Any}()
config[:seed]=$seed
config[:elite_frac]=$(elite_frac)
config[:w_new]=$(w_new)
config[:w_prior]=$(w_prior)
config[:maxsteps]=$maxsteps
symbolic_ce(; outdir=Pkg.dir("ExprSearch/results/Symbolic_CE_ParamStudy/seed$(seed)ef$(elite_frac)wnew$(w_new)wprior$(w_prior)maxsteps$maxsteps"), config...)
"""

notifydone() = sendifttt(;value1="$COMP done")

A = JuliaSource[]
for args in product(V_seeds, V_efs, V_wnews, V_wprior, V_maxsteps)
    push!(A, JuliaSource(template_ce(args...)))
end

#include("script.jl")
#fs=julia_process(A, np); notifydone()

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
# ---
# title: 1052. Grumpy Bookstore Owner
# id: problem1052
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Sliding Window
# link: <https://leetcode.com/problems/grumpy-bookstore-owner/description/>
# hidden: true
# ---
#
# Today, the bookstore owner has a store open for `customers.length` minutes.
# Every minute, some number of customers (`customers[i]`) enter the store, and
# all those customers leave after the end of that minute.
#
# On some minutes, the bookstore owner is grumpy.  If the bookstore owner is
# grumpy on the i-th minute, `grumpy[i] = 1`, otherwise `grumpy[i] = 0`.  When
# the bookstore owner is grumpy, the customers of that minute are not satisfied,
# otherwise they are satisfied.
#
# The bookstore owner knows a secret technique to keep themselves not grumpy for
# `X` minutes straight, but can only use it once.
#
# Return the maximum number of customers that can be satisfied throughout the
# day.
#
#
#
# **Example 1:**
#
#
#
#     Input: customers = [1,0,1,2,1,1,7,5], grumpy = [0,1,0,1,0,1,0,1], X = 3
#     Output: 16
#     Explanation:  The bookstore owner keeps themselves not grumpy for the last 3 minutes.
#     The maximum number of customers that can be satisfied = 1 + 1 + 1 + 1 + 7 + 5 = 16.
#
#
#
#
# **Note:**
#
#   * `1 <= X <= customers.length == grumpy.length <= 20000`
#   * `0 <= customers[i] <= 1000`
#   * `0 <= grumpy[i] <= 1`
#
#
## @lc code=start
using LeetCode

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


text:
"""
    function form_teams(
        method::String
        players::Array{Player}
    )::Array{Team}

Creates an array of teams from the given input of players.

Arguments
=========
- `method`:
    the strategy for forming teams.
    current strategies include:
        - "Next Best Player"
            assigns each player to be partners with the next best available player.
- `players`
    an array of players from which to form teams.
- `ranking_scheme`
    ranking scheme used to evaluate the players
"""

function form_teams(
    method::String,
    players::Array{Player},
    ranking_scheme::String,
)::Array{Team}
    if method == "Next Best Player"
        teams = []
        for i in 1:2:length(players)
            t = Team(players[i], players[i+1], ranking_scheme)
            push!(teams, t)
        end
        return teams
    else
        throw(
            DomainError(
                method,
                "please select from the following options:\n
                 1. Next Best Player",
            ),
        )
    end
end

max_stars_repo_name: Danial-Hussain/Tennis-Ranking

max_issues_repo_name: Danial-Hussain/Tennis-Ranking

max_forks_repo_name: Danial-Hussain/Tennis-Ranking


text:
# Thu 08 Nov 2012 05:07:44 PM EST
#
# NeHe Tut 6 - Rotate a textured cube
#
# Q - quit

# load necessary GLUT/OpenGL routines

global OpenGLver="1.0"
using OpenGL
using GLUT

### auxiliary functions

function cube(size)
  glBegin(GL_QUADS)
    # Front Face
    glTexCoord(0.0, 0.0)
    glVertex(-size, -size, size)
    glTexCoord(1.0, 0.0)
    glVertex(size, -size, size)
    glTexCoord(1.0, 1.0)
    glVertex(size, size, size)
    glTexCoord(0.0, 1.0)
    glVertex(-size, size, size)

    # Back Face
    glTexCoord(1.0, 0.0)
    glVertex(-size, -size, -size)
    glTexCoord(1.0, 1.0)
    glVertex(-size, size, -size)
    glTexCoord(0.0, 1.0)
    glVertex(size, size, -size)
    glTexCoord(0.0, 0.0)
    glVertex(size, -size, -size)

    # Top Face
    glTexCoord(0.0, 1.0)
    glVertex(-size, size, -size)
    glTexCoord(0.0, 0.0)
    glVertex(-size, size, size)
    glTexCoord(1.0, 0.0)
    glVertex(size, size, size)
    glTexCoord(1.0, 1.0)
    glVertex(size, size, -size)

    # Bottom Face
    glTexCoord(1.0, 1.0)
    glVertex(-size, -size, -size)
    glTexCoord(0.0, 1.0)
    glVertex(size, -size, -size)
    glTexCoord(0.0, 0.0)
    glVertex(size, -size, size)
    glTexCoord(1.0, 0.0)
    glVertex(-size, -size, size)

    # Right Face
    glTexCoord(1.0, 0.0)
    glVertex(size, -size, -size)
    glTexCoord(1.0, 1.0)
    glVertex(size, size, -size)
    glTexCoord(0.0, 1.0)
    glVertex(size, size, size)
    glTexCoord(0.0, 0.0)
    glVertex(size, -size, size)

    # Left Face
    glTexCoord(0.0, 0.0)
    glVertex(-size, -size, -size)
    glTexCoord(1.0, 0.0)
    glVertex(-size, -size, size)
    glTexCoord(1.0, 1.0)
    glVertex(-size, size, size)
    glTexCoord(0.0, 1.0)
    glVertex(-size, size, -size)
  glEnd()
end

### end of auxiliary functions

# initialize variables

global window

global xrot      = 0.0
global yrot      = 0.0
global zrot      = 0.0

global tex       = Array(Uint32,1) # generating 1 texture

global cube_size = 1.0

width            = 640
height           = 480

# load textures from images

function LoadGLTextures()
    global tex

    img, w, h = glimread(expanduser("~/.julia/GLUT/Examples/NeHe/tut6/NeHe.bmp"))

    glGenTextures(1,tex)
    glBindTexture(GL_TEXTURE_2D,tex[1])
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
    glTexImage2D(GL_TEXTURE_2D, 0, 3, w, h, 0, GL_RGB, GL_UNSIGNED_BYTE, img)
end

# function to init OpenGL context

function initGL(w::Integer,h::Integer)
    glViewport(0,0,w,h)
    LoadGLTextures()
    glClearColor(0.0, 0.0, 0.0, 0.0)
    glClearDepth(1.0)
    glDepthFunc(GL_LESS)
    glEnable(GL_DEPTH_TEST)
    glShadeModel(GL_SMOOTH)

    # enable texture mapping
    glEnable(GL_TEXTURE_2D)

    glMatrixMode(GL_PROJECTION)
    glLoadIdentity()

    gluPerspective(45.0,w/h,0.1,100.0)

    glMatrixMode(GL_MODELVIEW)
end

# prepare Julia equivalents of C callbacks that are typically used in GLUT code

function ReSizeGLScene(w::Int32,h::Int32)
    if h == 0
        h = 1
    end

    glViewport(0,0,w,h)

    glMatrixMode(GL_PROJECTION)
    glLoadIdentity()

    gluPerspective(45.0,w/h,0.1,100.0)

    glMatrixMode(GL_MODELVIEW)

    return nothing
end

_ReSizeGLScene = cfunction(ReSizeGLScene, Void, (Int32, Int32))

function DrawGLScene()
    global tex
    global xrot
    global yrot
    global zrot

    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glLoadIdentity()

    glTranslate(0.0, 0.0, -5.0)

    glRotate(xrot,1.0,0.0,0.0)
    glRotate(yrot,0.0,1.0,0.0)
    glRotate(zrot,0.0,0.0,1.0)

    glBindTexture(GL_TEXTURE_2D,tex[1])
    cube(cube_size)

    xrot +=0.2
    yrot +=0.3
    zrot +=0.4

    glutSwapBuffers()

    return nothing
end

_DrawGLScene = cfunction(DrawGLScene, Void, ())

function keyPressed(the_key::Char,x::Int32,y::Int32)
    if the_key == int('q')
        glutDestroyWindow(window)
    end

    return nothing # keyPressed returns "void" in C. this is a workaround for Julia's "automatically return the value of the last expression in a function" behavior.
end

_keyPressed = cfunction(keyPressed, Void, (Char, Int32, Int32))

# run GLUT routines

glutInit()
glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_ALPHA | GLUT_DEPTH)
glutInitWindowSize(width, height)
glutInitWindowPosition(0, 0)

window = glutCreateWindow("NeHe Tut 6")

glutDisplayFunc(_DrawGLScene)
glutFullScreen()

glutIdleFunc(_DrawGLScene)
glutReshapeFunc(_ReSizeGLScene)
glutKeyboardFunc(_keyPressed)

initGL(width, height)

glutMainLoop()