text:
#=
Base julia sizehint! + push! garbage compared to
resize! + setindex!

Remove when sizehinting in base julia isn't garbage
=#

mutable struct NestedPushVector{T, K<:PushVector{T}, V<:AbstractVector{K}} <: AbstractVector{K}
    "Vector used for storage."
    parent::V
    "Number of elements held by `parent`."
    len::Int
    "Inner sizehint"
    inner_sh::Int
    "Outer sizehint"
    outer_sh::Int
end

# Simplifications
const PV{T} = PushVector{T, Vector{T}} where {T}
const NPV{T} = NestedPushVector{T, PV{T}, Vector{PV{T}}} where {T}

function NestedPushVector{T}(inner_sh::Integer = 4, outer_sh::Integer = 4) where {T}
    inner_sh ≥ 0 && outer_sh ≥ 0 || throw(DomainError((inner_sh, outer_sh), "Invalid initial size."))
    NestedPushVector(
        [PushVector{T}(inner_sh) for _ in 1:outer_sh],
        0,
        inner_sh,
        outer_sh
    )
end

@inline Base.length(v::NestedPushVector) = v.len
@inline Base.size(v::NestedPushVector) = (v.len, )

function extend!(v::NestedPushVector{T}, n::Int) where {T}
    ish = v.inner_sh
    prev_l = length(v.parent)
    resize!(v.parent, prev_l + n)
    for i in 1:n
        v.parent[prev_l + i] = PushVector{T}(ish)
    end
    nothing
end

@inline function Base.getindex(v::NestedPushVector, i)
    @boundscheck checkbounds(v, i)
    @inbounds v.parent[i]
end

@inline function Base.setindex!(v::NestedPushVector, x::PushVector, i)
    @boundscheck checkbounds(v, i)
    @inbounds v.parent[i] = x
end

function freenext!(v::NestedPushVector)
    v.len += 1
    if v.len > length(v.parent)
        extend!(v, v.len)
    end
    nothing
end

function Base.push!(v::NestedPushVector, x::PushVector)
    v.len += 1
    if v.len > length(v.parent)
        extend!(v, v.len)
    end
    v.parent[v.len] = x
    v
end

function Base.empty!(v::NestedPushVector)
    resize!(v.parent, v.outer_sh)
    v.len = 0
    for pv in v.parent
        resize!(pv.parent, v.inner_sh)
        pv.len = 0
    end
    nothing
end


text:
#Demo 2: Chaining ABCD matrices
#-------------------------------------------------------------------------------

using CMDimCircuits
CMDimCircuits.@using_CData()

#Get a demo display:
include(CMDimCircuits.demoplotcfgscript)
#Normally use something like:
#CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
LBLAX_MAGDB = "Magnitude (dB)"
LBLAX_FREQ = "Frequency (Hz)"
color1 = cons(:a, line = set(color=:red, width=2))
color2 = cons(:a, line = set(color=:blue, width=2))

#==Input data
===============================================================================#
_logspace(start, stop, n) = 10 .^ range(start, stop=stop, length=n)
C = capacitance(2e-12)
R = 20.0
#f = (1:10)*1e9
f = _logspace(log10(1e3), log10(20e9), 100)
xrange = cons(:a, xyaxes = set(xmin=1e6, xmax=100e9))

#==Computations
===============================================================================#
#Use DataF1 (function of 1 argument) to augment expressivity:
f = DataF1(f, f)
_1 = ones(f)

#Compute ABCD matrix of series R:
sng = true #Compute using single ABCD matrix?
if sng
	#Compute single ABCD matrix:
	@show TR = series(:ABCD, impedance(R))
else
	#Compute ABCD{DataF1}:
	TR = series(:ABCD, impedance(R*_1))
end

#Compute ABCD matrix of shunt capacitors:
TC1 = shunt(:ABCD, admittance(C, f=f))
#TC2 = shunt(:ABCD, impedance(2C, f=f)) #Add asymmetry #TODO: fix multiplication of C
TC2 = shunt(:ABCD, impedance(capacitance(2*C.v), f=f)) #Add asymmetry

#Compute Pi network:
T = TC1*TR*TC2

S = Network(:S, T)
(s11, s12, s21, s22) = mx2elem(S)

#Intermediate network parameters (verify conversion routines):
nplist = [:Y, :Z, :ABCD, :H, :G, :T]

#==Generate plot
===============================================================================#
is11 = 1; is22 = 2; i𝛵 = 3; #Strip indices

plot = cons(:plot, nstrips=3,
	ystrip1 = set(axislabel=LBLAX_MAGDB, striplabel="Input Reflection (s11)"),
	ystrip2 = set(axislabel=LBLAX_MAGDB, striplabel="Output Reflection (s22)"),
	ystrip3 = set(axislabel=LBLAX_MAGDB, striplabel="Transmission Coefficients (𝛵)"),
	xaxis = set(scale=:log, label=LBLAX_FREQ),
)
push!(plot,
	cons(:wfrm, dB20(s11), color1, label="s11", strip=is11),
	cons(:wfrm, dB20(s22), color1, label="s22", strip=is22),
	cons(:wfrm, dB20(s12), color1, label="s12", strip=i𝛵),
	cons(:wfrm, dB20(s21), color2, label="s21", strip=i𝛵),
)

Sref = S

let S, s11, s12, s21, s22 #HIDEWARN_0.7
for np in nplist
	X = Network(np, Sref)
	S = Network(:S, X)
	(s11, s12, s21, s22) = mx2elem(S)
	push!(plot,
		cons(:wfrm, dB20(s11), color1, label="s11 ($np)", strip=is11),
		cons(:wfrm, dB20(s22), color1, label="s22 ($np)", strip=is22),
		cons(:wfrm, dB20(s12), color1, label="s12 ($np)", strip=i𝛵),
		cons(:wfrm, dB20(s21), color2, label="s21 ($np)", strip=i𝛵),
	)
end
end

pcoll = push!(cons(:plotcoll, title="Pi Network Cascade + Param Conversion"), plot)
	pcoll.ncolumns = 1

#==Display results in pcoll
===============================================================================#
EasyPlot.displaygui(pcoll)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll

ext: jl


text:
using Colors

include("helpers/page.jl")
include("helpers/doc.jl")

intro = md"""
$(title(2, "Behavior"))
$(vskip(1em))

# Functions
$(vskip(1em))

$(
   showdocs([hasstate, keypress, clickable, selectable, send, recv, wire])
)

$(vskip(4em))

"""

function main(window)
    docpage(intro)
end


text:
using Test
using RenderingGeometry

include("definitions.jl")
include("../test_logging.jl")
include("pair_constructor_test.jl")
include("pair_conversion_test.jl")
include("pair_operator_test.jl")

function run_pair_tests()
    @testset "Pair tests" begin
        test_all_pair_constructors()
        test_all_pair_conversions()
        test_all_pair_operators()
    end
end


text:
using Brillouin, PlotlyJS, Test
using Brillouin.KPaths: Bravais
using StaticArrays

# NB: For the plot functionality, there's not that many meaningful things to test against
#     in the output since it's a graphical thing. So, instead, we just check that things
#     run without erroring by checking that the return type is as expected

@testset "PlotlyJS" begin
    # ------------------------------------------------------------------------------------ #
    # Wigner-Seitz cells and k-space path visualization
    # 3D
    sgnum = 147
    Rs = convert(SVector{3,SVector{3,Float64}}, [[1.0,0,0], [-.5,√3/2,0], [0,0,1.25]])
    Gs = Bravais.reciprocalbasis(Rs)

    cell = wignerseitz(Gs)
    kp = irrfbz_path(sgnum, Rs)
    @test plot(cell) isa PlotlyJS.SyncPlot
    @test plot(kp) isa PlotlyJS.SyncPlot
    @test plot(cell, kp) isa PlotlyJS.SyncPlot

    # 2D
    sgnum = 5
    cntr = Bravais.centering(sgnum, 2) # c-type centered
    Rs = convert(SVector{2,SVector{2,Float64}}, [[1.0,0], [0,.8]])
    pGs = Bravais.primitivize(Bravais.reciprocalbasis(Rs), cntr)

    cell = wignerseitz(pGs)
    kp = irrfbz_path(sgnum, Rs)
    @test plot(cell) isa PlotlyJS.SyncPlot
    @test plot(kp) isa PlotlyJS.SyncPlot
    @test plot(cell, kp) isa PlotlyJS.SyncPlot

    # ------------------------------------------------------------------------------------ #
    # Band structure plots (example from docs)
    sgnum = 227
    Rs = [[1,0,0], [0,1,0], [0,0,1]]
    kp = irrfbz_path(sgnum, Rs)
    pGs = basis(kp)
    kpi = interpolate(kp, 100)

    function ϵ(k; γ::Real=1.0)
        kx = 2π*(-k[1]+k[2]+k[3])
        ky = 2π*(+k[1]-k[2]+k[3])
        kz = 2π*(+k[1]+k[2]-k[3])
        return 4γ * (cos(kx/2)*cos(ky/2) + cos(ky/2)*cos(kz/2) + cos(kz/2)*cos(kx/2))
    end

    band1 = ϵ.(kpi)
    band2 = 20 .- (1/2).*band1

    # basic functionality
    @test plot(kpi, [band1]) isa PlotlyJS.SyncPlot
    @test plot(kpi, reduce(hcat, [band1, band2])) isa PlotlyJS.SyncPlot # matrix input

    # annotations
    offset = 2maximum(band2)
    band3 = band1 .+ offset
    annotations = Dict(:X => [1:2 => "X, bands 1:2", 30.0 => "X, some frequency"],
                       :Γ => [1 => "Γ, band 1", 3 => "Γ, band 3"])
    @test plot(kpi, [band1, band2, band3]; annotations) isa PlotlyJS.SyncPlot

    # band highlights
    band_highlights = Dict(2:3 => attr(color=:black, width=3))
    @test plot(kpi, [band1, band2, band3]; band_highlights) isa PlotlyJS.SyncPlot

    # title/y-label
    @test plot(kpi, [band1]; title="hello world", ylabel="ωa/2πc") isa PlotlyJS.SyncPlot
end


text:
"""
    module FELinearAlgebra
Provides Finite Element linear algebra matrices
"""
module FELinearAlgebra

import LinearAlgebra
import Base: *, +, -, ==, \
import Base: getindex, setindex!, size

export
    # Types
    LocalMatrix,
    JacobianMatrix,
    SolutionVector,
    # Functions
    L2_error!,
    L2_norm,
    lp_error!,
    lp_norm,
    true_size

include("SolutionVector.jl")
include("LocalMatrix.jl")

end

max_stars_repo_name: NoseKnowsAll/DGToolkit

max_issues_repo_name: NoseKnowsAll/DGToolkit

max_forks_repo_name: NoseKnowsAll/DGToolkit


text:
"""
This is Tutorial_Flegrei.jl. It contains all the necessary commands to import the data files from
geophysical models at the volcano and represent them in cartesian and UTM coordinates.
The following steps are performed:
1. Export earthquake data with a time label in both formats.
2. Export a velocity models in both formats.
3. Scale and interpolate a shear-wave model and export in both formats.

You will need to download the zipped folder containing all files from:
[https://seafile.rlp.net/f/ff2c8424274c4d56b1f7/](https://ngdc.noaa.gov/mgg/global/global.html)
Remember to work in the downloaded directory.

"""
using DelimitedFiles, GeophysicalModelGenerator, Glob, GeoStats

# 1. define where the two files are located on your computer
data_80s            = readdlm("SeismicLocations/Seismicity_UTM_1983_1984.txt", '\t', skipstart=0, header=false);
data_00s            = readdlm("SeismicLocations/Seismicity_UTM_2005_2016.txt", ' ', skipstart=0, header=false);

# 1.1 create a single file and detine the timing and locations
data                = vcat(data_80s, data_00s)
time                = data[:,1];
WE                  = data[:,2];
SN                  = data[:,3];
depth               = data[:,4];

# 1.2 save in both formats
EQ_Data_Cart        = CartData(WE, SN, depth, (Depth = depth * m, Time = time * yr,));
Write_Paraview(EQ_Data_Cart, "CF_Earthquakes_Cartesian", PointsData=true)
EQ_Data_UTM         = UTMData(WE, SN, depth, 33, true, (Depth = depth * m, Time = time * yr,));
Data_set_UTM        =   convert(GeophysicalModelGenerator.GeoData, EQ_Data_UTM)
Write_Paraview(Data_set_UTM, "CF_Earthquakes_UTM", PointsData=true)

# 2 load the velocity model and assign the locations of nodes and velocities + Vp/Vs
data            =   readdlm("TravelTimeTomography/modvPS.dat", '\t', Float64, skipstart=0, header=false);
WE              =   data[:,1];
SN              =   data[:,2];
depth           =   data[:,3];
Vp              =   data[:,4];
Vs              =   data[:,5];
VpVs            =   data[:,6];

# 2 Reshape everything and save in both formats.
resolution      =   (length(unique(depth)),  length(unique(SN)), length(unique(WE)))
dim_perm        =   [3 2 1]
we              =   permutedims(reshape(WE, resolution), dim_perm);
sn              =   permutedims(reshape(SN, resolution), dim_perm);
depth           =   permutedims(reshape(depth, resolution), dim_perm);
Vp3d            =   permutedims(reshape(Vp, resolution), dim_perm);
Vs3d            =   permutedims(reshape(Vs, resolution), dim_perm);
Vp_Vs3d         =   permutedims(reshape(VpVs, resolution), dim_perm);
Data_set_Cartesian  =   CartData(we, sn, depth, (vp = Vp3d * (km / s), vs = Vs3d * (km / s), vpvs = Vp_Vs3d,))
Write_Paraview(Data_set_Cartesian, "CF_Velocity_Cartesian")
Data_set        =   UTMData(we, sn, depth, 33, true, (vp = Vp3d * (km / s), vs = Vs3d * (km / s), vpvs = Vp_Vs3d,))
Data_set_UTM    =   convert(GeophysicalModelGenerator.GeoData, Data_set)
Write_Paraview(Data_set_UTM, "CF_Velocity_UTM")

