text:
using Pkg
Pkg.activate(".")
using EquivalenceClassesConstructor
using Printf, DataStructures
using JLD2

include("./vertexIntTriple.jl")

const nBose = 150
const nFermi = 150
const shift = 0
path = "/scratch/usr/hhpstobb/grids/b$(nBose)_f$(nFermi)_s$(shift)"

# ===== Test with integers =====
println("Integer test")
const nB = UInt32(2*maximum((2*nFermi+shift*ceil(nBose/2),2*nBose+1))+1)
const nBh = floor(UInt32, nB/2)
(nB^2 > typemax(UInt32)) && throw(ArgumentError("nBose or nFermi too large for UInt32 operations."))
const nB2 = UInt32(nB*nB)
const nBm1 = UInt32(nB-1)
const oobConst = tripleToInt(nB,0,0,nB,nB2)

@fastmath @inline function reverseSymm(z::UInt32)
    r,k = divrem(z,nB)
    i,j = divrem(r,nB)
    return tripleToInt(UInt32(-i+nBm1), UInt32(-j+nBm1-1),UInt32(-k+nBm1-1),nB,nB2)
end

@fastmath @inline function symm_map_identity(z::UInt32)
    (z,)
end

@fastmath @inline function symm_map_int(z::UInt32)
    r,k = divrem(z,nB)
    i,j = divrem(r,nB)
    iu = UInt32(i)
    ju = UInt32(j)
    ku = UInt32(k)
    ni = UInt32(-i+nBm1)
    nj = UInt32(-j+nBm1)
    nk = UInt32(-k+nBm1)
    t1 = ku+nj-nBh #i and j both contain nBh, we need to subtract the double shift
    t2 = iu+ju-nBh
    t3 = iu+ku-nBh
    t4 = ju+nk-nBh
    r1 = tripleToInt(ni, nj-UInt32(1), nk-UInt(1),nB,nB2) # c.c.: Op 1 (*)
    r2 = tripleToInt(iu, ku, ju,nB,nB2)                   # time reversal: Op 2 (1)
    r3 = oobConst                                         # double Crossing: Op 3 (-)
    r4 = oobConst                                         # crossing c: Op 4 (-)
    r5 = oobConst                                         # crossing c^t: Op 5 (+)
    if t2 < nB
        if t1 < nB
            r3 = tripleToInt(t1, ju, t2,nB,nB2)
        end
        if t3 < nB
            r5 = tripleToInt(ni, t3, t2,nB,nB2)
        end
    end
    if (t3 < nB && t4 < nB)
        r4 = tripleToInt(t4, t3, ku,nB,nB2)
    end
    return r1,r2,r3,r4,r5
end

function uint_to_index(parents::Dict{UInt32,UInt32}, ops::Dict{UInt32,UInt32}, vl::Array{UInt32,1})
    parents_new = Array{Int64}(undef, length(vl))
    ops_new = Array{Int64}(undef, length(vl))
    lookup = Dict(zip(vl,1:length(vl)))
    for (i,el) in enumerate(vl)
        parents_new[i] = lookup[parents[el]]
        ops_new[i] = Int64(ops[el])
    end
    return parents_new, ops_new
end

println("Constructing Array")
freqList = [(i,j,k) for i in (-nBose:nBose) for j in (-nFermi:nFermi-1) .- trunc(Int64,shift*i/2) for k in (-nFermi:nFermi-1) .- trunc(Int64,shift*i/2)]
const freqList_int = map(x->tripleToInt(x..., nBh,nB,nB2), freqList)
const len_freq = (2*nBose+1)*(2*nFermi)^2

const mm_2 = Mapping(symm_map_int)
println("Starting Computation 3")
maxF = nBose + 2*nFermi + shift*ceil(Int,nBose/2) + 5
headerstr= @sprintf("%10d", maxF)

@time parents_int, ops_int = find_classes(mm_2.f, freqList_int, UInt32.([1, 2, 3, 4, 5]), vl_len=len_freq);
@time freqRed_map, freqList_min_int = minimal_set(parents_int, freqList_int)
freqList_min = intToTriple.(freqList_min_int)
@time parents,ops = uint_to_index(parents_int, ops_int, freqList_int)

mkpath(path)
write_fixed_width(path*"/freqList.dat", freqList_min, sorted=true, header_add=headerstr);
#write_JLD("freqList_2.jld", rm_2, expMap)
#save("freqList.jld", "ExpandMap", expMap, "ReduceMap", redMap, "base", nB, "nFermi", 2*nFermi, "nBose", 2*nBose+1, "shift", shift, "offset", nBh)
base = nB
offset = nBh
@save path*"/freqList.jld2" freqRed_map freqList freqList_min parents ops nFermi nBose shift base offset


text:
function solve_str8ts(
    grid,
    white;
    backtrack = true,
    all_solutions = false,
    keep_logs = false,
    logging = [:Info, :Table],
)
    straight_tables = Vector{Array{Int,2}}(undef, 9)
    for i in 1:9
        straight_tables[i] = get_straights(collect(1:9), i)
    end

    m = Model(optimizer_with_attributes(
        CS.Optimizer,
        "backtrack" => backtrack,
        "keep_logs" => keep_logs,
        "all_solutions" => all_solutions,
        "logging" => logging,
    ))
    @variable(m, 0 <= x[1:9, 1:9] <= 9, Int)

    # set variables
    for r in 1:9, c in 1:9
        if grid[r, c] != 0
            @constraint(m, x[r, c] == grid[r, c])
        elseif white[r, c] == 1
            @constraint(m, x[r, c] >= 1)
        else # black ones without a number
            @constraint(m, x[r, c] == 0)
        end
    end

    straights = []
    for r in 1:9
        found_straight = false
        vec = Vector{Tuple{Int,Int}}()
        for c in 1:9
            if white[r, c] == 1
                found_straight = true
                push!(vec, (r, c))
            elseif found_straight
                push!(straights, copy(vec))
                empty!(vec)
                found_straight = false
            end
        end
        if found_straight
            push!(straights, copy(vec))
        end
    end

    for c in 1:9
        found_straight = false
        vec = Vector{Tuple{Int,Int}}()
        for r in 1:9
            if white[r, c] == 1
                found_straight = true
                push!(vec, (r, c))
            elseif found_straight
                push!(straights, copy(vec))
                empty!(vec)
                found_straight = false
            end
        end
        if found_straight
            push!(straights, copy(vec))
        end
    end

    for r in 1:9
        variables = [x[r, c] for c = 1:9 if white[r, c] == 1 || grid[r, c] != 0]
        @constraint(m, variables in CS.AllDifferent())
    end
    for c in 1:9
        variables = [x[r, c] for r = 1:9 if white[r, c] == 1 || grid[r, c] != 0]
        @constraint(m, variables in CS.AllDifferent())
    end

    for straight in straights
        len = length(straight)
        variables = [x[s[1], s[2]] for s in straight]
        @constraint(m, variables in CS.TableSet(straight_tables[len]))
    end

    optimize!(m)
    status = JuMP.termination_status(m)
    return status, m, x
end

function get_straights(numbers, len)
    sort!(numbers)
    nrows = (length(numbers) - len + 1) * factorial(len)
    table = Array{Int64}(undef, (nrows, len))
    i = 1
    for j in 1:(length(numbers) - len + 1)
        l = numbers[j:(j + len - 1)]
        for row in permutations(l, len)
            table[i, :] = row
            i += 1
        end
    end
    return table
end

