text:
# Recreates the analysis from the original data files

using SMLMAssociationAnalysis_NCB
using FileIO
using StatsBase
using Statistics
using StatsPlots
using SimpleANOVA

### Load saved data

samplenames = ["A", "B", "C", "D"]

nreplicates = 3
nsamples = 4
ncells = 10

outputdir = "output"
datapath = joinpath(outputdir, "results.jld2")

experimentresults = load(datapath)["experimentresults"]

medianmeasurements = Array{Float64,4}(undef, ncells, nreplicates, 4, 2)
montecarlomeasurements = Array{Float64,4}(undef, ncells, nreplicates, 4, 2)
positivecontrolmontecarlomeasurements = Array{Float64,4}(undef, ncells, nreplicates, 4, 2)
negativecontrolmontecarlomeasurements = Array{Float64,4}(undef, ncells, nreplicates, 4, 2)

for k ∈ 1:2
    for i ∈ 1:nsamples
        samplemedianresults = Array{Float64,2}(undef, ncells, nreplicates)
        samplelessthan10results = Array{Float64,2}(undef, ncells, nreplicates)
        samplepositivecontrollessthan10results = Array{Float64,2}(undef, ncells, nreplicates)
        samplenegativecontrollessthan10results = Array{Float64,2}(undef, ncells, nreplicates)

        for j ∈ 1:nreplicates
            replicateresults = experimentresults[k][j][i]
            lessthanlimitreplicate = [(x.distances .< 200) .& (x.percentileranks .< 0.1) for x ∈ replicateresults]
            positivecontrollessthanlimitreplicate = [(x.positivecontrol_distances .< 200) .& (x.positivecontrol_percentileranks .< 0.1) for x ∈ replicateresults]
            negativecontrollessthanlimitreplicate = median.([count.([(x.negativecontrol_distances[m] .< 200) .& (x.negativecontrol_percentileranks[m] .< 0.1) for m ∈ eachindex(x.negativecontrol_distances)]) for x ∈ replicateresults])
            lessthan10 = count.(lessthanlimitreplicate) ./ length.(lessthanlimitreplicate)
            positivecontrollessthan10 = count.(positivecontrollessthanlimitreplicate) ./ length.(lessthanlimitreplicate)
            negativecontrollessthan10 = negativecontrollessthanlimitreplicate ./ length.(lessthanlimitreplicate)
            mediandistances = map(x -> x.mediandistance, replicateresults)
            samplemedianresults[:, j] = mediandistances
            samplelessthan10results[:, j] = lessthan10
            samplepositivecontrollessthan10results[:, j] = positivecontrollessthan10
            samplenegativecontrollessthan10results[:, j] = negativecontrollessthan10
        end
        medianmeasurements[:, :, i, k] = samplemedianresults
        montecarlomeasurements[:, :, i, k] = samplelessthan10results
        positivecontrolmontecarlomeasurements[:, :, i, k] = samplepositivecontrollessthan10results
        negativecontrolmontecarlomeasurements[:, :, i, k] = samplenegativecontrollessthan10results
    end
end

medianmeasurements = cat(
    cat(medianmeasurements[:, :, 1:2, 1], medianmeasurements[:, :, 3:4, 1], dims = 4),
    cat(medianmeasurements[:, :, 1:2, 2], medianmeasurements[:, :, 3:4, 2], dims = 4),
    dims = 5,
)
montecarlomeasurements = cat(
    cat(montecarlomeasurements[:, :, 1:2, 1], montecarlomeasurements[:, :, 3:4, 1], dims = 4),
    cat(montecarlomeasurements[:, :, 1:2, 2], montecarlomeasurements[:, :, 3:4, 2], dims = 4),
    dims = 5,
)
positivecontrolmontecarlomeasurements = cat(
    cat(positivecontrolmontecarlomeasurements[:, :, 1:2, 1], positivecontrolmontecarlomeasurements[:, :, 3:4, 1], dims = 4),
    cat(positivecontrolmontecarlomeasurements[:, :, 1:2, 2], positivecontrolmontecarlomeasurements[:, :, 3:4, 2], dims = 4),
    dims = 5,
)
negativecontrolmontecarlomeasurements = cat(
    cat(negativecontrolmontecarlomeasurements[:, :, 1:2, 1], negativecontrolmontecarlomeasurements[:, :, 3:4, 1], dims = 4),
    cat(negativecontrolmontecarlomeasurements[:, :, 1:2, 2], negativecontrolmontecarlomeasurements[:, :, 3:4, 2], dims = 4),
    dims = 5,
)

normalizedmontecarlomeasurements = (montecarlomeasurements .- negativecontrolmontecarlomeasurements) ./ (positivecontrolmontecarlomeasurements .- negativecontrolmontecarlomeasurements)

#### Analysis
show("Automatically executing this section isn't ideal, as results may not all print, although graphics will be saved. Intended to go through one-by-one as part of the analysis.")

### Median Exp 2