# 3 Create a list of text files containing sections and loop through them
list_files        = glob("AmbientNoiseTomography/*.txt");
li                = size(list_files, 1);
for i = 1:li
    nameFile        = list_files[i];
    name_vts        = nameFile[24:26];
    data            = readdlm(nameFile, '\t', Float64);
    WE              = data[:,1];
    SN              = data[:,2];
    depth           = data[:,3];
    Vs              = data[:,4];

# 3.1 Constrain them on a smaller grid:
    ind             = findall((WE .>= 419000) .& (WE .<= 435000) .& (SN .>= 4514000) .& (SN .<= 4528000));
    WE              = WE[ind];
    SN              = SN[ind];
    depth           = depth[ind];
    Vs              = Vs[ind];

  # Then create a 3D lat,long regular grid from UTM - measurements are all in the same time zone:
    l               = length(WE);

    n_WE            = minimum(WE):100:maximum(WE);
    n_SN            = minimum(SN):100:maximum(SN);
    we, sn, Depth   = XYZGrid(n_WE, n_SN, depth[1]);
    Vs_3D           = zeros(size(Depth));
    Cgrid           = GeoStats.CartesianGrid((size(we, 1), size(we, 2)), (minimum(we), minimum(sn)), (we[2,2,1] - we[1,1,1], sn[2,2,1] - sn[1,1,1]))
    coord           = PointSet([WE[:]'; SN[:]']);
    Geo             = georef((Vs = Vs[:],), coord);
    P               = EstimationProblem(Geo, Cgrid, :Vs);
    S               = IDW(:Vs => (;neighbors=2));
    sol             = solve(P, S);
    sol_Vs          = values(sol).Vs;
    Vs_2D           = reshape(sol_Vs, size(domain(sol)));
    Vs_3D[:,:,1]    = Vs_2D;
    Data_set_Cart   = CartData(we, sn, Depth, (Vs = Vs_3D  * (km / s),))
    Write_Paraview(Data_set_Cart, "CF_Noise" * name_vts * "_Cartesian")
    Data_set        = UTMData(we, sn, Depth, 33, true, (Vs = Vs_3D * (km / s),));
    Data_set_UTM    = convert(GeophysicalModelGenerator.GeoData, Data_set);
    Write_Paraview(Data_set_UTM, "CF_Noise" * name_vts * "_UTM")
end


text:
# ---
# title: 1556. Thousand Separator
# id: problem1556
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/thousand-separator/description/>
# hidden: true
# ---
#
# Given an integer `n`, add a dot (".") as the thousands separator and return it
# in string format.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 987
#     Output: "987"
#
#
# **Example 2:**
#
#
#
#     Input: n = 1234
#     Output: "1.234"
#
#
# **Example 3:**
#
#
#
#     Input: n = 123456789
#     Output: "123.456.789"
#
#
# **Example 4:**
#
#
#
#     Input: n = 0
#     Output: "0"
#
#
#
#
# **Constraints:**
#
#   * `0 <= n < 2^31`
#
#
## @lc code=start
using LeetCode

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


text:
using ITensors,
      Test

@test LatticeBond(1, 2) == LatticeBond(1, 2, 0.0, 0.0, 0.0, 0.0, "")
@testset "Square lattice" begin
    sL = squareLattice(3, 4)
    @test length(sL) == 17
end

@testset "Triangular lattice" begin
    tL = triangularLattice(3, 4)
    @test length(tL) == 23
end


text:
using MeshBasisFunctions
using Meshes
using Test

points = Point2[(0,0), (1,0), (0,1), (1,1),
                (0.25,0.5), (0.75,0.5), (0.5, 0), (0.5, 1.0)]
tris = connect.([(1,5,3), (1, 7, 5), (7, 6, 5), (7, 2, 6),
                 (4, 6, 2), (4, 8, 6), (8, 5, 6), (8, 3, 5)],
                Triangle)
mesh = SimpleMesh(points, tris)
md = meshdata(mesh, Dict(0 => (z = randn(8), )))

basis = PiecewiseLinearBasis{2}()

interp_pts = [Point(0.0, 0.0),
              Point(0.5, 0.5),
              Point(0.34, 0.98),
              Point(0.12, 0.84)]
interp_op = interpolation_operator(interp_pts,
                                   mesh,
                                   basis)
interp1 = interp_op * values(md, 0)[:z]
interp2 = [interpolate(pt, md, :z, basis) for pt in interp_pts]

@testset "MeshBasisFunctions.jl" begin
    @test all(interp1 .== interp2)
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

using .Convex
using .ECOS

function _kliep_coeffs(K_nu, K_de, dre::KLIEP, optlib::Type{ConvexLib})
  # retrieve parameters
  σ, b = dre.σ, size(K_de, 2)

  # number of numerator and denominator samples
  n_nu, n_de = size(K_nu, 1), size(K_de, 1)

  # constants for objective and constraints
  K = K_nu
  k = sum(K_de, dims=1)

  # objective function and constraints
  α = Convex.Variable(b); w = K*α
  objective   = sum(log(w[i]) for i in 1:n_nu)
  constraints = [α ≥ 0, dot(α,k) == n_de]

  # optimization problem
  problem = Convex.maximize(objective, constraints)

  # solve problem with ECOS solver
  Convex.solve!(problem, ECOS.Optimizer(verbose=false))

  # optimal coefficients
  vec(α.value)
end


text:
using FractionalDiffEq
using Plots, LaTeXStrings

s="\$ODE\\ Example\$"

T = 30
h = 0.05
tspan = collect(0.05:h:T)

f(x)=1/2*(-exp(-x)-sin(x)-cos(x)+2)
target=f.(tspan)

rightfun(x) = sin(x)
prob = MultiTermsFODEProblem([1, 1], [2, 1], rightfun, [0, 0], T)
sol = solve(prob, h, FODEMatrixDiscrete())

plot(sol, title=s, legend=:bottomright, label="ODE Numerical Solution!")

plot!(tspan, target, lw=3,ls=:dash,label="ODE Analytical Solution!")


text:
geneCols = Dict{String,Symbol}(
   "aliases" => :GeneSymbolAlias,
   "symbol" => :OfficialGeneSymbol,
   "region" => :GeneRegion,
   "regionStart" => :GeneStart,
   "regionEnd" => :GeneEnd,
   "id" => :EnsemblGeneId,
   "version" => :EnsemblGeneVersion,
   "biotype" => :EnsemblGeneBiotype,
   "name" => :GeneName
);
transcriptCols = Dict{String,Symbol}(
   "region" => :TranscriptRegion,
   "regionStart" => :TranscriptStart,
   "regionEnd" => :TranscriptEnd,
   "id" => :EnsemblTranscriptId,
   "version" => :EnsemblTranscriptVersion,
   "biotype" => :EnsemblTranscriptBiotype
);
proteinCols = Dict{String,Symbol}(
   "id" => :EnsemblProteinId,
   "version" => :EnsemblProteinVersion
);
elTypeDict = Dict{String, Union}(
   # "aliases" => Vector{String},
   "aliases" => Union{String, Missing},
   "symbol" => Union{String, Missing},
   "region" => Union{String, Missing},
   "regionStart" => Union{Int, Missing},
   "regionEnd" => Union{Int, Missing},
   "id" => Union{String, Missing},
   "version" => Union{Int, Missing},
   "biotype" => Union{String, Missing},
   "name" => Union{String, Missing}
);

function parseGngJson(json::Array{Any,1})

   #########################
   # Count number or rows and create columns
   #  In the case of the columns only look at the first entry
   #  First four columns are inputDB, inputID, targetID and targetDB
   nCols = 4;
   colNames = Symbol["InputId", "InputSourceDb", "TargetId", "TargetDb"];
   colNamesUnstack = Symbol["InputId", "InputSourceDb"];
   colElTypes = Union[Union{String, Missing}, Union{String, Missing}, Union{String, Missing}, Union{String, Missing}];
   # Gene column names and types
   map(i->haskey(geneCols,i) && (push!(colNames, geneCols[i]),
                                 push!(colNamesUnstack, geneCols[i]),
                                 push!(colElTypes, elTypeDict[i])), keys(json[1]["row"][1]["gene"]))
   # Transcript:
   # Dict{String,Any} with 4 entries:
   # "peptide"   => Dict{String,Any}(Pair{String,Any}("id", "ENSP00000381185"),Pair{String,Any}("targetIds", Any[Dict{String,Any}(Pair{String,Any}("dbName", "UniProtKB/TrEMBL"),Pair{String,Any}("id", "A0A024R2B3"))]),Pair{String,A…
   # "id"        => "ENST00000398117"
   # "targetIds" => Any[]
   # "region"    => "18"
   map(i->haskey(transcriptCols,i) && (push!(colNames, transcriptCols[i]),
                                       push!(colNamesUnstack, transcriptCols[i]),
                                       push!(colElTypes, elTypeDict[i])), keys(json[1]["row"][1]["gene"]["transcript"][1]))
   # number of cols minus 1 because of targetIds, e.g.:
   # Dict{String,Any} with 3 entries:
   # "id"        => "ENSP00000381185"
   # "targetIds" => Any[Dict{String,Any}(Pair{String,Any}("dbName", "UniProtKB/TrEMBL"),Pair{String,Any}("id", "A0A024R2B3"))]
   # "version"   => 1
   map(i->haskey(proteinCols,i) && (push!(colNames, proteinCols[i]),
                                    push!(colNamesUnstack, proteinCols[i]),
                                    push!(colElTypes, elTypeDict[i])), keys(json[1]["row"][1]["gene"]["transcript"][1]["peptide"]))
   ##################
   # Rows
   nRows = 0;
   for entry in json
      # For each gene entry
      # nRows += 1;
      # Number of external gene identifiers
      gene = entry["row"][1]["gene"];
      nRows += (length(gene["targetIds"]) * length(gene["transcript"]));
      for transcript in gene["transcript"]
         # Transcript:
         # Dict{String,Any} with 4 entries:
         # "peptide"   => Dict{String,Any}(Pair{String,Any}("id", "ENSP00000381185"),Pair{String,Any}("targetIds", Any[Dict{String,Any}(Pair{String,Any}("dbName", "UniProtKB/TrEMBL"),Pair{String,Any}("id", "A0A024R2B3"))]),Pair{String,A…
         # "id"        => "ENST00000398117"
         # "targetIds" => Any[]
         # "region"    => "18"
         nRows += length(transcript["targetIds"]);
         # Peptide:
         # Dict{String,Any} with 3 entries:
         # "id"        => "ENSP00000381185"
         # "targetIds" => Any[Dict{String,Any}(Pair{String,Any}("dbName", "UniProtKB/TrEMBL"),Pair{String,Any}("id", "A0A024R2B3"))]
         # "version"   => 1
         nRows += length(transcript["peptide"]["targetIds"]);
      end
      # println(entry["row"][1])
   end
   nRows

   x = DataFrames.DataFrame(
      colElTypes,
      colNames,
      nRows
      # 100
   );

   lowG = 1;
   highG = 1;
   lowT = 1;
   highT = 1;
   for entry in json
      # For each gene entry
      lowT = lowG;
      highT = lowT;
      highG = lowG;
      # Number of external gene identifiers
      gene = entry["row"][1]["gene"];
      externalIds = String[];
      externalDbs = String[];
      # In the case of external IDs add entries to arrays
      !isempty(gene["targetIds"]) && map(i->(push!(externalIds, i["id"]), push!(externalDbs, i["dbName"])), gene["targetIds"]);
      externalGeneIds = copy(externalIds);
      externalGeneDbs = copy(externalDbs);
      # nRows += length(gene["targetIds"])
      for transcript in gene["transcript"]
         # In the case of external IDs add entries to arrays
         !isempty(transcript["targetIds"]) && map(i->(push!(externalIds, i["id"]), push!(externalDbs, i["dbName"])), transcript["targetIds"]);
         peptide = transcript["peptide"];
         !isempty(peptide["targetIds"]) && map(i->(push!(externalIds, i["id"]), push!(externalDbs, i["dbName"])), peptide["targetIds"]);
         highT += (length(externalIds) - 1);
         # for
         # map(i->haskey(geneCols,i) && (push!(colNames, geneCols[i]), push!(colElTypes, elTypeDict[i])), keys(json[1]["row"][1]["gene"]))
         # for epEntry in transcript["peptide"]["targetIds"]

         map(key->haskey(transcriptCols,key) && (x[transcriptCols[key]][lowT:highT] = transcript[key]), keys(transcript))
         map(key->haskey(proteinCols,key) && (x[proteinCols[key]][lowT:highT] = peptide[key]), keys(peptide))

         # map(key->haskey(geneCols,key) && (println("col: $(geneCols[key]) = $(gene[key])")), keys(gene))
         # map(key->haskey(geneCols,key) && (println("$(collect(lowT:highT)) -> $(fill(gene[key], high))")), keys(gene))

         x[:TargetId][lowT:highT] = externalIds;
         x[:TargetDb][lowT:highT] = externalDbs;
         externalIds = copy(externalGeneIds);
         externalDbs = copy(externalGeneDbs);
         # x[low,:EnsemblGeneId] = gene["id"];
         # x[low,:EnsemblTranscriptId] = transcript["id"];
         # x[low,:EnsemblProteinId] = peptide["id"];
         # println("[$i] : $externalId, $externalDb")
         # println.("enspId: $(transcript["peptide"]["id"]) -> $(transcript["peptide"]["targetIds"])")
         lowT = highT + 1;
         highT = lowT;
      end
      highG = highT-1;
      # Fill in gene information
      for (key, value) in gene
         if haskey(geneCols, key) && !isempty(value)
            val = value;
            if key == "aliases"
               val = join(value, ",");
            end
            x[geneCols[key]][lowG:highG] = val;
         end
      end
      # map(key->haskey(geneCols,key) && (x[geneCols[key]][lowG:highG] = fill(join(gene[key], ","), (highG - lowG + 1))), keys(gene));
      x[:InputId][lowG:highG] = entry["row"][1]["inputId"];
      x[:InputSourceDb][lowG:highG] = entry["row"][1]["inputDb"];
      lowG = highG + 1;
   end
   x

   return x;