@testset "Str8ts no backtrack" begin
    grid = zeros(Int, (9, 9))
    grid[1, :] = [0 0 0 0 0 0 0 3 0]
    grid[2, :] = [0 0 0 0 8 0 0 0 0]
    grid[3, :] = [0 0 0 0 0 0 1 0 0]
    grid[4, :] = [0 0 0 0 0 6 0 0 0]
    grid[5, :] = [0 0 1 0 0 0 0 0 8]
    grid[6, :] = [0 0 0 4 0 0 9 0 0]
    grid[7, :] = [0 0 3 5 0 0 0 0 0]
    grid[8, :] = [0 0 0 0 0 0 0 0 2]
    grid[9, :] = [5 0 9 0 0 0 0 0 0]

    white = zeros(Int, (9, 9))
    white[1, :] = [1 1 1 1 0 0 1 1 0]
    white[2, :] = [0 1 1 1 0 1 1 1 0]
    white[3, :] = [0 1 1 0 1 1 0 1 1]
    white[4, :] = [1 1 0 1 1 0 1 1 1]
    white[5, :] = [1 1 1 1 1 1 1 1 1]
    white[6, :] = [1 1 1 0 1 1 0 1 1]
    white[7, :] = [1 1 0 1 1 0 1 1 0]
    white[8, :] = [0 1 1 1 0 1 1 1 0]
    white[9, :] = [0 1 1 0 0 1 1 1 1]

    status, m, x = solve_str8ts(grid, white; all_solutions = true, logging = [])
    com = CS.get_inner_model(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test MOI.get(m, MOI.ResultCount()) == 1
    @test convert.(Int, JuMP.value.(x[1, :])) == [7, 9, 6, 8, 0, 0, 2, 3, 0]
    @test convert.(Int, JuMP.value.(x[2, :])) == [0, 6, 5, 7, 8, 4, 3, 2, 0]
    @test convert.(Int, JuMP.value.(x[3, :])) == [0, 3, 4, 0, 6, 5, 1, 8, 9]
    @test convert.(Int, JuMP.value.(x[4, :])) == [4, 5, 0, 2, 3, 6, 8, 9, 7]
    @test convert.(Int, JuMP.value.(x[5, :])) == [2, 4, 1, 3, 5, 9, 7, 6, 8]
    @test convert.(Int, JuMP.value.(x[6, :])) == [3, 1, 2, 4, 7, 8, 9, 5, 6]
    @test convert.(Int, JuMP.value.(x[7, :])) == [1, 2, 3, 5, 4, 0, 6, 7, 0]
    @test convert.(Int, JuMP.value.(x[8, :])) == [0, 7, 8, 6, 0, 3, 5, 4, 2]
    @test convert.(Int, JuMP.value.(x[9, :])) == [5, 8, 9, 0, 0, 2, 4, 1, 3]
end

ext: jl


text:
"""
$(DocStringExtensions.README)
"""
module GalacticOptim

using DocStringExtensions
using Reexport
@reexport using SciMLBase
using Requires
using DiffResults
using Logging, ProgressLogging, Printf, ConsoleProgressMonitor, TerminalLoggers, LoggingExtras
using ArrayInterface, Base.Iterators

import SciMLBase: OptimizationProblem, OptimizationFunction, AbstractADType, __solve

include("solve/solve.jl")
include("function/function.jl")

function __init__()
    # Optimization backends
    @require BlackBoxOptim="a134a8b2-14d6-55f6-9291-3336d3ab0209" include("solve/blackboxoptim.jl")
    @require CMAEvolutionStrategy="8d3b24bd-414e-49e0-94fb-163cc3a3e411" include("solve/cmaevolutionstrategy.jl")
    @require Evolutionary="86b6b26d-c046-49b6-aa0b-5f0f74682bd6" include("solve/evolutionary.jl")
    @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" include("solve/flux.jl")
    @require MultistartOptimization="3933049c-43be-478e-a8bb-6e0f7fd53575" include("solve/multistartoptimization.jl")
    @require NLopt="76087f3c-5699-56af-9a33-bf431cd00edd" include("solve/nlopt.jl")
    @require Optim="429524aa-4258-5aef-a3af-852621145aeb" include("solve/optim.jl")
    @require QuadDIRECT="dae52e8d-d666-5120-a592-9e15c33b8d7a" include("solve/quaddirect.jl")

    # AD backends
    @require FiniteDiff="6a86dc24-6348-571c-b903-95158fe2bd41" include("function/finitediff.jl")
    @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" include("function/forwarddiff.jl")
    @require ReverseDiff="37e2e3b7-166d-5795-8a7a-e32c996b4267" include("function/reversediff.jl")
    @require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" include("function/tracker.jl")
    @require Zygote="e88e6eb3-aa80-5325-afca-941959d7151f" include("function/zygote.jl")
end

export solve

end # module


text:
using LasIO
using LazIO
using FileIO
using StaticArrays
using GeoArrays; const G = GeoArrays
using Colors
using PyPlot
using PyCall
using Plots
using ArchGDAL ; const AG = ArchGDAL ## make AG the shortcut to ArchGDAL, by doing this I call AG. rather than ArchGDAL.
using JSON

## ------------FUNCTIONS-------------- ##
"Store coordinates"
function coordinates(points,header)
    # coords_array = SArray{Tuple{4},Float64,1,4}[]
    coords_dict = Dict{Int64,Array{Float64,1}}()
    for (ind,p) in enumerate(points)
        x = xcoord(p,header)
        y = ycoord(p,header)
        z = zcoord(p,header)
        inten = intensity(p)
        ret = return_number(p)
        nret = number_of_returns(p)
        classificc = classification(p)
        gps = gps_time(p)
        coords_dict[ind] = [x,y,z,inten]
        # push!(coords_array,(ind, x, y, z))
    end
    return coords_dict
end

"Returns minimum boundary point"
function minboundrp(point::SArray{Tuple{3},Float64,1,3},voxelsize::SArray{Tuple{3},Float64,1,3})
    x,y,z = point
    zpx = x - voxelsize[1]
    zpy = y - voxelsize[2]
    zpz = z - voxelsize[3]
    newpoint = SArray{Tuple{3},Float64,1,3}(zpx,zpy,zpz)
    return newpoint
end

"Returns maximum boundary point"
function maxboundrp(point::SArray{Tuple{3},Float64,1,3},voxelsize::SArray{Tuple{3},Float64,1,3})
    x,y,z = point
    zpx = x + voxelsize[1]
    zpy = y + voxelsize[2]
    zpz = z + voxelsize[3]
    newpoint = SArray{Tuple{3},Float64,1,3}(zpx,zpy,zpz)
    return newpoint
end

"Local maxima, returns indexes"
function findlocalmaxima_indexes(signal)
    inds = Int[]
    if length(signal)>1
        if signal[1]>signal[2]
            push!(inds,1)
        end
        for i=2:length(signal)-1
            if signal[i-1]<signal[i]>signal[i+1]
                push!(inds,i)
            end
        end
        if signal[end]>signal[end-1]
            push!(inds,length(signal))
        end
    end
    inds
end

"Creates a histogram for every voxel"
function myhistogram(gr_points_to_voxels)
    alldepth = Float64[]; count_voxs = Int64[]
    allkeys = collect(keys(gr_points_to_voxels))
    for i in allkeys
        voxel = gr_points_to_voxels[i]
        len = length(voxel)
        # depth = Float64[] #vector with Floats
        for i in voxel
            z = coords_dict[i][3]
            push!(alldepth,z)
        end
        push!(count_voxs,len)
    end
    len = length(alldepth)
    minn = minimum(alldepth)
    maxx = maximum(alldepth)
    return minn,maxx,len,alldepth, allkeys
end

"Range of the z_start and z_stop of all the histogram bins of a voxel"
function ranges(rang1,rang2,voxel)
    ids =  Dict{Int64,Float64}()
    for i in voxel
        if rang1 <= coords_dict[i][3] <= rang2
            ids[i] = coords_dict[i][3]
        end
    end
    return ids
end

"Write confindence values as an attribute of a point"
function write_confidence(writeoutput, ds, confiden_vals, median_nden, median_nint, mean_ndis)
    confidenkeys = collect(keys(confiden_vals))
    LazIO.write(writeoutput, ds.header) do io
        for (i,p) in enumerate(ds) #I need to interate only to dataset, since it is Lazpoints
            if i in confidenkeys
                #1. density, 2. distance, 3. intensity
                if confiden_vals[i][1] > median_nden &&   confiden_vals[i][2] > mean_ndis && confiden_vals[i][3] > median_nint
                    p.classification = UInt(1)
                elseif confiden_vals[i][1] > median_nden && confiden_vals[i][2] <= mean_ndis  && confiden_vals[i][3] > median_nint
                    p.classification = UInt(2)
                elseif confiden_vals[i][1] > median_nden &&   confiden_vals[i][2] > mean_ndis && confiden_vals[i][3] <= median_nint
                    p.classification = UInt(3)
                elseif confiden_vals[i][1] <= median_nden &&  confiden_vals[i][2] > mean_ndis &&  confiden_vals[i][3] > median_nint
                    p.classification = UInt(4)
                elseif  confiden_vals[i][1] <= median_nden  &&  confiden_vals[i][2] > mean_ndis &&  confiden_vals[i][3] <= median_nint
                    p.classification = UInt(5)
                elseif confiden_vals[i][1] > median_nden && confiden_vals[i][2] <= mean_ndis  && confiden_vals[i][3] <= median_nint
                    p.classification = UInt(6)
                elseif confiden_vals[i][1] <= median_nden &&  confiden_vals[i][2] <= mean_ndis  &&confiden_vals[i][3] > median_nint
                    p.classification = UInt(7)
                elseif  confiden_vals[i][1] <= median_nden  &&  confiden_vals[i][2] <= mean_ndis   && confiden_vals[i][3] <= median_nint
                    p.classification = UInt(8)
                end
                LazIO.writepoint(io,p)
            end
        end
    end
end

function confidence(ds, confiden_vals, median_nden, median_nint, mean_ndis)
    confidenkeys = collect(keys(confiden_vals))
    conf = Dict()
    for (i,p) in enumerate(ds) #I need to interate only to dataset, since it is Lazpoints
        if i in confidenkeys
            #1. density, 2. distance, 3. intensity
            if confiden_vals[i][1] > median_nden &&   confiden_vals[i][2] > mean_ndis && confiden_vals[i][3] > median_nint
                confid = 1
            elseif confiden_vals[i][1] > median_nden && confiden_vals[i][2] <= mean_ndis  && confiden_vals[i][3] > median_nint
                confid = 2
            elseif confiden_vals[i][1] > median_nden &&   confiden_vals[i][2] > mean_ndis && confiden_vals[i][3] <= median_nint
                confid = 3
            elseif confiden_vals[i][1] <= median_nden &&  confiden_vals[i][2] > mean_ndis &&  confiden_vals[i][3] > median_nint
                confid = 4
            elseif  confiden_vals[i][1] <= median_nden  &&  confiden_vals[i][2] > mean_ndis &&  confiden_vals[i][3] <= median_nint
                confid = 5
            elseif confiden_vals[i][1] > median_nden && confiden_vals[i][2] <= mean_ndis  && confiden_vals[i][3] <= median_nint
                confid = 6
            elseif confiden_vals[i][1] <= median_nden &&  confiden_vals[i][2] <= mean_ndis  &&confiden_vals[i][3] > median_nint
                confid = 7
            elseif  confiden_vals[i][1] <= median_nden  &&  confiden_vals[i][2] <= mean_ndis   && confiden_vals[i][3] <= median_nint
                confid = 8
            end
            conf[i] = [confid]
        end
    end
    return conf
end

## ---------------MAIN-----------------##
"read - import dataset from a json file"
workdir = dirname(@__FILE__)
open(joinpath(workdir,"voxelization.json"), "r") do f
    global mydict = Dict(JSON.parse(read(f,String)))  # parse and transform data
end

#get the filename
filenn_NL1 = joinpath(workdir,mydict["dataname"] * ".laz")
filenn_NL1_las = File{format"LAZ_"}(filenn_NL1)

#read specific LAS file
header, points = LazIO.load(filenn_NL1_las) #point cloud as LAS points format
n = length(points) #number of points

#read dataset
ds = LazIO.open(filenn_NL1) #dataset
bb = LasIO.boundingbox(header) #boundary box of the pointcloud
coords_dict = coordinates(points,header) #coordinates

#### "VOXELIZATION_ALGORITHM, Nourian et.al." ####
#voxelsize in x,y dimensions
minb = bb[3]; maxb = bb[6]
voxel_z = round((maxb - minb),digits=1)
voxelsize = SArray{Tuple{3},Float64,1,3}(mydict["vx"], mydict["vy"], voxel_z)

#new boundary - lower left and upper right
bmin = SArray{Tuple{3},Float64,1,3}(bb.xmin, bb.ymin, bb.zmin)
bmax = SArray{Tuple{3},Float64,1,3}(bb.xmax, bb.ymax, bb.zmax)
nb_min = minboundrp(bmin,voxelsize) #new boundary- lower left
nb_max = maxboundrp(bmax,voxelsize) #new boundary- upper right

#dictionary all the points (indexes) --> i,j,k
vox_dict = Dict{Int64,Array{Int64,1}}()
for (ind,p) in enumerate(coords_dict)
    i = round(abs(coords_dict[ind][1] - nb_min[1])/voxelsize[1])
    j = round(abs(coords_dict[ind][2] - nb_min[2])/voxelsize[2])
    k = round(abs(coords_dict[ind][3] - nb_min[3])/voxelsize[3])
    vox_dict[ind] = [i,j,k]
end

#UNIQUE VALUES in the dictionary
#Store INDEXES of the points that correspond to EVERY Voxel
val = unique(values(vox_dict))
#dictionary all the voxels (keys) --> all the indexes of the points PER Voxel (values)
gr_points_to_voxels = Dict{Int64,Array{Int64,1}}()
for (i,every) in enumerate(val)
    result = filter((k,v)->v == every,vox_dict)
    gr_points_to_voxels[i] = collect(keys(result))
end

#PLOTS
#Return some characteristics of a histogram
minn, maxx, len, alldepth, allkeys = myhistogram(gr_points_to_voxels)

nplot = length(gr_points_to_voxels) # number of plots
num_plots = [i for i in range(1,stop = length(gr_points_to_voxels)-(length(gr_points_to_voxels)-nplot))] #selecting specific number of plots for plotting

#Create histogram plots and save the required points
c = 1
imin =  Dict{Int64,Array{Float64,1}}() # store the point id --> density, distance, intensity
watersurface =  Int64[] # just the watersurface points
for i in allkeys
    if c <= length(num_plots)
        voxel = gr_points_to_voxels[i] #points of a voxel
        npoints = length(voxel) #amount of points per voxel
        if length(voxel) > 1 #if there are points
            zvals = Float64[] #vector with Floats
            for i in voxel
                push!(zvals,coords_dict[i][3])
            end
            len = length(zvals); minn = minimum(zvals); maxx = maximum(zvals)
            dis = abs(maxx - minn);  width = 0.05  #5cm spacing
            nb = Int64(round(dis/width)) #number of bins per voxel

            if nb > 1
                #Create a histogram PER Voxel, with a specific Number of Bins (nb)
                plt.figure()
                plt.hist(zvals,bins=nb) #check it here
                plt.ylabel("Count of points")
                plt.xlabel("Z value (m)")
                plt.title("Voxel "*string(i)*": Count of points with respect to z value")

                #Extract info from a histogram
                counts, bins, bars = plt.hist(zvals, bins=nb)
                bins = [round(i,digits=2) for i in bins] #round them to 2 digits
                println("mycounts : ",counts)
                println("my bins : ",bins)
                println("Voxel: ",i)
                println()

                #1. MAX BIN
                # find all INDEXES the peaks from counts
                peaks = findlocalmaxima_indexes(counts)
                println("My peaks in the histogram: ",peaks)
                # just one peak
                if length(peaks) == 1
                    val = peaks[1] #extract that element

                    density = counts[val]/npoints #density
                    distance = -1.0 #distance

                    rang = bins[val:val+1] #find the values of Z range of that bin

                    if bins[val] >= mean(bins) #if the HIGHEST peak is on the RIGHT side of the Mean

                        # find the lowest point of the BIN
                        ids_above = ranges(rang[1],rang[2],voxel)
                        if length(ids_above) > 0
                            minpoint_high = findmin(ids_above)
                            push!(watersurface,minpoint_high[2]) # store it as a watersurface point
                        else
                            continue
                        end

                    else    #if it is BELOW the mean
                        ids_below = ranges(rang[1],rang[2],voxel)
                        if length(ids_below) > 0
                            minpoint_low = findmin(ids_below) #1. height, 2. index in the voxel
                            inten = coords_dict[minpoint_low[2]][4] #intensity
                            imin[minpoint_low[2]] = [minpoint_low[1], counts[val], distance, inten]
                            else
                            continue
                        end
                    end

                elseif length(peaks) > 1

                    println("I have more than 1 peak")
                    println("length of bins ",length(bins))
                    println(counts, bins)

                    values = sort(counts[peaks], rev=true) #access the values of peaks, sort from the high to the low
                    print(values)

                    #1.
                    first = values[1]
                    getind_first = findall(x->x==first, counts)

                    #if the first HIGHEST peak is below the mean, then, it's the minpoint
                    if bins[getind_first[1]] <=  mean(bins)

                        density = first/npoints #density
                        distance = -1.0 #distance

                        rang_first = bins[getind_first[1]:getind_first[1]+1] #range in z level o f a bin
                        ids_fir = ranges(rang_first[1],rang_first[2],voxel)

                        if length(ids_fir) > 0
                            minpoint_low = findmin(ids_fir)
                            inten = coords_dict[minpoint_low[2]][4]
                            imin[minpoint_low[2]] = [minpoint_low[1], first, distance, inten]
                        else
                            continue
                        end
                    else #the first HIGHEST is above the mean, so watersurface!

                        #1. watersurface point
                        rang_first = bins[getind_first[1]:getind_first[1]+1]
                        ids_fir = ranges(rang_first[1],rang_first[2],voxel)

                        if length(ids_fir) > 0
                            minpoint_high = findmin(ids_fir)
                            push!(watersurface,minpoint_high[2])
                        else
                            continue
                        end
                        #2.
                        sec = values[2]
                        getind_sec = findall(x->x==sec, counts)

                        density = sec/npoints #density
                        distance = abs(bins[getind_first[1]] - bins[getind_sec[1]]) #distance

                        rang_sec = bins[getind_sec[1]:getind_sec[1] + 1]
                        ids_sec =  ranges(rang_sec[1],rang_sec[2],voxel)

                        if length(ids_sec) > 0
                            minpoint_low = findmin(ids_sec)
                            inten = coords_dict[minpoint_low[2]][4]
                            imin[minpoint_low[2]] = [minpoint_low[1], sec, distance, inten]
                        else
                            continue
                        end
                    end
                else
                    continue
                end
                println()
                #To SAVE the plots
                path = "C:/Users/alexandr/OneDrive - Stichting Deltares/Desktop/Thesis/P3_processing/with_greenstowa/Working_to/199NL3_Voxels/"
                path = joinpath(path,"nbins_" * string(length(bins)-1) * "_voxel_" * string(i))
                # plt.savefig(path,dpi=300)
                plt.close()
                global c+=1
            else
                imin[i] = [0, 0, 0, 0, 0] #if the number of bins is lower than 2, then I initialize all zeros
            end
        end
    end
end

#Get MAXIMUM values of density, distance, intensity for normalization procedure
den=[]; dis=[]; int=[];
for id in collect(values(imin))
    push!(den,id[2]); push!(dis,id[3]); push!(int,id[4])
end
maxden = maximum(den); maxdis = maximum(dis); maxint = maximum(int)

#Compute the Confidence values
confiden_vals = Dict{Int64,Array{Float64,1}}()
for i in keys(imin)
    if imin[i][3] == -1
        confiden_vals[i] = [imin[i][2]/maxden, 0, imin[i][4]/maxint]
    else
        confiden_vals[i] = [imin[i][2]/maxden, imin[i][3]/maxdis, imin[i][4]/maxint]
    end
end

#Define the threshold values of density, intensity and distance
nden=[]; ndis=[]; nint=[];
for id in collect(values(confiden_vals))
    push!(nden,id[1]); push!(ndis,id[2]); push!(nint,id[3])
end
#Mean and Median thresholds for the density, intensity and distance
mean_nden = mean(nden); mean_ndis = mean(ndis) ; mean_nint = mean(nint)
median_nden = median(nden); median_dis = median(ndis); median_nint = median(nint)

#write confidence values PER point in a las file
writeoutput = joinpath(workdir, mydict["dataname"] * mydict["outputname"] * ".laz")
write_confidence(writeoutput, ds, confiden_vals, median_nden, median_nint, mean_ndis)

## ------Write a RASTER --------- ##
np = pyimport("numpy")
rasterio = pyimport("rasterio")

#Create the extend of the array
keyes = collect(keys(vox_dict))
xxtend = Int64(round((nb_max[1] - nb_min[1])/voxelsize[1]))
yxtend = Int64(round((nb_max[2] - nb_min[2])/voxelsize[1]))
bands = 5 #number of bands
arr = np.zeros((xxtend,yxtend,bands)) #array in the bound box extend

#get confidence values per point
myconf = confidence(ds, confiden_vals, median_nden, median_nint, mean_ndis)

#store data into a multidimensional array
all = []
for key in keyes
    if key in collect(keys(imin))
        println("key: ",key)
        row, col, dontneedit = vox_dict[key]
        if row!=0 && col!=0 && dontneedit!=0
            arr[row,col,1] = imin[key][1] #z
            arr[row,col,2] = imin[key][2] #density
            arr[row,col,3] = imin[key][3] #distance
            arr[row,col,4] = imin[key][4] #intensity
            arr[row,col,5] = myconf[key][5] #confidence
            push!(all,(row,col))
        else
            continue
        end
    end
end

#Write raster output with the multidimensional array
rot = np.rot90(arr, k=1, axes=(0, 1)) #rotate the matrix
outtiff = joinpath(workdir,mydict["dataname"] * mydict["rasteroutput"])
f = rasterio.open(outtiff,
                  "w+",
                  driver= "GTiff",
                  width= xxtend,
                  height= yxtend,
                  count= bands,
                  dtype= rasterio.float64,
                  crs = "EPSG:28992",
                  nodata= -9999,)
f.write(rot[:,:,1],1)
f.write(rot[:,:,2],2)
f.write(rot[:,:,3],3)
f.write(rot[:,:,4],4)
f.write(rot[:,:,5],5)
f.close()


text:
# Information about mathematical functions needed for expression parsing

#################################################
#
#    Additional math functions
#
#################################################

sqr(x::Number) = x * x
rcp(x::Number) = one(x) / x
rsqrt(x::Number) = one(x) / sqrt(x)
rcbrt(x::Real) = one(x) / cbrt(x)

logit(x::Real) = log(x / (one(x) - x))
xlogx(x::Real) = x > zero(x) ? x * log(x) : zero(x)
xlogy(x::Real, y::Real) = x > zero(x) ? x * log(y) : zero(x)

logistic(x::Real) = one(x) / (one(x) + exp(-x))
invlogistic(y::Real) = -log(one(y) / y - one(y))

softplus(x::Real) = log(one(x) + exp(x))
invsoftplus(x::Real) = log(exp(x) - one(x))

logsumexp{T<:Real}(x::T, y::T) = x > y ? x + log(one(T) + exp(y - x)) : y + log(one(T) + exp(x - y))

@vectorize_1arg Number sqr
@vectorize_1arg Number rcp
@vectorize_1arg Real rsqrt
@vectorize_1arg Real rcbrt

@vectorize_1arg Real logit
@vectorize_1arg Real logistic


text:
module Exceptions

import Revise
using SearchLight

export NotPersistedException, MissingDatabaseConfigurationException, DatabaseAdapterException
export UnretrievedModelException, InvalidModelException

struct NotPersistedException <: Exception
  model
end

Base.showerror(io::IO, ex::NotPersistedException) = print(io, "NotPersistedException: Model is not persisted \n$(ex.model)")

struct MissingDatabaseConfigurationException <: Exception
  msg::String
end
MissingDatabaseConfigurationException() = MissingDatabaseConfigurationException("The database configuration can not be found")

struct DatabaseAdapterException <: Exception
  msg::String
end

struct UnretrievedModelException <: Exception
  model
  id
end

Base.showerror(io::IO, ex::UnretrievedModelException) = print(io, "UnretrievedModelException: Model can not be retrieved for id $(ex.id)")

struct InvalidModelException <: Exception
  model
  errors
end

Base.showerror(io::IO, ex::InvalidModelException) =
  print(io, "Validation errors for $(typeof(ex.model)): $(join( map(e -> "$(e.field) $(e.error_message)", ex.errors), ", "))")

end


text:
module Day1

import DelimitedFiles

# https://adventofcode.com/2021/day/1
# To do this, count the number of times a depth measurement increases from the previous measurement.
# (There is no measurement before the first measurement.)
# In the example above, the changes are as follows:

# 199 (N/A - no previous measurement)
# 200 (increased)
# 208 (increased)
# 210 (increased)
# 200 (decreased)
# 207 (increased)
# 240 (increased)
# 269 (increased)
# 260 (decreased)
# 263 (increased)
# In this example, there are 7 measurements that are larger than the previous measurement.

# How many measurements are larger than the previous measurement?

function num_increases(x::Vector{Int64})
    return sum(x[i] > x[i-1] for i=2:length(x))
end

function num_3_window_increases(x::Vector{Int64})
    count = 0
    for i = 2:(length(x)-2)
        previous_window = x[i-1:i+1]
        this_window = x[i:i+2]
        if sum(this_window) > sum(previous_window)
            count += 1
        end
    end
    return count
end

function read_input()
    input_filepath = joinpath(@__DIR__, "input.txt")
    return DelimitedFiles.readdlm(input_filepath, ',', Int64, '\n')[:,1]
end

function part1_answer()
    return num_increases(read_input())
end

function part2_answer()
    return num_3_window_increases(read_input())
end

end


text:
###
  Name     : c10_30_QuantEcon_simulated.jl
  Book     : Hands-on Data Science with Anaconda )
  Publisher: Packt Publishing Ltd.
  Author   : Yuxing Yan and James Yan
  Date     : 4/24/2018
  email    : yany@canisius.edu
             paulyxy@hotmail.com