# Check for unusual cases
p1 = boxplot(medianmeasurements[:, :, 1, 1, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Median distance (nm)"))
p2 = boxplot(medianmeasurements[:, :, 2, 1, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Median distance (nm)"))
p3 = boxplot(medianmeasurements[:, :, 1, 2, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Median distance (nm)"))
p4 = boxplot(medianmeasurements[:, :, 2, 2, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Median distance (nm)"))

plot(p1, p2, p3, p4, layout = grid(2, 2), legend = :none, plot_title = "MDM2-p53")
savefig(joinpath(outputdir, "Mdm2-p53 boxplots.png"))

#ZResid/ZPred plot and Levene's test
mediansflat = [medianmeasurements[:, :, 1, 1, 1] medianmeasurements[:, :, 2, 1, 1] medianmeasurements[:, :, 1, 2, 1] medianmeasurements[:, :, 2, 2, 1]]
z = zscore(mediansflat)
zpred = repeat(mean(z, dims = 1), 10)
zresid = z .- zpred
scatter(
    zresid,
    zpred,
    xaxis = ("Standardized Residual (ZResid)"),
    yaxis = ("Standardized Predicted Value (ZPred)"),
    legend = :none,
)
savefig(joinpath(outputdir, "Mdm2-p53 zresid-zpred.png"))

levene(mediansflat)

# qqnorm, skewness, kurtosis
p = [qqnorm(zresid[:, i]) for i ∈ 1:12]
plot(p..., layout = grid(4, 3), legend = :none, plot_title = "MDM2-p53 qqnorm")
savefig(joinpath(outputdir, "Mdm2-p53 qqnorm.png"))

[skewness(mediansflat[:, i]) for i ∈ 1:12]
[kurtosis(mediansflat[:, i]) for i ∈ 1:12]

# anova

medianresult = anova(
    medianmeasurements[:, :, :, :, 1],
    [nested],
    factornames = ["Replicate", "Doxycycline", "Nutlin-3a"],
)

plot(medianresult)
savefig(joinpath(outputdir, "Mdm2-p53 interactionplot.png"))

### Monte Carlo Exp 2

# Check for unusual cases
p1 = boxplot(montecarlomeasurements[:, :, 1, 1, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))
p2 = boxplot(montecarlomeasurements[:, :, 2, 1, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))
p3 = boxplot(montecarlomeasurements[:, :, 1, 2, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))
p4 = boxplot(montecarlomeasurements[:, :, 2, 2, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))

plot(p1, p2, p3, p4, layout = grid(2, 2), legend = :none, plot_title = "Mdm2-p53")
savefig(joinpath(outputdir, "Mdm2-p53 boxplots montecarlo.png"))

#ZResid/ZPred plot and Levene's test
montecarloflat = [montecarlomeasurements[:, :, 1, 1, 1] montecarlomeasurements[:, :, 2, 1, 1] montecarlomeasurements[:, :, 1, 2, 1] montecarlomeasurements[:, :, 2, 2, 1]]
z = zscore(montecarloflat)
zpred = repeat(mean(z, dims = 1), 10)
zresid = z .- zpred
scatter(
    zresid,
    zpred,
    xaxis = ("Standardized Residual (ZResid)"),
    yaxis = ("Standardized Predicted Value (ZPred)"),
    legend = :none,
)
savefig(joinpath(outputdir, "Mdm2-p53 zresid-zpred montecarlo.png"))

levene(montecarloflat)

# qqnorm, skewness, kurtosis
p = [qqnorm(zresid[:, i]) for i ∈ 1:12]
plot(p..., layout = grid(4, 3), legend = :none, plot_title = "MDM2-p53 qqnorm montecarlo")
savefig(joinpath(outputdir, "Mdm2-p53 qqnorm montecarlo.png"))

[skewness(montecarloflat[:, i]) for i ∈ 1:12]
[kurtosis(montecarloflat[:, i]) for i ∈ 1:12]

# anova

montecarloresult = anova(
    montecarlomeasurements[:, :, :, :, 1],
    [nested],
    factornames = ["Replicate", "Doxycycline", "Nutlin-3a"],
)

plot(montecarloresult)
savefig(joinpath(outputdir, "Mdm2-p53 interactionplot montecarlo.png"))

montecarloresultNutMinus = anova(
    montecarlomeasurements[:, :, :, 1, 2],
    [nested],
    factornames = ["Replicate", "Doxycycline"],
)
montecarloresultNutPlus = anova(
    montecarlomeasurements[:, :, :, 2, 2],
    [nested],
    factornames = ["Replicate", "Doxycycline"],
)

### Paper figure

nutdoxgroups = repeat([2, 1, 4, 3], inner=30)
p53_mdm2_median = [medianmeasurements[:,:,1,1,1] |> vec; medianmeasurements[:,:,2,1,1] |> vec; medianmeasurements[:,:,1,2,1] |> vec; medianmeasurements[:,:,2,2,1] |> vec]
boxplot(nutdoxgroups, p53_mdm2_median, outliers=false,
        label=["- Dox", "+ Dox"],
        guidefontsize=52,
        tickfontsize=48,
        legend=:none,
        top_margin=10mm,
        bottom_margin=-25mm,
        left_margin=-3mm,
        seriescolor=[:white :lightgray],
        line=(6, 1.0),
        gridopacity=0.3,
        xgrid=:none,
        size=(1024,2048),
        xaxis=("Condition", (1:4, ["-Dox -Nut", "+Dox -Nut", "-Dox +Nut", "+Dox +Nut"])),
        yaxis=("Median exclusive pairwise distance (nm)", (0,500)))
dotplot!(nutdoxgroups, p53_mdm2_median, mode = :density, label="", marker=(8, repeat([:orange, :darkblue, :darkred], inner=10), stroke(0)))
savefig( "p53_mdm2_median_boxplot.png")

p53_mdm2_montecarlo = [montecarlomeasurements[:,:,1,1,1] |> vec; montecarlomeasurements[:,:,2,1,1] |> vec; montecarlomeasurements[:,:,1,2,1] |> vec; montecarlomeasurements[:,:,2,2,1] |> vec]
groupedboxplot(nutdoxgroups, p53_mdm2_montecarlo, outliers=false,
        label=["- Dox", "+ Dox"],
        guidefontsize=12,
        tickfontsize=72,
        legend=:none,
        left_margin=20mm,
        top_margin=5mm,
        bottom_margin=5mm,
        seriescolor=[:white :lightgray],
        line=(6, 1.0),
        size=(1024,2048),
        xaxis=("Condition", (1:4, ["-Dox -Nut", "+Dox -Nut", "-Dox +Nut", "+Dox +Nut"])),
        yaxis=("Fraction bound", (0,0.15)))
groupeddotplot!(nutdoxgroups, p53_mdm2_montecarlo, mode = :density, label="", marker=(8, repeat([:orange, :darkblue, :darkred], inner=10), stroke(0)))
savefig(joinpath(outputdir, "p53_mdm2_montecarlo_boxplot.png"))

p53_mdm2_montecarlo = [normalizedmontecarlomeasurements[:,:,1,1,1] |> vec; normalizedmontecarlomeasurements[:,:,2,1,1] |> vec; normalizedmontecarlomeasurements[:,:,1,2,1] |> vec; normalizedmontecarlomeasurements[:,:,2,2,1] |> vec]
groupedboxplot(nutdoxgroups, p53_mdm2_montecarlo, outliers=false,
        label=["- Dox", "+ Dox"],
        guidefontsize=52,
        tickfontsize=48,
        legend=:none,
        top_margin=10mm,
        bottom_margin=-25mm,
        left_margin=5mm,
        seriescolor=[:white :lightgray],
        line=(6, 1.0),
        gridopacity=0.3,
        xgrid=:none,
        size=(1024,2048),
        xaxis=("Condition", (1:4, ["-Dox -Nut", "+Dox -Nut", "-Dox +Nut", "+Dox +Nut"])),
        yaxis=("Fraction bound", (-0.04,0.1), -0.04:0.02:0.1))
groupeddotplot!(nutdoxgroups, p53_mdm2_montecarlo, mode = :density, label="",
        marker=(8, repeat([:orange, :darkblue, :darkred], inner=10), stroke(0)))
savefig("p53_mdm2_normalizedmontecarlo_boxplot.png")

### Median Exp 3

# Check for unusual cases
p1 = boxplot(medianmeasurements[:, :, 1, 1, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Median distance (nm)"))
p2 = boxplot(medianmeasurements[:, :, 2, 1, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Median distance (nm)"))
p3 = boxplot(medianmeasurements[:, :, 1, 2, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Median distance (nm)"))
p4 = boxplot(medianmeasurements[:, :, 2, 2, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Median distance (nm)"))

plot(p1, p2, p3, p4, layout = grid(2, 2), legend = :none, plot_title = "MEG3-p53")
savefig(joinpath(outputdir, "MEG3-p53 boxplots.png"))

#ZResid/ZPred plot and Levene's test
mediansflat = [medianmeasurements[:, :, 1, 1, 2] medianmeasurements[:, :, 2, 1, 2] medianmeasurements[:, :, 1, 2, 2] medianmeasurements[:, :, 2, 2, 2]]
z = zscore(mediansflat)
zpred = repeat(mean(z, dims = 1), 10)
zresid = z .- zpred
scatter(
    zresid,
    zpred,
    xaxis = ("Standardized Residual (ZResid)"),
    yaxis = ("Standardized Predicted Value (ZPred)"),
    legend = :none,
)
savefig(joinpath(outputdir, "MEG3-p53 zresid-zpred.png"))

levene(mediansflat)

# qqnorm, skewness, kurtosis
p = [qqnorm(zresid[:, i]) for i ∈ 1:12]
plot(p..., layout = grid(4, 3), legend = :none, plot_title = "MEG3-p53 qqnorm")
savefig(joinpath(outputdir, "MEG3-p53 qqnorm.png"))

[skewness(mediansflat[:, i]) for i ∈ 1:12]
[kurtosis(mediansflat[:, i]) for i ∈ 1:12]

# anova

medianresult = anova(medianmeasurements[:, :, :, :, 2], [nested], factornames = ["Replicate", "Doxycycline", "RNA"])

plot(medianresult)
savefig(joinpath(outputdir, "MEG3-p53 interactionplot.png"))

### Monte Carlo Exp 3

# Check for unusual cases
p1 = boxplot(montecarlomeasurements[:, :, 1, 1, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))
p2 = boxplot(montecarlomeasurements[:, :, 2, 1, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))
p3 = boxplot(montecarlomeasurements[:, :, 1, 2, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))
p4 = boxplot(montecarlomeasurements[:, :, 2, 2, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))

plot(p1, p2, p3, p4, layout = grid(2, 2), legend = :none, plot_title = "MEG3-p53")
savefig(joinpath(outputdir, "MEG3-p53 boxplots montecarlo.png"))

#ZResid/ZPred plot and Levene's test
montecarloflat = [montecarlomeasurements[:, :, 1, 1, 2] montecarlomeasurements[:, :, 2, 1, 2] montecarlomeasurements[:, :, 1, 2, 2] montecarlomeasurements[:, :, 2, 2, 2]]
z = zscore(montecarloflat)
zpred = repeat(mean(z, dims = 1), 10)
zresid = z .- zpred
scatter(
    zresid,
    zpred,
    xaxis = ("Standardized Residual (ZResid)"),
    yaxis = ("Standardized Predicted Value (ZPred)"),
    legend = :none,
)
savefig(joinpath(outputdir, "MEG3-p53 zresid-zpred montecarlo.png"))

levene(montecarloflat)

# qqnorm, skewness, kurtosis
p = [qqnorm(zresid[:, i]) for i ∈ 1:12]
plot(p..., layout = grid(4, 3), legend = :none, plot_title = "MEG3-p53 qqnorm montecarlo")
savefig(joinpath(outputdir, "MEG3-p53 qqnorm montecarlo.png"))

[skewness(montecarloflat[:, i]) for i ∈ 1:12]
[kurtosis(montecarloflat[:, i]) for i ∈ 1:12]

# anova

montecarloresult = anova(
    montecarlomeasurements[:, :, :, :, 2],
    [nested],
    factornames = ["Replicate", "Doxycycline", "RNA"],
)

plot(montecarloresult)
savefig(joinpath(outputdir, "MEG3-p53 interactionplot montecarlo.png"))

montecarloresultMEG3 = anova(
    montecarlomeasurements[:, :, :, 1, 2],
    [nested],
    factornames = ["Replicate", "Doxycycline"],
)
montecarloresultGAPDH = anova(
    montecarlomeasurements[:, :, :, 2, 2],
    [nested],
    factornames = ["Replicate", "Doxycycline"],
)

### Paper figure

import Plots.mm
rnagroups = repeat([1,2], inner = 60)
doxgroups = repeat([2, 1], inner = 30, outer = 2)
p53_meg3_median = [medianmeasurements[:,:,1,1,2] |> vec; medianmeasurements[:,:,2,1,2] |> vec; medianmeasurements[:,:,1,2,2] |> vec; medianmeasurements[:,:,2,2,2] |> vec]
groupedboxplot(rnagroups, p53_meg3_median, group = doxgroups, outliers=false,
        label=["- Dox", "+ Dox"],
        guidefontsize=52,
        tickfontsize=48,
        legend=:none,
        top_margin=10mm,
        bottom_margin=-25mm,
        seriescolor=[:white :lightgray],
        line=(6, 1.0),
        gridopacity=0.3,
        xgrid=:none,
        xaxis=("RNA", (1:2, ["MEG3", "GAPDH"])),
        size=(1024,2048),
        yaxis=("Median exclusive pairwise distance (nm)", (0,1200)))
groupeddotplot!(rnagroups, p53_meg3_median, group = doxgroups, mode = :density, label="", marker=(8, repeat([:orange, :darkblue, :darkred], inner=10), stroke(0)))
savefig("p53_meg3_median_boxplot.png")

p53_meg3_montecarlo = [montecarlomeasurements[:,:,1,1,2] |> vec; montecarlomeasurements[:,:,2,1,2] |> vec; montecarlomeasurements[:,:,1,2,2] |> vec; montecarlomeasurements[:,:,2,2,2] |> vec]
groupedboxplot(rnagroups, p53_meg3_montecarlo, group = doxgroups, outliers=false,
        label=["- Dox", "+ Dox"],
        guidefontsize=12,
        tickfontsize=72,
        legend=:none,
        left_margin=20mm,
        top_margin=5mm,
        bottom_margin=5mm,
        seriescolor=[:white :lightgray],
        line=(6, 1.0),
        size=(1024,2048),
        xaxis=("RNA", (1:2, ["MEG3", "GAPDH"])),
        yaxis=("Fraction bound", (0,0.2)))
groupeddotplot!(rnagroups, p53_meg3_montecarlo, group = doxgroups, mode = :density, label="", marker=(8, repeat([:orange, :darkblue, :darkred], inner=10), stroke(0)))
savefig(joinpath(outputdir, "p53_meg3_montecarlo_boxplot.png"))

p53_meg3_montecarlo = [normalizedmontecarlomeasurements[:,:,1,1,2] |> vec; normalizedmontecarlomeasurements[:,:,2,1,2] |> vec; normalizedmontecarlomeasurements[:,:,1,2,2] |> vec; normalizedmontecarlomeasurements[:,:,2,2,2] |> vec]
groupedboxplot(rnagroups, p53_meg3_montecarlo, group = doxgroups, outliers=false,
        label=["- Dox", "+ Dox"],
        guidefontsize=52,
        tickfontsize=48,
        legend=:none,
        top_margin=10mm,
        bottom_margin=-25mm,
        seriescolor=[:white :lightgray],
        line=(6, 1.0),
        gridopacity=0.3,
        xgrid=:none, #######################
        size=(1024,2048),
        xaxis=("RNA", (1:2, ["MEG3", "GAPDH"])),
        yaxis=("Fraction bound", (-0.05,0.2), -0.05:0.05:0.2))
groupeddotplot!(rnagroups, p53_meg3_montecarlo, group = doxgroups, mode = :density, label="",
        marker=(8, repeat([:orange, :darkblue, :darkred], inner=10), stroke(0)))
savefig(joinpath(outputdir, "p53_meg3_normalizedmontecarlo_boxplot.png"))


text:
immutable CacheChunkVector{T,A,N} <: AbstractVector{T}
  cache::CacheStack{CacheVector{T,A}}
  data::ChunkedVector{T,CacheVector{T,A},N}
end

function (::Type{CacheChunkVector{T,A,N}}){T,A,N}()
  cache = CacheStack{CacheVector{T,A}}()
  data = ChunkedVector{T,CacheVector{T,A},N}(CacheVector{T,A}[CacheVector(cache, A())])
  CacheChunkVector{T,A,N}(cache, data)
end

@forward CacheChunkVector.data Base.size, Base.getindex, Base.setindex!, Base.push!,
  ChunkIter, chunks

@iter xs::CacheChunkVector -> ChunkIter(xs)

function Base.start(xs::CacheChunkVector)
  it = ChunkIter(xs)
  touch!(chunks(xs)[1])
  return SubIter(it, start(it))
end

@inline function Base.next(xs::CacheChunkVector, sub::SubIter)
  i, j = sub.state
  if j == chunksize(sub.iter)
    i, j = i+1, 1
    @inbounds touch!(chunks(xs)[i])
  else
    i, j = i, j+1
  end
  @inbounds return chunks(xs)[i].view[j], SubIter(sub.iter, (i, j))
end

# Convenience Alias

export BigVector

immutable BigVector{T,N} <: AVector{T}
  data::CacheChunkVector{T,DiskVector{T},N}
end

(::Type{BigVector{T}}){T}(chunk::Integer = prevpow2(150*1024^2÷sizeof(T))) =
  BigVector{T,chunk}(CacheChunkVector{T,DiskVector{T},chunk}())

function BigVector(xs)
  v = BigVector{eltype(xs)}()
  for x in xs push!(v, x) end
  return v
end

@forward BigVector.data Base.size, Base.getindex, Base.setindex!, Base.push!, chunks

@iter xs::BigVector -> xs.data

# Merge sort

function merge{T}(xs::AVector{T}, ys::AVector{T})
  v = BigVector{T}()
  ix, nx, iy, ny = 1, length(xs), 1, length(ys)
  while true
    if ix ≤ nx && (iy > ny || xs[ix] ≤ ys[iy])
      push!(v, xs[ix])
      ix += 1
    elseif iy ≤ ny && (ix > nx || ys[iy] ≤ xs[ix])
      push!(v, ys[iy])
      iy += 1
    else
      break
    end
  end
  return v
end

function sort_mem(xs)
  ys = BigVector(sort!(copy(xs)))
  gc()
  return ys
end

function Base.sort(xs::BigVector)
  n = 4*1024^3÷sizeof(eltype(xs))
  if length(xs) ≤ n
    return sort_mem(xs)
  else
    left, right = slice(xs, 1:length(xs)÷2), slice(xs, length(xs)÷2+1:length(xs))
    if length(left) ≤ n
      return merge(sort_mem(left), sort_mem(right))
    else
      return merge(sort(BigVector(left)), sort(BigVector(right)))
    end
  end
end

ext: jl


text:
###
### writemime
###
# only mime writeable to PNG if 2D (used by IJulia for example)
mimewritable(::MIME"image/png", img::AbstractImage) = sdims(img) == 2 && timedim(img) == 0
mimewritable{C<:Colorant}(::MIME"image/png", img::AbstractArray{C}) = sdims(img) == 2 && timedim(img) == 0
# Colors.jl turns on SVG display of colors, which leads to poor
# performance and weird spacing if you're displaying images. We need
# to disable that here.
# See https://github.com/JuliaLang/IJulia.jl/issues/229 and Images #548
mimewritable(::MIME"image/svg+xml", img::AbstractImage)  = false
mimewritable{C<:Color}(::MIME"image/svg+xml", img::AbstractMatrix{C}) = false

# This is used for output by IJulia. Really large images can make
# display very slow, so we shrink big images.  Conversely, tiny images
# don't show up well, so in such cases we repeat pixels.
@compat function Base.show(io::IO, mime::MIME"image/png", img::AbstractImage; mapi=mapinfo_writemime(img), minpixels=10^4, maxpixels=10^6)
    assert2d(img)
    A = data(img)
    nc = ncolorelem(img)
    npix = length(A)/nc
    while npix > maxpixels
        # Big images
        A = restrict(A, coords_spatial(img))
        npix = length(A)/nc
    end
    if npix < minpixels
        # Tiny images
        fac = ceil(Int, sqrt(minpixels/npix))
        r = ones(Int, ndims(img))
        r[coords_spatial(img)] = fac
        A = repeat(A, inner=r)
    end
    imgcopy = shareproperties(img, A)
    save(Stream(format"PNG", io), imgcopy)
end

@compat function Base.show(stream::IO, mime::MIME"image/png", img::AbstractImageIndexed; kwargs...)
    @compat show(stream, mime, convert(Image, img); kwargs...)
end
@compat function Base.show{C<:Colorant}(stream::IO, mime::MIME"image/png", img::AbstractMatrix{C}; kwargs...)
    @compat show(stream, mime, Image(img, spatialorder=["y","x"]); kwargs...)
end

function mapinfo_writemime(img; maxpixels=10^6)
    if length(img) <= maxpixels
        return mapinfo_writemime_(img)
    end
    mapinfo_writemime_restricted(img)
end

to_native_color{T<:Colorant}(::Type{T}) = base_color_type(T){UFixed8}
to_native_color{T<:Color}(::Type{T}) = RGB{UFixed8}
to_native_color{T<:TransparentColor}(::Type{T}) = RGBA{UFixed8}

mapinfo_writemime_{T <:Colorant}(img::AbstractImage{T}) = Images.mapinfo(to_native_color(T), img)
mapinfo_writemime_(img::AbstractImage) = Images.mapinfo(UFixed8,img)

mapinfo_writemime_restricted{T<:Colorant}(img::AbstractImage{T}) = ClampMinMax(to_native_color(T), 0.0, 1.0)
mapinfo_writemime_restricted(img::AbstractImage) = Images.mapinfo(UFixed8, img)


text:
module Applications
using  ...Destructuring.Function
using  ...Helper
import Base.Meta: isexpr, quot
export @letds, @macrods, @anonds, @funds

DA = Applications

#----------------------------------------------------------------------------
# Utilities
#----------------------------------------------------------------------------

function destructure_all(patterns, values, body)
  destructure(:($(patterns...),), :(Expr(:tuple, $(values...))), body)
end

#----------------------------------------------------------------------------
# Destructuring let
#----------------------------------------------------------------------------

macro letds(lets...)
  esc(letds(lets[1:end-1], lets[end]))
end

function letds(lets, body)
  patterns, values = unzip(map(x->x.args, lets))
  destructure_all(patterns, values, body)
end

#----------------------------------------------------------------------------
# Destructuring macro
#----------------------------------------------------------------------------

macro macrods(signature, body)
  name      = signature.args[1]
  patterns  = Expr(:tuple, signature.args[2:end]...)

  esc(:(macro $name(args...)
          $DA.@letds $patterns = Expr(:tuple, args...) $body end))
end

#----------------------------------------------------------------------------
# Destructuring anonymous function
#----------------------------------------------------------------------------

@macrods anonds(patterns -> body) begin
  esc(anonds(patterns, body))
end

anonds(p::Symbol, body) = anonds_impl(:($p,), body)
anonds(p::Expr,  body)  = anonds_impl(p.head == :tuple ? p : :($p,), body)

function anonds_impl(patterns, body)
  :((args...) ->
      $DA.@letds $patterns = Expr(:tuple, args...) $body)
end

#----------------------------------------------------------------------------
# Destructuring named function
#----------------------------------------------------------------------------

@macrods funds(name(*{patterns}) = body) begin
  esc(funds(name, patterns, body))
end

function funds(name, patterns, body)
  :(function $name(args...)
      $DA.@letds $(Expr(:tuple, patterns...)) = Expr(:tuple, args...) $body end)
end

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

end


text:
@testset "1005.maximize-sum-of-array-after-k-negations.jl" begin
    @test largest_sum_after_k_negations([4,2,3], 1) == 5
    @test largest_sum_after_k_negations([3,-1,0,2], 3) == 6
    @test largest_sum_after_k_negations([2,-3,-1,5,-4], 2) == 13
end


text:
mutable struct REINFORCE{T,TP,TV} <: AbstractPolicyAgent where{T<:Real,TP<:AbstractPolicy,TV<:AbstractFuncApprox}
    π::TP
    vf::TV
    αₚ::T
    αᵥ::T
    γ::T
    avgGrad::Bool
    discountGrad::Bool
    states::Vector
    actions::Vector
    rewards::Array{T, 1}
    function REINFORCE(π::AbstractPolicy, vf::AbstractFuncApprox, αₚ::T, αᵥ::T, γ::T, avgGrad::Bool=false, discountGrad::Bool=true) where {T<:Real}
        new{T,typeof(π),typeof(vf)}(clone(π), clone(vf), αₚ, αᵥ, γ, avgGrad, discountGrad, [], [], T[])
    end
end

function act!(agent::REINFORCE, env::AbstractEnvironment, rng::AbstractRNG)
    push!(agent.states, env.state)

    get_action!(agent.π, env.state, rng)
    push!(agent.actions, agent.π.action)

    reward = step!(env, agent.π.action, rng)
    push!(agent.rewards, reward)

    if is_terminal(env)
        Gs = compute_returns(agent.rewards, agent.γ)
        update_withBaseline!(agent, Gs)
        new_episode!(agent, rng)
    end
end

function compute_returns(rewards::Array{T, 1}, γ::T)::Array{T, 1} where {T<:Real}
    G = 0.
    Gs = zeros(T, length(rewards))
    for t in length(rewards):-1:1
        G = rewards[t] + γ * G
        Gs[t] = G
    end
    return Gs
end

function update_withBaseline!(agent::REINFORCE, Gs::Array{T, 1}) where {T<:Real}
    γₜ = 1.0
    gₚ = zeros(get_num_params(agent.π))
    gᵥ = zeros(get_num_params(agent.vf))
    ∇ = zeros(get_num_params(agent.π))
    for t in 1:length(agent.states)
        v = call_gradient!(gᵥ, agent.vf, agent.states[t], 1)
        δ = Gs[t] - v
        @. gᵥ *= δ * agent.αᵥ
        add_to_params!(agent.vf, gᵥ)

        gradient_logp!(gₚ, agent.π, agent.states[t], agent.actions[t])
        @. ∇ += δ * agent.αₚ * γₜ * gₚ

        if agent.discountGrad
            γₜ *= agent.γ
        end
    end
    if agent.avgGrad
        ∇ ./= length(agent.states)
    end
    add_to_params!(agent.π, ∇)
end

function new_episode!(agent::REINFORCE, rng::AbstractRNG)
    empty!(agent.states)
    empty!(agent.actions)
    empty!(agent.rewards)
end

function clone(agent::REINFORCE)::REINFORCE
    a = REINFORCE(agent.π, agent.vf, agent.αₚ, agent.αᵥ, agent.γ, agent.avgGrad, agent.discountGrad)
    a.states = deepcopy(agent.states)
    a.actions = deepcopy(agent.actions)
    a.rewards = deepcopy(agent.rewards)
    return a
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AvailableProvidersListParameters <: SwaggerModel
    azureLocations::Any # spec type: Union{ Nothing, Vector{String} } # spec name: azureLocations
    country::Any # spec type: Union{ Nothing, String } # spec name: country
    state::Any # spec type: Union{ Nothing, String } # spec name: state
    city::Any # spec type: Union{ Nothing, String } # spec name: city

    function AvailableProvidersListParameters(;azureLocations=nothing, country=nothing, state=nothing, city=nothing)
        o = new()
        validate_property(AvailableProvidersListParameters, Symbol("azureLocations"), azureLocations)
        setfield!(o, Symbol("azureLocations"), azureLocations)
        validate_property(AvailableProvidersListParameters, Symbol("country"), country)
        setfield!(o, Symbol("country"), country)
        validate_property(AvailableProvidersListParameters, Symbol("state"), state)
        setfield!(o, Symbol("state"), state)
        validate_property(AvailableProvidersListParameters, Symbol("city"), city)
        setfield!(o, Symbol("city"), city)
        o
    end
end # type AvailableProvidersListParameters

const _property_map_AvailableProvidersListParameters = Dict{Symbol,Symbol}(Symbol("azureLocations")=>Symbol("azureLocations"), Symbol("country")=>Symbol("country"), Symbol("state")=>Symbol("state"), Symbol("city")=>Symbol("city"))
const _property_types_AvailableProvidersListParameters = Dict{Symbol,String}(Symbol("azureLocations")=>"Vector{String}", Symbol("country")=>"String", Symbol("state")=>"String", Symbol("city")=>"String")
Base.propertynames(::Type{ AvailableProvidersListParameters }) = collect(keys(_property_map_AvailableProvidersListParameters))
Swagger.property_type(::Type{ AvailableProvidersListParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_AvailableProvidersListParameters[name]))}
Swagger.field_name(::Type{ AvailableProvidersListParameters }, property_name::Symbol) =  _property_map_AvailableProvidersListParameters[property_name]

function check_required(o::AvailableProvidersListParameters)
    true
end

function validate_property(::Type{ AvailableProvidersListParameters }, name::Symbol, val)
end


text:
function is_palindrome(n)
  reversed, temp = 0, n

  while temp != 0
    reversed = 10 * reversed + temp % 10;
    temp = div(temp, 10)
  end

  reversed == n
end

function solution()
  max = 0

  for a in 100:999
    for b in 100:a
      p = a * b

      if p > max && is_palindrome(p)
        max = p
      end
    end
  end

  max
end

if length(ARGS) == 1 && ARGS[1] == "-a"
  println(solution())
  exit(0)
end

for line in eachline(STDIN)
  iters = int(line)

  start = time_ns()
  for i in 1:iters
    solution()
  end
  end_ = time_ns()

  println(end_ - start)
end


text:
using LazySets
"""
    PolytopeComplement

The complement to a given set. Note that in general, a `PolytopeComplement` is not necessarily a convex set.
Also note that `PolytopeComplement`s are open by definition.

### Examples
```julia
julia> H = Hyperrectangle([0,0], [1,1])
Hyperrectangle{Int64}([0, 0], [1, 1])

julia> PC = complement(H)
PolytopeComplement of:
  Hyperrectangle{Int64}([0, 0], [1, 1])

julia> center(H) ∈ PC
false

julia> high(H).+[1,1] ∈ PC
true
```
"""
struct PolytopeComplement{S<:LazySet}
    P::S
end

Base.show(io::IO, PC::PolytopeComplement) = (println(io, "PolytopeComplement of:"), println(io, "  ", PC.P))
LazySets.issubset(s, PC::PolytopeComplement) = LazySets.is_intersection_empty(s, PC.P)
LazySets.is_intersection_empty(s, PC::PolytopeComplement) = LazySets.issubset(s, PC.P)
LazySets.tohrep(PC::PolytopeComplement) = PolytopeComplement(convert(HPolytope, PC.P))
Base.in(pt, PC::PolytopeComplement) = pt ∉ PC.P
complement(PC::PolytopeComplement)  = PC.P
complement(P::LazySet) = PolytopeComplement(P)
# etc.

"""
    Problem{P, Q}(network::Network, input::P, output::Q)

Problem definition for neural verification.

The verification problem consists of: for all  points in the input set,
the corresponding output of the network must belong to the output set.
"""
struct Problem{P, Q}
    network::Network
    input::P
    output::Q
end

"""
    Result
Supertype of all result types.

See also: [`BasicResult`](@ref), [`CounterExampleResult`](@ref), [`AdversarialResult`](@ref), [`ReachabilityResult`](@ref)
"""
abstract type Result end

status(result::Result) = result.status

function validate_status(st::Symbol)
    @assert st ∈ (:holds, :violated, :unknown) "unexpected status code: `:$st`.\nOnly (:holds, :violated, :unknown) are accepted"
    return st
end

"""
    BasicResult(status::Symbol)

Result type that captures whether the input-output constraint is satisfied.
Possible status values:\n
    :holds (io constraint is satisfied always)\n
    :violated (io constraint is violated)\n
    :unknown (could not be determined)
"""
struct BasicResult <: Result
    status::Symbol
end

"""
    CounterExampleResult(status, counter_example)

Like `BasicResult`, but also returns a `counter_example` if one is found (if status = :violated).
The `counter_example` is a point in the input set that, after the NN, lies outside the output set.
"""
struct CounterExampleResult <: Result
    status::Symbol
    counter_example::Vector{Float64}
    CounterExampleResult(s, ce) = new(validate_status(s), ce)
end

"""
    AdversarialResult(status, max_disturbance)

Like `BasicResult`, but also returns the maximum allowable disturbance in the input (if status = :violated).
"""
struct AdversarialResult <: Result
	status::Symbol
	max_disturbance::Float64
    AdversarialResult(s, md) = new(validate_status(s), md)
end

"""
    ReachabilityResult(status, reachable)

Like `BasicResult`, but also returns the output reachable set given the input constraint (if status = :violated).
"""
struct ReachabilityResult <: Result
	status::Symbol
	reachable::Vector{<:AbstractPolytope}
    ReachabilityResult(s, r) = new(validate_status(s), r)
end

# Additional constructors:
CounterExampleResult(s) = CounterExampleResult(s, Float64[])
AdversarialResult(s)    = AdversarialResult(s, -1.0)
ReachabilityResult(s)   = ReachabilityResult(s, AbstractPolytope[])


text:
"""
`module  JAC.InternalConversion`
    ... a submodel of JAC that contains all methods for computing photo-excitation-autoionization cross sections and
        rates.
"""
module InternalConversion

    using Printf, ..AutoIonization, ..Basics, ..ManyElectron, ..Radial, ..PhotoEmission, ..TableStrings

    """
    `struct  InternalConversion.Settings`
        ... defines a type for the details and parameters of computing photon-impact excitation-autoionization
            pathways |i(N)>  --> |m(N)>  --> |f(N-1)>.

        + multipoles              ::Array{EmMultipole,1}               ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                  ::Array{UseGauge,1}                  ... Specifies the gauges to be included into the computations.
        + printBefore  ::Bool                               ... True, if all energies and lines are printed before their evaluation.
        + selectPathways          ::Bool                               ... True if particular pathways are selected for the computations.
        + selectedPathways        ::Array{Tuple{Int64,Int64,Int64},1}  ... List of list of pathways, given by tupels (inital, inmediate, final).
        + maxKappa                ::Int64                              ... Maximum kappa value of partial waves to be included.
    """
    struct Settings
        multipoles                ::Array{EmMultipole,1}
        gauges                    ::Array{UseGauge,1}
        printBefore    ::Bool
        selectPathways            ::Bool
        selectedPathways          ::Array{Tuple{Int64,Int64,Int64},1}
        maxKappa                  ::Int64
    end

    """
    `InternalConversion.Settings()`  ... constructor for the default values of photon-impact excitation-autoionizaton settings.
    """
    function Settings()
        Settings( Basics.EmMultipole[], Basics.UseGauge[], false,  false, Tuple{Int64,Int64,Int64}[], 0)
    end

    # `Base.show(io::IO, settings::InternalConversion.Settings)`
    #   ... prepares a proper printout of the variable settings::InternalConversion.Settings.
    function Base.show(io::IO, settings::InternalConversion.Settings)
        println(io, "multipoles:              $(settings.multipoles)  ")
        println(io, "gauges:                  $(settings.gauges)  ")
        println(io, "printBefore:  $(settings.printBefore)  ")
        println(io, "selectPathways:          $(settings.selectPathways)  ")
        println(io, "selectedPathways:        $(settings.selectedPathways)  ")
        println(io, "maxKappa:                $(settings.maxKappa)  ")
    end

    """
    `struct  InternalConversion.Channel`
        ... defines a type for a photon-impact excitaton & autoionization channel that specifies all quantum numbers,
            phases and amplitudes.

        + excitationChannel  ::PhotoEmission.Channel     ... Channel that describes the photon-impact excitation process.
        + augerChannel       ::AutoIonization.Channel    ... Channel that describes the subsequent Auger/autoionization process.
    """
    struct  Channel
        excitationChannel    ::PhotoEmission.Channel
        augerChannel         ::AutoIonization.Channel
    end

    """
    `struct  InternalConversion.Pathway`
        ... defines a type for a photon-impact excitation pathway that may include the definition of different excitation and
            autoionization channels and their corresponding amplitudes.

        + initialLevel        ::Level                  ... initial-(state) level
        + intermediateLevel   ::Level                  ... intermediate-(state) level
        + finalLevel          ::Level                  ... final-(state) level
        + photonEnergy        ::Float64                 ... energy of the (incoming) electron
        + electronEnergy      ::Float64                 ... energy of the (finally outgoing, scattered) electron
        + crossSection        ::EmProperty              ... total cross section of this pathway
        + hasChannels         ::Bool                    ... Determines whether the individual excitation and autoionization channels are defined
                                                            in terms of their multipole, gauge, free-electron kappa, phases and the total
                                                            angular momentum/parity as well as the amplitude, or not.
        + channels            ::Array{InternalConversion.Channel,1}  ... List of channels of this pathway.
    """
    struct  Pathway
        initialLevel          ::Level
        intermediateLevel     ::Level
        finalLevel            ::Level
        photonEnergy          ::Float64
        electronEnergy        ::Float64
        crossSection          ::EmProperty
        hasChannels           ::Bool
        channels              ::Array{InternalConversion.Channel,1}
    end

    """
    `InternalConversion.Pathway()`
        ... 'empty' constructor for an photon-impact excitation-autoionization pathway between a specified initial, intermediate and
            final level.
    """
    function Pathway()
        Pathway(Level(), Level(), Level(), 0., 0., 0., false, InternalConversion.Channel[] )
    end

    # `Base.show(io::IO, pathway::InternalConversion.Pathway)`
    # ... prepares a proper printout of the variable pathway::InternalConversion.Pathway.
    function Base.show(io::IO, pathway::InternalConversion.Pathway)
        println(io, "initialLevel:               $(pathway.initialLevel)  ")
        println(io, "intermediateLevel:          $(pathway.intermediateLevel)  ")
        println(io, "finalLevel:                 $(pathway.finalLevel)  ")
        println(io, "photonEnergy                $(pathway.photonEnergy)  ")
        println(io, "electronEnergy              $(pathway.electronEnergy)  ")
        println(io, "crossSection:               $(pathway.crossSection)  ")
        println(io, "hasChannels:                $(pathway.hasChannels)  ")
        println(io, "channels:                   $(pathway.channels)  ")
    end

    """
    `InternalConversion.computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                                     settings::InternalConversion.Settings; output=true)`
        ... to compute the multiphoton transition amplitudes and all properties as requested by the given settings. A list of
            lines::Array{InternalConversion.Lines} is returned.
    """
    function  computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                           settings::InternalConversion.Settings; output=true)
        println("")
        printstyled("InternalConversion.computeLines(): The computation of internal conversion amplitudes starts now ... \n", color=:light_green)
        printstyled("---------------------------------------------------- ---------------------------------------------- \n", color=:light_green)
        println("")
        #
        pathways = "Not yet implemented !"
        #
        if    output    return( pathways )
        else            return( nothing )
        end
    end

end # module


text:
function plot_obstacles(circles)
    for  circle in circles
        x,y,r = circle
        plot_circle!((x,y),r,color=:red,label="")
    end
end

function plot_trajectory!(X;kwargs...)
    plot!(X[1,:],X[2,:];kwargs...)
end

function plot_solution(X;kwargs...)
    p = plot(aspect_ratio=:equal,xlim=[0,12],ylim=[0,12])
    plot_obstacles(circles)
    plot_trajectory!(X; kwargs...)
end

### Solver Options ###
dt = 0.01
opts = TrajectoryOptimization.SolverOptions()
opts.square_root = false
opts.verbose = true
# opts.c1 = 1e-4
# opts.c2 = 2.0
opts.cost_tolerance_intermediate = 1e-1
opts.constraint_tolerance = 1e-6
opts.cost_tolerance = 1e-6
opts.iterations_outerloop = 100
opts.iterations = 1000
# opts.iterations_linesearch = 50
opts.constraint_decrease_ratio = 0.1
opts.outer_loop_update_type = :default
######################

### Set up model, objective, and solver ###
# Model
model, obj_uncon = TrajectoryOptimization.Dynamics.dubinscar!
n,m = model.n,model.m

# Objective
# -Circle obstacles, state, and control constraints
obj_uncon_obs = copy(obj_uncon)
obj_uncon_obs.x0[:] = [0;0;0]
obj_uncon_obs.xf[:] = [12;12;pi/2]
obj_uncon_obs.tf = 2.5
obj_uncon_obs.Qf[:,:] = 100.0*Diagonal(I,n)
obj_uncon_obs.R[:,:] = (1e-2)*Diagonal(I,m)

u_min = [-100; -100]
u_max = [100; 100]
x_min = [-20; -20; -100]
x_max = [20; 20; 100]

n_circles = 15
c = zeros(n_circles)
cI_obstacles, circles = generate_random_circle_obstacle_field(n_circles)

p = plot(aspect_ratio=:equal,xlim=[0,10],ylim=[0,10])
plot_obstacles(circles)
display(p)

# Solve with iLQR
obj_con_obs = TrajectoryOptimization.ConstrainedObjective(obj_uncon_obs,cI=cI_obstacles)# u_min=u_min, u_max=u_max,x_min=x_min, x_max=x_max,cI=cI) # constrained objective
solver = Solver(model, obj_con_obs, integration=:rk3, dt=dt, opts=opts)
solver.opts.verbose = false
solver.opts.resolve_feasible = false
n,m,N = get_sizes(solver)
X0 = line_trajectory(solver)
U = ones(m,N)
sol, stats = TrajectoryOptimization.solve(solver,X0,U)

plot_solution(to_array(sol.X),width=2,color=:black,label="iLQR")

# Solve with DIRCOL
res_d, stat_d = solve_dircol(solver,X0,U)

plot_trajectory!(res_d.X,width=2,color=:blue,linestyle=:dash,label="DIRCOL")

solver_uncon = Solver(model,obj_uncon_obs,dt=dt)
solver_uncon.opts.verbose = true
res_d, stat_d = solve_dircol(solver_uncon,X0,U)
plot(res_d.X[1,:],res_d.X[2,:])

plot_trajectory!

### Results ###
if opts.verbose
    println("Final state (foh): $(sol_foh_con_obs.X[:,end])")
    println("Final state (zoh): $(sol_zoh_con_obs.X[:,end])")

    println("Final cost (foh): $(sol_foh_con_obs.cost[sol_foh_con_obs.termination_index])")
    println("Final cost (zoh): $(sol_zoh_con_obs.cost[sol_zoh_con_obs.termination_index])")

    println("Termination index\n foh: $(sol_foh_con_obs.termination_index)\n zoh: $(sol_foh_con_obs.termination_index)")

    # plot(sol_foh_con_obs.U')
    # plot!(sol_zoh_con_obs.U')
    #
    # ## Plot obstacle field and trajectory
    # plot((solver_zoh_con_obs.obj.x0[1],solver_zoh_con_obs.obj.x0[2]),marker=(:circle,"red"),label="x0",xlim=(-1.1,11.1),ylim=(-1.1,11.1))
    # plot!((solver_zoh_con_obs.obj.xf[1],solver_zoh_con_obs.obj.xf[2]),marker=(:circle,"green"),label="xf")
    #
    # theta = linspace(0,2*pi,100)
    # for k = 1:n_circles
    #     x_circle = circles[3][k]*cos.(theta)
    #     y_circle = circles[3][k]*sin.(theta)
    #     plot!(x_circle+circles[1][k],y_circle+circles[2][k],color="red",width=2,fill=(100),legend=:none)
    # end
    #
    # plot!(sol_foh_con_obs.X[1,:],sol_foh_con_obs.X[2,:])
    # plot!(sol_zoh_con_obs.X[1,:],sol_zoh_con_obs.X[2,:])
end
###############


text:
module Utils

#
# This contains general purpose, messy stuff not easily fitting elsewhere.
# A fair chunk is no longer used :(
#
using ArgCheck
using Base: Integer, String, Bool
using DataFrames
using Dates
using Base.Unicode
using CSV
using BudgetConstraints
using Printf

export
   @exported_enum,
   BR_DIGITS,
   ≈,
   addsysnotoname,
   age_in_years,
   age_then,
   basiccensor,
   coarse_match,
   date_string,
   diff_between,
   df_diff,
   eq_nearest_p,
   extract_digits,
   get_if_set,
   get_project_path,
   has_non_z,
   haskeys,
   index_of_field,
   is_zero_or_missing,
   isapprox,
   isordered,
   loadtoframe,
   make_start_stops,
   mult_dict!,
   mult,
   md_format,
   nearest,
   nearesti,
   nearz,
   not_zero_or_missing,
   operate_on_struct!,
   pretty,
   qstrtodict,
   riskyhash,
   to_md_table,
   todays_date,
   uprate_struct!

"""
Make a new dataframe with the difference between the fields between
start_col and end_col and other fields just copied
diff is difference df2 - df1, Frames should have identical other cols.
"""
function df_diff( df1, df2 :: DataFrame, start_col::Int, end_col :: Int  ) :: DataFrame
   @argcheck size( df1 ) == size( df2 )
   ## maybe check that the non diffed fields are all the same too..
   d = copy(df1)
   d[:,start_col:end_col] = df2[:,start_col:end_col] .- df1[:,start_col:end_col]
   return d
end

"""
col number of given field name in a dataframe
"""
function index_of_field( d :: DataFrame, name :: String ) :: Int
   n = names(d)
   return findfirst( n .== name)
end

"""
For threading households
Dup of function in SurveyDataWeighting
"""
function make_start_stops( nrows::Int, num_threads::Int )::Tuple
   start = zeros(Int, num_threads)
   stop = zeros(Int, num_threads)
   chunk_size = Int(trunc(nrows/num_threads))-1;
   p = 1;
   for i in 1:num_threads
         start[i] = p
         p = p + chunk_size;
         if i < num_threads;
            stop[i] = p;
         else;
            stop[i] = nrows; # possibly different number on last thread
         end;
         p = p + 1;
   end;
   return (start, stop)
end

function extract_digits( b :: Integer, r :: UnitRange ) :: Integer
   parse(Int,string(b)[r])
end

"""
get some digits from a Digits, from the left, viewed as decimal.
e.g extract_digits( 7890, 2,3 ) = 89
extract_digits( 78901, 1,2 ) = 78
"""
function extract_digits( b :: Integer, fd::Int, ld::Int ) :: Integer
   # crude but effective
   extract_digits(b, [fd:ld])
end

"""
The index of the DataFrame row with a date field nearest to Date `d`, assuming dates are
sorted low->high. Should really be a binary search but isn't quite.
"""
function nearest( d::Date, df :: DataFrame, col::Symbol=:Date ) :: Integer
   pos = 0
   n = size(df)[1];
   p = -1;
   mid = n÷2
   mind = d - df[mid,col]
   if mind == Day(0)
       return mid
   elseif mind < Day(0)
       direction = -1
       las = 1
       mid -= 1
       if mid < 1 # idiot check for v short lists
           return 1
       end
   else
       direction = 1
       las = n
       mid += 1
       if mid > n # idiot check for v short lists
           return n
       end
   end
   mind = abs(mind)
   pmind = mind
   # println( "p=$p mid=$mid direction=$direction las=$las mind=$mind")
   for i in mid:direction:las
       dd = df[i,col]
       diff = dd - d
       diff = abs( diff ); # print( diff )
       if diff < mind
         mind = diff
         p = i
       end
       if pmind < mind # we've gone past
         # println( "returning pmind=$pmind mind=$mind")
         return p
       end
       pmind = mind
  end
  return p
end

function not_zero_or_missing( thing :: Union{Missing, Number }) :: Bool
   if ismissing( thing )
      return false
   end
   return thing != 0
end

function is_zero_or_missing( thing :: Union{Missing, Number }) :: Bool
   if ismissing( thing )
      return true
   end
   return thing == 0
end

"""
return the index of element of comps that's closest to x
"""
function nearesti( x :: Real , comps... )
   n = length( comps )
   dist = zeros( n )
   for i in 1:n
       dist[i] = abs( x - comps[i])
   end
   i = findmin( dist )
   return i[1]
end

"""
return the element of comps that's closest to x
"""
function nearest( x :: Real , comps... )
   return comps[ nearesti( x, comps )]
end

"""
The thing (comps...) that x is closest to, but returning zero even
if 0 isn't in comps..
"""
function nearz( x :: Real, comps ... ) :: Real
    if x ≈ 0
        return 0.0
    end
    return nearest( x, comps ... )
end

#
# this has a higher top income than the BC default
#
const DSETT = BudgetConstraints.DEFAULT_SETTINGS
const BC_SETTINGS = BCSettings(0.0,20_000.0,DSETT.increment,DSETT.tolerance,true,DSETT.maxdepth)

"""
finds the matches for a single recipient tuple `recip` in a data set `donor`.

each of the recip and donor should be structured as follows

firstvar_1, firstvar_2, firstvar_2 <- progressively coarsened first variable with the `_1` needed exactly as is;
then secondvar_1 .. thirdvar_1 .. _2 and so on. Variables can actually be in any order in the frame.

`vars` list of `firstvar`, `secondvar` and so on, in the order you want them coarsened
`max_coarsens` stop after _2, _3 etc. coarsened variables.

returns a named tuple:
     matches->indexes of rows that match
     quality->numerical indication of quality of each match: +1 for each first match, +4 for each 2nd +9 3rd, etc.
     results_matrix->nobs x nvars matrix with each cell indicating match quality (1..max_coarsens, or -9 for no match).
"""
function coarse_match(
    recip :: DataFrameRow,
    donor :: DataFrame,
    vars  :: Vector{Symbol},
    max_coarsens:: Int ) :: NamedTuple
    nobs = size( donor )[1]
    nvars = size( vars )[1]
    matches = fill( true, nobs )
    results = fill(-99, nobs,nvars)
    # println( size( results ))
    quality = fill( 0, nobs )
    for row in 1:nobs
        for col in 1:nvars
            for cs in 1:max_coarsens
                sym = Symbol("$(String(vars[col]))_$(cs)")
                if(donor[row,sym] == recip[sym])
                    results[row,col] = cs
                    break;
                end
            end # coarsens
        end # match vars
    end # rows
    quality = fill( 0, nobs )
    for row in 1:nobs
        for col in 1:nvars
            if results[row,col] == -99 # no match
                matches[row] = false
                quality[row] = -9
                break;
            else
                quality[row] += results[row,col]^2 # maybe
            end
        end # no match for some matching var break
    end # check each row
    return (matches=matches,results_matrix=results,quality=quality)
end

date_string() = lowercase(Dates.format( todays_date(), "d-u-Y"))

function todays_date() :: Date
    Date( now())
end

"""
Age in years now of someone born on some date.
Annoyingly difficult to get exactly right...
"""
function age_now( born :: Date ) :: Integer
   ago = todays_date()-born
   return Int(trunc(Dates.value( todays_date()-born )/365.25))
end

"""
These are trivial but annoying to get right.
"""
function age_then( age :: Int, when :: Int )::Int
    ago = Dates.year(todays_date()) - when
    age - ago
end

function age_then( age :: Int, when :: Date = todays_date() ) :: Int
    age_then( age, Dates.year( when ))
end

function eq_nearest_p( a :: Real, b :: Real )
   if round( a, digits=2) != round( b, digits=2 )
      println( "NOT Equal a = $a b= $b ")
      return false
   end
   return true
end

@generated function operate_on_struct!( rec::Rec, x::NR, f::Function ) where Rec where NR<:Number
     assignments = [
         :( rec.$name = f(rec.$name, x) ) for name in fieldnames(Rec)
     ]
     quote $(assignments...) end
end

function mult( v :: AbstractFloat, x :: NR )::NR where NR <: Number
   v*x
end

function mult( v :: AbstractArray, x :: NR )::AbstractArray where NR <: Number
   mult.( v, x )
end

function mult( v::Any, x :: NR )::Any where NR <: Number
   v
end

function m( v :: Any, x :: NR ) where NR <: Number
   mult( v, x )
end

function uprate_struct!( rec::Rec, x::NR ) where Rec where NR<:Number
   operate_on_struct!( rec, x, m )
end

function addsysnotoname(names, sysno)::Array{Symbol,1}
   a = Array{Symbol,1}(undef, 0)
   for n in names
      push!(a, Symbol(String(n) * "_$sysno"))
   end
   a
    # Symbol.(String.( names ).*"_sys_$sysno")
end

"""
Does the ordering here follow stye guide
"""
function get_if_set(dict::Dict, key, default; operation=nothing)
   v = default
   if haskey(dict, key)
      v = dict[key]
      if operation !== nothing
         v = operation( v )
      end
   end
   v
end

import Base.isapprox
import Base.≈

"""
true if the keys are the same and all the elements
   compare approx equal, else false
"""
function isapprox( d1::Dict, d2::Dict ) :: Bool
   k1 = Set( collect( keys(d1)))
   k2 = Set( collect( keys(d2)))
   if k1 != k2
      # println(" resturn ne $(k1) !== $(k2)")
      return false
   end
   for k in k1
      if ! (d1[k] ≈ d2[k])
         # println( "return ne $k $(d1[k]) ≈ $(d2[k])")
         return false
      end
   end
   return true
end

#const ≈ = isapprox

"""
return a dict with all the elements in common where m2[k]-m1[k] is possible.
"""
function diff_between(m2::Dict, m1::Dict)::Dict
   out = Dict()
   keybs = intersect(keys(m1), keys(m2))
   for k in keybs
      v1 = m1[k]
      v2 = m2[k]
      try
         d = v2 - v1
         out[k] = d
      catch

      end # exception
   end # loop[]
   out
end

"""
This is for multiplying incomes where there
   are some (raw) from the data and some,
   possibly the same elements, from earlier
   calculations (net)
"""
function mult(; data::Dict{K,T}, calculated :: Dict{K,T}, included :: Dict{K,T}) :: T where T<:Number where K
   s = zero(T)
   kg = keys( data )
   kn = keys( calculated )
   ki = keys( included )
   ka = union( kg, kn )
   ka = intersect( ka, ki )
   for k in ka
      # choose the calculated value if there is one, otherwise the data one
      i = k in kn ? calculated[k] : data[k]
      s += i*included[k]
   end
   return s;
end # mult

function has_non_z( d::AbstractDict, k ) :: Bool
   if haskey( d, k )
      if typeof(d[k]) <: Number
         return d[k] != 0
      end
   end
   return false
end

function mult_dict!( m :: Dict, n :: Number )
   for (k,v) in m
      try
         m[k] = v*n
      catch
         ;
      end
   end
end

"""
parse an html query string like "sex=male&joe=22&bill=21342&z=12.20"
into a dict. If the value looks like a number, it's parsed into either an Int64 or Float64
"""
function qstrtodict(query_string::AbstractString)::Dict{AbstractString,Any}
   d = Dict{AbstractString,Any}()
   strip(query_string)
   if (query_string == "")
      return d
   end
   as = split(query_string, "&")
   for a in as
      try
         aa = split(a, "=")
         k = aa[1]
         v = aa[2]
         try
            v = parse(Int64, v)
         catch
            try
               v = parse(Float64, v)
            catch
            end
         end
         if v in ["false","False","FALSE"]
            v = false
         elseif v in ["true","True","TRUE"]
            v = true
         end
         d[k] = v
      catch

      end
   end
   d
end

function haskeys( d :: AbstractDict, keys ... ) :: Bool
   for k in keys
      if haskey(d, k)
         return true
      end
   end
   return false
end

function haskeys( d :: AbstractDict, keys :: AbstractSet ) :: Bool
   for k in keys
      if haskey(d, k)
         return true
      end
   end
   return false
end

function haskeys( d :: AbstractDict, keys :: AbstractArray ) :: Bool
   for k in keys
      if haskey(d, k)
         return true
      end
   end
   return false
end

function isordered( v :: AbstractArray )
   n = size(v)[1]
   for i in 2:n
      if v[i] < v[i]-1
         return false
      end
   end
   return true
end

"""
returns the string converted to a form suitable to be used as (e.g.) a Symbol,
with leading/trailing blanks removed, forced to lowercase, and with various
characters replaced with '_' (at most '_' in a run).
"""
function basiccensor(s::AbstractString)::AbstractString
   s = strip(lowercase(s))
   s = replace(s, r"[ \-,\t–]" => "_")
   s = replace(s, r"[=\:\)\('’‘]" => "")
   s = replace(s, r"[\";:\.\?\*”“]" => "")
   s = replace(s, r"_$" => "")
   s = replace(s, r"^_" => "")
   s = replace(s, r"^_" => "")
   s = replace(s, r"\/" => "_or_")
   s = replace(s, r"\&" => "_and_")
   s = replace(s, r"\+" => "_plus_")
   s = replace(s, r"_\$+$" => "")
   if occursin(r"^[\d].*", s)
      s = string("v_", s) # leading digit
   end
   s = replace(s, r"__" => "_")
   s = replace(s, r"__" => "_")
   s = replace(s, r"__" => "_") # fixme neater way?
   s = replace(s, r"^_" => "")
   s = replace(s, r"_$" => "")
   return s
end

"""
a_string_or_symbol_like_this => "A String Or Symbol Like This"
"""
function pretty(a)
   s = string(a)
   s = strip(lowercase(s))
   s = replace(s, r"[_]" => " ")
   Unicode.titlecase(s)
end

"""
 macro to define an enum and automatically
 add export statements for its elements
 see: https://discourse.julialang.org/t/export-enum/5396
"""
macro exported_enum(name, args...)
   esc(quote
      @enum($name, $(args...))
      export $name
      for a in $args
         local av = string(a)
         :(export $av)
      end
   end)
end

"""
load a file into a dataframe and force all the identifiers into
lower case
"""
function loadtoframe(filename::AbstractString)::DataFrame
   println( "loading $filename")
    df = CSV.File(filename, delim = '\t') |> DataFrame #
    lcnames = Symbol.(lowercase.(string.(names(df))))
    rename!(df, lcnames)
    df
end

"""
Age now for someone with this birthday.
If today is 2020-01-27 then:

   * dob(1958-01-26) = 62
   * dob(1958-01-27) = 62
   * dob(1958-01-28) = 61

And so on.
"""
function age_in_years(
   dob :: Dates.TimeType,
   to_date :: Dates.TimeType = Dates.now() ) :: Integer
   @argcheck dob <= to_date
   y_to = year(to_date)
   m_to = month(to_date)
   d_to = day(to_date)
   y_dob = year(dob)
   m_dob = month(dob)
   d_dob = day(dob)
   age = y_to - y_dob
   if m_dob > m_to # check if you've not yet had your birthday ..
      age -= 1
   elseif (m_dob == m_to ) && (d_dob > d_to )
      age -= 1
   end
   age
end

"""
1234.456 => 1,234.45
"""
function format_delimited( n :: Number; groupsize :: Int = 3, delim :: Char = ',', decimal_delim :: Char = '.', prec :: Int = 2 ) :: String

end

# FIXME!! THIS TS TERRIBLE..
function get_project_path()
   path = splitpath(pwd())
   n = size(path)[1]
   if path[end] == "test"
      n -= 1
   end
   join( path[1:n],"/")*"/"
end

function md_format( a :: Union{AbstractArray,Tuple} )::String
   s = ""
   for i in eachindex(a)
      vs = md_format(a[i])
      s *= "[$i = $(vs)]"
      if i != lastindex(a)
         s *= ", "
      end
   end
   return s
end

function md_format( a :: AbstractDict )::String
   s = ""
   n = length(a)
   i = 0
   for (k,v) in a
      i += 1
      vs = md_format(v)
      s *= "[$k = $vs]"
      if i != n
         s *= ", "
      end
   end
   return s
end

function md_format( a :: AbstractFloat )::String
   return @sprintf( "%0.2f", a)
end

function md_format( a )
   "$a"
end

function is_a_struct( T::Type )::Bool
   if ! isstructtype( T )
      return false
   end
   return ! (( T <: AbstractArray)||(T<:AbstractDict)||(T<:Real)||(T<:AbstractString)||(T<:Symbol))
end

"""
Crude but more-or-less effective thing that prints out a struct (which may contain other structs) as
a markdown table.
"""
function to_md_table( f; exclude=[], depth=0 ) :: String
    F = typeof(f)
    @assert isstructtype( F )
    names = fieldnames(F)
    prinames = []
    structnames = []

    for n in names
        v = getfield(f,n)
        T = typeof(v)
        if n in exclude
            ;
        elseif is_a_struct( T )
            push!(structnames, n )
        else
            push!(prinames, n )
        end
    end
    s = """

    |            |              |
    |:-----------|-------------:|
    """
    for n in prinames
      v = getfield(f,n)
      pn = pretty(n)
      vs = md_format(v)
      s *= "|**$(pn)**|$vs|\n"
    end
    s *= "\n\n"
    depth += 1
    for n in structnames
        v = getfield(f,n)
        s *= "#"*repeat( "#", depth ) * " " * pretty(n) *"\n"
        s *= to_md_table( v, exclude=exclude, depth=depth )
    end
    return s;
end

"""
Hash of all (??) the elements of the things in the struct `things`
"""
function riskyhash( f, h :: UInt = UInt(0) ) :: UInt
      F = typeof(f)
      @assert isstructtype( F )
      names = fieldnames(F)
      prinames = []
      structnames = []

      for n in names
         v = getfield(f,n)
         T = typeof(v)
         if is_a_struct( T )
            push!(structnames, n )
         else
            push!(prinames, n )
         end
      end
      for n in prinames
         v = getfield(f,n)
         h = hash( v, h )
      end
      for n in structnames
         v = getfield(f,n)
         h = riskyhash( v, h )
      end
   return h
end # riskyhash

"""
Hash value of a collection of things, each of which should be
a struct of some sort. I'm sure a Proper Developer would decide that
this is A Bad Thing, but it's useful for caching in some visualisation apps.
"""
function riskyhash( things :: AbstractVector )
   h = UInt(0)
   for f in things
      h = riskyhash( f, h )
   end
   h
end

end # module

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
using PopGen
using Test

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


text:
module Util

using Distributions
import Random
using StatsFuns

include("SkewT.jl")
include("misc.jl")

end # module

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
module TestAnchors

using Reagents
using Reagents.Internal: withanchor, nextanchor!, finish!, Offer, Reaction
using Test
using ..TestCancellableContainers: blocking_treiberstack

function module_context(f)
    Reagents.Internal.use_anchors() && return f()
    Reagents.Internal.enable_anchors()
    try
        Base.invokelatest(f)
    finally
        Reagents.Internal.disable_anchors()
    end
end

function can_test()
    if Reagents.Internal.use_anchors()
        return true
    else
        @test_skip nothing
        return false
    end
end

function test_simple_commit()
    can_test() || return
    ref1 = Reagents.Ref{Int}(0)
    ref2 = Reagents.Ref{Int}(0)
    # Using two CASes to avoid the single-CAS optimization:
    reagent = Reagents.CAS(ref1, 0, 111) ⨟ Reagents.CAS(ref2, 0, 222)
    anc = withanchor() do
        reagent()
    end
    k, v = nextanchor!(anc)
    @test k === :commit
    @test v.offer === nothing
    @test v.rx isa Reaction
    @test ref1[] == ref2[] == 0
    finish!(anc)
    @test ref1[] == 111
    @test ref2[] == 222
end

function test_simple_race()
    can_test() || return
    ref1 = Reagents.Ref{Int}(0)
    ref2 = Reagents.Ref{Int}(0)
    # Using two CASes to avoid the single-CAS optimization:
    r1 = Reagents.CAS(ref1, 0, 111) ⨟ Reagents.CAS(ref2, 0, 222)
    r2 = Reagents.CAS(ref1, 0, 333) ⨟ Reagents.CAS(ref2, 0, 444)
    a1 = withanchor() do
        Reagents.try(r1)
    end
    a2 = withanchor() do
        Reagents.try(r2)
    end
    k1, _ = nextanchor!(a1)
    k2, _ = nextanchor!(a2)
    @test k1 == k2 === :commit
    @test ref1[] == ref2[] == 0
    @test finish!(a2) !== nothing
    @test ref1[] == 333
    @test ref2[] == 444
    @test finish!(a1) === nothing
end

function test_blocking_treiberstack_late_take_win()
    can_test() || return
    b = blocking_treiberstack(Int)
    a1 = withanchor(==(:commit)) do
        put!(b, 111)
    end
    a2 = withanchor(==(:commit)) do
        take!(b)
    end
    nextanchor!(a1)
    nextanchor!(a2)
    @test b.data.head[] === nothing
    @test finish!(a2) == 111
    @test b.data.head[] === nothing
    @test finish!(a1) === nothing
    @test b.data.head[] === nothing
end

end  # module


text:
### A Pluto.jl notebook ###
# v0.14.2

using Markdown
using InteractiveUtils

# ╔═╡ 44c12b60-97e4-11eb-23f7-eb2bf0847809
using Plots, SSpline, SpectralVis

# ╔═╡ d5fbb184-9eee-41db-8180-1cf5c5b42db7
plotlyjs()

# ╔═╡ 26ffc9f9-e1ec-4187-9eb6-4b0c3fc04223
begin
	x=[10,20,30,40,50]
	y=[19,30,30,15,10]
end

# ╔═╡ 1bb1f1e1-b366-4b31-873b-e4e78d1dca43
marsspec=mars_sky_opportunity()

# ╔═╡ 4f6b10bb-810d-49c4-ab9e-97c935e94504
marsspline=extrap(spline3(marsspec.λ,marsspec.l,:natural),380:830,:linear)

# ╔═╡ e2292549-9f32-427d-aa5b-1b816374caa1
marssky= interp(marsspline,380:1:830)

# ╔═╡ bc9a854d-4fc5-47da-8f09-1e21afe70522
begin
	scatter(marsspec.λ,marsspec.l,markersize=3,yrange=0.0:1.0)
	plot!(marssky[1],marssky[2])
end

# ╔═╡ ae0407ff-283a-4bf8-81ea-f0ac027aab8f
begin
	ms=LSpec(marssky[1],marssky[2])
	mr=RSpec(marssky[1],marssky[2])
end

# ╔═╡ fceaba8b-3c12-4883-8140-0930a75a590e
begin
	two=ms*mr
	three=ms*mr*mr*mr*mr*mr*mr*mr*mr
	plot(ms.λ,ms.l)
	plot!(two.λ,two.l)
	plot!(three.λ,three.l)
end

# ╔═╡ 208ed350-6c4e-40b2-9b3a-0a821b18a778

# ╔═╡ c1485749-eb3b-4f51-8814-416a39e7f5e0
begin
	lin=spline1(x,y)
	linint=interp(lin,10:.5:50)
	cub=spline3(x,y,:natural)
	cubint=interp(cub,10:.5:50)
	cubext=extrap(cub,0:60,:linear)
	extint=interp(cubext,0:.5:60)
end

# ╔═╡ b1299c43-c5a7-434e-8543-1715e23dc33f
begin
	scatter(x,y,color=:black,markersize=3,label="points")
	plot!(linint[1],linint[2],color=:grey,style=:dash,label="linear")
	plot!(cubint[1],cubint[2],color=:cyan,label="cubic nat.")
	plot!(extint[1],extint[2],color=:orange,linestyle=:dash,label="lin. extrap.")
end

# ╔═╡ f536734c-b7c0-4589-9be6-f4abea07faa8
function extrapo(Spl,xrange,extr::Symbol)
    envx=collect(xrange)
    xs=0.0 # start of extrapolatin range
    xe=0.0 # end of extrapolation range
    envx[1] < Spl.x[1] ? xs=envx[1] : xs=Spl.x[1]
    envx[1] < Spl.x[1] ? xs=envx[1] : xs=Spl.x[1]
    envx[end] > Spl.x[end] ? xe=envx[end] : xe=Spl.x[end]

    x=copy(Spl.x)
    a=copy(Spl.a)
	nf=nextfloat(Float64(x[end]))

    if typeof(Spl) in (LinearSpline, QuadraticSpline, CubicSpline)
        b=copy(Spl.b)
    end
    if typeof(Spl) in (QuadraticSpline, CubicSpline)
        c=copy(Spl.c)
    end
    if typeof(Spl) == CubicSpline
        d=copy(Spl.d)
    end
	if extr == :zero
		if xs < x[1]
			pushfirst!(x,xs)
			pushfirst!(a, 0.0)
			if typeof(Spl) in (LinearSpline, QuadraticSpline, CubicSpline)
				pushfirst!(b,0.0)
			end
			if typeof(Spl) in (QuadraticSpline,CubicSpline)
				pushfirst!(c,0.0)
			end
			if typeof(Spl) == CubicSpline
				pushfirst!(d,0.0)
			end
		end
		if xe ≥ x[end]
			if typeof(Spl) in (LinearSpline, QuadraticSpline, CubicSpline)
				x[end]=nf
				push!(x,xe)
				push!(a,0.0)
				push!(b,0.0)
			end
			if typeof(Spl) in (QuadraticSpline,CubicSpline)
				push!(c,0.0)
			end
			if typeof(Spl) == CubicSpline
				push!(d,0.0)
			end
		end
	elseif extr == :boundary
		de=x[end]-x[end-1]
		if xs < x[1]
			pushfirst!(x,xs)
			pushfirst!(a, a[1])
			if typeof(Spl) in (LinearSpline, QuadraticSpline, CubicSpline)
				pushfirst!(b,0.0)
			end
			if typeof(Spl) in (QuadraticSpline,CubicSpline)
				pushfirst!(c,0.0)
			end
			if typeof(Spl) == CubicSpline
				pushfirst!(d,0.0)
			end
		end
		if xe ≥ x[end]
			if typeof(Spl) == LinearSpline
				push!(a,Spl.b[end]*de+Spl.a[end])
				push!(b,0.0)
				x[end]=nf
			end
			if typeof(Spl) == QuadraticSpline
				push!(a,c[end]*de^2 + b[end]*de + a[end])
				push!(b,0.0)
				push!(c,0.0)
				x[end]=nf
			end
			if typeof(Spl) == CubicSpline
				push!(a,d[end]*de^3 + c[end]*de^2 + b[end]*de + a[end])
				push!(b,0.0)
				push!(c,0.0)
				push!(d,0.0)
				x[end]=nf
			end
		end
		push!(x,xe)
	elseif extr == :linear
		if xs < x[1]
			pushfirst!(x,xs)
			#
			pushfirst!(a,a[1]+b[1]*(x[1]-x[2]))
			if typeof(Spl) in (LinearSpline, QuadraticSpline, CubicSpline)
				pushfirst!(b,b[1])
			end
			if typeof(Spl) in (QuadraticSpline,CubicSpline)
				pushfirst!(c,0.0)
			end
			if typeof(Spl) == CubicSpline
				pushfirst!(d,0.0)
			end
		end
		if xe ≥ x[end]
			if typeof(Spl) == LinearSpline
				push!(a,a[end]+b[end]*(x[end]-x[end-1]))
				push!(b,b[end])
			end
			if typeof(Spl) == QuadraticSpline
				push!(c,0.0)
			end
			if typeof(Spl) == CubicSpline
				de=x[end]-x[end-1]
				push!(a,d[end]*de^3 + c[end]*de^2 + b[end]*de + a[end])
				push!(b,3*d[end]*de^2 + 2*c[end]*de + b[end])
				push!(c,0.0)
				push!(d,0.0)
				x[end]=nf
			end
		end
		push!(x,xe)
    elseif extr == :quadratic
        if xs < x[1]
            pushfirst!(x,xs)
            if typeof(Spl) in (LinearSpline, QuadraticSpline, CubicSpline)
                pushfirst!(b,x[1]*((a[3]-a[2])+x[2]*(a[1]-b[3])+x[3]*(a[2]-a[1]))/((x[1]-x[2])*(x[1]-x[3])*(x[2]-x[3])))
                pushfirst!(a,a[1]/(b[1]*(x[1]-xs)))
            end
            if typeof(Spl) in (QuadraticSpline, CubicSpline)
                pushfirst!(c,0.0)
            end
            if typeof(Spl) == CubicSpline
                pushfirst!(d,0.0)
            end
        end
        if xe ≥ x[end]
            push!(x,xe)
            if typeof(Spl) in (LinearSpline, QuadraticSpline, CubicSpline)
                push!(b,x[end-2]*((a[end]-a[end-1])+x[end-1]*(a[end-2]-b[end])+x[end]*(a[end-1]-a[end-2]))/((x[end-2]-x[end-1])*(x[end-2]-x[end])*(x[end-1]-x[end])))
                push!(a,a[end]/(b[end]*(x[end]-xs)))
            end
            if typeof(Spl) in (QuadraticSpline, CubicSpline)
                push!(c,0.0)
            end
            if typeof(Spl) == CubicSpline
                push!(d,0.0)
            end
        end
    else error("Extrapolation type $extr does not exist.")
    end
    if typeof(Spl) == NearestNeighborSpline
        NearestNeighborSpline(a,x)
    elseif typeof(Spl) == LinearSpline
        LinearSpline(b,a,x)
    elseif typeof(Spl) == QuadraticSpline
        QuadraticSpline(c,b,a,x)
    else
        CubicSpline(d,c,b,a,x)
    end
end

# ╔═╡ Cell order:
# ╠═44c12b60-97e4-11eb-23f7-eb2bf0847809
# ╠═d5fbb184-9eee-41db-8180-1cf5c5b42db7
# ╠═26ffc9f9-e1ec-4187-9eb6-4b0c3fc04223
# ╠═1bb1f1e1-b366-4b31-873b-e4e78d1dca43
# ╠═4f6b10bb-810d-49c4-ab9e-97c935e94504
# ╠═e2292549-9f32-427d-aa5b-1b816374caa1
# ╠═bc9a854d-4fc5-47da-8f09-1e21afe70522
# ╠═ae0407ff-283a-4bf8-81ea-f0ac027aab8f
# ╠═fceaba8b-3c12-4883-8140-0930a75a590e
# ╠═208ed350-6c4e-40b2-9b3a-0a821b18a778
# ╠═c1485749-eb3b-4f51-8814-416a39e7f5e0
# ╠═b1299c43-c5a7-434e-8543-1715e23dc33f
# ╠═f536734c-b7c0-4589-9be6-f4abea07faa8


text:
consensus(u, v) = (isequal(u.x,v.x) || v.x == nothing ? u.x : (u.x == nothing ? v.x : nothing))

function Base.:+(u::Poly, v::Poly)
    x = consensus(u, v)
    x === nothing && return poly(sym(u) + sym(v), nothing)
    q = Poly(u.T, x)

    for k = 0:max(u.n, v.n)
        q[k] = u[k] + v[k]
    end
    q
end

Base.:+(u, v::Poly) = polyas(v, u) + v
Base.:+(u::Poly, v) = u + polyas(u, v)

function Base.:-(u::Poly, v::Poly)
    x = consensus(u, v)
    x === nothing && return poly(sym(u) - sym(v), nothing)
    q = Poly(u.T, x)

    for k = 0:max(u.n, v.n)
        q[k] = u[k] - v[k]
    end
    q
end

Base.:-(u::Poly, v) = u - polyas(u, v)
Base.:-(u, v::Poly) = polyas(v, u) - v

function Base.:*(u::Poly, v::Poly)
    x = consensus(u, v)
    x === nothing && return poly(sym(u) * sym(v), nothing)
    q = Poly(u.T, x)

    for i = 0:u.n
        for j = 0:v.n
            q[i+j] += u[i] * v[j]
        end
    end
    q
end

Base.:*(u::Poly, v) = u * polyas(u, v)
Base.:*(u, v::Poly) = polyas(v, u) * v

function divide_poly(u::Poly, v::Poly; op=/)
    x = consensus(u, v)
    x === nothing && return poly(sym(u) / sym(v), nothing)
    q = Poly(u.T, x)
    r = copy(u)
    m, n = u.n, v.n

    for i = m-n:-1:0
        q[i] = op(r[i+n], v[n])
        for j = n+i:-1:i
            r[j] = r[j] - q[i] * v[j-i]
        end
    end

    update_order!(r)

    q, r
end

function Base.:/(u::Poly, v::Poly)
    q, r = divide_poly(u, v)
    q
end

Base.:/(u::Poly, v) = u / polyas(u, v)
Base.:/(u, v::Poly) = polyas(v, u) / v

function Base.:%(u::Poly, v::Poly)
    q, r = divide_poly(u, v)
    r
end

Base.:%(u::Poly, v) = u % polyas(u, v)
Base.:%(u, v::Poly) = polyas(v, u) % v

function pseudo_divide_poly(u::Poly, v::Poly)
    m, n = u.n, v.n
    α = v[n]^(m-n+1)
    q, r = divide_poly(u*α, v; op=÷)
    return q, r, α
end

Base.:^(p::Poly, k) = poly(sym(p)^k, p.x)

##############################################################################

function gcd_naive(u::Poly, v::Poly)
    while degree(v) > 0
        println(u, '\t', v)
        u, v = v, u % v
    end
    return iszero(v) ? u : 1
end

function gcd_monic(u::Poly, v::Poly)
    if degree(u) == 0 || degree(v) == 0
        return gcd(cont(u), cont(v))
    end
    u = u / leading(u)
    while degree(v) > 0
        v = v / leading(v)
        u, v = v, u % v
    end
    return iszero(v) ? u : 1
end

function gcd_extended(u::Poly, v::Poly; op=/)
    if degree(u) == 0 || degree(v) == 0
        return gcdx(cont(u), cont(v))
    end

    sᵤ = one(u.T)
    tᵤ = zero(u.T)
    sᵥ = zero(v.T)
    tᵥ = one(v.T)

    while degree(v) > 0
        q, r = divide_poly(u, v; op)
        s, t = sᵤ - q*sᵥ, tᵤ - q*tᵥ
        u, sᵤ, tᵤ = v, sᵥ, tᵥ
        v, sᵥ, tᵥ = r, s, t
    end

    l = leading(u)
    if iszero(l)
        return u, sᵤ, tᵤ
    else
        return u / l, sᵤ / l, tᵤ / l
    end
end

function Base.gcd(u::Poly, v::Poly)
    if u.T <: Rational && u.T <: Rational
        return gcd_monic(u, v)
    else
        return gcd_naive(u, v)
    end
end

Base.gcd(u::Poly, v) = gcd(u, polyas(u, v))
Base.gcd(u, v::Poly) = gcd(polyas(v, u), v)

Base.gcdx(u::Poly, v::Poly) = gcd_extended(rationalize(u), rationalize(v))
Base.gcdx(u::Poly, v) = gcdx(u, polyas(u, v))
Base.gcdx(u, v::Poly) = gcdx(polyas(v, u), v)

#############################################################################

Symbolics.derivative(p::Poly, x) = isequal(p.x, x) ? derivative(p) : derivative(sym(p), x)

function Symbolics.derivative(p::Poly)
    q = Poly(p.T, p.x)
    for i = 1:p.n
        q[i-1] = i * p[i]
    end
    q
end

function integrate(p::Poly)
    q = Poly(p.T, p.x)
    for i = 0:p.n
        q[i+1] = p[i] / (i+1)
    end
    q
end

###############################################################################

function rationalize(p::Poly)
    p.T <: Rational && return p
    try
        q = Poly(Rational{BigInt}, p.x, p.n, Rational{BigInt}(p.a₀), Dict())
        for k = 1:p.n
            q[k] = Rational{BigInt}(p[k])
        end
        return q
    catch
        return nothing
    end
end

function to_monic(p::Poly)
    p = rationalize(p)
    c = leading(p)
    n = degree(p)
    q = Poly(p.T, p.x)
    for i = 0:n
        q[i] = p[i] * c ^ (n-1-i)
    end
    c, q
end

function from_monic(p::Poly, c)
    q = Poly(p.T, p.x)
    n = degree(p)
    for i = 0:n
        q[i] = p[i] / c ^ (n-1-i)
    end
    q
end

function integer_poly(p::Poly)
    p = rationalize(p)
    l = lcm(denominator.(extract_coef(p))...)
    return 1//l, p*l
end


text:
using Distributed
using Observables

@enum Status st_waiting st_working st_done st_error

status_str(st::Status) = (" ", "→", "✓", "✗")[Int(st)+1]
indent(str, n) = join((" "^n * line for line in split(str,r"(\r)?\n")), "\n")
default_summarize( els, statuses, messages, ids, indentlevel=8 ) = sprint() do io
    for (i,(el,st,msg,id)) in enumerate(zip(els,statuses,messages,ids))
        print(io, status_str(st), " ", join(el,","))
        isnothing(id) || print(io, "\t(on worker $id)")
        println(io)
        println(io, indent(msg, indentlevel))
    end
end

"""
    summary, task = opmap(f, [::AbstractWorkerPool], c...;
                          schedule_now=true, summarize=default_summarize,
                          kwargs... )

Transform collection(s) `c` by applying `f` to each element using available
workers and tasks, by calling `Distributed.pmap`. The function `f` receives as
its first argument a callable `setmessage`, and may call `setmessaage(msg)` to
set a status message describing the worker's current status.
The returned `summary` is an `Observable` containing a summary string of all
workers' statuses including their status messages.
The summary may be customized by setting the `summarize` keyword argument.

The application of `f` across available workers is performed asynchronously in
the returned `task::Task`. Setting `schedule_now=false` prevents `opmap` from
scheduling the task.

The `ObservablePmap` module must be loaded on all workers before calling
`opmap`, e.g. by `@everywhere using ObservablePmap`.

Other positional and keyword arguments to `opmap` are the same as for `pmap`.

Example:
```jldoctest
@everywhere using ObservablePmap
using Observables # for `map` method below
summ, task = opmap(1:5; on_error=identity) do setmessage, x
    for i=1:5
        setmessage("\$i @ \$x")
        (10i+x==33) && error("Error at \$i")
        sleep(rand())
    end
end
html = map(x -> HTML("<pre>\$x</pre>"), summ)
```
"""
opmap( f, c...; kwargs... ) = opmap( f, default_worker_pool(), c...; kwargs... )
function opmap( f, p::AbstractWorkerPool, c...;
               schedule_now=true, summarize=default_summarize, kwargs... )
    els = collect(zip(c...))
    workers_ch = RemoteChannel(()->Channel(Inf))

    pmap_task = Task(()->pmap(p, eachindex(els); kwargs...) do i
        setmessage(msg) = put!( workers_ch, (myid(), i, msg) )
        setmessage(st_working)
        try
            f( setmessage, els[i]... )
        catch e
            setmessage(sprint(showerror,e,catch_backtrace()))
            setmessage(st_error)
            rethrow()
        end
        setmessage(st_done)
    end)

    messages = Any["" for el in els]
    worker_ids = Union{Int,Nothing}[nothing for el in els]
    statuses = [st_waiting for el in els]
    active = Set(eachindex(els))
    summ = Observable(summarize(els,statuses,messages,worker_ids))

    summary_task = Task(()->begin
        while !isempty(active)
            id, i, msg = take!(workers_ch)
            worker_ids[i] = id
            if msg isa Status
                statuses[i] = msg
                if msg ∈ (st_done, st_error)
                    pop!(active, i)
                end
            else
                messages[i] = msg
            end
            # `invokelatest` as workaround to world-age issues (`Observables` issue #50)
            Base.invokelatest(setindex!, summ, summarize(els,statuses,messages,worker_ids))
            #summ[] = summarize(els,statuses,messages,worker_ids)
        end
    end)

    schedule_now && schedule(pmap_task)
    schedule(summary_task)

    summ, pmap_task, summary_task
end


text:
#=
RandomLessonsLearned.jl
Jacob Child
April 23, 2022
=#

#quick test function to understand iterating
function test_iterator(x)
    for i in x
        if i <= .3
            println(i+7)
        elseif i > .3
            println(i)
        end

    end

end

# returning more than one thing is
# return a , b , c etc... gives a tuple
zbar, a = zbar_calculations(c, p, x)

# ***look into broadcasting, it looks like you can iterate through
#an array using a broadcast function without every term in the equation having a dot operator?


text:
@testset "Kolmogorov-Smirnov test" begin
    x = [
        392,585,410,414,384,423,456,549,455,406,470,402,416,431,435,413,468,435,
        432,438,425,392,410,498,475,431,460,440,420,422,409,417,370,500,501,478,
        454,444,447,393,446,454,500,433,444,427,412,433,460,434
    ]
    res_kstest = ks_test(x, Normal(mean(x), std(x)))
    @test res_kstest.statistic ≈ 0.119089045870544
    @test res_kstest.n_obs == 50
    @test res_kstest.pval ≈ 0.477419559223632
end


text:
srand(1)

# --- generate a random (n,ne,np) network.
# --- d = matrix diagonal
# --- ne = number of edges
# --- np = number of positive edges
function randomnkp (d,ne,np)
    n = length(d)                       # int matrix diagonal - never change
    nn = ne-np                          # negative edges
    nz = n*(n-1)-ne                     # number of zero-edges
    intlist = [ones(Int64,np);-ones(Int64,nn);zeros(Int64,nz)]
    rect = reshape(shuffle(intlist),n,n-1)
    rrect = [rect zeros(Int64,n)]       # now insert back the diagonal
    lrect = [zeros(Int64,n) rect]
    interrand = tril(rrect - diagm(diag(rrect))) + triu(lrect-diagm(diag(lrect)))
    return interrand
end

# --- return a (n,ne,np) preserved shuffled version of a network
function shufflenet (intr,slf)
    n = length(slf)
    nw = similar(intr)                  # this will be the randomized version
    shf = shuffle(slf)                  # keep number of self degredating genes
    tl = tril(intr)                     # remove diag and join upper/lower triang
    tu = triu(intr)
    rect = tu[:,2:end]+tl[:,1:end-1]    # shuffle this
    rect = reshape(shuffle(reshape(rect,length(rect),1)[:,1]),n,n-1)
    rrect = [rect zeros(Int64,n)]       # now insert back the shuffled diagonal
    lrect = [zeros(Int64,n) rect]
    nw = tril(rrect - diagm(diag(rrect))) + triu(lrect-diagm(diag(lrect)))
    return (nw,shf)
end

# --- generate the symmetric adjacency matrix with elements 0,1 ---
function symadj(m)
    n = length(m[:,1])
    ad = zeros(Int64,n,n)
    for (i=1:n) for (j=i+1:n)       # adjacency matrix - without diagonals
        ad[i,j] = m[i,j]==0 && m[j,i]==0 ? 0 : 1
        ad[j,i] = ad[i,j]
    end end
    return ad
end

## --- check connectivity -----------
function isconnected (adj)
    n = length(adj[:,1])                # number of nodes
    nsq = n*n
    dum = adj
    (sumprev,sumdum) = (0,0)
    for (i=1:n)
        dum = sign(dum + dum*adj)
        sumdum = sum(dum)
        if (sumdum == sumprev)
            break
        end
        sumprev = sumdum
    end
    return (sumdum==nsq)
end

## --- node update function (integrator): yeast's majority rule -----------
function upd(s,intr,slf,n)
    sss = similar(s)
    inp = intr*s
    for (i=1:n)
        sss[i] = (inp[i]>0)?1:(inp[i]<0)?0:(slf[i]==1)?0:s[i]
    end
    return sss
end

## ---- output is the attractor list + basin sizes ---
## ---- calculated from a subset of initial conditions ---
function getcyclessample(intr,slf,ninit)

    n = length(slf)
    attrlist = zeros(Int64,n)           # s=0 is always an attractor
    keyzeros = join(zeros(Int64,n))
    attdict = [keyzeros=>0]    # values are the basin sizes
    for (i=1:ninit)       # calculate basins using ninit initial states
        s = rand(0:1,n)
        skey = join(s)
        pathdict = Dict()
        while (!haskey(attdict,skey) && !haskey(pathdict,skey))
            pathdict[skey] = "-"
            s = upd(s,intr,slf,n)
            skey=join(s)
        end
        if (haskey(attdict,skey)) attdict[skey] += 1
        else                        # this is a new cycle or fp
            attdict[skey] = 1
            attrlist = [attrlist s]
        end
    end
    basinlist = [attdict[join(attrlist[:,i])]/ninit for i in 1:length(attrlist[1,:])]
    return attrlist,basinlist
end

## ---- find number of coop nodes in a state ---
function ncoop (s,intr)
    n = length(s)
    act = intr*diagm(s) # keep only the columns for active nodes
    nc=0
    for (i=1:n)
#        if (abs(sum(act[i,:])) == sum(abs(act[i,:])))
        if (!in(1,act[i,:]) || !in(-1,act[i,:]))
            nc += 1
        end
    end
    return nc
end

## ---- find coop coef for sampled cycles of a network ---
function coopattr (intr,slf,ninit)
    n = length(slf)
    (listatt,listb) = getcyclessample(intr,slf,ninit)
    listperiod=Int64[];
    listcop=Float64[];
    listact=Float64[];
    npointattr = 0;
    for (i=1:length(listatt[1,:]))        # go over each attractor
        s = listatt[:,i]
#        if (sum(s) != 0)
        ncp = ncoop(s,intr)
#            nac = sum(abs(intr)*s)
        nac = sum(s)
        period = 1
        s=upd(s,intr,slf,n)
        while (s != listatt[:,i])
            ncp += ncoop(s,intr)
#            nac += sum(abs(intr)*s)
            nac += sum(s)
            period += 1
            s=upd(s,intr,slf,n)
        end
        push!(listperiod,period)
        push!(listcop,ncp/period)
        push!(listact,nac/period)
        if (period==1) npointattr += 1 end
    #        end
    end
#    if (length(listperiod)==1)  return -1,-1,-1,-1 # no fixed point other than (0,..0)
#    else                    # also returns the fracion of active nodes
#        return (1.0*mean(listcop)/n,mean(listact)/n,length(listc[1,:]),npointattr)
    return listatt,listb,listperiod,listcop,listact
#    end
end

# --- calculate average cooperativity over random networks
n = 11
ntotal = 100000;            # ensemble size
numinit = 2^11              # num of initial cond for fp determination

for ne=29:29
    for np=15:15
        dirname = "$n.$ne.$np"
        if (!isdir(dirname))
            run(`mkdir $dirname`)
        end
        nself = 5
        selfsh = [ones(Int64,nself),zeros(Int64,n-nself)]
        listc=zeros(Float64,ntotal)        # list of coherence values
        lista=zeros(Float64,ntotal)        # list of active node fractions
        listl=zeros(Float64,ntotal)        # list of fraction of point attractors
        intersh = randomnkp(selfsh,ne,np) # generate one sample with given (nkp)
        i=0
        while (i<ntotal)
            (intersh,selfsh) = shufflenet(intersh,selfsh) # shuffle
            if (isconnected(symadj(intersh))) # check if connected
#                (ncop,actratio,natt,npointattr) = coopattr(intersh,selfsh,numinit)
                (attlist,basinlist,listperiod,listcop,listact) = coopattr(intersh,selfsh,numinit)
                if (length(attlist) > 1)
                    i += 1              # increase network count
                    edges = zeros(Int64,ne+sum(selfsh),3)
                    idx=1
                    for j=1:n
                        if selfsh[j] == 1
                            edges[idx,:] = [j j -1] # self degrading genes
                            idx += 1
                        end
                        for k=1:n
                            if intersh[j,k] != 0
                                edges[idx,:] = [j k intersh[j,k]]
                                idx += 1
                            end end end
                    writedlm("$dirname/$i.dat",edges,' ')
                    writedlm("$dirname/$i.att",[attlist' basinlist listperiod listcop/n listact],' ')
                    if (mean(listcop/n)==1)
                        print("$i\n")
                    end
                end
            end
        end
    end
end

max_stars_repo_name: akabakcioglu/network-motifs

max_issues_repo_name: akabakcioglu/network-motifs

max_forks_repo_name: akabakcioglu/network-motifs


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to compute the Celestial Intermediate Origin (CIO).
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
#   [2] Vallado, D. A (06-Feb-2018). Consolidated Errata of Fundamentals of
#       Astrodynamics and Applications 4th Ed.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export cio_iau2006

"""
    cio_iau2006(JD_TT::Number)

Compute the coordinates `X` and `Y` of the Celestial Intermediate Pole (CIP)
with respect to the Geocentric Celestial Reference Frame (GCRF), and the CIO
locator `s`. The algorithm is based on the IAU-2006 theory.

The CIO locator `s` provides the position of the CIO on the Equator of the CIP
corresponding to the kinematical definition of the non-rotation origin in the
GCRS when the CIP is moving with respect to the GCRS between the reference epoch
and the epoch due to precession and nutation **[1]**(p. 214).

# Returns

* The coordinate `X` of the CIP w.r.t. the GCRF.
* The coordinate `Y` of the CIP w.r.t. the GCRF.
* The CIO locator `s`.

# References

- **[1]**: Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
    Microcosm Press, Hawthorn, CA, USA.
"""
function cio_iau2006(JD_TT::Number)
    # Compute the Julian Centuries from `JD_TT`.
    T_TT = (JD_TT - JD_J2000) / 36525

    # Auxiliary variables
    # ===================

    a2d = 1 / 3600
    d2r = π / 180
    a2r = a2d * d2r

    # Fundamental arguments
    # =====================

    # Luni-solar part.
    M_s, M_m, u_Mm, D_s, Ω_m = luni_solar_args_iau2006(JD_TT)

    # Planetary part.
    λ_M☿, λ_M♀, λ_Me, λ_M♂, λ_M♃, λ_M♄, λ_M⛢, λ_M♆, p_λ = planetary_args_iau2006(JD_TT)

    # X position of the CIP
    # ==========================================================================

    ΔX = _iau2006_sum(
        (
            coefs_iau2006_CIP_x0,
            coefs_iau2006_CIP_x1,
            coefs_iau2006_CIP_x2,
            coefs_iau2006_CIP_x3,
            coefs_iau2006_CIP_x4,
        ),
        T_TT,
        M_s,
        M_m,
        u_Mm,
        D_s,
        Ω_m,
        λ_M☿,
        λ_M♀,
        λ_Me,
        λ_M♂,
        λ_M♃,
        λ_M♄,
        λ_M⛢,
        λ_M♆,
        p_λ
    )

    X = @evalpoly(
        T_TT,
        -0.016_617,
        +2004.191_898,
        -0.429_782_9,
        -0.198_618_34,
        +0.000_007_578,
        +0.000_005_928_5
    ) + ΔX

    # Convert to [rad].
    X *= a2r

    # Y position of the CIP
    # ==========================================================================

    ΔY = _iau2006_sum(
        (
            coefs_iau2006_CIP_y0,
            coefs_iau2006_CIP_y1,
            coefs_iau2006_CIP_y2,
            coefs_iau2006_CIP_y3,
            coefs_iau2006_CIP_y4
        ),
        T_TT,
        M_s,
        M_m,
        u_Mm,
        D_s,
        Ω_m,
        λ_M☿,
        λ_M♀,
        λ_Me,
        λ_M♂,
        λ_M♃,
        λ_M♄,
        λ_M⛢,
        λ_M♆,
        p_λ
    )

    Y = @evalpoly(
        T_TT,
        -0.006_951,
        -0.025_896,
        -22.407_274_7,
        +0.001_900_59,
        +0.001_112_526,
        +0.000_000_135_8
    ) + ΔY

    # Convert to [rad].
    Y *= a2r

    # Parameter `s` (CIO locator)
    # ==========================================================================
    #
    # The value `s` provides the position of the CIO on the Equator of the CIP
    # corresponding to the kinematical definition of the non-rotation origin in
    # the GCRS when the CIP is moving with respect to the GCRS between the
    # reference epoch and the epoch due to precession and nutation [1, p. 214].

    Δs = _iau2006_sum(
        (
            coefs_iau2006_CIO_s0,
            coefs_iau2006_CIO_s1,
            coefs_iau2006_CIO_s2,
            coefs_iau2006_CIO_s3,
            coefs_iau2006_CIO_s4
        ),
        T_TT,
        M_s,
        M_m,
        u_Mm,
        D_s,
        Ω_m,
        λ_M☿,
        λ_M♀,
        λ_Me,
        λ_M♂,
        λ_M♃,
        λ_M♄,
        λ_M⛢,
        λ_M♆,
        p_λ
    )

    s = @evalpoly(
        T_TT,
        # We must convert this term to [arcsec] to match the units.
        #    ||
        # |------|
        (-X * Y / 2) / a2r + 0.000_094,
        +0.003_808_65,
        -0.000_122_68,
        -0.072_574_11,
        +0.000_027_98,
        +0.000_015_65
    ) + Δs

    # Convert to [rad].
    s *= a2r

    return X, Y, s
end


text:
using Test
using XPORTA: IEQ

@testset "./src/game_conversions.jl" begin

using BellScenario

@testset "Conversions to BellGames" begin
    @testset  "Vector{Int64} -> BellGame, LocalSignaling" begin
        norm_facet = [2,1,0,0,1,2,-2,-1,-2,2]
        gen_facet = [2,1,0,0,0,1,2,0,0,1,0,2,4]
        scenario = LocalSignaling(3,4,2)

        bell_game = convert(BellGame, norm_facet, scenario, rep="normalized")
        bell_game2 = convert(BellGame, gen_facet, scenario, rep="generalized")

        @test bell_game isa BellGame
        @test bell_game2 isa BellGame
        @test bell_game == [2 0 0;1 1 1;0 2 0;0 0 2]
        @test bell_game.β == 4
        @test bell_game2 == [2 0 0;1 1 1;0 2 0;0 0 2]
        @test bell_game2.β == 4

        @test_throws DomainError convert(BellGame, norm_facet, scenario, rep="non-signaling")
    end

    @testset "IEQ -> BellGames, BlackBox" begin
        gen_ieq = IEQ(inequalities = [1 1 0 0 0 1 2;1 1 0 0 0 0 1] )
        scenario = BlackBox(3,2)

        bell_games = convert(Vector{BellGame}, gen_ieq, scenario, rep="generalized")

        @test bell_games isa Vector{BellGame}
        @test bell_games == BellGame.([[1 0;1 0;0 1],[1 0;1 0;0 0]],[2,1])

        norm_ieq = IEQ(inequalities = [1 1 -1 -1 1;1 1 0 0 1])

        bell_games = convert(Vector{BellGame}, norm_ieq, scenario, rep="normalized")
        bell_games2 = convert(Vector{BellGame}, norm_ieq, scenario)

        @test bell_games isa Vector{BellGame}
        @test bell_games2 isa Vector{BellGame}
        @test bell_games == BellGame.([[1 0;1 0;0 1],[1 0;1 0;0 0]],[2,1])
        @test bell_games2 == BellGame.([[1 0;1 0;0 1],[1 0;1 0;0 0]],[2,1])
    end

    @testset "Vector{Int64} -> BellGame, BipartiteNonSignaling" begin
        scenario = BipartiteNonSignaling(2,2,2,2)

        ns_facet = [1, 0, 1, 0, -1, -1, -1, 1, 1]
        bell_game = convert(BellGame, ns_facet, scenario)

        @test bell_game isa BellGame
        @test bell_game == [1 0 0 1;1 1 1 0;1 1 1 0;0 1 1 0]
        @test bell_game.β == 3

        norm_facet = [0, 0, 1, 0, 1, -1, 0, 0, 0, 1, 0, 1, 1]
        bell_game = convert(BellGame, norm_facet, scenario, rep="normalized")

        @test bell_game isa BellGame
        @test bell_game.β == 2
        @test bell_game == [0 1 0 1;0 2 0 0;1 0 0 1;0 1 0 0]
    end
end

@testset "Conversions from BellGames" begin

    @testset "BellGame -> Vector{Int64}" begin
        BG = BellGame([2 0 0;1 1 1;0 2 0;0 0 2], 4)

        gen_facet = convert(Vector{Int64}, BG, rep="generalized")

        @test gen_facet isa Vector{Int64}
        @test gen_facet == [2,1,0,0,0,1,2,0,0,1,0,2,4]

        norm_facet = convert(Vector{Int64}, BG, rep="normalized")

        @test norm_facet isa Vector{Int64}
        @test norm_facet == [2,1,0,0,1,2,-2,-1,-2,2]

        @test_throws DomainError convert(Vector{Int64}, BG, rep="non-signaling")
    end

    @testset "BellGames -> IEQ" begin
        bell_games = BellGame.([[1 0;0 1],[1 0;0 0]],[2,1])
        gen_ieq = convert(IEQ, bell_games, rep="generalized")

        @test gen_ieq isa IEQ
        @test gen_ieq.inequalities == [1 0 0 1 2;1 0 0 0 1]

        norm_ieq = convert(IEQ, bell_games, rep="normalized")

        @test norm_ieq isa IEQ
        @test norm_ieq.inequalities == [1 -1 1;1 0 1]

        norm_ieq2 = convert(IEQ, bell_games)

        @test norm_ieq2 isa IEQ
        @test norm_ieq2.inequalities == [1 -1 1;1 0 1]
    end

    @testset "BellGame -> Facet BipartiteNonSignaling" begin
        @testset "chsh spot check" begin
            BG = BellGame([1 0 0 1;1 1 1 0;1 1 1 0;0 1 1 0], 3)
            scenario = BipartiteNonSignaling(2,2,2,2)

            facet = convert(Vector{Int64}, BG, scenario)

            @test facet isa Vector{Int64}
            @test length(facet) == LocalPolytope.vertex_dims(scenario, "non-signaling") + 1
            @test facet == [1,0,1,0,-1,-1,-1,1,1]
        end
    end
end

@testset "identity conversions" begin
    @testset "BellGames -> IEQ -> BellGames" begin
        bell_games = BellGame.([[2 0 0 0;0 2 0 0;0 0 2 0;1 1 1 0],[1 0 0 0;1 0 0 0;1 0 0 0;0 0 0 0]],[4,1])
        scenario = LocalSignaling(4,4,2)

        ieq = convert(IEQ, bell_games)
        converted_bell_games = convert(Vector{BellGame}, ieq, scenario)

        @test converted_bell_games == bell_games
    end

    @testset "BipartiteNonSignaling facet -> BellGame -> facet" begin
        @testset "simple scenarios full check" for scenario_params in [
            (2,2,2,2),(3,2,2,2),(2,3,2,2),(2,2,3,2),(2,2,2,3),
            (3,3,2,2),(3,2,3,2),(3,2,2,3),(2,3,3,2),(2,3,2,3),
            (2,2,3,3)
        ]
            scenario = BipartiteNonSignaling(scenario_params...)

            vertices = LocalPolytope.vertices(scenario)
            facet_dict = LocalPolytope.facets(vertices)

            for facet in facet_dict["facets"]
                BG = convert(BellGame, facet, scenario)

                @test BG isa BellGame

                facet2 = convert(Vector{Int64}, BG, scenario)

                @test facet2 isa Vector{Int64}
                @test facet2 == facet
            end
        end
    end
end

end


text:
"""
Update the toeplitz matrix A using vector x
"""
function toeplitz!(A, x)
   n = length(x)
   for i in 1:n
	   for j in i:n
	       A[i,j] = x[j-i+1]
	       A[j,i] = x[j-i+1]
       end
   end
   return A
end


text:
module CartesianCutCell

using LinearAlgebra
using SparseArrays
using FillArrays
using LazyArrays

export dirichlet, periodic
export mac, stencil

export mesh
export restrict

export laplacian
export forwarddiff, backwarddiff
export forwardinterp, backwardinterp
export integrate

include("boundary.jl")
include("stencil.jl")
include("mesh.jl")
include("differentiation.jl")
include("interpolation.jl")
include("restrict.jl")
include("capacity.jl")

end