end

function unstackDf(x::D, keyCols::S = Vector{Symbol}(), keyColsElType::T = Vector{DataType}()) where {
   D<:DataFrames.DataFrame,
   S<:Vector{Symbol},
   T<:Vector{DataType}}

   if isempty(keyCols)
      keyCols = names(x)[map(i->!in(i,["TargetId", "TargetDb"]), collect(String,names(x)))];
      keyColsElType = DataFrames.eltypes(x)[map(i->!in(i,["TargetId", "TargetDb"]), collect(String,names(x)))];
   elseif length(keyCols != length(keyColsElType))
      error("keyCols and keyColsElType have to be of same length.!");
   end
   colNames = copy(keyCols);
   colElTypes = copy(keyColsElType);

   append!(colNames, DataFrames.identifier.(collect(unique(x[:TargetDb]))));
   append!(colElTypes, fill(Union{String, Missing}, (length(colNames) - length(keyCols))));
   nRows = DataFrames.nrow(unique(x[keyCols]));

   xDf = DataFrames.DataFrame(
      colElTypes,
      colNames,
      nRows
   );
   sortByCols = copy(keyCols);
   push!(sortByCols, :TargetDb);
   sort!(x, cols = sortByCols);
   idx = 0;
   key = "";
   for row in DataFrames.eachrow(x)
      keyCurr = join(map(i->i[2], row[keyCols]));
      if key != keyCurr
         idx += 1;
         key = keyCurr;
         map(i->xDf[idx,i[1]] = i[2], row[keyCols]);
      end
      colId = DataFrames.identifier(row[:TargetDb])
      xDf[idx, colId] = DataFrames.ismissing(xDf[idx, colId]) ?
                           row[:TargetId] :
                           xDf[idx, colId] * "," * row[:TargetId];
   end
   return xDf;
end

max_stars_repo_name: StefanHaunsberger/GeneNameGenieJ

max_issues_repo_name: StefanHaunsberger/GeneNameGenieJ

max_forks_repo_name: StefanHaunsberger/GeneNameGenieJ


text:
@testset "NonOverlappingOrthotopes2DisjunctionLP: $(fct_type), orthotopes $(n_ortho), dimension $(dim), $(T)" for fct_type in ["vector of variables", "vector affine function"], n_ortho in [2, 3], dim in [2, 3], T in [Int, Float64]
    mock = MOIU.MockOptimizer(DisjunctionModel{T}())
    model = COIB.NonOverlappingOrthotopes2DisjunctionLP{T}(mock)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    if fct_type == "vector of variables"
        @test MOI.supports_constraint(
            model,
            MOI.VectorOfVariables,
            CP.Disjunction{NTuple{4, MOI.LessThan{T}}},
        )
        @test MOI.supports_constraint(
            model,
            MOI.VectorOfVariables,
            CP.Disjunction{NTuple{6, MOI.LessThan{T}}},
        )
        @test MOIB.supports_bridging_constraint(
            model,
            MOI.VectorOfVariables,
            CP.NonOverlappingOrthotopes{CP.UNCONDITIONAL_NONVERLAPPING_ORTHOTOPES},
        )
    elseif fct_type == "vector affine function"
        @test MOI.supports_constraint(
            model,
            MOI.VectorAffineFunction{T},
            CP.Disjunction{NTuple{4, MOI.LessThan{T}}},
        )
        @test MOI.supports_constraint(
            model,
            MOI.VectorAffineFunction{T},
            CP.Disjunction{NTuple{6, MOI.LessThan{T}}},
        )
        @test MOIB.supports_bridging_constraint(
            model,
            MOI.VectorAffineFunction{T},
            CP.NonOverlappingOrthotopes{CP.UNCONDITIONAL_NONVERLAPPING_ORTHOTOPES},
        )
    else
        @assert false
    end

    if T == Int
        x_pos, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:(dim * n_ortho)])
        x_sze, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:(dim * n_ortho)])
        x_end, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:(dim * n_ortho)])
    elseif T == Float64
        x_pos = MOI.add_variables(model, dim * n_ortho)
        x_sze = MOI.add_variables(model, dim * n_ortho)
        x_end = MOI.add_variables(model, dim * n_ortho)
    end
    x_ortho = [
        [
            x_pos[(1 + (i - 1) * dim) : (i * dim)],
            x_sze[(1 + (i - 1) * dim) : (i * dim)],
            x_end[(1 + (i - 1) * dim) : (i * dim)],
        ]
        for i in 1:n_ortho
    ]
    x = vcat([vcat(x_ortho[i]...) for i in 1:n_ortho]...)

    # Test the construction of x before going on (no loops to ease checking
    # for correctness).
    @test length(x_pos) == dim * n_ortho
    @test length(x_sze) == dim * n_ortho
    @test length(x_end) == dim * n_ortho
    @test length(x) == 3 * dim * n_ortho
    if n_ortho >= 1
        @test x[(0 * dim + 1):(1 * dim)] == x_pos[(0 * dim + 1):(1 * dim)]
        @test x[(1 * dim + 1):(2 * dim)] == x_sze[(0 * dim + 1):(1 * dim)]
        @test x[(2 * dim + 1):(3 * dim)] == x_end[(0 * dim + 1):(1 * dim)]
    end
    if n_ortho >= 2
        @test x[(3 * dim + 1):(4 * dim)] == x_pos[(1 * dim + 1):(2 * dim)]
        @test x[(4 * dim + 1):(5 * dim)] == x_sze[(1 * dim + 1):(2 * dim)]
        @test x[(5 * dim + 1):(6 * dim)] == x_end[(1 * dim + 1):(2 * dim)]
    end
    if n_ortho >= 3
        @test x[(6 * dim + 1):(7 * dim)] == x_pos[(2 * dim + 1):(3 * dim)]
        @test x[(7 * dim + 1):(8 * dim)] == x_sze[(2 * dim + 1):(3 * dim)]
        @test x[(8 * dim + 1):(9 * dim)] == x_end[(2 * dim + 1):(3 * dim)]
    end
    if n_ortho >= 4
        @assert false
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables(x)
    elseif fct_type == "vector affine function"
        MOIU.vectorize(x)
    else
        @assert false
    end

    c = MOI.add_constraint(model, fct, CP.NonOverlappingOrthotopes(n_ortho, dim))

    for i in 1:(3 * dim *  n_ortho)
        @test MOI.is_valid(model, x[i])
    end
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.VectorOfVariables, CP.NonOverlappingOrthotopes{CP.UNCONDITIONAL_NONVERLAPPING_ORTHOTOPES}}(-1)]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.NonOverlappingOrthotopes{CP.UNCONDITIONAL_NONVERLAPPING_ORTHOTOPES}) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == Tuple{Type}[]
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
            (MOI.VectorAffineFunction{T}, CP.Disjunction{NTuple{n, MOI.LessThan{T}} where n}),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == 0
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}()) == n_ortho * dim
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, CP.Disjunction{NTuple{n, MOI.LessThan{T}} where n}}()) == Int(n_ortho * (n_ortho - 1) / 2)

        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}()) == collect(values(bridge.cons_ends))
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, CP.Disjunction{NTuple{n, MOI.LessThan{T}} where n}}()) == collect(values(bridge.cons_disjunction))
    end

    @testset "End-point constraints" begin
        @test length(bridge.cons_ends) == n_ortho * dim
        for i in 1:n_ortho
            for d in 1:dim
                @test MOI.is_valid(model, bridge.cons_ends[i, d])
                @test MOI.get(model, MOI.ConstraintSet(), bridge.cons_ends[i, d]) == MOI.EqualTo(zero(T))
                f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_ends[i, d])

                @test length(f.terms) == 3

                t1 = f.terms[1]
                @test t1.coefficient === one(T)
                @test t1.variable === x_pos[(i - 1) * dim + d]

                t2 = f.terms[2]
                @test t2.coefficient === one(T)
                @test t2.variable === x_sze[(i - 1) * dim + d]

                t3 = f.terms[3]
                @test t3.coefficient === -one(T)
                @test t3.variable === x_end[(i - 1) * dim + d]
            end
        end
    end

    @testset "Disjunction constraints" begin
        @test length(bridge.cons_disjunction) == Int(n_ortho * (n_ortho - 1) / 2)
        for i in 1:n_ortho
            for j in 1:n_ortho
                if i < j
                    @test MOI.is_valid(model, bridge.cons_disjunction[i, j])
                    @test MOI.get(model, MOI.ConstraintSet(), bridge.cons_disjunction[i, j]) isa CP.Disjunction

                    f_vector = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_disjunction[i, j])
                    @test length(f_vector.terms) == 3 * 2 * dim

                    f_scalars = MOIU.scalarize(f_vector)
                    @test length(f_scalars) == 2 * dim

                    for d in 1:dim
                        f1 = f_scalars[(d - 1) * 2 + 1]
                        f2 = f_scalars[(d - 1) * 2 + 2]

                        # Doesn't work, because == requires that the order
                        # of the terms is the same.
                        # @test f1 == one(T) * x_pos[(i - 1) * dim + d] + one(T) * x_sze[(i - 1) * dim + d] - one(T) * x_pos[(j - 1) * dim + d]

                        t = f1.terms[1]
                        @test t.coefficient === one(T)
                        @test t.variable === x_pos[(i - 1) * dim + d]

                        t = f1.terms[2]
                        @test t.coefficient === -one(T)
                        @test t.variable === x_pos[(j - 1) * dim + d]

                        t = f1.terms[3]
                        @test t.coefficient === one(T)
                        @test t.variable === x_sze[(i - 1) * dim + d]

                        t = f2.terms[1]
                        @test t.coefficient === -one(T)
                        @test t.variable === x_pos[(i - 1) * dim + d]

                        t = f2.terms[2]
                        @test t.coefficient === one(T)
                        @test t.variable === x_pos[(j - 1) * dim + d]

                        t = f2.terms[3]
                        @test t.coefficient === one(T)
                        @test t.variable === x_sze[(j - 1) * dim + d]
                    end
                end
            end
        end
    end
end


text:
using IcanProj
using PyPlot
using DataFrames
using PyCall
using JFSM2
using CSV
using Statistics

import IcanProj.project_results

function project_results(df::DataFrame, on::Symbol)

    map_values = fill(NaN, 1550, 1195)

    map_values[df[:ind_julia]] = df[on]

    return map_values

end

function plot_map(df, variable, cb_label, cb_unit, title_label, figpath)

    resmap = project_results(df, variable)

    file = joinpath(figpath, String(variable) * ".png")

    cb_text =  cb_label * " " * cb_unit

    py"""

    import matplotlib.pyplot as plt

    fig, ax = plt.subplots()

    #fig.figsize = (6, 12)

    ax.set_title($(title_label))

    img_plot = ax.imshow($(resmap), cmap = 'jet', vmin = -0.1, vmax = 1.1)

    cbar = fig.colorbar(img_plot)

    cbar.set_label($(cb_text))

    #cbar.set_clim(0, 1)

    ax.grid(linestyle='dotted')

    plt.tick_params(axis='both', left=False, top=False, right=False, bottom=False,
                    labelleft=False, labeltop=False, labelright=False, labelbottom=False)

    plt.tight_layout()

    fig.savefig($(file), dpi = 200, bbox_inches='tight')

    plt.close('all')

    """

    @info "finished plotting"

end

# Load results

path = dirname(pathof(IcanProj))

df_res = CSV.File(joinpath(path, "..", "data", "table_errors_50km.txt")) |> DataFrame

df_links = CSV.read(joinpath(path, "..", "data", "df_links.csv")) |> DataFrame

df_all = join(df_links, df_res, on = :ind_50km)

cfgs = cfg_table()

# Plot average error for all cfgs with exchng = 0

for variable in ["swe", "latmo", "hatmo", "melt"]

    cfg_subset = cfgs[convert(Array{Bool}, cfgs[:exchng] .== 0), :]

    colnames = Symbol.("nse_$(variable)_cfg" .* string.(cfg_subset[:cfg]))

    df_all[Symbol("nse_$(variable)_exchng0")] =  mean(convert(Array, df_all[colnames]), dims = 2)[:]

end

# Plot average error for all cfgs with exchng = 1

for variable in ["swe", "latmo", "hatmo", "melt"]

    cfg_subset = cfgs[convert(Array{Bool}, cfgs[:exchng] .== 1), :]

    colnames = Symbol.("nse_$(variable)_cfg" .* string.(cfg_subset[:cfg]))

    df_all[Symbol("nse_$(variable)_exchng1")] =  mean(convert(Array, df_all[colnames]), dims = 2)[:]