###

using QuantEcon
P = [0.4 0.6; 0.2 0.8];
mc = MarkovChain(P)
x = simulate(mc, 100000);
mean(x .== 1)
#
mc2 = MarkovChain(P, ["employed", "unemployed"])
simulate(mc2, 4)


text:
mutable struct Dendrite <: AbstractDendrite
    model::Model.ModelData
    soma::AbstractSoma

    # These fields are not used yet.
    taoEff::Float64
    length::Float64

    # Minimum value. Typically 0.0
    min_psp::Float64

    # Contains Compartments
    compartments::Array{AbstractCompartment,1}

    # Average weight over time
    average::Float64
    synapses::Int64

    function Dendrite(soma::AbstractSoma, model::Model.ModelData)
        o = new()
        o.model = model
        soma.dendrite = o
        o.compartments = Array{AbstractCompartment,1}()
        o.soma = soma
        o
    end
end

function initialize!(dendrite::AbstractDendrite)
    # Set properties based on model. These drive the other properties.
    dendrite.taoEff = Model.tao_eff(dendrite.model)
    dendrite.length = Model.dendrite_length(dendrite.model)
    dendrite.min_psp = Model.dendrite_min_psp(dendrite.model)
    dendrite.synapses = Model.synapses(dendrite.model)

    println("___ Dendrite properties ___")
    println("| taoEff: ", dendrite.taoEff)
    println("| length: ", dendrite.length)
    println("| min_psp: ", dendrite.min_psp)
    println("---------------------------")

    for compartment in dendrite.compartments
        initialize!(compartment)
    end
end

function add_compartment!(dendrite::AbstractDendrite, compartment::Compartment)
    push!(dendrite.compartments, compartment)
end

function reset!(dendrite::AbstractDendrite)
    for compartment in dendrite.compartments
        reset!(compartment)
    end
end

function AP_efficacy(dendrite::AbstractDendrite, distance::Float64)
    if distance < dendrite.length
        return 1.0
    end

    exp(-(dendrite.length - distance) / dendrite.taoEff)
end

function integrate!(dendrite::AbstractDendrite, span_t::Int64, t::Int64)
    psp = 0.0
    total_weight = 0.0

    for compartment in dendrite.compartments
        (sum, total) = integrate!(compartment, span_t, t)
        psp += sum
        total_weight += total
    end

    psp = max(psp, dendrite.min_psp)

    dendrite.average = total_weight / Float64(dendrite.synapses)

    # Collecting is centralized in streams.jl for consistency.
    collect_dendrite!(dendrite.soma.samples, dendrite, span_t)

    psp
end


text:
import Images
"""
    load_imgs(path, dims, numerical)

Reads all the volume images of an entire directory. Data must be stored as follows:
The required input "path" is the path to a directory containing one sub-directory
for every volume that has to be read. Each sub-directory contains the images making
up one volume in .tif format. if the image names contain alphabetical characters
(e.g. image001.tif, image002.tif ect.) the variable "numerical" has to be set to false, if
the image names are just numerical values (e.g 1.tif, 2.tif, 3.tif, ect), it has to
be set to true. Any file not ending in .tif will not be read.

Data will be stored in a 4D array [a x b c x n] with n being the number of
sub-directories, a and b the resolution of the .tif images and c the number of
images. The the resolution of the images has to be the same over all the volumes.
"""
function load_imgs(path, dims, numerical)
    inds = 0
    for (root, dirs, files) in walkdir(path)
        if length(dirs) != 0
            inds = length(dirs)
        end
    end
    image_matrix = zeros(dims[1], dims[2], dims[3], inds)
    vol_names = []
    for (root, dirs, files) in walkdir(path)
        for ind in 1:length(dirs)
          push!(vol_names, dirs[ind])
            ind_path = joinpath(path, dirs[ind])
            items = readdir(ind_path)
            if numerical
                sort!(items, by=x->parse(Int,split(x, ".")[1]))
            end
            for i in items
                if i[end-3:end] != ".tif"
                    deleteat!(items, findall(x->x==i, items))
                end
            end
            for img in 1:length(items)
                im_path = joinpath(ind_path, items[img])
                image = Images.load(im_path)
                image = Images.Gray.(image)
                mat = convert(Array{Float32}, image)
                if size(mat, 1) != 128
                end
                image_matrix[:, :, img, ind] = mat
            end
        end
    end
    return image_matrix, vol_names
end