end

# Plot averaged results

info = [(:nse_swe_exchng0,   "NSE", "(-)", "SWE exchnge = 0"),
        (:nse_latmo_exchng0, "NSE", "(-)", "Latmo exchnge = 0"),
        (:nse_hatmo_exchng0, "NSE", "(-)", "Hatmo exchnge = 0"),
        (:nse_melt_exchng0,  "NSE", "(-)", "Melt exchnge = 0"),
        (:nse_swe_exchng1,   "NSE", "(-)", "SWE exchnge = 1"),
        (:nse_latmo_exchng1, "NSE", "(-)", "Latmo exchnge = 1"),
        (:nse_hatmo_exchng1, "NSE", "(-)", "Hatmo exchnge = 1"),
        (:nse_melt_exchng1,  "NSE", "(-)", "Melt exchnge = 1")]

figpath = joinpath(path, "..", "plots", "error_maps_1")

for (v, cl, cu, t) in info

    plot_map(df_all, v, cl, cu, t, figpath)

end

#variable = :nse_rnet_cfg1

#cb_label = "Correlation coeffcient"

#cb_unit = "(-)"

#title_label = "SWE"

#figpath = "dummy"

#plot_map(df_all, variable, cb_label, cb_unit, title_label, figpath)

#mapres = project_results(df, :r_swe_cfg1)

#figure()
#imshow(mapres)


text:
using Statistics
x = rand(10^7);
data = 1 ./ (x .^ 2);
a,b,c,d = Statistics.quantile!(data, [0, .001, .999, 1]);
#used_bits(d) # Naive required bit-width
println("min=",a,", 0.1%=", b, ", 99.9%=",c,", max=",d)


text:
s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true)
@testset "test ac polar opf" begin
    @testset "3-bus case" begin
        result = run_acdcopf("../test/data/case3.m", ACPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 5907; atol = 1e0)
    end
    @testset "5-bus asymmetric case" begin
        result = run_acdcopf("../test/data/case5_asym.m", ACPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 17551; atol = 1e0)
    end
    @testset "5-bus ac dc case" begin
        result = run_acdcopf("../test/data/case5_acdc.m", ACPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 194.14; atol = 1e0)
    end
    @testset "5-bus ac dc case with 2 seperate ac grids" begin
        result = run_acdcopf("../test/data/case5_2grids.m", ACPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 397.36; atol = 1e0)
    end
    @testset "24-bus rts ac dc case with three zones" begin
        result = run_acdcopf("../test/data/case24_3zones_acdc.m", ACPPowerModel, ipopt_solver; setting = s)

    @test isapprox(result["objective"], 150228.15; atol = 1e0)
    end
end

@testset "test dc opf" begin
    @testset "3-bus case" begin
        result = run_acdcopf("../test/data/case3.m", DCPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 5638.97; atol = 1e0)
    end
    @testset "5-bus asymmetric case" begin
        result = run_acdcopf("../test/data/case5_asym.m", DCPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 17479; atol = 1e0)
    end
    @testset "5-bus ac dc case" begin
        result = run_acdcopf("../test/data/case5_acdc.m", DCPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 171.69; atol = 1e0)
    end
    @testset "5-bus ac dc case with 2 seperate ac grids" begin
        result = run_acdcopf("../test/data/case5_2grids.m", DCPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 371.01; atol = 1e0)
    end

    @testset "24-bus rts ac dc case with three zones" begin
       result = run_acdcopf("../test/data/case24_3zones_acdc.m", DCPPowerModel, ipopt_solver; setting = s)

       @test result["termination_status"] == LOCALLY_SOLVED
       @test isapprox(result["objective"], 143714.75; atol = 1e0)
    end
end

@testset "test soc (BIM) opf" begin
    @testset "3-bus case" begin
        result = run_acdcopf("../test/data/case3.m", SOCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 5746.7; atol = 1e0)
    end
    @testset "5-bus asymmetric case" begin
        result = run_acdcopf("../test/data/case5_asym.m", SOCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 14999; atol = 1e0)
    end
    @testset "5-bus ac dc case" begin
        result = run_acdcopf("../test/data/case5_acdc.m", SOCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 183.76; atol = 1e0)
    end
    @testset "5-bus ac dc case with 2 seperate ac grids" begin
        result = run_acdcopf("../test/data/case5_2grids.m", SOCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 363.50; atol = 1e0)
    end
    @testset "24-bus rts ac dc case with three zones" begin
        result = run_acdcopf("../test/data/case24_3zones_acdc.m", SOCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 150156.24; atol = 1e0)
    end
end

@testset "test soc distflow opf_bf" begin
    @testset "3-bus case" begin
        result = run_acdcopf_bf("../test/data/case3.m", SOCBFPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 5746.7; atol = 1e0)
    end
    @testset "5-bus asymmetric case" begin
        result = run_acdcopf_bf("../test/data/case5_asym.m", SOCBFPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 14999; atol = 1e0)
    end
    @testset "5-bus ac dc case" begin
        result = run_acdcopf_bf("../test/data/case5_acdc.m", SOCBFPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 183.91; atol = 1e0)
    end
    @testset "5-bus ac dc case with 2 seperate ac grids" begin
        result = run_acdcopf_bf("../test/data/case5_2grids.m", SOCBFPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 363.50; atol = 1e0)
    end
    @testset "24-bus rts ac dc case with three zones" begin
        result = run_acdcopf_bf("../test/data/case24_3zones_acdc.m", SOCBFPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 150156.26; atol = 1e0)
    end
end

@testset "test qc opf" begin
    @testset "3-bus case" begin
        result = run_acdcopf("../test/data/case3.m", QCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 5780; atol = 1e0)
    end
    @testset "5-bus asymmetric case" begin
        result = run_acdcopf("../test/data/case5_asym.m", QCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 15921; atol = 1e0)
    end
    @testset "5-bus ac dc case" begin
        result = run_acdcopf("../test/data/case5_acdc.m", QCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 183.76; atol = 1e0)
    end
    @testset "24-bus rts ac dc case with three zones" begin
        result = run_acdcopf("../test/data/case24_3zones_acdc.m", QCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 150156.25; atol = 1e0)
    end
end

@testset "test qc opf with trilinear convexhull relaxation" begin
    @testset "3-bus case" begin
        result = run_acdcopf("../test/data/case3.m", QCWRTriPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 5817.58; atol = 1e0)
    end
    # @testset "5-bus asymmetric case" begin
    #     result = run_acdcopf("../test/data/case5_asym.m", QCWRTriPowerModel, ipopt_solver)
    #
    #     @test result["termination_status"] == LOCALLY_SOLVED
    #     @test isapprox(result["objective"], 15816.9; atol = 1e0)
    # end
    @testset "5-bus ac dc case" begin
        result = run_acdcopf("../test/data/case5_acdc.m", QCWRTriPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 183.76; atol = 1e0)
    end
    @testset "24-bus rts ac dc case with three zones" begin
        result = run_acdcopf("../test/data/case24_3zones_acdc.m", QCWRTriPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"],  150156.25; atol = 1e0)
    end
end

# @testset "test sdp opf" begin
#     @testset "3-bus case" begin
#         result = run_acdcopf("../test/data/case3.m", SDPWRMPowerModel, scs_solver; setting = s)
#
#         @test result["termination_status"] == :Optimal
#         @test isapprox(result["objective"], 5851.3; atol = 2e0)
#     end
#     # @testset "5-bus ac dc case" begin
#     #     result = run_acdcopf("../test/data/case5_acdc.m", SDPWRMPowerModel, mosek)
#     #
#     #     @test result["termination_status"] == :Optimal
#     #     @test isapprox(result["objective"], 194.12; atol = 1e0)
#     # end
#     # TODO replace this with smaller case, way too slow for unit testing
#     # @testset "24-bus rts ac dc case with three zones" begin
#     #    result = run_acdcopf("../test/data/case24_3zones_acdc.m", SDPWRMPowerModel, scs_solver)
#     #
#     #    @test result["termination_status"] == :Optimal
#     #    @test isapprox(result["objective"], 150179.96; atol = 1e0)
#     # end
# end


text:
"""
    info(registry)

Show information about a registry, specifically its name, description
and fields.

## Examples

{cell}
```julia
using FeatureRegistries: exampleregistry, info
registry = exampleregistry()
info(registry)
```
"""
function info end

struct RegistryInfo
    registry::Registry
end

function Base.show(io::IO, registryinfo::RegistryInfo)
    registry = registryinfo.registry
    println(io,
        "Information on registry ",
        crayon"underline", getfield(registry, :name), crayon"reset", "\n\n",
        "$(crayon"bold")Fields$(crayon"reset"): ")

    fs = registry.fields
    data = reduce(vcat, [
            reshape([
                f.name,
                AnsiTextCell(k in getfield(registry, :id) ? "$(crayon"red bold"):$k (ID)" : ":$k"),
                AnsiTextCell("$(crayon"dark_gray")$(_fieldtype(f))"),
                f.description],
            1, :)
            for (k, f) in pairs(fs)
        ])
    pretty_table(
        io,
        data,
        alignment=:l,
        vlines = :all,
        hlines = :none,
        vcrop_mode=:middle,
        tf=PrettyTables.tf_borderless,
        header = ["Name", "Column", "Type", "Description"]
    )

    if !isempty(registry.description)
        println(io, "\n$(crayon"bold")Description$(crayon"reset"): ")
        show(io, RichCell(Markdown.parse(registry.description)))
    end
end

_fieldtype(::Field{T}) where T = T
info(registry::Registry) = RegistryInfo(registry)


text:
module DocTestPkg

push!(LOAD_PATH, joinpath(Pkg.dir("DocTestPkg"), "src"))
using SupportMod

# DocTestPkg
export Equation, func, calcIt

"""
  Type holding all data for an instance of an equation
"""
type Equation{Tsol, Tres}
  q::Array{Tsol, 3}
  res::Array{Tres, 3}
  t::Float64
end

"""
  Outer constructor for [`Equation`](@ref)
"""
function Equation(numDofPerNode::Integer, numNodesPerElement::Integer,
                 numEl::Integer)
  q = zeros(numDofPerNode, numNodesPerElement, numEl)
  res = zeros(numDofPerNode, numNodesPerElement, numEl)
  t = 0.0

  return Equation{Float64, Float64}(q, res, t)
end

"""
Returns double the number `x` plus `1`
"""
func(y) = 2y + 1

"""
  This function calculates something
"""
function calcIt(x::Integer)
  return x + 2
end

end # module


text:
using QuantumStateBase
using Test
using LinearAlgebra
using ClassicalOrthogonalPolynomials

const QSB = QuantumStateBase
const DIM = 100

@testset "QuantumStateBase.jl" begin
    # polynomial
    include("polynomial.jl")

    # state
    include("operator.jl")
    include("state.jl")

    # wigner
    include("wigner_util.jl")
    include("wigner.jl")
end


text:
using ExaTron
using LinearAlgebra
using SparseArrays
using Random

using Profile
using PProf

Random.seed!(0)
n = 8
A = ExaTron.TronDenseMatrix{Array}(n)
L = ExaTron.TronDenseMatrix{Array}(n)
for j=1:n, i=j:n
    v = rand(1)[1]
    L.vals[i,j] = v
end
A.vals .= (L.vals * transpose(L.vals))

g = ones(n)
x = zeros(n)
xc = zeros(n)
xl = -Inf*ones(n)
xu = Inf*ones(n)
eval_f_cb(x) = 0.5*(transpose(x)*A.vals*x) + transpose(g)*x
function eval_grad_f_cb(x, grad_f)
    grad_f .= A.vals*x .+ g
end
function eval_h_cb(x, mode, rows, cols, scale, lambda, values)
    if mode == :Structure
        nz = 1
        for j=1:n,i=j:n
            rows[nz] = i
            cols[nz] = j
            nz += 1
        end
    else
        nz = 1
        for j=1:n,i=j:n
            values[nz] = A.vals[i,j]
            nz += 1
        end
    end
end

prob = ExaTron.createProblem(n, xl, xu, Int(floor((n*(n+1))/2)),
                             eval_f_cb, eval_grad_f_cb, eval_h_cb)#; :matrix_type=>:Dense)

howmany = 10
for i=1:howmany
    fill!(prob.x, 0.0)
    @profile ExaTron.solveProblem(prob)
end
pprof()


text:
using NPZ, Plots, NipponColors
cd(@__DIR__)

D = dict_color()

# sampler
gmv = npzread("generated_macroscopic_variables_normalized.npy")
∇gmv = npzread("generated_macroscopic_variables_normalized_gradients.npy")

# solver generator
smv = npzread("sod_macroscopic_variables_normalized.npy")
∇smv = npzread("sod_macroscopic_variables_gradients.npy")

# plot
p1 = scatter(gmv[:, 1], gmv[:, 2], marker = (:circle, 3, 0.6, D["ukon"], stroke(1, 0.1, :gray27, :dot)), legend=:none, xlims=(-1, 1), ylims=(0, 1), xlabel="U", ylabel="T")
savefig(p1, "illu_generator_u.pdf")

p2 = scatter(smv[:, 1], smv[:, 2], marker = (:circle, 3, 0.6, D["ukon"], stroke(1, 0.1, :gray27, :dot)), legend=:none, xlims=(-1, 1), ylims=(0, 1), xlabel="U", ylabel="T")
savefig(p2, "illu_solver_u.pdf")

p3 = scatter(∇gmv[:, 1], ∇gmv[:, 2], marker = (:circle, 3, 0.6, D["ukon"], stroke(1, 0.1, :gray27, :dot)), legend=:none, xlims=(-10, 10), ylims=(-10, 10), xlabel="∇U", ylabel="∇T")
savefig(p3, "illu_generator_du.pdf")

p4 = scatter(∇smv[:, 1], ∇smv[:, 2], marker = (:circle, 3, 0.6, D["ukon"], stroke(1, 0.1, :gray27, :dot)), legend=:none, xlims=(-10, 10), ylims=(-10, 10), xlabel="∇U", ylabel="∇T")
savefig(p4, "illu_solver_du.pdf")

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
module Surrogates
using LinearAlgebra
using Distributions

abstract type AbstractSurrogate <: Function end
include("utils.jl")
include("Radials.jl")
include("Kriging.jl")
include("Sampling.jl")
include("Optimization.jl")
include("Lobachesky.jl")
include("LinearSurrogate.jl")
include("InverseDistanceSurrogate.jl")
include("SecondOrderPolynomialSurrogate.jl")
include("SthenoKriging.jl")
include("RandomForestSurrogate.jl")
include("NeuralSurrogate.jl")
include("Wendland.jl")
include("MOE.jl")
include("VariableFidelity.jl")
include("PolynomialChaos.jl")
include("Earth.jl")

current_surrogates = ["Kriging","LinearSurrogate","LobacheskySurrogate","NeuralSurrogate",
                          "RadialBasis","RandomForestSurrogate","SecondOrderPolynomialSurrogate","Wendland"]
export current_surrogates_MOE
export AbstractSurrogate, SamplingAlgorithm
export Kriging, RadialBasis, add_point!, current_estimate, std_error_at_point
export linearRadial,cubicRadial,multiquadricRadial,thinplateRadial
export sample, GridSample, UniformSample, SobolSample, LatinHypercubeSample, LowDiscrepancySample
export RandomSample, KroneckerSample, GoldenSample
export SRBF,LCBS,EI,DYCORS,SOP,surrogate_optimize
export LobacheskySurrogate, lobachesky_integral, lobachesky_integrate_dimension
export LinearSurrogate
export RandomForestSurrogate
export SVMSurrogate
export NeuralSurrogate
export InverseDistanceSurrogate
export SecondOrderPolynomialSurrogate
export SthenoKriging
export Wendland
export RadialBasisStructure, KrigingStructure, LinearStructure, InverseDistanceStructure
export LobacheskyStructure, NeuralStructure, RandomForestStructure, SecondOrderPolynomialStructure
export WendlandStructure
export MOE
export VariableFidelitySurrogate
export PolynomialChaosSurrogate
export EarthSurrogate

end


text:
module ScalarWrapper

using JavaCall

const BaseMultiStartMultivariateOptimizer = @jimport org.hipparchus.optim.BaseMultiStartMultivariateOptimizer
const BaseMultivariateOptimizer = @jimport org.hipparchus.optim.BaseMultivariateOptimizer
const BaseOptimizer = @jimport org.hipparchus.optim.BaseOptimizer
const Class = @jimport java.lang.Class
const ConvergenceChecker = @jimport org.hipparchus.optim.ConvergenceChecker
const Enum = @jimport java.lang.Enum
const GoalType = @jimport org.hipparchus.optim.nonlinear.scalar.GoalType
const GradientMultivariateOptimizer = @jimport org.hipparchus.optim.nonlinear.scalar.GradientMultivariateOptimizer
const JString = @jimport java.lang.String
const LeastSquaresConverter = @jimport org.hipparchus.optim.nonlinear.scalar.LeastSquaresConverter
const LineSearch = @jimport org.hipparchus.optim.nonlinear.scalar.LineSearch
const MultiStartMultivariateOptimizer = @jimport org.hipparchus.optim.nonlinear.scalar.MultiStartMultivariateOptimizer
const MultivariateFunction = @jimport org.hipparchus.analysis.MultivariateFunction
const MultivariateFunctionMappingAdapter = @jimport org.hipparchus.optim.nonlinear.scalar.MultivariateFunctionMappingAdapter
const MultivariateFunctionPenaltyAdapter = @jimport org.hipparchus.optim.nonlinear.scalar.MultivariateFunctionPenaltyAdapter
const MultivariateOptimizer = @jimport org.hipparchus.optim.nonlinear.scalar.MultivariateOptimizer
const MultivariateVectorFunction = @jimport org.hipparchus.analysis.MultivariateVectorFunction
const Object = @jimport java.lang.Object
const ObjectiveFunction = @jimport org.hipparchus.optim.nonlinear.scalar.ObjectiveFunction
const ObjectiveFunctionGradient = @jimport org.hipparchus.optim.nonlinear.scalar.ObjectiveFunctionGradient
const OptimizationData = @jimport org.hipparchus.optim.OptimizationData
const Optional = @jimport java.util.Optional
const PointValuePair = @jimport org.hipparchus.optim.PointValuePair
const RandomVectorGenerator = @jimport org.hipparchus.random.RandomVectorGenerator
const RealMatrix = @jimport org.hipparchus.linear.RealMatrix
const UnivariatePointValuePair = @jimport org.hipparchus.optim.univariate.UnivariatePointValuePair

export GoalType
export GradientMultivariateOptimizer
export LeastSquaresConverter
export LineSearch
export MultiStartMultivariateOptimizer
export MultivariateFunctionMappingAdapter
export MultivariateFunctionPenaltyAdapter
export MultivariateOptimizer
export ObjectiveFunction
export ObjectiveFunctionGradient
export bounded_to_unbounded
export compare_to
export compute_objective_value
export describe_constable
export equals
export get_class
export get_convergence_checker
export get_declaring_class
export get_evaluations
export get_goal_type
export get_iterations
export get_lower_bound
export get_max_evaluations
export get_max_iterations
export get_objective_function
export get_objective_function_gradient
export get_optima
export get_start_point
export get_upper_bound
export hash_code
export name
export notify
export notify_all
export optimize
export ordinal
export search
export to_string
export unbounded_to_bounded
export value
export value_of
export values
export wait

include("GradientWrapper/GradientWrapper.jl")
include("NoderivWrapper/NoderivWrapper.jl")
include("goal_type.jl")
include("gradient_multivariate_optimizer.jl")
include("least_squares_converter.jl")
include("line_search.jl")
include("multi_start_multivariate_optimizer.jl")
include("multivariate_function_mapping_adapter.jl")
include("multivariate_function_penalty_adapter.jl")
include("multivariate_optimizer.jl")
include("objective_function.jl")
include("objective_function_gradient.jl")

end


text:
module NOVAS

using DataDeps, LinearAlgebra, Memoize, StaticArrays, DelimitedFiles, FortranFiles,
      Polynomials
import Base.show

include("data.jl")
include("constants.jl")
include("utils.jl")
include("nutation.jl")
include("daf.jl")
include("naif_codes.jl")
include("spk.jl")
include("solarsystem.jl")
include("main.jl")
# Thin C wrapper
include("LibNOVAS.jl")

end


text:
module DeduplicateFiles

using CRC

export DeduplicationFile, same_file, identical_files,
       delete_duplicate_file, indexfiles, deduplicate_files

"A structure with data about files to keep/delete."
type DeduplicationFile
  start::String
  realstart::String
  relpath::String
  realpath::String
  dirname::String
  dirinode::UInt64
  basename::String
  stat::Base.Filesystem.StatStruct
end

"A hash that is unique for the file"
file_id(x::DeduplicationFile) = hash(x.basename, hash(x.dirinode, hash(x.stat.device)))

"Determine if two paths are hard-linked to the same inode"
function is_hardlink(a,b)
  sa = stat(a)
  sb = stat(b)
  return sa.device == sb.device && sa.inode == sb.inode
end
is_hardlink(a::DeduplicationFile,b::DeduplicationFile) =
  a.stat.device == b.stat.device && a.stat.inode == b.stat.inode

"""
Determine if two paths point to the same directory entry.

This function returns `true` iff `a` and `b` both refer to the same directory
entry. This can happen if directories on either path are symbolic or hard
links to directories on the other path.

If `a` and `b` point to the same directory entry, then deleting `a` will also
delete `b`, so it would be bad to mistake `b` for a backup copy of `a`.

If `a` and/or `b` are symbolic links, then this function checks whether the
*links* are stored in the same directory entry. It does *not* check if they
point to the same or different files.

NOTE: This has *never been tested* on a file-system with hard-linked
directories. (Why would you use hard-linked directories?)
"""
function same_file(a::AbstractString, b::AbstractString)
  basename(a) == basename(b) || return false
  dira = dirname(abspath(a))
  dirb = dirname(abspath(b))
  realpath(dira) != realpath(dirb) || return true
  sa = stat(dira)
  sb = stat(dirb)
  return sa.device == sb.device && sa.inode == sb.inode
end

"Determine if two files are byte-for-byte identical"
function identical_files(a::IO, b::IO; buflen::Int=1048576)
  stat(a).size == stat(b).size || return false
  # TODO: Option to pass buffers as arguments and avoid repeated allocation.
  buf_a = Array{UInt8}(buflen)
  buf_b = Array{UInt8}(buflen)
  while (n_a = readbytes!(a, buf_a)) > 0
    n_b = readbytes!(b, buf_b)
    n_a == n_b || return false
    view(buf_a, 1:n_a) == view(buf_b, 1:n_b) || return false
  end
  return true
end
function identical_files(a::AbstractString, b::AbstractString; kwargs...)
  !same_file(a,b) || throw(ArgumentError("Cannot compare a file to itself."))
  result = false
  open(a, "r") do af
    open(b, "r") do bf
      result = identical_files(af, bf; kwargs...)
    end
  end
  return result
end

"""
   `delete_duplicate_file(; delete=file_to_delete, keep=file_to_keep)`

Delete a file if, and only if, it is an identical copy of another file.

The files to delete and keep are compared byte-for-byte. If they are
identical, then `file_to_delete` is deleted, and the function returns `true`.

If the files differ, or if both file paths refer to the same file (e.g. via
symbolic links) then no file is removed, and the function returns `false`.

The `file_to_delete` may not be a symbolic link. (The `file_to_keep` may be a
symbolic link, but if it points to the `file_to_delete`, then, obviously, that
file will not be deleted.)

This function cannot be called as `delete_duplicate_file("file1", "file2")`
because there might be confusion as to which file should be kept.

`delete_duplicate(..., dry_run=true)` does not delete files, but only returns
`true` or `false` depending on if a file would have been deleted.

`delete_duplicate(..., replace_with=:symlink)` replaces the deleted file with
a symbolic link to the `file_to_keep`.

`delete_duplicate(..., replace_with=:hardlink)` replaces the deleted file with
a hard link. If the files are not on the same file system, then an error is
thrown and no file is deleted. (Requires `ln` to be available via system call.)

`delete_duplicate(..., delete_hardlinks=true)` deletes the `file-to-delete`
even if it is a hardlink to `file-to-keep` (and therefore didn't occupy any
extra disk space).
"""
function delete_duplicate_file(; delete::AbstractString="", keep::AbstractString="",
                               dry_run::Bool=false, replace_with::Symbol=:nothing,
                               delete_hardlinks::Bool=false, verbose::Bool=false)
  (!isempty(delete) && !isempty(keep)) || throw(ArgumentError("Missing file name."))
  !islink(delete) || error("The file to delete must not be as symbolic link.")
  replace_with in [:nothing, :symlink, :hardlink] ||
      throw(ArgumentError("Illegal value for replace_with: ", string(replace_with)))

  # Follow symlinks now, so that we're not caught by a symlink that changes
  # between comparing the files and deleting.
  delete = realpath(delete)
  keep = realpath(keep)

  replace_with != :hardlink || stat(delete).device == stat(keep).device ||
      error("Cannot create hard-links across file systems.")

  !same_file(delete, keep) || return false

  hl = is_hardlink(delete, keep)

  hl || identical_files(delete, keep) || return false

  verbose && println(dry_run?"Would delete ":"Deleting ", delete, " as duplicate of ", keep, " (", stat(delete).size, " byters)")

  if !dry_run && (!hl || ( hl && delete_hardlinks && replace_with != :hardlink ))
    rm(delete)
    if replace_with == :symlink
      symlink(keep, delete)
    elseif replace_with == :hardlink
      run(`ln $keep $delete`)
    end
  end

  # If deleting the `delete` file somehow removed the `keep` file,
  # then there's a serious problem. Stop immediately!
  # (This should never happen.)
  isfile(keep) || error("Deleting ", delete, " somehow removed ", keep)

  return true
end

"""
   `indexfiles(start; idx, follow_symlinks=false)`

Create an index of all files found when walking the directory tree.
Returns an iterator over `DeduplicationFile` objects.
"""
function indexfiles(start::AbstractString; idx=Dict{UInt64,DeduplicationFile}(),
      follow_symlinks=false)
  realstart = realpath(start)
  for (root, dirs, files) in walkdir(realstart, follow_symlinks=follow_symlinks)
    ds = stat(root)
    for file in files
      ff = joinpath(root, file)
      if !islink(ff)
        df = DeduplicationFile(start, realstart, relpath(ff, realstart),
            realpath(ff), dirname(realpath(ff)), ds.inode, file, stat(ff))
        idx[file_id(df)] = df
      end
    end
  end
  return values(idx)
end

function indexfiles(starts::AbstractArray; idx=Dict{UInt64,DeduplicationFile}(),
    follow_symlinks=false)
  for start in starts
    indexfiles(start, idx=idx, follow_symlinks=follow_symlinks)
  end
  return values(idx)
end

"""
   `processdups(by, proc, list)`

Apply the `by` function to each item in `list`, for every set of two or more
items that generated identical `by` values, call `proc` on that set.
The `proc` function should return a list, and processdups returns a
concatenation of all returned lists.
"""
function processdups(by::Function, proc::Function, list::AbstractArray; quickpair::Bool=false)
  if quickpair && length(list)==2
    return(proc(list))
  end

  ret = Array{Any,1}[]

  length(list) >= 2 || return ret

  bylist = map(i->(i, by(list[i])), 1:length(list));

  # First sort the list...
  sort!(bylist, by=i->i[2])

  # ...then go through the list. Duplicates will be next to each other.
  ci = 1;
  cmp = bylist[ci][2]
  for i=1:length(bylist)
    if i==length(bylist) || bylist[i+1][2] != cmp
      if i>ci
        sublist = map(j->list[bylist[j][1]], ci:i)
        ret = vcat(ret, proc(sublist))
      end
      ci = i+1;
      if i<length(bylist)
        cmp = bylist[ci][2];
      end
    end
  end
  return ret
end

"""
   `deldups(list, dfun; dry_run)`

Determine which files from a list of suspected duplicates that should be
deleted according to the decision function `dfun`, and delete them if they
actually are duplicates.
"""
function deldups(list, dfun; kwargs...)
  deleted = Array{Tuple{DeduplicationFile, DeduplicationFile}, 1}();
  sort!(list, lt=dfun)
  lmax = length(list)
  i = 1
  while i<lmax
    for j=length(list):-1:i+1
      if dfun(list[i], list[j]) && delete_duplicate_file(delete=list[i].realpath,
             keep=list[j].realpath; kwargs...)
        push!(deleted, (list[i], list[j]))
        lmax = j
        break
      end
    end
    i = i+1
  end
  return deleted
end

const crc32 = crc(CRC_32)

"Run crc32 on the first `buflen` bytes of a file"
function partialcrc32(a::IO; buflen::Int=1048576)
  buf_a = Array{UInt8}(buflen)
  n_a = readbytes!(a, buf_a)
  if n_a != length(buf_a)
    buf_a = buf_a[1:n_a]
  end
  return crc32(buf_a)
end

"Divide list by crc, and run deldups on each sub-list"
crcdeldups(list, dfun; kwargs...) =
  processdups(x->open(crc32, x.realpath), lst->deldups(lst, dfun; kwargs...), list, quickpair=true)

"Divide list by partial crc, and run crcdeldups on each sub-list"
pcrcdeldups(list, dfun; kwargs...) =
  processdups(x->open(partialcrc32, x.realpath), lst->crcdeldups(lst, dfun; kwargs...), list)

"Divide list by size, then run pcrcdeldups/crcdeldups on each sub-sub-list"
sizcrcdeldups(list, dfun; kwargs...) =
  processdups(x->x.stat.size, lst->lst[1].stat.size>1048576?pcrcdeldups(lst, dfun; kwargs...):crcdeldups(lst, dfun; kwargs...), list)

"""
   `list = deduplicate_files(startdirs, dfun)`

Scan `startdirs` recursively for files, call dfun on pairs of suspected
duplicates, and delete identical copies when dfun returns true.
"""
function deduplicate_files(startdirs, dfun; verbose::Bool=false, kwargs...)
  list = indexfiles(startdirs)
  if verbose
    println("Indexed ", length(list), " files.")
  end
  push!(kwargs, (:verbose, verbose))
  sizcrcdeldups(collect(list), dfun; kwargs...)
end

# Define a `show` method for tuples of deduplication files, to make it easier to
# read the returned array from `deduplicate_files`.
import Base.show
show(io::IO, t::Tuple{DeduplicationFile, DeduplicationFile}) = print(io,
  joinpath(t[1].start, t[1].relpath), " duplicate of ", joinpath(t[2].start, t[2].relpath))

"The total size (in bytes) of an array of `DeduplicationFile`s"
totalsize(x::DeduplicationFile) = x.stat.size
totalsize(x::Tuple{DeduplicationFile, DeduplicationFile}) = totalsize(x[1])
totalsize(x::AbstractArray) = sum(totalsize.(x))

# TODO: Check with dfun that at least one duplicate might be deleted
#       before doing CRC.

end # module


text:
using CSV, DataFrames, Feather, Gadfly, IJulia, JSON, RDatasets, StatsBase

# data = readdlm("iris_dataset.csv", ',')

iris = CSV.read("iris_dataset.csv", header=true, delim=',', copycols = true)

print(iris)

print(describe(iris))

# print(describe(iris, stats=[:q25, :q75, :first, :last]))

# print(describe(iris, stats=:all))

# Pkg.add("Statistics")

using Statistics

print(mean(iris[:sepallength]))

for x in names(iris)[1:end-1]
    for y in names(iris)[1:end-1]
        println("$x /t $y /t $(cor(iris[x], iris[y]))")
    end
    println("--------------------------------------------------")
end

# columns to array
sepallength = Array(iris[:sepallength])
# print(sepallength)

# DataFrame to matrix
irisarr = convert(Array, iris[:,:])

iris[1:3, [:petallength, :petalwidth]]

# CSV.write("test_iris.csv", iris)


text:
@testset "Kriging" begin
    d, n = 3, 10
    x = 2π * rand(d, n)
    y1 = [sum(sin, view(x, :, i)) / d for i in 1:n]
    y2 = [sum(cos, view(x, :, i)) / d for i in 1:n]
    y = vcat(transpose(y1), transpose(y2))
    mZero = MeanZero()
    kern = SE(0.0, 0.0)
    ntest = 5
    xtest = randn(d, ntest)
    yt1 = [sum(sin, view(xtest, :, i)) / d for i in 1:ntest]
    yt2 = [sum(cos, view(xtest, :, i)) / d for i in 1:ntest]
    ytest = vcat(transpose(yt1), transpose(yt2))
    krg = Kriging(x, y, mZero, kern)
    @testset "predict samples" begin
        y_pred, sig2 = predict_full(krg, x)
        @test maximum(abs, y - y_pred) ≈ 0.0 atol = 0.1
    end
    @testset "predict tests" begin
        y_pred, sig2 = predict_full(krg, xtest)
    end
    @testset "tune local" begin
        tune!(krg, ([-10.0, -10.0], [10.0, 10.0]), nothing; verbose = true)
    end
    @testset "tune" begin
        tune!(
            krg,
            ([-10.0, -10.0], [10.0, 10.0]),
            nothing;
            global_eval = 10000,
            verbose = true,
        )
    end
    @testset "append" begin
        append!(krg, xtest, ytest)
    end
    @testset "get" begin
        x1 = getx(krg, 1)
        y1 = gety(krg, 1)
        @test x1 ≈ x[:, 1] atol = 1e-6
        @test y1 ≈ y[:, 1] atol = 1e-6
    end
    @testset "get all" begin
        x1 = getx(krg)
        y1 = gety(krg)
        @test x1 ≈ hcat(x,xtest) atol = 1e-6
        @test y1 ≈ hcat(y,ytest) atol = 1e-6
    end
end


text:
module FastLevenshtein

    function fastlevenshtein(v1::Vector, v2::Vector, insert_cost::Real=1.0, delete_cost::Real=1.0, update_cost::Real=1.0)::Real
        """
        Efficiently Compute the Levenshtein (edit) distance between two vectors using dynamic programming.
        :param v1: The first vector is a Mx1 vector.
        :param v2: The second vector is a Nx1 vector.
        :return: The minimum number of edits required to make the vectors the same.
        """
        M::Int64, N::Int64 = length(v1), length(v2)

        # This table holds the subproblems we have to solve to get an answer
        table::Array{Real} = zeros(Real, (M+1, N+1))

        for i::Int64 ∈ 0:M
            for j::Int64 ∈ 0:N
                # if the length of the 1st vector is 0 then we need at least j edits
                if i == 0
                    table[i+1, j+1] = convert(AbstractFloat, j)
                # if the length of the 2nd vector is 0 then we need at least i edits
                elseif j == 0
                    table[i+1, j+1] = convert(AbstractFloat, i)
                # if the values are the same then the minumum number of edits hasn't changed
                elseif v1[i] == v2[j]
                    table[i+1, j+1] = table[i, j]
                # if the value are different then the minimum edits has increased by 1
                else
                    table[i+1, j+1] = min(
                        table[i+1, j] + insert_cost, # The minimum edits required if i > j (insert value @ current index)
                        table[i, j+1] + delete_cost, # The minimum edits required if i < j (remove value @ current index)
                        table[i, j]   + update_cost  # The minimum edits required if i = j (update value @ current index)
                    )
                end
            end
        end

        return table[M+1, N+1]
    end

    function fastlevenshtein(s1::AbstractString, s2::AbstractString, insert_cost::Real=1.0, delete_cost::Real=1.0, update_cost::Real=1.0)::Real
        """
        Efficiently Compute the Levenshtein (edit) distance between two strings using dynamic programming.
        :param s1: The first string contains M characters.
        :param s2: The second string contains N characters.
        :return: The minimum number of edits required to make the strings the same.
        """
        return fastlevenshtein(Vector{Char}(s1), Vector{Char}(s2), insert_cost, delete_cost, update_cost)
    end

    export fastlevenshtein
end


text:
export  fejer,
        fejer2,
        clenshaw_curtis,
        quadgp,
        golubwelsch,
        gauss,
        radau,
        lobatto
"""
    fejer(N::Int)
Fejer's first quadrature rule.
"""
function fejer(N::Int)
    @assert N >= 1 "N has to be positive"
    N == 1 && return zeros(1), [2.]
    θ = map(x->(2x-1)*pi/(2N),1:N)
    M = N ÷ 2
    return cos.(θ), map(x->2/N*(1-2*sum(cos(2*n*x)/(4*n^2-1) for n=1:M)),θ)
end

"""
    fejer2(n::Int)
Fejer's second quadrature rule according to [Waldvogel, J. Bit Numer Math (2006) 46: 195](https://doi.org/10.1007/s10543-006-0045-4).
"""
function fejer2(n::Int)
    @assert n >= 2
    N = 1:2:n-1
    m = n - length(N)
    v0 = push!(push!(map(x->2/x/(x-2.),N),1/last(N)), zeros(m)...)
    @inbounds v2 = -v0[1:n] - v0[n+1:-1:2];
    wf2 = real(ifft(v2))
    return map(x->cos(x*pi/n),0:n), push!(wf2, first(wf2))
end

"""
    clenshaw_curtis(n::Int)
Clenshaw-Curtis quadrature according to [Waldvogel, J. Bit Numer Math (2006) 46: 195](https://doi.org/10.1007/s10543-006-0045-4).
"""
function clenshaw_curtis(n::Int)
    @assert n >= 2
    N = 1:2:n-1
    l = length(N)
    m = n - l
    v0 = push!(push!(map(x->2/x/(x-2.),N),1/last(N)), zeros(m)...)
    @inbounds v2 = -v0[1:n] - v0[n+1:-1:2]
    g0= -ones(Float64,n)
    g0[1+l] += n
    g0[1+m] += n
    g = 1 / (n^2 - 1 + mod(n,2) ) * g0
    wcc = real(ifft(v2 + g));
    return map(i->cos(i/n*pi),0:n), push!(wcc,first(wcc))
end

"""
    quadgp(weight::Function,lb::Real,ub::Real,N::Int=10;quadrature::Function=clenshaw_curtis,bnd::Float64=Inf)
general purpose quadrature based on Gautschi, "Orthogonal Polynomials: Computation and Approximation", Section 2.2.2, pp. 93-95

Compute the `N`-point quadrature rule for `weight` with support (`lb`, `ub`).
The quadrature rule can be specified by the keyword `quadrature`.
The keyword `bnd` sets the numerical value for infinity.
"""
function quadgp(weight::Function,lb::Real,ub::Real,N::Int=10;quadrature::Function=clenshaw_curtis,bnd::Float64=Inf)
    @assert lb < ub "inconsistent interval bounds"
    t_fej,w_fej = quadrature(N)
    t, w, d = zeros(Float64,N),zeros(Float64,N),zeros(Float64,N)
    # transformation
    if (-bnd<lb<bnd) && (-bnd<ub<bnd)
        t = 0.5*((ub-lb)*t_fej .+ (ub+lb))
        d = 0.5*(ub-lb)*ones(Float64,size(t_fej))
    elseif -bnd>=lb && (-bnd<ub<bnd)
        t = ub .- (1 .- t_fej)./(1 .+ t_fej)
        d = 2 ./(1 .+ t_fej).^2
    elseif (-bnd<lb<bnd) && ub>=bnd
        t = lb .+ (1 .+ t_fej)./(1 .- t_fej)
        d = 2 ./(1 .- t_fej).^2
    elseif -bnd>=lb && ub>=bnd
        t = t_fej./(1 .- t_fej.^2)
        d = (1 .+ t_fej.^2)./( 1 .- t_fej.^2).^2
    end
    w = w_fej.*weight.(t).*d
    return t, w
end

####################################################
####################################################
####################################################
# Quadrature rules based on recurrence coefficients

function golubwelsch(α::AbstractVector{<:Real}, β::AbstractVector{<:Real}, maxiter::Int=30)
    N = length(α) - 1
    a, β0 = copy(α[1:N]), β[1]
    w = zero(a)
    special_eigenproblem!(a,copy(sqrt.(β)),w,maxiter)
    idx = sortperm(a)
    a[idx], (β0*w.^2)[idx]
end

golubwelsch(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly}) = golubwelsch(op.α,op.β)