"""
    read_landmarks(path, num_landmarks, group)

Specific function to read avizo landmark data into an array.
Reads all the files ending in .Ascii in the directory specified
as "path". Data will be stored in a 2D array [c x n] with c being the
number of individual 3D coordinates (30 coordinates in the case of 10 landmarks)
and n being the number of landmark files read.

The variable Group specifies the group of landmarks in the Avizo file that has
to be read (e.g "@1" in the case of group 1)
"""
function read_landmarks(path, num_landmarks, group)
    whole_dir = readdir(path)
    lm_files = []
    for item in whole_dir
        if '.' in item
            if item[end-4:end] == "Ascii"
                push!(lm_files, joinpath(path, item))
            end
        else
            continue
        end
    end
    coordinates = zeros(length(lm_files), num_landmarks * 3)
    for file in 1:length(lm_files)
        f = open(lm_files[file])
        fil = readlines(f)
        close(f)
        coords = []
        counter = 0
        for i in fil
            if i == group
                counter = 1
                continue
            end
            if 1<=counter<=num_landmarks
                push!(coords, i)
                counter+=1
            end
        end
        for point in 1:length(coords)
            single_co = split(coords[point], " ")
            single_co = single_co[1:3]
            for co in 1:3
                if parse(Float64, single_co[co][1:end-5]) * 10.0^parse(Float64, single_co[co][end-1:end]) <= 256
                    coordinates[file, point*3 - (3-co)] = parse(Float64, single_co[co][1:end-5]) * 10.0^parse(Float64, single_co[co][end-1:end])
                else
                    coordinates[file, point*3 - (3-co)] = parse(Float64, single_co[co][1:end-5]) * 10.0^-parse(Float64, single_co[co][end-1:end])
                end
            end
        end
    end
    return convert(Array, coordinates')
end

max_stars_repo_name: imi-diener/AutoLM

max_issues_repo_name: imi-diener/AutoLM

max_forks_repo_name: imi-diener/AutoLM


text:
struct GravityParameters
    μ::Float64
    J2::Float64
    R::Float64
    function GravityParameters()
        GM_MARS = 4.2828375816e13
        J2_MARS = 1960.45e-6
        R_MARS = 3386200.0
        new(GM_MARS,J2_MARS,R_MARS)
    end
end

@inline function gravity(g::GravityParameters,r::SVector{3, T}) where T

    nr = norm(r)

    x,y,z = r

    # precompute repeated stuff
    Re_r_sqr = 1.5*g.J2*(g.R/nr)^2
    five_z_sqr = 5*z^2/nr^2

    return  (-g.μ/nr^3)*SA[x*(1 - Re_r_sqr*(five_z_sqr - 1)),
                           y*(1 - Re_r_sqr*(five_z_sqr - 1)),
                           z*(1 - Re_r_sqr*(five_z_sqr - 3))]
end


text:
# -----------------------------------------------------------------------------
# Integration with AtomsBase
# -----------------------------------------------------------------------------

# Internal type used to represent a MassBody and retain the element metadata
struct ElementMassBody{cType<:Real,mType<:Real} <: Body
    r::SVector{3,cType} # in LENGTH_UNIT
    v::SVector{3,cType} # in VELOCITY_UNIT
    m::mType            # in MASS_UNIT
    s::Symbol
    n::Int
end
function ElementMassBody(r::SVector{3,<:Unitful.Length}, v::SVector{3,<:Unitful.Velocity}, e::Element)
    ElementMassBody{Float64,Float64}(austrip.(r), austrip.(v), austrip(e.atomic_mass), Symbol(e.symbol), e.number)
end

# Convert AtomsBase Atom to NBodySimulator body
function ElementMassBody(atom::Atom)
    ElementMassBody{Float64,Float64}(austrip.(position(atom)), austrip.(velocity(atom)), austrip(atomic_mass(atom)), AtomsBase.atomic_symbol(atom), atomic_number(atom))
end
# Convert NBodySimulator body to AtomsBase Atom
# TODO: support more boundary conditions
function AtomsBase.Atom(b::ElementMassBody, boundary_conditions::CubicPeriodicBoundaryConditions)
    Atom(b.s, mod.(b.r, boundary_conditions.L) .* LENGTH_UNIT, b.v .* VELOCITY_UNIT)
end

# Convert AtomsBase boundary conditions to NBodySimulator boundary conditions
function nbs_boundary_conditions(system::AbstractSystem)
    # TODO: support more boundary conditions
    box = bounding_box(system)
    @assert box[1][1] == box[2][2] == box[3][3]
    @assert boundary_conditions(system) == [Periodic(), Periodic(), Periodic()]
    CubicPeriodicBoundaryConditions(austrip(box[1][1]))
end
# Convert NBodySimulator boundary conditions to AtomsBase boundary conditions
# TODO: support more boundary conditions
ab_boundary_conditions(::CubicPeriodicBoundaryConditions) = [Periodic(), Periodic(), Periodic()]

# Convert NBodySimulator boundary conditions to AtomsBase bounding box
# TODO: support more boundary conditions
function ab_bounding_box(boundary_conditions::CubicPeriodicBoundaryConditions)
    box_size = boundary_conditions.L * LENGTH_UNIT
    z = zero(typeof(box_size))
    [[box_size, z, z], [z, box_size, z], [z, z, box_size]]
end

# Convert AtomsBase AbstractSystem to Vector of NBodySimulator bodies
function bodies(system::AbstractSystem)
    ElementMassBody.(system)
end
# Convert Vector of NBodySimulator bodies to AtomsBase FlexibleSystem
function AtomsBase.FlexibleSystem(bodies::AbstractVector{<:ElementMassBody}, boundary_conditions::BoundaryConditions)
    particles = Fix2(Atom, boundary_conditions).(bodies)
    FlexibleSystem(particles, ab_bounding_box(boundary_conditions), ab_boundary_conditions(boundary_conditions))
end
# Convert Vector of NBodySimulator bodies to AtomsBase FastSystem
function AtomsBase.FastSystem(bodies::AbstractVector{<:ElementMassBody}, boundary_conditions::BoundaryConditions)
    particles = Fix2(Atom, boundary_conditions).(bodies)
    FastSystem(particles, ab_bounding_box(boundary_conditions), ab_boundary_conditions(boundary_conditions))
end

# Extract the atomic symbol of each body in the system
get_atomic_symbols(system::PotentialNBodySystem{ElementMassBody}) = [b.s for b ∈ system.bodies]
# Extract the atomic number of each body in the system
get_atomic_numbers(system::PotentialNBodySystem{ElementMassBody}) = [b.n for b ∈ system.bodies]

# -----------------------------------------------------------------------------
# Convenience methods for generating starting configurations with element data
# -----------------------------------------------------------------------------

function NBodySimulator.generate_bodies_in_cell_nodes(n::Integer, symbol::Union{Integer,AbstractString,Symbol}, L::Unitful.Length, reference_temp::Unitful.Temperature; rng = MersenneTwister(n))
    e = elements[symbol]
    average_velocity = √(u"k" * reference_temp / e.atomic_mass)
    generate_bodies_in_cell_nodes(n, symbol, average_velocity, L, rng = rng)
end
function NBodySimulator.generate_bodies_in_cell_nodes(n::Integer, symbol::Union{Integer,AbstractString,Symbol}, average_velocity::Unitful.Velocity, L::Unitful.Length; rng = MersenneTwister(n))
    velocities = average_velocity * randn(rng, Float64, (3, n))
    e = elements[symbol]
    bodies = ElementMassBody[]

    count = 1
    dL = L / (ceil(n^(1 / 3)))
    for x ∈ dL/2:dL:L, y ∈ dL/2:dL:L, z ∈ dL/2:dL:L
        if count > n
            break
        end
        push!(bodies, ElementMassBody(SVector(x, y, z), SVector{3}(velocities[:, count]), e))
        count += 1
    end
    return bodies
end


text:
"""
A dense table stores the data as a vector of row-tuples.
"""
immutable DenseTable{Index,DataTypes} <: AbstractTable{Index,DefaultKey()}
    data::Vector{Row{Index,DataTypes}}

    function DenseTable(data_in)
        check_row(Index,DataTypes)
        new(data_in)
    end
end

DenseTable{Index}(::Index) = DenseTable{Index,eltypes(Index)}(Vector{Row{Index,eltypes(Index)}}())
DenseTable{Index,DataTypes}(::Index,data::Vector{DataTypes}) = DenseTable{Index,DataTypes}(Vector{Row{Index,DataTypes}}(data))
# Constructor that takes index and seperate values
# Conversion between Dense and normal Tables??
DenseTable{Index,StorageTypes}(table::Table{Index,StorageTypes}) = DenseTable(Index,collect(zip(table.data...)))

Base.length{Index,DataTypes}(table::DenseTable{Index,DataTypes}) = length(table.data)
Base.size{Index,DataTypes}(table::DenseTable{Index,DataTypes}) = (length(Index),length(table.data))
Base.size{Index,DataTypes}(table::DenseTable{Index,DataTypes},i::Int) = i == 1 ? length(Index) : (i == 2 ? length(table.data) : error("Tables are two-dimensional"))
Base.eltype{Index,DataTypes}(table::DenseTable{Index,DataTypes}) = eltype(table.data)

Base.push!{Index,DataTypes}(table::DenseTable{Index,DataTypes},row::Row{Index,DataTypes}) = push!(table.data,row)
Base.push!{Index,DataTypes}(table::DenseTable{Index,DataTypes},data_in::DataTypes) = push!(table.data,Row{Index,DataTypes}(data_in))
Base.push!{Index,DataTypes}(table::DenseTable{Index,DataTypes},data_in...) = push!(table,Row{Index,DataTypes}((data_in...))) # TODO probably slow
Base.append!{Index,DataTypes}(table::DenseTable{Index,DataTypes},table_in::DenseTable{Index,DataTypes}) = append!(table.data,table_in.data)
Base.append!{Index,DataTypes}(table::DenseTable{Index,DataTypes},rows::Vector{Row{Index,DataTypes}}) = append!(table.data,rows)

# Some output
function Base.summary{Index,DataTypes}(io::IO,table::DenseTable{Index,DataTypes})
    println(io,"$(length(table.data[1]))-row DenseTable with columns $Index")
end

function Base.show{Index,DataTypes}(io::IO,table::DenseTable{Index,DataTypes})
    summary(io,table)
    show(io,table.data)
end


text:
# module Operators
#
# export
#     δx_c2f,
#     δx_f2c,
#     δx_e2f,
#     δx_f2e,
#     δy_c2f,
#     δy_f2c,
#     δy_e2f,
#     δy_f2e,
#     δz_c2f,
#     δz_f2c,
#     δz_e2f,
#     δz_f2e,
#     avgx_c2f,
#     avgx_f2c,
#     avgx_f2e,
#     avgy_c2f,
#     avgy_f2c,
#     avgy_f2e,
#     avgz_c2f,
#     avgz_f2c,
#     avgz_f2e,
#     div_f2c,
#     div_c2f,
#     δx_f2c_ab̄ˣ,
#     δy_f2c_ab̄ʸ,
#     δz_f2c_ab̄ᶻ,
#     div_flux,
#     δx_c2f_ūˣūˣ,
#     δy_e2f_v̄ˣūʸ,
#     δz_e2f_w̄ˣūᶻ,
#     u∇u,
#     δx_e2f_ūʸv̄ˣ,
#     δy_c2f_v̄ʸv̄ʸ,
#     δz_e2f_w̄ʸv̄ᶻ,
#     u∇v,
#     δx_e2f_ūᶻw̄ˣ,
#     δy_e2f_v̄ᶻw̄ʸ,
#     δz_c2f_w̄ᶻw̄ᶻ,
#     u∇w,
#     fv,
#     fu,
#     δx²_c2f2c,
#     δy²_c2f2c,
#     δz²_c2f2c,
#     κ∇²,
#     δx²_f2c2f,
#     δy²_f2c2f,
#     δx²_f2e2f,
#     δy²_f2e2f,
#     δz²_f2e2f,
#     δz²_f2c2f,
#     𝜈∇²u,
#     𝜈∇²v,
#     𝜈∇²w,
#     ∇²_ppn,
#     ∇h_u

include("ops_regular_cartesian_grid.jl")

# end


text:
# This file contains code that was a part of Julia
# License is MIT: see LICENSE.md

export Unicode
module Unicode
export normalize, graphemes, isassigned
const normalize  = Base.UTF8proc.normalize_string
const graphemes  = Base.UTF8proc.graphemes
const isassigned = Base.UTF8proc.is_assigned_char
end

## Start of code from operators.jl =================================================
##
## It is used to support the new string searching syntax on v0.6.2

"""
    Fix2(f, x)

A type representing a partially-applied version of function `f`, with the second
argument fixed to the value "x".
In other words, `Fix2(f, x)` behaves similarly to `y->f(y, x)`.
"""
struct Fix2{F,T} <: Function
    f::F
    x::T

    Fix2(f::F, x::T) where {F,T} = new{F,T}(f, x)
    Fix2(f::Type{F}, x::T) where {F,T} = new{Type{F},T}(f, x)
end

(f::Fix2)(y) = f.f(y, f.x)

"""
    isequal(x)

Create a function that compares its argument to `x` using [`isequal`](@ref), i.e.
a function equivalent to `y -> isequal(y, x)`.

The returned function is of type `Base.Fix2{typeof(isequal)}`, which can be
used to implement specialized methods.
"""
isequal(x) = Fix2(isequal, x)

const EqualTo = Fix2{typeof(isequal)}

"""
    ==(x)

Create a function that compares its argument to `x` using [`==`](@ref), i.e.
a function equivalent to `y -> y == x`.

The returned function is of type `Base.Fix2{typeof(==)}`, which can be
used to implement specialized methods.
"""
==(x) = Fix2(==, x)

"""
    in(x)

Create a function that checks whether its argument is [`in`](@ref) `x`, i.e.
a function equivalent to `y -> y in x`.

The returned function is of type `Base.Fix2{typeof(in)}`, which can be
used to implement specialized methods.
"""
in(x) = Fix2(in, x)

const OccursIn = Fix2{typeof(in)}

## end of code from operators.jl =================================================

## Start of codeunits support from basic.jl ======================================
##
##It is used for CodeUnit support in pre v0.7 versions of Julia

## code unit access ##

codeunit(str::AbstractString) = UInt8
codeunit(::Type{<:AbstractString}) = UInt8
ncodeunits(str::AbstractString) = sizeof(str)

"""
    CodeUnits(s::AbstractString)

Wrap a string (without copying) in an immutable vector-like object that accesses the code units
of the string's representation.
"""
struct CodeUnits{T,S<:AbstractString} <: DenseVector{T}
    s::S
    CodeUnits(s::S) where {S<:AbstractString} = new{codeunit(s),S}(s)
end

length(s::CodeUnits) = ncodeunits(s.s)
sizeof(s::CodeUnits{T}) where {T} = ncodeunits(s.s) * sizeof(T)
size(s::CodeUnits) = (length(s),)
strides(s::CodeUnits) = (1,)
@propagate_inbounds getindex(s::CodeUnits, i::Int) = codeunit(s.s, i)
IndexStyle(::Type{<:CodeUnits}) = IndexLinear()
start(s::CodeUnits) = 1
next(s::CodeUnits, i) = (@_propagate_inbounds_meta; (s[i], i+1))
done(s::CodeUnits, i) = (@_inline_meta; i == length(s)+1)

write(io::IO, s::CodeUnits) = write(io, s.s)

unsafe_convert(::Type{Ptr{T}},    s::CodeUnits{T}) where {T} = unsafe_convert(Ptr{T}, s.s)
unsafe_convert(::Type{Ptr{Int8}}, s::CodeUnits{UInt8}) = unsafe_convert(Ptr{Int8}, s.s)

"""
    codeunits(s::AbstractString)

Obtain a vector-like object containing the code units of a string.
Returns a `CodeUnits` wrapper by default, but `codeunits` may optionally be defined
for new string types if necessary.
"""
codeunits(s::AbstractString) = CodeUnits(s)

## end of codeunits support ============================================================

size(cp::Chr) = ()
ndims(cp::Chr) = 0
ndims(::Type{<:Chr}) = 0
length(cp::Chr) = 1
lastindex(cp::Chr) = 1
getindex(cp::Chr) = cp
first(cp::Chr) = cp
last(cp::Chr) = cp
start(cp::Chr) = false
next(cp::Chr, state) = (cp, true)
done(cp::Chr, state) = state
isempty(cp::Chr) = false
in(x::AbsChar, y::Chr) = x == y
in(x::Chr, y::AbsChar) = x == y
in(x::Chr, y::Chr) = x == y
-(x::AbsChar, y::Chr) = Int(x) - Int(y)
-(x::Chr, y::AbsChar) = Int(x) - Int(y)
-(x::Chr, y::Chr) = Int(x) - Int(y)
-(x::Chr, y::Integer) = Chr((Int32(x) - Int32(y))%UInt32)
+(x::Chr, y::Integer) = Chr((Int32(x) + Int32(y))%UInt32)
+(x::Integer, y::Chr) = y + x
show(io::IO, cp::Chr)  = show(io, Char(cp))
print(io::IO, cp::Chr) = print(io, Char(cp))

codepoint(v::Char) = v%UInt32

unsafe_crc32c(a, n, crc) = ccall(:jl_crc32c, UInt32, (UInt32, Ptr{UInt8}, Csize_t), crc, a, n)

function sizeof(str::SubString{T}) where {T<:Str}
    is_multi(str) || return str.endof
    str.endof == 0 && return 0
    _nextind(CodeUnitMulti(), str.string, str.offset + str.endof) - str.offset - 1
end

occurs_in(str::String, hay::String) = contains(hay, str)
occurs_in(chr::Char,   hay::String) = contains(hay, string(chr))
Base.contains(hay::AbstractString, str::Str)     = occurs_in(str, hay)
Base.contains(hay::Str, str::AbstractString)     = occurs_in(str, hay)
Base.contains(hay::Str, str::Str)                = occurs_in(str, hay)
Base.contains(hay::AbstractString, chr::AbsChar) = occurs_in(chr, hay)
Base.contains(hay::AbstractString, pat::Regex)   = occurs_in(pat, hay)

const utf8crc = Base.crc32c

import Base: isalnum, isgraph, islower, isupper
const is_alphanumeric = isalnum
const is_graphic      = isgraph
const is_lowercase    = islower
const is_uppercase    = isupper

# These are deprecated in v0.7
for sym in (:bin, :oct, :dec, :hex)
    @eval import Base:$sym
    @eval ($sym)(x::Chr, p::Int) = ($sym)(codepoint(x), p, false)
    @eval ($sym)(x::Chr)         = ($sym)(codepoint(x), 1, false)
end

function repeat(ch::Char, cnt::Integer)
    cnt > 1 && return String(_repeat(CodeUnitMulti(), UTF8CSE, ch%UInt32, cnt))
    cnt < 0 && repeaterr(cnt)
    cnt == 0 ? empty_string : string(Char(ch%UInt32))
end

@noinline index_error(s::AbstractString, i::Integer) =
    throw(UnicodeError(UTF_ERR_INVALID_INDEX, Int(i), codeunit(s, i)))

function thisind(str::String, pos::Integer)
    @boundscheck 0 < pos <= ncodeunits(str) || boundserr(str, pos)
    pnt = pointer(str) + pos - 1
    pos - (checkcont(pnt) ? (checkcont(pnt - 1) ? (checkcont(pnt - 2) ? 3 : 2) : 1) : 0)
end

typemin(::Type{String}) = ""

macro preserve(args...)
    syms = args[1:end-1]
    for x in syms
        isa(x, Symbol) || error("Preserved variable must be a symbol")
    end
    #=
    s, r = gensym(), gensym()
    esc(quote
        $s = $(Expr(:gc_preserve_begin, syms...))
        $r = $(args[end])
        $(Expr(:gc_preserve_end, s))
        $r
        $(args[end])
    end)
    =#
    esc(quote ; $(args[end]) ; end)
end

Base.SubString(str::AbstractString, rng::UnitRange) = SubString(str, first(rng), last(rng))

Base.checkbounds(::Type{Bool}, s::AbstractString, i::Integer) = 1 <= i <= ncodeunits(s)

const _Chars = Union{<:Chr,Tuple{Vararg{<:Chr}},AbstractVector{<:Chr},Set{<:Chr},Base.Chars}

starts_with(str::AbstractString, chars::_Chars) = !is_empty(str) && first(str) in chars

ends_with(str::AbstractString, chars::_Chars) = !is_empty(str) && last(str) in chars

function Base.lstrip(s::AbstractString, chars::_Chars)
    e = lastindex(s)
    for (i, c) in pairs(s)
        c in chars || return SubString(s, i, e)
    end
    SubString(s, e+1, e)
end

function Base.rstrip(s::AbstractString, chars::_Chars)
    r = RevString(s)
    i = start(r)
    while !done(r,i)
        c, j = next(r,i)
        c in chars || return s[1:end-i+1]
        i = j
    end
    s[1:0]
end

Base.strip(s::AbstractString, chars::_Chars) = lstrip(rstrip(s, chars), chars)


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Test
using Behavior.Gherkin: Given, When, ScenarioStep, Scenario, Background
using Behavior: FromMacroStepDefinitionMatcher, findstepdefinition

@testset "Parameters           " begin
    @testset "Matching against parameters; Definition has one String parameter; has value bar" begin
        given = Given("some bar")
        stepdef_matcher = FromMacroStepDefinitionMatcher("""
            using Behavior: @given

            @given("some {String}") do context, v end
        """)

        stepdefinitionmatch = findstepdefinition(stepdef_matcher, given)

        @test stepdefinitionmatch.variables == ["bar"]
    end

    @testset "Matching against parameters; Definition has one empty String parameter; has value bar" begin
        given = Given("some bar")
        stepdef_matcher = FromMacroStepDefinitionMatcher("""
            using Behavior: @given

            @given("some {}") do context, v end
        """)

        stepdefinitionmatch = findstepdefinition(stepdef_matcher, given)

        @test stepdefinitionmatch.variables == ["bar"]
    end

    @testset "Matching against parameters; Definition has one String parameter; has value baz" begin
        given = Given("some baz")
        stepdef_matcher = FromMacroStepDefinitionMatcher("""
            using Behavior: @given

            @given("some {String}") do context, v end
        """)

        stepdefinitionmatch = findstepdefinition(stepdef_matcher, given)

        @test stepdefinitionmatch.variables == ["baz"]
    end

    @testset "Matching against parameters; Definition has two String parameters;  has values bar and fnord" begin
        given = Given("some bar and fnord")
        stepdef_matcher = FromMacroStepDefinitionMatcher("""
            using Behavior: @given

            @given("some {String} and {String}") do context, v1, v2 end
        """)

        stepdefinitionmatch = findstepdefinition(stepdef_matcher, given)

        @test stepdefinitionmatch.variables == ["bar", "fnord"]
    end

    @testset "Scenario step has a variable foo; Args has :foo => bar" begin
        stepdefmatcher = FromMacroStepDefinitionMatcher("""
            using Behavior: @given

            @given("some value {String}") do context, v
                @expect v == "bar"
            end
        """)
        executor = Behavior.Executor(stepdefmatcher)

        given = Given("some value bar")
        scenario = Scenario("Description", String[], ScenarioStep[given])

        scenarioresult = Behavior.executescenario(executor, Background(), scenario)

        @test isa(scenarioresult.steps[1], Behavior.SuccessfulStepExecution)
    end

    @testset "Scenario step has two String parameters; Arguments are bar, fnord" begin
        stepdefmatcher = FromMacroStepDefinitionMatcher("""
            using Behavior: @given

            @given("some values {} and {}") do context, v1, v2
                @expect v1 == "bar"
                @expect v2 == "fnord"
            end
        """)
        executor = Behavior.Executor(stepdefmatcher)

        given = Given("some values bar and fnord")
        scenario = Scenario("Description", String[], ScenarioStep[given])

        scenarioresult = Behavior.executescenario(executor, Background(), scenario)

        @test isa(scenarioresult.steps[1], Behavior.SuccessfulStepExecution)
    end

    @testset "Typed parameters" begin
        @testset "Definition has one Int parameter; has value 123" begin
            given = Given("some value 123")
            stepdef_matcher = FromMacroStepDefinitionMatcher("""
                using Behavior: @given

                @given("some value {Int}") do context, v end
            """)

            stepdefinitionmatch = findstepdefinition(stepdef_matcher, given)

            @test stepdefinitionmatch.variables == [123]
        end

        @testset "Definition has one Float64 parameter; has value 234.0" begin
            given = Given("some value 234.0")
            stepdef_matcher = FromMacroStepDefinitionMatcher("""
                using Behavior: @given

                @given("some value {Float64}") do context, v end
            """)

            stepdefinitionmatch = findstepdefinition(stepdef_matcher, given)

            @test stepdefinitionmatch.variables == [234.0]
        end

        @testset "Definition has one Bool parameter; has value true" begin
            given = Given("some value true")
            stepdef_matcher = FromMacroStepDefinitionMatcher("""
                using Behavior: @given

                @given("some value {Bool}") do context, v end
            """)

            stepdefinitionmatch = findstepdefinition(stepdef_matcher, given)

            @test stepdefinitionmatch.variables == [true]
        end

        @testset "Definition has one Bool parameter; has value false" begin
            given = Given("some value false")
            stepdef_matcher = FromMacroStepDefinitionMatcher("""
                using Behavior: @given

                @given("some value {Bool}") do context, v end
            """)

            stepdefinitionmatch = findstepdefinition(stepdef_matcher, given)

            @test stepdefinitionmatch.variables == [false]
        end

        @testset "Definition has one Bool parameter; has value false" begin
            given = Given("some value false")
            stepdef_matcher = FromMacroStepDefinitionMatcher("""
                using Behavior: @given

                @given("some value {Bool}") do context, v end
            """)

            stepdefinitionmatch = findstepdefinition(stepdef_matcher, given)

            @test stepdefinitionmatch.variables == [false]
        end
    end
end


text:
module SampleData

using Celeste: Model, DeterministicVI
import Celeste: Synthetic
import Celeste.SDSSIO: RunCamcolField, load_field_images, load_field_catalog,
                       SDSSDataSet, SDSSBackground, SDSSPSFMap

using Distributions
using StaticArrays
import WCS
import FITSIO
import DataFrames

export empty_model_params,
       sample_ce, perturb_params,
       sample_star_fluxes, sample_galaxy_fluxes,
       gen_sample_star_dataset, gen_sample_galaxy_dataset,
       gen_two_body_dataset, gen_three_body_dataset, gen_n_body_dataset,
       make_elbo_args, true_star_init

const sample_star_fluxes = [
    4.451805E+03,1.491065E+03,2.264545E+03,2.027004E+03,1.846822E+04]
const sample_galaxy_fluxes = [
    1.377666E+01, 5.635334E+01, 1.258656E+02,
    1.884264E+02, 2.351820E+02] * 100  # 1x wasn't bright enough

# A world coordinate system where the world and pixel coordinates are the same.
const wcs_id = WCS.WCSTransform(2,
                    cd = Float64[1 0; 0 1],
                    ctype = ["none", "none"],
                    crpix = Float64[1, 1],
                    crval = Float64[1, 1]);

# globals to hold already-loaded test images and catalogs
const DATADIR = joinpath(Pkg.dir("Celeste"), "test", "data")
const DATASET = SDSSDataSet(DATADIR)
const SDSS_FIELD_IMAGES =
    Dict{RunCamcolField,Vector{Image{SDSSBackground,SDSSPSFMap}}}()
const SDSS_FIELD_CATALOGS = Dict{RunCamcolField, Vector{CatalogEntry}}()

# Fetch sdss data for field `rcf` if not already present. See DATADIR/Makefile
# for possible make targets.
function fetch_sdss_data(rcf, make_target)
    # ensure images and catalog are downloaded
    wd = pwd()
    cd(DATADIR)
    make_output = readstring(`make $(make_target) RUN=$(rcf.run) CAMCOL=$(rcf.camcol) FIELD=$(rcf.field)`)

    # only print something if we actually ran a download
    if !startswith(make_output, "make: Nothing to be done for")
        print(make_output)
    end

    cd(wd)
end
fetch_sdss_data(run, camcol, field, make_target) =
    fetch_sdss_data(RunCamcolField(run, camcol, field), make_target)

"""
    get_sdss_images(run, camcol, field) -> Vector{<:Image}

Return lazily-loaded images from the given SDSS field.
"""
function get_sdss_images(run, camcol, field)
    rcf = RunCamcolField(run, camcol, field)
    if !haskey(SDSS_FIELD_IMAGES, rcf)
        fetch_sdss_data(rcf, "all")
        SDSS_FIELD_IMAGES[rcf] = load_field_images(DATASET, rcf)
    end
    return SDSS_FIELD_IMAGES[rcf]
end

"""
    get_sdss_catalog(run, camcol, field) -> Vector{CatalogEntry}

Return lazily-loaded catalog for the given SDSS field.
"""
function get_sdss_catalog(run, camcol, field)
    rcf = RunCamcolField(run, camcol, field)
    if !haskey(SDSS_FIELD_CATALOGS, rcf)
        fetch_sdss_data(rcf, "photoobj")
        SDSS_FIELD_CATALOGS[rcf] = load_field_catalog(DATASET, rcf)
    end
    return SDSS_FIELD_CATALOGS[rcf]
end

"""
Turn a images and vector of catalog entries into elbo arguments
that can be used with Celeste.
"""
function make_elbo_args(images::Vector{<:Image},
                        catalog::Vector{CatalogEntry};
                        active_source=-1,
                        patch_radius_pix::Float64=NaN,
                        include_kl=true)
    patches = Model.get_sky_patches(images,
                                    catalog,
                                    radius_override_pix=patch_radius_pix)
    S = length(catalog)
    active_sources = active_source > 0 ? [active_source] :
                                          S <= 3 ? collect(1:S) : [1,2,3]
    ElboArgs(images, patches, active_sources; include_kl=include_kl)
end

function empty_model_params(S::Int)
    vp = [DeterministicVI.generic_init_source([ 0., 0. ]) for s in 1:S]
    ElboArgs(Image[],
             vp,
             Matrix{ImagePatch}(S, 0),
             collect(1:S))
end

function sample_ce(pos, is_star::Bool)
    CatalogEntry(pos, is_star, sample_star_fluxes, sample_galaxy_fluxes,
        0.1, 0.7, pi/4, 4.0)
end

# for testing away from the truth, where derivatives != 0
function perturb_params(vp)
    for vs in vp
        vs[ids.is_star] = [ 0.4, 0.6 ]
        vs[ids.pos[1]] += .8
        vs[ids.pos[2]] -= .7
        vs[ids.flux_loc] -= log(10)
        vs[ids.flux_scale] *= 25.
        vs[ids.gal_frac_dev] += 0.05
        vs[ids.gal_axis_ratio] += 0.05
        vs[ids.gal_angle] += pi/10
        vs[ids.gal_radius_px] *= 1.2
        vs[ids.color_mean] += 0.5
        vs[ids.color_var] =  1e-1
    end
end

function cropped_sample_images(new_H, new_W)
    srand(1)
    sample_images = get_sdss_images(3900, 6, 269)
    images = deepcopy(sample_images)

    for img in images
        img.H = new_H
        img.W = new_W
        img.pixels = zeros(Float32, new_H, new_W)
        img.nelec_per_nmgy = img.nelec_per_nmgy[1:new_H]
        img.wcs = wcs_id
    end

    return images
end

function gen_sample_star_dataset(; perturb=true)
    images = cropped_sample_images(20, 23)
    catalog = [sample_ce([10.1, 12.2], true),]
    Synthetic.gen_images!(images, catalog)
    ea = make_elbo_args(images, catalog)

    vp = Vector{Float64}[DeterministicVI.catalog_init_source(ce) for ce in catalog]
    if perturb
        perturb_params(vp)
    end

    ea, vp, catalog
end

function gen_sample_galaxy_dataset(; perturb=true, include_kl=true)
    images = cropped_sample_images(20, 23)
    catalog = [sample_ce([8.5, 9.6], false),]
    Synthetic.gen_images!(images, catalog)
    ea = make_elbo_args(images, catalog; include_kl=include_kl)

    vp = Vector{Float64}[DeterministicVI.catalog_init_source(ce) for ce in catalog]
    if perturb
        perturb_params(vp)
    end

    ea, vp, catalog
end

# A small two-body dataset for quick unit testing.  These objects
# will be too close to be identifiable.
function gen_two_body_dataset(; perturb=true)
    images = cropped_sample_images(20, 23)
    catalog = [
        sample_ce([4.5, 3.6], false),
        sample_ce([10.1, 12.1], true)
    ]
    Synthetic.gen_images!(images, catalog)
    ea = make_elbo_args(images, catalog)

    vp = Vector{Float64}[DeterministicVI.catalog_init_source(ce) for ce in catalog]
    if perturb
        perturb_params(vp)
    end

    ea, vp, catalog
end

function gen_three_body_dataset(; perturb=true)
    images = cropped_sample_images(112, 238)
    catalog = [
        sample_ce([4.5, 3.6], false),
        sample_ce([60.1, 82.2], true),
        sample_ce([71.3, 100.4], false),
    ];
    Synthetic.gen_images!(images, catalog);
    ea = make_elbo_args(images, catalog);

    vp = Vector{Float64}[DeterministicVI.catalog_init_source(ce) for ce in catalog]
    if perturb
        perturb_params(vp)
    end

    ea, vp, catalog
end

"""
Generate a large dataset with S randomly placed bodies and non-constant
background.
"""
function gen_n_body_dataset(S::Int; patch_pixel_radius=20., perturb=true)
    images = cropped_sample_images(900, 1000)

    fluxes = [4.451805E+03,1.491065E+03,2.264545E+03,2.027004E+03,1.846822E+04]
    locations = rand(2, S) .* [900, 1000]
    world_locations = WCS.pix_to_world(images[3].wcs, locations)
    catalog = CatalogEntry[CatalogEntry(world_locations[:, s], true,
            fluxes, fluxes, 0.1, 0.7, pi/4, 4.0) for s in 1:S];

    Synthetic.gen_images!(images, catalog);

    # Make non-constant background.
    for b=1:5
        images[b].nelec_per_nmgy = fill(images[b].nelec_per_nmgy[1], images[b].H)
        images[b].sky = SDSSBackground(fill(images[b].sky[1,1], images[b].H, images[b].W),
                                       collect(1:images[b].H), collect(1:images[b].W),
                                       ones(images[b].H))
    end

    ea = make_elbo_args(
        images, catalog, patch_radius_pix=patch_pixel_radius)

    vp = Vector{Float64}[DeterministicVI.catalog_init_source(ce) for ce in catalog]
    if perturb
        perturb_params(vp)
    end

    ea, vp, catalog
end

function true_star_init()
    ea, vp, catalog = gen_sample_star_dataset(perturb=false)

    vp[1][ids.is_star] = [ 1.0 - 1e-4, 1e-4 ]
    vp[1][ids.flux_scale] = 1e-4
    vp[1][ids.flux_loc] = log(sample_star_fluxes[3]) - 0.5 * vp[1][ids.flux_scale]
    vp[1][ids.color_var] = 1e-4

    ea, vp, catalog
end

end # End module


text:
## LEARNING CURVES

"""
    curve = learning_curve(mach; resolution=30,
                                 resampling=Holdout(),
                                 repeats=1,
                                 measure=default_measure(machine.model),
                                 rows=nothing,
                                 weights=nothing,
                                 operation=nothing,
                                 range=nothing,
                                 acceleration=default_resource(),
                                 acceleration_grid=CPU1(),
                                 rngs=nothing,
                                 rng_name=nothing)

Given a supervised machine `mach`, returns a named tuple of objects
suitable for generating a plot of performance estimates, as a function
of the single hyperparameter specified in `range`. The tuple `curve`
has the following keys: `:parameter_name`, `:parameter_scale`,
`:parameter_values`, `:measurements`.

To generate multiple curves for a `model` with a random number
generator (RNG) as a hyperparameter, specify the name, `rng_name`, of
the (possibly nested) RNG field, and a vector `rngs` of RNG's, one for
each curve. Alternatively, set `rngs` to the number of curves desired,
in which case RNG's are automatically generated. The individual curve
computations can be distributed across multiple processes using
`acceleration=CPUProcesses()` or `acceleration=CPUThreads()`. See the
second example below for a demonstration.

```julia
X, y = @load_boston;
atom = @load RidgeRegressor pkg=MultivariateStats
ensemble = EnsembleModel(atom=atom, n=1000)
mach = machine(ensemble, X, y)
r_lambda = range(ensemble, :(atom.lambda), lower=10, upper=500, scale=:log10)
curve = learning_curve(mach; range=r_lambda, resampling=CV(), measure=mav)
using Plots
plot(curve.parameter_values,
     curve.measurements,
     xlab=curve.parameter_name,
     xscale=curve.parameter_scale,
     ylab = "CV estimate of RMS error")
```

If using a `Holdout()` `resampling` strategy (with no shuffling) and
if the specified hyperparameter is the number of iterations in some
iterative model (and that model has an appropriately overloaded
`MLJModelInterface.update` method) then training is not restarted from scratch
for each increment of the parameter, ie the model is trained
progressively.

```julia
atom.lambda=200
r_n = range(ensemble, :n, lower=1, upper=250)
curves = learning_curve(mach; range=r_n, verbosity=0, rng_name=:rng, rngs=3)
plot!(curves.parameter_values,
     curves.measurements,
     xlab=curves.parameter_name,
     ylab="Holdout estimate of RMS error")

```
    learning_curve(model::Supervised, X, y; kwargs...)
    learning_curve(model::Supervised, X, y, w; kwargs...)

Plot a learning curve (or curves) directly, without first constructing
a machine.

### Summary of key-word options

- `resolution` - number of points generated from `range` (number model
  evaluations); default is `30`

- `acceleration` - parallelization option for passing to `evaluate!`;
  an instance of `CPU1`, `CPUProcesses` or `CPUThreads` from the
  `ComputationalResources.jl`; default is `default_resource()`

- `acceleration_grid` - parallelization option for distributing each
  performancde evaluation

- `rngs` - for specifying random number generator(s) to be passed to
  the model (see above)

- `rng_name` - name of the model hyper-parameter representing a random
  number generator (see above); possibly nested

Other key-word options are documented at [`TunedModel`](@ref).

"""
learning_curve(mach::Machine{<:Supervised}; kwargs...) =
    learning_curve(mach.model, mach.args...; kwargs...)

function learning_curve(model::Supervised, args...;
                        resolution=30,
                        resampling=Holdout(),
                        weights=nothing,
                        measures=nothing,
                        measure=measures,
                        rows=nothing,
                        operation=nothing,
                        ranges::Union{Nothing,ParamRange}=nothing,
                        range::Union{Nothing,ParamRange},
                        repeats=1,
                        acceleration=default_resource(),
                        acceleration_grid=CPU1(),
                        verbosity=1,
                        rngs=nothing,
                        rng_name=nothing,
                        check_measure=true)

    range !== nothing || error("No param range specified. Use range=... ")

    if rngs != nothing
        rng_name == nothing &&
            error("Having specified `rngs=...`, you must specify "*
                  "`rng_name=...` also. ")
        if rngs isa Integer
            rngs = MersenneTwister.(1:rngs)
        elseif rngs isa AbstractRNG
            rngs = [rngs, ]
        elseif !(rngs isa AbstractVector{<:AbstractRNG})
            error("`rng` must have type `Integer` , `AbstractRNG` or "*
                  "`AbstractVector{<:AbstractRNG}`. ")
        end
    end

    if (acceleration isa CPUProcesses &&
        acceleration_grid isa CPUProcesses)
        message =
        "The combination acceleration=$(acceleration) and"*
        " acceleration_grid=$(acceleration_grid) is"*
        "  not generally optimal. You may want to consider setting"*
        " `acceleration = CPUProcesses()` and"*
        " `acceleration_grid = CPUThreads()`."
       @warn message
     end
    if (acceleration isa CPUThreads &&
        acceleration_grid isa CPUProcesses)
        message =
        "The combination acceleration=$(acceleration) and"*
        " acceleration_grid=$(acceleration_grid) isn't supported. \n"*
        "Resetting to"*
        " `acceleration = CPUProcesses()` and"*
        " `acceleration_grid = CPUThreads()`."
        @warn message
        acceleration = CPUProcesses()
        acceleration_grid = CPUThreads()
     end
   _acceleration = _process_accel_settings(acceleration)
    tuned_model = TunedModel(model=model,
                             range=range,
                             tuning=Grid(resolution=resolution,
                                         shuffle=false),
                             resampling=resampling,
                             operation=operation,
                             measure=measure,
                             train_best=false,
                             weights=weights,
                             repeats=repeats,
                             acceleration=acceleration_grid)

    tuned = machine(tuned_model, args...)

    results =
        _tuning_results(rngs, _acceleration, tuned, rows, rng_name, verbosity)

    parameter_name=results.parameter_names[1]
    parameter_scale=results.parameter_scales[1]
    parameter_values=[results.parameter_values[:, 1]...]
    measurements = results.measurements

    return (parameter_name=parameter_name,
            parameter_scale=parameter_scale,
            parameter_values=parameter_values,
            measurements=measurements)
end

_collate(plotting1, plotting2) =
    merge(plotting1,
          (measurements=hcat(plotting1.measurements,
                             plotting2.measurements),))

# fallback:
#_tuning_results(rngs, acceleration, tuned, rows, rngs_name, verbosity) =
#    error("acceleration=$acceleration unsupported. ")

# single curve:
_tuning_results(rngs::Nothing,
                acceleration,
                tuned,
                rows,
                rngs_name,
                verbosity) = _single_curve(tuned, rows, verbosity)

function _single_curve(tuned, rows, verbosity)
    fit!(tuned, rows=rows, verbosity=verbosity, force=true)
    tuned.report.plotting
end

# CPU1:
function _tuning_results(rngs::AbstractVector,
                         acceleration::CPU1,
                         tuned,
                         rows,
                         rng_name,
                         verbosity)

    old_rng = recursive_getproperty(tuned.model.model, rng_name)
    n_rngs = length(rngs)

    p = Progress(n_rngs,
         dt = 0,
         desc = "Evaluating Learning curve with $(n_rngs) rngs: ",
         barglyphs = BarGlyphs("[=> ]"),
         barlen = 18,
         color = :yellow)

    verbosity < 1 ||  update!(p,0)

    ret = mapreduce(_collate, rngs) do rng
              recursive_setproperty!(tuned.model.model, rng_name, rng)
              fit!(tuned, rows=rows, verbosity=verbosity-1, force=true)
              r =tuned.report.plotting
              verbosity < 1 || begin
                      p.counter += 1
                      ProgressMeter.updateProgress!(p)
                    end
              r
         end

    recursive_setproperty!(tuned.model.model, rng_name, old_rng)

    return ret
end

# CPUProcesses:
function _tuning_results(rngs::AbstractVector,
                         acceleration::CPUProcesses,
                         tuned,
                         rows,
                         rng_name,
                         verbosity)

    old_rng = recursive_getproperty(tuned.model.model, rng_name)
    n_rngs = length(rngs)

    ret = @sync begin

    p = Progress(n_rngs,
         dt = 0,
         desc = "Evaluating Learning curve with $(n_rngs) rngs: ",
         barglyphs = BarGlyphs("[=> ]"),
         barlen = 18,
         color = :yellow)

    channel = RemoteChannel(()->Channel{Bool}(min(1000, n_rngs)), 1)

    # printing the progress bar
    verbosity < 1 || begin
                update!(p,0)
                @async while take!(channel)
                    p.counter +=1
                    ProgressMeter.updateProgress!(p)
                end
             end

    ret_ = @distributed (_collate) for rng in rngs
            recursive_setproperty!(tuned.model.model, rng_name, rng)
            fit!(tuned, rows=rows, verbosity=verbosity-1, force=true)
            r=tuned.report.plotting
            verbosity < 1 || put!(channel, true)
            r
        end

     verbosity < 1 || put!(channel, false)
     ret_
 end
    recursive_setproperty!(tuned.model.model, rng_name, old_rng)
    return ret
end

# CPUThreads:
@static if VERSION >= v"1.3.0-DEV.573"
    function _tuning_results(rngs::AbstractVector,
                             acceleration::CPUThreads,
                             tuned,
                             rows,
                             rng_name,
                             verbosity)

    n_threads = Threads.nthreads()
    if n_threads == 1
        return _tuning_results(rngs,
                               CPU1(),
                               tuned,
                               rows,
                               rng_name,
                               verbosity)
    end

    old_rng = recursive_getproperty(tuned.model.model, rng_name)
    n_rngs = length(rngs)
    ntasks = acceleration.settings
    partitions = MLJBase.chunks(1:n_rngs, ntasks)

     p = Progress(n_rngs,
         dt = 0,
         desc = "Evaluating Learning curve with $(n_rngs) rngs: ",
         barglyphs = BarGlyphs("[=> ]"),
         barlen = 18,
         color = :yellow)

    ch = Channel{Bool}(length(partitions))

    ret_ = Vector(undef, length(partitions))

    @sync begin
        verbosity < 1 || begin
                          update!(p,0)
                          @async while take!(ch)
                            p.counter +=1
                            ProgressMeter.updateProgress!(p)
                          end
                        end

        # One t_tuned per task
        ## deepcopy of model is because other threads can still change the state
        ## of tuned.model.model
        tmachs = [tuned,
                  [machine(TunedModel(model = deepcopy(tuned.model.model),
                                      range=deepcopy(tuned.model.range),
                                      tuning=tuned.model.tuning,
                                      resampling=tuned.model.resampling,
                                      operation=tuned.model.operation,
                                      measure=tuned.model.measure,
                                      train_best=tuned.model.train_best,
                                      weights=tuned.model.weights,
                                      repeats=tuned.model.repeats,
                                      acceleration=tuned.model.acceleration),
                           tuned.args...) for _ in 2:length(partitions)]...]
        @sync for (i,rng_part) in enumerate(partitions)
            Threads.@spawn begin
                ret_[i] = mapreduce(_collate, rng_part) do k
                    recursive_setproperty!(tmachs[i].model.model,
                                           rng_name, rngs[k])
                    fit!(tmachs[i], rows=rows,
                         verbosity=verbosity-1, force=true)
                    verbosity < 1 || put!(ch, true)
                    tmachs[i].report.plotting
                end
            end
        end
        verbosity < 1 || put!(ch, false)
    end

        ret =  reduce(_collate, ret_)
        recursive_setproperty!(tuned.model.model, rng_name, old_rng)
        return ret
    end

end

max_stars_repo_name: alan-turing-institute/MLJTuning

max_issues_repo_name: alan-turing-institute/MLJTuning

max_forks_repo_name: alan-turing-institute/MLJTuning


text:
""
function run_opf_bf_dvorkin(file, model_type::Type{T}, optimizer; kwargs...) where T <: _PM.AbstractBFModel
    return _PM.run_model(file, model_type, optimizer, build_opf_bf_dvorkin; kwargs...)
end

""
function build_opf_bf_dvorkin(pm::_PM.AbstractPowerModel)
    _PM.variable_bus_voltage(pm)
    _PM.variable_gen_power(pm)
    _PM.variable_branch_power(pm)
    _PM.variable_branch_current(pm)
    # variable_dcline_power(pm)

    _PM.objective_min_fuel_and_flow_cost(pm)

    _PM.constraint_model_current(pm)

    for i in _PM.ids(pm, :ref_buses)
        _PM.constraint_theta_ref(pm, i)
    end

    for i in _PM.ids(pm, :bus)
        _PM.constraint_power_balance(pm, i)
    end

    for i in _PM.ids(pm, :branch)
        _PM.constraint_power_losses(pm, i)
        _PM.constraint_voltage_magnitude_difference(pm, i)

        _PM.constraint_voltage_angle_difference(pm, i)

        _PM.constraint_thermal_limit_from(pm, i)
        _PM.constraint_thermal_limit_to(pm, i)
    end

    # for i in ids(pm, :dcline)
    #     constraint_dcline_power_losses(pm, i)
    # end
end

""
function run_opf_bf_dvorkin_cc(file, model_type::Type{T}, optimizer; kwargs...) where T <: _PM.AbstractBFModel
    return _PM.run_model(file, model_type, optimizer, build_opf_bf_dvorkin_cc; kwargs...)
end

""
function build_opf_bf_dvorkin_cc(pm::_PM.AbstractPowerModel)
    _PM.variable_bus_voltage(pm, bounded=false)
    _PM.variable_gen_power(pm, bounded=true)
    _PM.variable_branch_power(pm, bounded=false)
    _PM.variable_branch_current(pm)
    variable_gen_power_response(pm)
    # variable_dcline_power(pm)

    _PM.objective_min_fuel_and_flow_cost(pm)

    _PM.constraint_model_current(pm)

    for i in _PM.ids(pm, :ref_buses)
        _PM.constraint_theta_ref(pm, i)
    end

    for i in _PM.ids(pm, :bus)
        _PM.constraint_power_balance(pm, i)
    end
    for i in _PM.ids(pm, :gen)
        constraint_gen_bounds_cc(pm,i)
    end
    for i in _PM.ids(pm, :branch)
        _PM.constraint_power_losses(pm, i)
        _PM.constraint_voltage_magnitude_difference(pm, i)

        # _PM.constraint_voltage_angle_difference(pm, i)

        # _PM.constraint_thermal_limit_from(pm, i)
        # _PM.constraint_thermal_limit_to(pm, i)
    end

    # for i in ids(pm, :dcline)
    #     constraint_dcline_power_losses(pm, i)
    # end
end


text:
module RunSettings
    #
    # This module contains things needed to control one run e.g. the output destination, number of households to use andd so on.
    #
    using Parameters

    export
        Settings,
        DEFAULT_SETTINGS,
        MT_Routing,
        uc_full,
        lmt_full,
        modelled_phase_in

    @enum MT_Routing uc_full lmt_full modelled_phase_in
    @with_kw mutable struct Settings
        uid :: Int = 1 # placeholder for maybe a user somewhere
        run_name :: String = "default_run"
        start_year :: Int = 2015
        end_year :: Int = 2018
        scotland_full :: Bool = true
        weighted :: Bool = false
        household_name = "model_households_scotland"
        people_name    = "model_people_scotland"

        num_households :: Int = 0
        num_people :: Int = 0
        prices_file = "indexes_sep_1_2021.tab"
        to_y :: Int = 2021
        to_q :: Int = 2
        output_dir :: String = joinpath(tempdir(),"output")
        # ... and so on
        means_tested_routing :: MT_Routing = uc_full
    end

    const DEFAULT_SETTINGS = Settings()

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Pre-defined formats for LaTeX tables.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export tf_latex_default, tf_latex_simple, tf_latex_modern, tf_latex_booktabs

const tf_latex_default = LatexTableFormat()

const tf_latex_simple = LatexTableFormat(
    top_line       = "\\hline",
    header_line    = "\\hline",
    mid_line       = "\\hline",
    bottom_line    = "\\hline",
    left_vline     = "|",
    mid_vline      = "|",
    right_vline    = "|",
    header_envs    = ["textbf"],
    subheader_envs = ["texttt"]
   )

const tf_latex_modern = LatexTableFormat(
    top_line       = "\\noalign{\\hrule height 2pt}",
    header_line    = "\\noalign{\\hrule height 2pt}",
    mid_line       = "\\noalign{\\hrule height 1pt}",
    bottom_line    = "\\noalign{\\hrule height 2pt}",
    left_vline     = "!{\\vrule width 2pt}",
    mid_vline      = "!{\\vrule width 1pt}",
    right_vline    = "!{\\vrule width 2pt}",
    header_envs    = ["textbf"],
    subheader_envs = ["texttt"]
   )

const tf_latex_booktabs = LatexTableFormat(
    top_line       = "\\toprule",
    header_line    = "\\midrule",
    mid_line       = "\\midrule",
    bottom_line    = "\\bottomrule",
    left_vline     = "",
    mid_vline      = "",
    right_vline    = "",
    header_envs    = ["textbf"],
    subheader_envs = ["texttt"],
   )


text:
# From Program 10
using LinearAlgebra

function meshgrid(x, y)
    x = reshape(x, 1, length(x))
    y = reshape(y, length(y), 1)
    X = repeat(x, length(y), 1)
    Y = repeat(y, 1, length(x))
    return X, Y
end


text:
struct DualSolution{dev, T, TVCPU <: AbstractVector{T}, TVGPU <: Union{Nothing, CuVector{T}}}
    cpu::TVCPU
    gpu::TVGPU
end
GPUUtils.whichdevice(::DualSolution{:CPU}) = CPU()
GPUUtils.whichdevice(::DualSolution{:GPU}) = GPU()

DualSolution(x) = DualSolution{:CPU}(x)
DualSolution{:CPU}(x) = DualSolution{:CPU, eltype(x), typeof(x), Nothing}(x, nothing)
function DualSolution{:GPU}(x)
    cux = CuArray(x)
    DualSolution{:GPU, eltype(x), typeof(x), typeof(cux)}(x, cux)
end
Base.getindex(d::DualSolution, i...) = d.cpu[i...]
togpu!(d::DualSolution{:GPU}) = copyto!(d.gpu, d.cpu)
tocpu!(d::DualSolution{:GPU}) = copyto!(d.cpu, d.gpu)
togpu!(d::DualSolution{:CPU}) = nothing
tocpu!(d::DualSolution{:CPU}) = nothing

function getdualterm(p, q, ρi, σ, x1, x, λi, ji::Tuple)
    j, i = ji
    σj = σ[j]
    Lj, Uj = minus_plus(x1[j], σj)
    xj = x[j]
    Ujxj = Uj - xj
    xjLj = xj - Lj
    Δ = ρi*σj/4
    return λi*((p[j,i] + Δ)/Ujxj + (q[j,i] + Δ)/xjLj)
end
function getdualterm(p0, q0, σ, x1, x, j::Int)
    σj = σ[j]
    Lj, Uj = minus_plus(x1[j], σj)
    return p0[j]/(Uj-x[j]) + q0[j]/(x[j]-Lj)
end
struct DualObjVal{TPD, TXUpdater, TDualSol}
    pd::TPD
    x_updater::TXUpdater
    λ::TDualSol
end
function (dobj::DualObjVal)(λ)
    @unpack pd = dobj
    @unpack r, r0, p0, q0, σ, x1, x, p, q, ρ = pd
    T = eltype(p)
    # Updates x to the Lagrangian minimizer for the input λ
    #println("Updating x")
    dobj.λ.cpu .= λ
    dobj.x_updater(dobj.λ)
    nv, nc = size(p)
    φ0 = r0 + dot(λ, r)
    φ = computeφ(φ0, p0, q0, p, q, ρ, σ, dobj.λ, x1, x)
    return -φ
end
function computeφ(φ0, p0::AbstractVector{T}, q0, p, q, ρ, σ, λ, x1, x) where {T}
    nv, nc = size(p)
    φ = φ0
    φ += tmapreduce(+, 1:nv, init = zero(T)) do j
        getdualterm(p0, q0, σ, x1, x, j)
    end
    for i in 1:nc
        λi = λ.cpu[i]
        ρi = ρ[i]
        φ += tmapreduce(+, 1:nv, init = zero(T)) do j
            getdualterm(p, q, ρi, σ, x1, x, λi, (j, i))
        end
    end
    return φ
end

function computeφ(φ0, p0::CuVector{T}, q0, p, q, ρ, σ, λ, x1, x, ::Val{blocksize} = Val(80), ::Val{threads} = Val(256)) where {T, blocksize, threads}
    nv, nc = size(p)

    φ = φ0
    result = similar(x, T, (blocksize,))
    args = (x, x1, p0, q0, σ, result, Val(threads))
    @cuda blocks = blocksize threads = threads computeφ_kernel1(args...)
    CUDAnative.synchronize()
    φ += sum(Array(result))

    for i in 1:nc
        args = (x, x1, p, q, ρ[i], σ, λ.cpu[i], i, result, Val(threads))
        @cuda blocks = blocksize threads = threads computeφ_kernel2(args...)
        CUDAnative.synchronize()
        φ += sum(Array(result))
    end

    return φ
end

function computeφ_kernel1(x::AbstractVector{T}, x1, p0, q0, σ, result, ::Val{LMEM}) where {T, LMEM}
    @mapreduce_block(j, length(x), +, T, LMEM, result, begin
        getdualterm(p0, q0, σ, x1, x, j)
    end)

    return
end

function computeφ_kernel2(x::AbstractVector{T}, x1, p, q, ρi, σ, λi, i, result, ::Val{LMEM}) where {T, LMEM}
    @mapreduce_block(j, length(x), +, T, LMEM, result, begin
        getdualterm(p, q, ρi, σ, x1, x, λi, (j, i))
    end)

    return
end

struct DualObjGrad{TPD, TXUpdater, TDualSol}
    pd::TPD
    x_updater::TXUpdater
    λ::TDualSol
end
function (dgrad::DualObjGrad)(∇φ::AbstractVector{T}, λ) where {T}
    compute_grad!(whichdevice(dgrad.pd), dgrad, ∇φ, λ)
end
function compute_grad!(::CPU, dgrad, ∇φ, λ)
    T = eltype(λ)
    @unpack pd = dgrad
    @unpack r, r0, p, q, σ, x1, x, ρ = pd
    # Updates x to the Lagrangian minimizer for the input λ
    dgrad.λ.cpu .= λ
    dgrad.x_updater(dgrad.λ)
    nv, nc = size(p)
    # Negate since we have a maximization problem
    for i in 1:nc
        ∇φ[i] = -r[i]
        ρi = ρ[i]
        ∇φ[i] -= tmapreduce(+, 1:nv, init = zero(T)) do j
            getgradterm(x, x1, p, q, ρi, σ, (j, i))
        end
    end
    return ∇φ
end
function compute_grad!(::GPU, dgrad, ∇φ::AbstractVector{T}, λ, ::Val{blocksize} = Val(80), ::Val{threads} = Val(256)) where {T, blocksize, threads}
    @unpack pd = dgrad
    @unpack r, r0, p, q, σ, x1, x, ρ = pd
    # Updates x to the Lagrangian minimizer for the input λ
    dgrad.λ.cpu .= λ
    dgrad.x_updater(dgrad.λ)
    nv, nc = size(p)
    # Negate since we have a maximization problem
    for i in 1:nc
        ∇φ[i] = -r[i]
        result = similar(x, T, (blocksize,))
        args = (x, x1, p, q, ρ[i], i, σ, result, Val(threads))
        @cuda blocks = blocksize threads = threads compute_grad_kernel(args...)
        CUDAnative.synchronize()
        ∇φ[i] -= sum(Array(result))
    end
    return ∇φ
end
function compute_grad_kernel(x::AbstractVector{T}, x1, p, q, ρi, i, σ, result, ::Val{LMEM}) where {T, LMEM}
    @mapreduce_block(j, length(x), +, T, LMEM, result, begin
        getgradterm(x, x1, p, q, ρi, σ, (j, i))
    end)

    return
end
function getgradterm(x, x1, p, q, ρi, σ, ji::Tuple)
    j, i = ji
    σj = σ[j]
    Lj, Uj = minus_plus(x1[j], σj)
    Ujxj = Uj - x[j]
    xjLj = x[j] - Lj
    Δ = ρi*σj/4
    return (p[j,i] + Δ)/Ujxj + (q[j,i] + Δ)/xjLj
end


text:
markproperties_spec =
[
	(:enter, VegaMarkPropertySet, nothing),
	(:exit, VegaMarkPropertySet, nothing),
	(:update, VegaMarkPropertySet, nothing),
	(:hover, VegaMarkPropertySet, nothing)
]

primitivefactory(:VegaMarkProperties, markproperties_spec)

max_stars_repo_name: jpthor/vega


text:
using ChromNet
using Memoize

export project_groupgm, edge_groupgm

"""
This projects a group graphical model down to a network among individual variables.

Projection is necessary in order to compare with standard modeling methods that
only connect individual datasets. The group filter determines which groups should
be considered for the projection (this may exclude large groups), and then the
edge merge function determines how to merge two edge values.
"""
function project_groupgm(G::AbstractMatrix, header::AbstractArray, groups::AbstractArray, groupFilter=nothing, edgeMerge::Function=(x,y)->abs(x) > abs(y) ? x : y)

    # by default only allow groups targeting the same factor
    @memoize function same_factor_filter(x)
        parts = split(x[2], ' ')
        referenceUid = ChromNetPaper.id2uniprot(parts[1])
        for i in 2:length(parts)
            if referenceUid != ChromNetPaper.id2uniprot(parts[i])
                return false
            end
        end
        true
    end
    thresholdVal = 0.0
    if groupFilter == nothing
        groupFilter = same_factor_filter
    elseif typeof(groupFilter) <: AbstractFloat
        thresholdVal = groupFilter
        groupFilter = x->begin
            same_factor_filter(x) && x[1] <= thresholdVal
        end
    end

    # create a map from header ids to indexes
    indexMap = Dict{ASCIIString,Int64}()
    for i in 1:length(header)
        indexMap[header[i]] = i
    end

    # fill in our projected matrix
    X = NaN*zeros(length(header),length(header))
    indexGroups = [[indexMap[x] for x in split(g[2])] for g in groups]
    for i in 1:length(indexGroups)
        if !groupFilter(groups[i]) continue end

        # fill in diagonal
        if length(indexGroups[i]) == 1
            ind = indexGroups[i][1]
            X[ind,ind] = G[i,i]
        end

        for j in i+1:length(indexGroups)
            if !groupFilter(groups[j]) || issubset(indexGroups[i], indexGroups[j]) || issubset(indexGroups[j], indexGroups[i])
                continue
            end

            for gi in indexGroups[i], gj in indexGroups[j]
                X[gj,gi] = X[gi,gj] = isnan(X[gj,gi]) ? G[j,i] : edgeMerge(X[gj,gi], G[j,i])
            end
        end
    end

    X
end

function edge_groupgm(C, ids, invFunc=inv; checkGroup=nothing)
    groups = build_groups(C, ids)
    G,headerG = build_groupgm(invFunc(C), ids, groups)
    project_groupgm(G, ids, groups, checkGroup)
end


text:
# ---
# title: 897. Increasing Order Search Tree
# id: problem897
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Tree, Depth-first Search, Recursion
# link: <https://leetcode.com/problems/increasing-order-search-tree/description/>
# hidden: true
# ---
#
# Given the `root` of a binary search tree, rearrange the tree in **in-order**
# so that the leftmost node in the tree is now the root of the tree, and every
# node has no left child and only one right child.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/17/ex1.jpg)
#
#
#
#     Input: root = [5,3,6,2,4,null,8,1,null,null,null,7,9]
#     Output: [1,null,2,null,3,null,4,null,5,null,6,null,7,null,8,null,9]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/17/ex2.jpg)
#
#
#
#     Input: root = [5,1,7]
#     Output: [1,null,5,null,7]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the given tree will be in the range `[1, 100]`.
#   * `0 <= Node.val <= 1000`
#
#
## @lc code=start
using LeetCode

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