"""
    gauss(N::Int,α::AbstractVector{<:Real},β::AbstractVector{<:Real})
    gauss(α::AbstractVector{<:Real},β::AbstractVector{<:Real})
    gauss(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly})
    gauss(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly})
Gauss quadrature rule, also known as Golub-Welsch algorithm

`gauss()` generates the `N` Gauss quadrature nodes and weights for a given weight function.
The weight function is represented by the `N` recurrence coefficients for the monic polynomials orthogonal
with respect to the weight function.

!!! note
    The function `gauss` accepts at most `N = length(α) - 1` quadrature points,
        hence providing at most an `(length(α) - 1)`-point quadrature rule.

!!! note
    If no `N` is provided, then `N = length(α) - 1`.
"""
function gauss(N::Int,α::AbstractVector{<:Real},β::AbstractVector{<:Real})
    N += 1
    @assert N > 0 "only positive N allowed"
    @assert length(α) == length(β) "inconsistent number of recurrence coefficients"
    N0 = length(α)
    @assert N0 >= N  "not enough recurrence coefficients"
    @inbounds golubwelsch(α[1:N],β[1:N])
end

gauss(α::AbstractVector{<:Real},β::AbstractVector{<:Real}) = gauss(length(α)-1,α,β)

gauss(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly}) = gauss(N::Int,op.α,op.β)

gauss(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly}) = gauss(op.α,op.β)

"""
    radau(N::Int,α::AbstractVector{<:Real},β::AbstractVector{<:Real},end0::Real)
    radau(α::AbstractVector{<:Real},β::AbstractVector{<:Real},end0::Real)
    radau(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},end0::Real)
    radau(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},end0::Real)
Gauss-Radau quadrature rule.
Given a weight function encoded by the recurrence coefficients `(α,β)`for the associated
orthogonal polynomials, the function generates the
nodes and weights `(N+1)`-point Gauss-Radau
quadrature rule for the weight function having a prescribed
node `end0` (typically at one of the end points of the support
interval of w, or outside thereof).

!!! note
    The function `radau` accepts at most `N = length(α) - 2` as an input,
        hence providing at most an `(length(α) - 1)`-point quadrature rule.

!!! note
    Reference: OPQ: A MATLAB SUITE OF PROGRAMS FOR GENERATING ORTHOGONAL POLYNOMIALS AND RELATED QUADRATURE RULES by Walter Gautschi
"""
function radau(N::Int,α::AbstractVector{<:Real},β::AbstractVector{<:Real},end0::Real)
    α_ = copy(α)
    @assert N > 0 "only positive N allowed"
    @assert length(α_) == length(β) > 0 "inconsistent number of recurrence coefficients"
    @assert length(α_) >= N + 1 "not enough recurrence coefficients"
    p0 = 0.
    p1 = 1.
    for n in Base.OneTo(N)
      pm1 = p0
      p0 = p1;
      @inbounds p1 = (end0 - α_[n])*p0 - β[n]*pm1;
    end
    @inbounds α_[N+1] = end0 - β[N+1]*p0/p1
    gauss(N+1 , α_, β)
end
radau(α::AbstractVector{<:Real},β::AbstractVector{<:Real},end0::Real) = radau(length(α)-2,α,β,end0)
radau(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},end0::Real) = radau(N,op.α,op.β,end0)
radau(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},end0::Real) = radau(op.α,op.β,end0::Real)

"""
    lobatto(N::Int,α::AbstractVector{<:Real},β::AbstractVector{<:Real},endl::Real,endr::Real)
    lobatto(α::AbstractVector{<:Real},β::AbstractVector{<:Real},endl::Real,endr::Real)
    lobatto(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},endl::Real,endr::Real)
    lobatto(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},endl::Real,endr::Real)
Gauss-Lobatto quadrature rule.
Given a weight function encoded by the recurrence coefficients for the associated
orthogonal polynomials, the function generates
the nodes and weights of the `(N+2)`-point Gauss-Lobatto
quadrature rule for the weight function, having two
prescribed nodes `endl`, `endr` (typically the left and right end
points of the support interval, or points to the left
resp. to the right thereof).

!!! note
    The function `radau` accepts at most `N = length(α) - 3` as an input,
    hence providing at most an `(length(α) - 1)`-point quadrature rule.

!!! note
    Reference: OPQ: A MATLAB SUITE OF PROGRAMS FOR GENERATING ORTHOGONAL POLYNOMIALS AND RELATED QUADRATURE RULES by Walter Gautschi
"""
function lobatto(N::Int,α_::AbstractVector{<:Real},β_::AbstractVector{<:Real},endl::Real,endr::Real)
    α, β = copy(α_), copy(β_)
    @assert N > 0 "only positive N allowed"
    @assert length(α) == length(β) > 0 "inconsistent number of recurrence coefficients"
    @assert endl < endr "inconsistent end points"
    @assert length(α) >= N+2 "not enough recurrence coefficients"
    p0l = 0.; p0r = 0.; p1l = 1.; p1r = 1.;
    for n=1:N+1
        pm1l=p0l
        p0l=p1l
        pm1r=p0r
        p0r=p1r
        @inbounds p1l=(endl-α[n])*p0l-β[n]*pm1l
        @inbounds p1r=(endr-α[n])*p0r-β[n]*pm1r
    end
    det = p1l * p0r - p1r * p0l
    @inbounds α[N+2] = (endl * p1l * p0r - endr * p1r * p0l) / det
    @inbounds β[N+2] = (endr - endl) * p1l * p1r / det;
    gauss(N+2,α,β)
end
lobatto(α::AbstractVector{<:Real},β::AbstractVector{<:Real},endl::Real,endr::Real) = lobatto(length(α)-3,α,β,endl,endr)
lobatto(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},endl::Real,endr::Real) = lobatto(N,op.α,op.β,endl,endr)
lobatto(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},endl::Real,endr::Real) = lobatto(op.α,op.β,endl,endr)


text:
export rule

@rule BIFMHelper(:in, Marginalisation) (m_out::Any, ) = m_out


text:
using Pkg
using CUDA
using ScatterNNlib
using DataFrames
using CSV
using BenchmarkTools
using BenchmarkTools: Trial, TrialEstimate, median, mean

Pkg.status("ScatterNNlib")
Pkg.status("CUDA")

ENV["JULIA_NUM_THREADS"] = 1

d = 50
nbins = 20
getinfo(te::TrialEstimate) = te.time, te.gctime, te.memory
getstats(t::Trial) = [getinfo(minimum(t)), getinfo(mean(t)), getinfo(maximum(t))]

metadata = DataFrame(device=String[], dim=Int[], sample=Int[], bins=Int[])
mintime = DataFrame(min_time=Float64[], min_gc=Float64[], min_mem=Int[])
meantime = DataFrame(mean_time=Float64[], mean_gc=Float64[], mean_mem=Int[])
maxtime = DataFrame(max_time=Float64[], max_gc=Float64[], max_mem=Int[])

for l = [2^5, 2^10, 2^15, 2^20]
    hist = zeros(Float32, d, nbins)
    δ = rand(Float32, d, l)
    idx = rand(1:nbins, l)

    hist_gpu = CuArray(hist)
    δ_gpu = CuArray(δ)
    idx_gpu = CuArray(idx)

    scatter_add!(hist, δ, idx)
    scatter_add!(hist_gpu, δ_gpu, idx_gpu)

    b_cpu = @benchmark scatter_add!($hist, $δ, $idx)
    b_gpu = @benchmark scatter_add!($hist_gpu, $δ_gpu, $idx_gpu)
    s_cpu = getstats(b_cpu)
    s_gpu = getstats(b_gpu)

    push!(metadata, ("cpu", d, l, nbins))
    push!(mintime, s_cpu[1])
    push!(meantime, s_cpu[2])
    push!(maxtime, s_cpu[3])

    push!(metadata, ("gpu", d, l, nbins))
    push!(mintime, s_gpu[1])
    push!(meantime, s_gpu[2])
    push!(maxtime, s_gpu[3])
end

data = hcat(metadata, mintime, meantime, maxtime)
CSV.write("benchmark/scatter_julia.tsv", data; delim="\t")

## Benchmark
# @benchmark scatter_add!($hist, $δ, $idx)
# CuArrays.@time scatter_add!(hist_gpu, δ_gpu, idx_gpu)


text:
# ```@setup slurmgen
# # all examples are run on an example global.toml and downsized maps.
# ARGS = ["example.toml"]
# ```

configfile = ARGS[1]

# #
# # [SLURM Commands for Spectra (slurmgen.jl)](@id slurmgen)
# This command generates SLURM commands that executes [rawspectra.jl](@ref rawspectra)
# on all the pairs of maps in the config.
## this file just prints out the SLURM commands required to compute the spectra
using TOML
config = TOML.parsefile(configfile)

# Let's generate the commands we need for likelihood spectra.
mapids = [k for k in keys(config["map"])]
cmd = "sbatch scripts/8core2hr.cmd"
for i in 1:length(mapids)
    for j in i:length(mapids)
        println("$(cmd) \"julia src/rawspectra.jl $(configfile) $(mapids[i]) $(mapids[j])\"")
    end
end

# I just paste these in and wait a few hours. The resulting spectra are deposited in the
# `scratch` dir in the config.

# ## Fit Noise Model
# These SLURM commands fit the raw spectra with the camspec noise model.
## loop over freqs and noise channels
println()
for freq in ("100", "143", "217")
    for spec in ("TT", "EE")
        println("$(cmd) \"julia src/fitnoisemodel.jl $(configfile) $(freq) $(spec)\"")
    end
end

# ## Generating Signal-Only Simulations
# These SLURM commands correct the covariances for the insufficiently-apodized point source
# holes.
## loop over freqs and noise channels
println()
cmd = "sbatch scripts/8core2hr.cmd"
freqs = ("100", "143", "217")
for i in 1:3
    for j in i:3
        for spec in ("TT", "TE", "ET", "EE")
            freq1, freq2 = freqs[i], freqs[j]
            println("$(cmd) \"julia src/signalsim.jl $(configfile) $(freq1) $(freq2) $(spec) 5000\"")
        end
    end
end

# ## Fit Signal-Only Simulation Corrections
# These SLURM commands correct the covariances for the insufficiently-apodized point source
# holes.
## loop over freqs and noise channels
println()
freqs = ("100", "143", "217")
for i in 1:3
    for j in i:3
        for spec in ("TT", "TE", "ET", "EE")
            freq1, freq2 = freqs[i], freqs[j]
            println("$(cmd) \"julia src/corrections.jl $(configfile) $(freq1) $(freq2) $(spec)\"")
        end
    end
end

##
function unique_splits(X, Y, f1, f2)
    if (X == Y) && (f1 == f2)
        return [(1,2)]
    elseif (X==Y) && (f1 != f2)
        return [(1,2), (2,1)]
    else
        return [(1,1), (1,2), (2,1), (2,2)]
    end
end

using Test
using Base.Iterators

@test unique_splits("T", "T", "100", "100") == [(1,2)]
@test unique_splits("T", "E", "100", "100") == [(1,1), (1,2), (2,1), (2,2)]
@test unique_splits("T", "T", "100", "143") == [(1,2), (2,1)]
@test unique_splits("E", "E", "100", "143") == [(1,2), (2,1)]
@test unique_splits("T", "E", "100", "143") == [(1,1), (1,2), (2,1), (2,2)]

include("util.jl")
specs = (
    (:TT,"100","100"), (:TT,"143","143"), (:TT,"143","217"), (:TT,"217","217"),
    (:EE,"100","100"), (:EE,"100","143"), (:EE,"100","217"), (:EE,"143","143"),
    (:EE,"143","217"), (:EE,"217","217"),
    (:TE,"100","100"), (:TE,"100","143"), (:TE,"100","217"), (:TE,"143","143"),
    (:TE,"143","217"), (:TE,"217","217"),
    (:ET,"100","143"), (:ET,"100","217"), (:ET,"143","217")  # combined with TE in plic
)

constituents = []

for spec in specs
    AB, f1, f2 = spec
    A, B = string(AB)
    for (s1, s2) in unique_splits(A, B, f1, f2)
        push!(constituents, (A, B, f1, f2, s1, s2))
    end
end

run_name = config["general"]["name"]
cmd = "sbatch scripts/8core2hr.cmd"
scriptpath = joinpath(config["scratch"], "scripts")
mkpath(scriptpath)
open(joinpath(scriptpath, "$(run_name)_gen_covmats.sh"), "w") do f
    nspecs = length(constituents)

    ## generate the upper triangle
    for i1 in 1:nspecs
        for i2 in i1:nspecs
            A, B, f1, f2, s1, s2 = constituents[i1]
            C, D, f3, f4, s3, s4 = constituents[i2]
            write(f, "$(cmd) \"julia src/covmat.jl $(configfile) $f1 $f2 $f3 $f4 $A $B $C $D $s1 $s2 $s3 $s4\"\n")
        end
    end
end

max_forks_repo_name: sgiardie/PSpipe


text:
using PyPlot

vars=[]
for i in 1:length(emu_out.population)
    m, v = gp_regression(emu_out.population[i], emu_out.emulators[i])
    push!(vars, v)
end

ion()
for i in 1:length(emu_out.threshold_schedule)
    epsilon = emu_out.threshold_schedule[i]
    figure()
    scatter3D(emu_out.population[i][:, 3], emu_out.population[i][:, 2], vars[i])
    title("Variance #$(i); ϵ=$(epsilon)")
    zlabel("Variance")
    xlabel("Parameter 3")
    ylabel("Parameter 2")
    figure(figsize=(13,4.8))
    subplot(121)
    scatter(emu_out.population[i][:, 3], emu_out.population[i][:, 2], marker=".")
    xlabel("Parameter 3")
    ylabel("Parameter 2")
    title("Accepted population #$(i); ϵ=$(epsilon)")
    subplot(122)
    scatter(emu_out.emulators[i].gp_training_x[:, 3], emu_out.emulators[i].gp_training_x[:, 2], marker="x")
    xlabel("Parameter 3")
    ylabel("Parameter 2")
    title("Emulator #$(i) training data")
end


text:
#   Adapted from
#   https://github.com/dabeaz/bitey/blob/master/test/test_M.py
#
#   Copyright (C) 2012
#   David M. Beazley (Dabeaz LLC), http://www.dabeaz.com
#   BSD 3-clause "New" or "Revised" License

using Base.Test
using BitcodeRunner

M = bitcode_library("ctest.bc")

@test M.add_char(4,5) == Int8(9)
@test M.add_short(4,5) == Int16(9)
@test M.add_long(4,5) == Int32(9)
@test M.add_longlong(4,5) == Int64(9)

@test M.add_float(4,5) == 9f0
@test M.add_double(4,5) == 9.0

# Test passing Ref's and Arrays
a = Ref(Int16(2))
M.mutate_short(a)
@test a[] == Int16(4)
a = Ref(Int32(2))
M.mutate_long(a)
@test a[] == Int32(4)
a = Ref(Int64(2))
M.mutate_longlong(a)
@test a[] == Int64(4)
a = Ref(2f0)
M.mutate_float(a)
@test a[] == 4f0
a = Ref(2.0)
M.mutate_double(a)
@test a[] == 4.0
a = [2.0]
M.mutate_double(a)
@test a[1] == 4.0

# Array tests
@test M.arr_sum_int(collect(UInt32(1):UInt32(4))) == Int32(10)
@test M.arr_sum_double(collect(1.0:4.0)) == 10.0

# Struct test
p = Ref(M.Point(3, 4))
q = Ref(M.Point(6, 8))
@test M.distance(p, q) == 5.0


text:
for N in [Float64, Rational{Int}, Float32]
    B1 = BallInf(N[0, 0], N(1))
    B2 = BallInf(N[4, -4], N(1))
    B3 = BallInf(N[1, -1], N(1))
    C = Complement(B1)

    # double-complement is the identity
    @test Complement(C) == B1

    # dimension
    @test dim(C) == 2

    # membership
    @test N[1, 1] ∉ C && N[2, 2] ∈ C

    # emptiness
    @test !isempty(C) && isempty(Complement(Universe(2)))

    # inclusion
    subset, point = ⊆(B2, C, true)
    @test B2 ⊆ C && subset && point == N[]
    for X in [B1, B3]
        subset, point = ⊆(X, C, true)
        @test X ⊈ C && !subset && point ∈ X && point ∉ C
    end

    # isdisjoint
    for X in [B2, B3]
        res, w = isdisjoint(X, C, true)
        @test !isdisjoint(X, C) && !res && w ∈ X && w ∈ C
        res, w = isdisjoint(C, X, true)
        @test !isdisjoint(C, X) && !res && w ∈ X && w ∈ C
    end
    res, w = isdisjoint(B1, C, true)
    @test isdisjoint(B1, C) && res && w == N[]
    res, w = isdisjoint(C, B1, true)
    @test isdisjoint(C, B1) && res && w == N[]

    # test convexity from the type
    @test isconvextype(typeof(Complement(Universe{N}(2))))
    @test isconvextype(typeof(Complement(EmptySet{N}(2))))
    @test isconvextype(typeof(Complement(HalfSpace(N[1], N(0)))))

    # concrete complement
    H = HalfSpace(N[1, 0], N(1)) # x <= 1
    @test complement(H) == HalfSpace(N[-1, 0], N(-1))
    @test constraints_list(Complement(H)) == [HalfSpace(N[-1, 0], N(-1))]
    # test fallback for <: LazySet
    c = complement(N[0 1; 1 0] * H)
    @test c isa UnionSetArray && length(array(c)) == 1
    @test first(array(c)) == HalfSpace(N[0, -1], N(-1)) # complement of y <= 1 is y >= 1

    # boundedness
    @test isboundedtype(typeof(Complement(Universe{N}(2))))
    @test !isboundedtype(typeof(Complement(EmptySet{N}(2))))
end


text:
using GeometryTypes, Colors, GLVisualize, GLWindow, FileIO, GLAbstraction, ModernGL

function read_ifs(filename)
    open(filename) do io
        function str()
            n = read(io, UInt32)
            String(read(io, UInt8, n))
        end
        ifs = str()
        zero = read(io, Float32)
        name = str()
        vertices = str()
        nverts = read(io, UInt32)
        verts = Vector{Point3f0}(nverts)
        for i = 1:nverts
            verts[i] = Point3f0(
                read(io, Float32),
                read(io, Float32),
                read(io, Float32)
            )
        end
        tris = str()
        nfaces = read(io, UInt32)
        faces = Vector{GLTriangle}(nfaces)
        for i = 1:nfaces
            faces[i] = GLTriangle(
                read(io, UInt32),
                read(io, UInt32),
                read(io, UInt32)
            )
        end
        (verts, faces)
    end
end

window = glscreen()

function display_data(folder)
    meshpaths = filter(x-> endswith(x, ".ifs"), readdir(folder))[1:1024]
    v0 = (Point3f0(typemax(Float32)), Point3f0(typemin(Float32)))
    faces = GLTriangle[]
    verts = Point3f0[]
    norms = Normal{3, Float32}[]
    lastvidx = 0
    for (i, meshpath) in enumerate(meshpaths)
        vs, fs = read_ifs(joinpath(folder, meshpath))
        mini, maxi = extrema(vs)
        ns = normals(vs, fs)
        x, y = ind2sub((32, 32), i)
        s = maximum(maxi .- mini)
        append!(verts, vs .+ Point3f0(x, y, 0f0))
        append!(norms, ns)
        append!(faces, fs .+ lastvidx)
        lastvidx += length(fs) * 3
    end
    _view(visualize(
        GLNormalMesh(faces = faces, vertices = verts, normals = norms),
        color = RGBA{Float32}(rand(RGB{Float32}), 1f0)
    ))
end
display_data(homedir() * "/3dstuff/models")

function rendloop(window, frame_times)
    while isopen(window)
        tic()
        glFinish()
        render_frame(window)
        GLWindow.swapbuffers(window)
        GLWindow.poll_glfw()
        GLWindow.poll_reactive()
        push!(frame_times, toq())
    end
    GLFW.DestroyWindow(window)
    window.handle = C_NULL
    frame_times
end
timsies = Float64[]
rendloop(window, timsies)

open("/home/s/.julia/v0.6/Visualize/test2.csv", "w") do io
    for t in timsies
        println(io, t)
    end
end
println(pwd())


text:
@testset "memory" begin

let
    a,b = Mem.info()
    # NOTE: actually testing this is pretty fragile on CI
    #=@test a == =# CUDAdrv.available_memory()
    #=@test b == =# CUDAdrv.total_memory()
end

# dummy data
T = UInt32
N = 5
data = rand(T, N)
nb = sizeof(data)

# buffers are untyped, so we use a convenience function to get a typed pointer
# we prefer to return a device pointer (for managed buffers) to maximize CUDAdrv coverage
typed_pointer(buf::Union{Mem.Device, Mem.Unified}, T) = convert(CuPtr{T}, buf)
typed_pointer(buf::Mem.Host, T)                       = convert(Ptr{T},   buf)

# allocations and copies
for srcTy in [Mem.Device, Mem.Host, Mem.Unified],
    dstTy in [Mem.Device, Mem.Host, Mem.Unified]

    dummy = Mem.alloc(srcTy, 0)
    Mem.free(dummy)

    src = Mem.alloc(srcTy, nb)
    unsafe_copyto!(typed_pointer(src, T), pointer(data), N)

    dst = Mem.alloc(dstTy, nb)
    unsafe_copyto!(typed_pointer(dst, T), typed_pointer(src, T), N)

    ref = Array{T}(undef, N)
    unsafe_copyto!(pointer(ref), typed_pointer(dst, T), N)

    @test data == ref

    if isa(src, Mem.Device) || isa(src, Mem.Unified)
        Mem.set!(typed_pointer(src, T), zero(T), N)
    end

    Mem.free(src)
    Mem.free(dst)
end

# asynchronous operations
let
    src = Mem.alloc(Mem.Device, nb)

    @test_throws ArgumentError unsafe_copyto!(typed_pointer(src, T), pointer(data), N; async=true)
    unsafe_copyto!(typed_pointer(src, T), pointer(data), N; async=true, stream=CuDefaultStream())

    Mem.set!(typed_pointer(src, T), zero(T), N; async=true, stream=CuDefaultStream())

    Mem.free(src)
end

# pinned memory
let
    # can only get GPU pointer if the pinned buffer is mapped
    src = Mem.alloc(Mem.Host, nb)
    @test_throws ArgumentError convert(CuPtr{T}, src)
    Mem.free(src)

    # create a pinned and mapped buffer
    src = Mem.alloc(Mem.Host, nb, Mem.HOSTALLOC_DEVICEMAP)

    # get the CPU address and copy some data
    cpu_ptr = convert(Ptr{T}, src)
    unsafe_copyto!(cpu_ptr, pointer(data), N)

    # get the GPU address and construct a fake device buffer
    gpu_ptr = convert(CuPtr{Cvoid}, src)
    gpu_obj = Mem.alloc(Mem.Device, nb)
    dst = similar(gpu_obj, gpu_ptr)
    Mem.free(gpu_obj)

    # copy data back from the GPU and compare
    ref = Array{T}(undef, N)
    unsafe_copyto!(pointer(ref), typed_pointer(dst, T), N)
    @test ref == data

    Mem.free(src)
    # NOTE: don't free dst, it's just a mapped pointer
end

# pinned memory with existing memory
if attribute(dev, CUDAdrv.DEVICE_ATTRIBUTE_HOST_REGISTER_SUPPORTED) != 0
    # can only get GPU pointer if the pinned buffer is mapped
    src = Mem.register(Mem.Host, pointer(data), nb)
    @test_throws ArgumentError convert(CuPtr{T}, src)
    Mem.unregister(src)

    # register a pinned and mapped buffer
    src = Mem.register(Mem.Host, pointer(data), nb, Mem.HOSTREGISTER_DEVICEMAP)

    # get the GPU address and construct a fake device buffer
    gpu_ptr = convert(CuPtr{Cvoid}, src)
    gpu_obj = Mem.alloc(Mem.Device, nb)
    dst = similar(gpu_obj, gpu_ptr)
    Mem.free(gpu_obj)

    # copy data back from the GPU and compare
    ref = Array{T}(undef, N)
    unsafe_copyto!(pointer(ref), typed_pointer(dst, T), N)
    @test ref == data

    Mem.unregister(src)
    # NOTE: don't unregister dst, it's just a mapped pointer
end

# unified memory
let
    src = Mem.alloc(Mem.Unified, nb)

    #Mem.prefetch(src, nb; device=CUDAdrv.DEVICE_CPU)
    Mem.advise(src, Mem.ADVISE_SET_READ_MOSTLY)

    # get the CPU address and copy some data
    cpu_ptr = convert(Ptr{T}, src)
    unsafe_copyto!(cpu_ptr, pointer(data), N)

    # get the GPU address and construct a fake device buffer
    gpu_ptr = convert(CuPtr{Cvoid}, src)
    gpu_obj = Mem.alloc(Mem.Device, nb)
    dst = similar(gpu_obj, gpu_ptr)
    Mem.free(gpu_obj)

    # copy data back from the GPU and compare
    ref = Array{T}(undef, N)
    unsafe_copyto!(pointer(ref), typed_pointer(dst, T), N)
    @test ref == data

    Mem.free(src)
end

end


text:
export push, pop, shift, unshift, insert, deleteat

@generated function push(vec::StaticVector, x)
    newtype = similar_type(vec, (length(vec) + 1 ,))
    exprs = vcat([:(vec[$i]) for i = 1:length(vec)], :x)
    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

@generated function unshift(vec::StaticVector, x)
    newtype = similar_type(vec, (length(vec) + 1 ,))
    exprs = vcat(:x, [:(vec[$i]) for i = 1:length(vec)])
    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

# unfortunately index can't be used directly by a @generated function... this code is pretty slow
@generated function insert(vec::StaticVector, index, x)
    newtype = similar_type(vec, (length(vec) + 1 ,))
    expr = :((Tuple(vec[1:index-1])..., x, Tuple(vec[index:end])...))
    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, expr))
    end
end

@generated function pop(vec::StaticVector)
    newtype = similar_type(vec, (length(vec) - 1 ,))
    exprs = [:(vec[$i]) for i = 1:length(vec)-1]
    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

@generated function shift(vec::StaticVector)
    newtype = similar_type(vec, (length(vec) - 1 ,))
    exprs = [:(vec[$i]) for i = 2:length(vec)]
    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

# unfortunately index can't be used directly by a @generated function... this code is pretty slow
@generated function deleteat(vec::StaticVector, index)
    newtype = similar_type(vec, (length(vec) - 1 ,))
    expr = :((Tuple(vec[1:index-1])..., Tuple(vec[index+1:end])...))
    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, expr))
    end
end

# TODO consider prepend, append (can use vcat, but eltype might change), and
# maybe splice (a bit hard to get statically sized)