text:
abstract type AbstractEoSComponent end

abstract type AbstractEoSMixture end

const NothingOrT{T} = Union{Nothing,T}

struct BrusilovskyEoSComponent{T<:Number} <: AbstractEoSComponent
    # meta information
    name::String

    # physical parameters
    Pc::T  # critical pressure
    acentric_factor::T
    RTc::T   # R * critical temperature
    molar_mass::T  # [kg mol⁻¹] molar mass
    carbon_number::Int  # [dimless] number of carbons

    # eos parameters
    ac::T     # explicit coefficient of the eos a_c
    b::T      # explicit coefficient of the eos b
    c::T      # explicit coefficient of the eos c
    d::T      # explicit coefficient of the eos d
    Psi::T    # primary coefficient of the eos - \Psi

    function BrusilovskyEoSComponent{T}(
        ;
        name::AbstractString="No Name",
        critical_pressure::Number=NaN,
        critical_temperature::Number=NaN,
        acentric_factor::Number=NaN,
        Omegac::Number=NaN,
        Zc::Number=NaN,
        Psi::Number=NaN,
        molar_mass::Number=NaN,
        carbon_number::Integer=0,
        kw...
    ) where {T}
        alpha = Omegac^3
        beta = Zc + Omegac - 1.0
        ds = sqrt(Omegac - 0.75)
        sigma = -Zc + Omegac * (0.5 + ds)
        delta = -Zc + Omegac * (0.5 - ds)
        RTc = GAS_CONSTANT_SI * critical_temperature
        rtp = RTc / critical_pressure
        ac = alpha * RTc * rtp
        b = beta * rtp
        c = sigma * rtp
        d = delta * rtp
        new{T}(
            name,
            critical_pressure,
            acentric_factor,
            RTc,
            molar_mass,
            carbon_number,
            ac, b, c, d, Psi
        )
    end
end
BrusilovskyEoSComponent(; x...) = BrusilovskyEoSComponent{Float64}(; x...)

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

for func in (:molar_mass, :name, :carbon_number)
    expr = :($(func)(c::BrusilovskyEoSComponent) = getfield(c, $(QuoteNode(func))))
    eval(expr)
    eval(:(export $func))
end

#=
Mixture
=#

struct BrusilovskyEoSMixture{T} <: AbstractEoSMixture
    components::Vector{BrusilovskyEoSComponent{T}}

    eij::Matrix{T} # constant  thermal binary interaction coefficient
    gij::Matrix{T} # linear    thermal binary interaction coefficient
    hij::Matrix{T} # quadratic thermal binary interaction coefficient

    function BrusilovskyEoSMixture(
        ;
        components::AbstractVector{BrusilovskyEoSComponent{T}},
        constant::AbstractMatrix,
        linear::AbstractMatrix,
        quadratic::AbstractMatrix,
        kw...
    ) where {T}
        new{T}(components, constant, linear, quadratic)
    end
end

@inline Base.@propagate_inbounds function Base.getindex(
    mix::BrusilovskyEoSMixture,
    i::Integer
)
    return mix.components[i]
end

struct BrusilovskyThermoBuffer{T}
    matr::Matrix{T}
    vec1::Vector{T}
    vec2::Vector{T}
end

function BrusilovskyThermoBuffer{T}(n::Integer) where {T}
    matr = Matrix{T}(undef, n, n)
    vec1 = similar(matr, (n,))
    vec2 = similar(vec1)
    return BrusilovskyThermoBuffer{T}(matr, vec1, vec2)
end

BrusilovskyThermoBuffer(n::Integer) = BrusilovskyThermoBuffer{Float64}(n)

function BrusilovskyThermoBuffer(mix::BrusilovskyEoSMixture{T}) where {T}
    nc = ncomponents(mix)
    return BrusilovskyThermoBuffer{T}(nc)
end

"""
    thermo_buffer(mix)

Create a buffer for intermediate calculations of mixture thermodynamic properties.

See also: [`pressure`](@ref), [`log_c_activity`](@ref), [`log_c_activity!`](@ref),
[`log_c_activity_wj`](@ref), [`log_c_activity_wj!`](@ref)
"""
thermo_buffer(mix::BrusilovskyEoSMixture) = BrusilovskyThermoBuffer(mix)

max_stars_repo_name: vvpisarev/CubicEoS

max_issues_repo_name: vvpisarev/CubicEoS

max_forks_repo_name: vvpisarev/CubicEoS


text:
using CUDA: CUDA

function cpu_gpu_agree(f, data; atol=1e-7)
    CUDA.functional() && @test f(data) == to_cpu(f(to_gpu(data)))
end

function cpu_gpu_agree_approx(f, data; atol=1e-7)
    CUDA.functional() && @test f(data) ≈ to_cpu(f(to_gpu(data))) atol=atol
end

ext: jl

max_stars_repo_name: Juice-jl/DiscriminativeCircuits

max_issues_repo_name: Juice-jl/DiscriminativeCircuits

max_forks_repo_name: Juice-jl/DiscriminativeCircuits


text:
const assets_dir = joinpath(@__DIR__, "..", "assets")
const conversions = Dict{Symbol,Float64}(:bohr2ang => 0.52917721092)
conversions[:ang2bohr] = 1/conversions[:bohr2ang]
# include("abinit/constants.jl")
include("qe/constants.jl")
include("wannier90/constants.jl")
include("elk/constants.jl")


text:
using Replay

repl_script = """
using PythonCall
jlv = [1,2,3]
pyv = Py(jlv)
pyv.append(4)
pyv
jlv
"""

replay(
    repl_script,
    stdout,
    julia_project=@__DIR__,
    use_ghostwriter=true,
    cmd="--color=yes",
)


text:
function fstep(state::Union{SimState, Nothing}, agent::Agent)
    #print(string(agent.id,"\n"))
end

a1 = Agent(fstep,nothing)
a2 = Agent(fstep,nothing)

f = Field2D(100.0,100.0,0.7,true)

setObjectLocation!(f,a1,Real2D(10.5,10.0))

@test getObjectLocation(f,a1) == Real2D(10.5,10.0)

setObjectLocation!(f,a1,Real2D(10.0,10.0))

@test getObjectLocation(f,a1) == Real2D(10.0,10.0)

setObjectLocation!(f,a2,Real2D(11.0,11.0))

@test getObjectLocation(f,a2) == Real2D(11.0,11.0)

@test length(getAllObjects(f)) == 2

clean!(f)

@test length(getAllObjects(f)) == 0

setObjectLocation!(f,a1,Real2D(1.0,1.0))
setObjectLocation!(f,a2,Real2D(1.0,1.399999))

@test length(getObjectsAtLocation(f,Real2D(1.0,1.0))) == 1

setObjectLocation!(f,a1, getObjectLocation(f,a2))

@test getObjectLocation(f,a1) == Real2D(1.0,1.399999)

@test length(getObjectsAtLocation(f, getObjectLocation(f,a1))) == 2

@test numObjectsAtLocation(f, getObjectLocation(f,a1)) == 2

f2 = Field2D(10.0,10.0,0.5,false)

pos = Real2D(5.2,5.2)
__distance = 3.0

setObjectLocation!(f2,Patch(nothing),Real2D(5.5,5.5))
setObjectLocation!(f2,Patch(nothing),Real2D(4.0,4.0))
setObjectLocation!(f2,Patch(nothing),Real2D(5.2,5.2))
setObjectLocation!(f2,Patch(nothing),Real2D(5.2,2.2))
setObjectLocation!(f2,Patch(nothing),Real2D(5.2,2.1))

@test length(getNeighborsWithinDistance(f2, pos , __distance)) == 4

setObjectLocation!(f2,Patch(nothing),Real2D(0.1,0.1))
__distance = 5.0

@test length(getNeighborsWithinDistance(f2, pos , __distance)) == 5

f3 = Field2D(10.0,10.0,0.5,true)

pos = Real2D(8.1,8.1)
__distance = 3.0

setObjectLocation!(f3,Patch(nothing),Real2D(7.0,8.1))
setObjectLocation!(f3,Patch(nothing),Real2D(8.1,1.0))
setObjectLocation!(f3,Patch(nothing),Real2D(8.1,1.3))
setObjectLocation!(f3,Patch(nothing),Real2D(8.1,9.9))

@test length(getNeighborsWithinDistance(f3, pos , __distance)) == 3

f4 = Field2D(10.0,10.0,0.6,true)

setObjectLocation!(f4,Patch(nothing),Real2D(7.0,8.1))
setObjectLocation!(f4,Patch(nothing),Real2D(8.1,1.0))
setObjectLocation!(f4,Patch(nothing),Real2D(8.1,1.3))
setObjectLocation!(f4,Patch(nothing),Real2D(8.1,9.9))

@test length(getNeighborsWithinDistance(f4, pos , __distance)) == 3

setObjectLocation!(f4,Patch(nothing),Real2D(1.0,9.7))

@test length(getNeighborsWithinDistance(f4, pos , __distance)) == 3

setObjectLocation!(f4,Patch(nothing),Real2D(1.1,8.1))

@test length(getNeighborsWithinDistance(f4, pos , __distance)) == 4

@test length(getNeighborsWithinDistance(f4,  Real2D(0.1,5.0) , __distance)) == 0

setObjectLocation!(f4,Patch(nothing),Real2D(9.9,2.1))
setObjectLocation!(f4,Patch(nothing),Real2D(2.9,4.5))

@test length(getNeighborsWithinDistance(f4,  Real2D(0.1,5.0) , __distance)) == 2

setObjectLocation!(f4,Patch(nothing),Real2D(2.9,8.0))

@test length(getNeighborsWithinDistance(f4,  Real2D(0.1,5.0) , __distance)) == 2


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule PCRE_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("PCRE")
JLLWrappers.@generate_main_file("PCRE", UUID("2f80f16e-611a-54ab-bc61-aa92de5b98fc"))
end  # module PCRE_jll


text:
# find hyperparameter values using gird search

include(pwd()*"/dlmodel.jl")

p = [0]

lambda_1 = collect(0:0.02:0.5)
lambda_2 = collect(0:0.02:0.5)
lambda_3 = [0]

model_fit_regularization = zeros(length(p),length(lambda_1),length(lambda_2),length(lambda_3),3+3+1+1)
model_fit_train = zeros(length(p),length(lambda_1),length(lambda_2),length(lambda_3),2000)
model_fit_val = zeros(length(p),length(lambda_1),length(lambda_2),length(lambda_3),2000)

# find hyperparameters using gridsearch
@time for i in 1:length(p)
    for j in 1:length(lambda_1)
        @printf "Percentage done %.2f\n" (j-1)/length(lambda_1)
        for k in 1:length(lambda_2)
            for l in 1:length(lambda_3)

                Random.seed!(1234) # set random numbers

                w = Any[ xavier(n_hidden_1,n_input), zeros(n_hidden_1,1),
                         #xavier(n_hidden_2,n_hidden_1), zeros(n_hidden_2,1),
                         xavier(n_out,n_hidden_1), zeros(n_out,1)]

                optim = optimizers(w, Adam)

                dropout_percentage = p[i]
                lambda = [lambda_1[j], lambda_2[k]]
                loss_train_vec, loss_val_vec = train(2000, dropout_percentage, lambda, w, optim)

                model_fit_train[i,j,k,l,:] = loss_train_vec
                model_fit_val[i,j,k,l,:] = loss_val_vec

                # calc predictions
                y_pred_train_proc = predict(w, x_train, 0)
                y_pred_val_proc = predict(w, x_val, 0)

                y_pred_train = class_nll(y_pred_train_proc,[1,2])
                y_pred_val = class_nll(y_pred_val_proc,[1,2])

                # loss
                model_fit_regularization[i,j,k,l,end-1] = loss(w,x_train, y_train, 0, lambda)/length(y_train)
                model_fit_regularization[i,j,k,l,end] = loss(w,x_val, y_val, 0, lambda)/length(y_val)

                # classification results
                model_fit_regularization[i,j,k,l,1:3] = classification_results(y_pred_train, y_train, [1,2])
                model_fit_regularization[i,j,k,l,4:6] = classification_results(y_pred_val, y_val, [1,2])

            end
        end
    end
end

# results without regularization
training_loss_no_reg = model_fit_regularization[1,1,1,1,end-1]
test_loss_no_reg = model_fit_regularization[1,1,1,1,end]
training_acc_no_reg = model_fit_regularization[1,1,1,1,1]
test_acc_no_reg = model_fit_regularization[1,1,1,1,4]

# find best hyperparameter settings
(~, idx_best_training_loss) = findmin(model_fit_regularization[:,:,:,:,end-1])
(~, idx_best_val_loss) = findmin(model_fit_regularization[:,:,:,:,end])

(~, idx_best_training_tp) = findmax(model_fit_regularization[:,:,:,:,2])
(~, idx_best_val_tp) = findmax(model_fit_regularization[:,:,:,:,5])

(~, idx_best_training_tn) = findmax(model_fit_regularization[:,:,:,:,3])
(~, idx_best_val_tn) = findmax(model_fit_regularization[:,:,:,:,6])

(~, idx_best_training_acc) = findmax(model_fit_regularization[:,:,:,:,1])
(~, idx_best_val_acc) = findmax(model_fit_regularization[:,:,:,:,4])

# 5,1,1,5
# 3,1,8,1

@printf "Best dropout percentage %.2f\n" p[1]
@printf "Best lambda_1 %.10f\n" lambda_1[1]
<<<<<<< HEAD
@printf "Best lambda_2 %.10f\n" lambda_2[4]
@printf "Best lambda_3 %.10f\n" lambda_3[1]
=======
@printf "Best lambda_2 %.10f\n" lambda_2[1]
@printf "Best lambda_3 %.10f\n" lambda_3[11]
>>>>>>> 2layernetwork

max_stars_repo_name: SamuelWiqvist/dl-compute-project

max_issues_repo_name: SamuelWiqvist/dl-compute-project

max_forks_repo_name: SamuelWiqvist/dl-compute-project


text:
type Parameter
    algorithm
    budget

    # common parameters that all algorithm should accept
    init_sample
    time_limit

    # for racos optimization
    precision
    uncertain_bits
    train_size
    positive_size
    negative_size
    probability
    replace_strategy

    # for tcp with python
    evaluation_server_num
    ip_port
    control_server_ip
    control_server_port
    objective_file
    func
    constraint
    show_x
    output_file

    # for pareto optimization
    isolationfunc

    # result
    positive_data
    negative_data

    function Parameter(; algorithm="asracos", budget=0, init_sample=Nullable(),
        time_limit=Nullable{Int64}(), precision=Nullable(),
        uncertain_bits=Nullable{Int64}(), train_size=0, positive_size=0, negative_size=0,
        probability=0.99, replace_strategy="WR", evaluation_server_num = 1,
        control_server_ip_port=Nullable{String}(), objective_file=Nullable{String}(),
        func="target_func", constraint=Nullable{String}(), show_x=false, output_file=Nullable{String}(),
        isolationfunc=x->0, autoset=true)

        temp = split(control_server_ip_port, ":")
        control_server_ip = temp[1]
        control_server_port = parse(Int64, temp[2])
        parameter = new(algorithm, budget, init_sample, time_limit,
        precision, uncertain_bits, train_size, positive_size, negative_size,
        probability, replace_strategy, evaluation_server_num, Nullable(), control_server_ip,
        control_server_port, objective_file, func, constraint, show_x, output_file,
        isolationfunc, [], [])
        if budget != 0 && autoset == true
            autoset!(parameter)
        end
        return parameter
    end
end

function autoset!(parameter)
    if parameter.budget < 3
        zoolog("parameter.jl: budget too small")
    elseif parameter.budget <= 50
        parameter.train_size = 4
        parameter.positive_size = 1
    elseif parameter.budget <= 100
        parameter.train_size = 6
        parameter.positive_size = 1
    elseif parameter.budget <= 1000
        parameter.train_size = 12
        parameter.positive_size = 2
    else
        parameter.train_size = 22
        parameter.positive_size = 2
    end
    parameter.negative_size = parameter.train_size - parameter.positive_size
end

function get_positive_data(parameter)
    result = take!(parameter.positive_data)
    return result
end

function get_negative_data(parameter)
    result = take!(parameter.negative_data)
    return result
end

function clear_parameter(parameter)
    parameter.control
end


text:
function lab_16(r::Robot)
    path = go_to_west_south_corner_and_return_path!(r; go_around_barriers = true)

    side = Ost
    while !isborder(r,Nord)
        go_to_border_and_return_path!(r, side; markers = true, go_around_barriers = true)
        side = inverse_side(side)
        go!(r, Nord)
    end

    go_to_border_and_return_path!(r, side; markers = true, go_around_barriers = true)

    go_to_west_south_corner_and_return_path!(r, go_around_barriers = true)
    go_by_path!(r, path)
end


text:
#==========================================================
This routine initializes an optimization problem to recover
a CNOT gate on a single qudit with 4 energy levels (and 2
guard states). The drift Hamiltonian in the rotating frame
is
              H_0 = - 0.5*ξ_a(a^†a^†aa),
where a is the annihilation operator for the qudit.
Here the control Hamiltonian includes the usual symmetric
and anti-symmetric terms
     H_{sym} = p(t)(a + a^†),    H_{asym} = q(t)(a - a^†)
which come from the rotating frame approximation and hence
we refer to these as "coupled" controls.
The problem parameters are:
                ω_a =  2π × 4.10595   Grad/s,
                ξ_a =  2π × 2(0.1099) Grad/s.
We use Bsplines with carrier waves with frequencies
0, ξ_a, 2ξ_a Grad/s.
==========================================================#
using LinearAlgebra
using Plots
pyplot()
using FFTW
using DelimitedFiles
using Printf
using Ipopt
using Random

Base.show(io::IO, f::Float64) = @printf(io, "%20.13e", f)

using Juqbox # quantum control module

Nosc = 1 # Number of coupled sub-systems = oscillators
N = 4 # Number of essential energy levels

Nguard = 2 # Number of extra guard levels
Ntot = N + Nguard # Total number of energy levels

T = 100.0 # Duration of gate.

# frequencies (in GHz, will be multiplied by 2*pi to get angular frequencies in the Hamiltonian matrix)
fa = 4.10336
xa = 0.2198
rot_freq = [fa] # Used to calculate the lab frame ctrl function

# setup drift Hamiltonian
number = Diagonal(collect(0:Ntot-1))

H0 = -0.5*(2*pi)*xa* (number*number - number) # xa is in GHz

# lowering matrix
amat = Bidiagonal(zeros(Ntot),sqrt.(collect(1:Ntot-1)),:U) # standard lowering matrix
adag = Array(transpose(amat));
Hsym_ops=[Array(amat + adag)]
Hanti_ops=[Array(amat - adag)]
H0 = Array(H0)

# Estimate time step
maxctrl = 0.001*2*pi * 8.5 #  9, 10.5, 12, 15 MHz

nsteps = calculate_timestep(T, H0, Hsym_ops, Hanti_ops, [maxctrl])
println("# time steps: ", nsteps)

use_bcarrier = true # Use carrier waves in the control pulses?

if use_bcarrier
    Nfreq = 3 # number of carrier frequencies
else
    Nfreq = 1
end

Ncoupled = length(Hsym_ops) # Here, Ncoupled = 1
om = zeros(Ncoupled,Nfreq)

if use_bcarrier
    om[1:Ncoupled,2] .= -2.0*pi *xa       # Note negative sign
    om[1:Ncoupled,3] .= -2.0*pi* 2.0*xa
end
println("Carrier frequencies [GHz]: ", om[1,:]./(2*pi))

maxamp = zeros(Nfreq)

if Nfreq >= 3
    const_fact = 0.45
    maxamp[1] = maxctrl*const_fact
    maxamp[2:Nfreq] .= maxctrl*(1.0-const_fact)/(Nfreq-1) # split the remainder equally
else
    # same threshold for all frequencies
    maxamp .= maxctrl/Nfreq
end

maxpar = maximum(maxamp)

Ident = Matrix{Float64}(I, Ntot, Ntot)

# CNOT trarget
utarget = Ident[1:Ntot,1:N]
utarget[:,3] = Ident[:,4]
utarget[:,4] = Ident[:,3]

omega1 = Juqbox.setup_rotmatrices([N], [Nguard], [fa])

# Compute Ra*utarget
rot1 = Diagonal(exp.(im*omega1*T))

# target in the rotating frame
vtarget = rot1*utarget

# Initial conditions
U0 = Ident[1:Ntot,1:N]

params = Juqbox.objparams([N], [Nguard],
                          T, nsteps, Uinit=U0, Utarget=vtarget,
                          Cfreq=om, Rfreq=rot_freq,
                          Hconst=H0, Hsym_ops=Hsym_ops, Hanti_ops=Hanti_ops)

# initial parameter guess

# D1 smaller than 5 does not work
D1 = 10 # Number of B-spline coefficients per segment
nCoeff = 2*Ncoupled*Nfreq*D1 # factor '2' is for sin/cos

# Random.seed!(12456)

startFromScratch = true
startFile="cnot-pcof-opt.dat"

# initial parameter guess
if startFromScratch
    pcof0 = maxpar*0.01 * rand(nCoeff)
    println("*** Starting from random pcof with amplitude ", maxpar*0.01)
else
    # use if you want to read the initial coefficients from file
    pcof0 = vec(readdlm(startFile))
    println("*** Starting from B-spline coefficients in file: ", startFile)
    nCoeff = length(pcof0)
    D1 = div(nCoeff, 2*Ncoupled*Nfreq)  # number of B-spline coeff per control function
end

# min and max coefficient values
useBarrier = true
minCoeff, maxCoeff = Juqbox.assign_thresholds_freq(maxamp, Ncoupled, Nfreq, D1)

samplerate = 32 # only used for plotting
casename = "cnot1" # base file name (used in optimize-once.jl)

maxIter = 200 # 0  # optional argument
lbfgsMax = 250 # optional argument
ipTol = 1e-5   # optional argument
acceptTol = ipTol # 1e-4 # acceptable tolerance
acceptIter = 15

println("*** Settings ***")
println("System Hamiltonian coefficients [GHz]: (fa, xa) =  ", fa, xa)
println("Total number of states, Ntot = ", Ntot, " Total number of guard states, Nguard = ", Nguard)
if use_bcarrier
  println("Using B-spline basis functions with carrier wave, # freq = ", Nfreq)
else
  println("Using regular B-spline basis functions")
end
println("Number of coefficients per spline = ", D1, " Total number of parameters = ", nCoeff)
for q=1:Nfreq
    println("Carrier frequency: ", om[q]/(2*pi), " GHz, max parameter amplitude = ", 1000*maxamp[q]/(2*pi), " MHz")
end
println("Tikhonov coefficients: tik0 = ", params.tik0)

# Estimate number of terms in Neumann series for time stepping (Default 3)
tol = eps(1.0); # machine precision
Juqbox.estimate_Neumann!(tol, params, [maxpar])

wa = Juqbox.Working_Arrays(params,nCoeff)

# prob = Juqbox.setup_ipopt_problem(params, wa, nCoeff, minCoeff, maxCoeff, maxIter, lbfgsMax, startFromScratch)
# println("Initial coefficient vector stored in 'pcof0'")

pcof0 = maxpar*0.5*(rand(nCoeff).-0.5)
pbjfv, U, mfidelityrot = Juqbox.traceobjgrad(pcof0,params,wa,true,false)

W = U[:,:,end]
W = W[:]
ndata = 10
nin = length(pcof0)
nout = 2*length(W)

function generate_data(nin,nout,ndata,params,wa)

    input_data = zeros(nin,ndata)
    output_data = zeros(nout,ndata)
    nouthalf = Int64(nout/2)

    for i = 1:ndata
        pcof0 = maxpar*0.5*(rand(nin).-0.5)
        objfv, U, mfidelityrot = Juqbox.traceobjgrad(pcof0,params,wa,true,false)
        W = U[:,:,end]
        W = W[:]
        input_data[:,i] = pcof0
        output_data[1:nouthalf,i] = real.(W)
        output_data[nouthalf+1:nout,i] = imag.(W)
    end
    return input_data, output_data
end

in_data, out_data = generate_data(nin,nout,ndata,params,wa)

using Flux
model = Chain(Dense(nin, 5, relu),
              Dense(5, nout))

loss(x, y) = Flux.Losses.mse(model(x), y)

opt = Descent()
parameters = Flux.params(model)
data = [(in_data, out_data)]

loss(in_data,out_data)
for epoch in 1:200
    Flux.train!(loss, parameters, data, opt)
end
loss(in_data,out_data)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule MathGL_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("MathGL")
JLLWrappers.@generate_main_file("MathGL", UUID("6834ddeb-87f2-5bbb-bfa4-c37572f854d4"))
end  # module MathGL_jll


text:
"""
`scatterplot(x, y; nargs...)` → `Plot`

Description
============

Draws the given points on a new canvas.
It uses the first parameter `x`
(which should be a vector or range) to denote
the horizontal position of each point,
and the second parameter `y`
(which should also be a vector or range)
as their vertical position.
This means that the two vectors have to have the same length.

Usage
======

    scatterplot(x, y; title = "", name = "", width = 40, height = 15, border = :solid, xlim = [0, 0], ylim = [0, 0], margin = 3, padding = 1, color = :auto, labels = true, canvas = BrailleCanvas, grid = true)

Arguments
==========

- **`x`** : The horizontal dimension for each point.

- **`y`** : The vertical dimension for each point.

- **`title`** : Text to display on the top of the plot.

- **`name`** : Annotation of the current drawing to displayed on the right

- **`width`** : Number of characters per row that should be used for plotting.

- **`height`** : Number of rows that should be used for plotting. Not applicable to `barplot`.

- **`border`** : The style of the bounding box of the plot. Supports `:solid`, `:bold`, `:dashed`, `:dotted`, `:ascii`, and `:none`.

- **`xlim`** : Plotting range for the x coordinate. `[0, 0]` stands for automatic.

- **`ylim`** : Plotting range for the y coordinate. `[0, 0]` stands for automatic.

- **`margin`** : Number of empty characters to the left of the whole plot.

- **`padding`** : Space of the left and right of the plot between the labels and the canvas.

- **`color`** : Color of the drawing. Can be any of `:blue`, `:red`, `:yellow`

- **`labels`** : Can be used to hide the labels by setting `labels=false`.

- **`canvas`** : The type of canvas that should be used for drawing.

- **`grid`** : Can be used to hide the gridlines at the origin

Returns
========

A plot object of type `Plot{T<:Canvas}`

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
=========

    julia> scatterplot(randn(50), randn(50), title = "My Scatterplot", color = :red)

                        My Scatterplot
           ┌────────────────────────────────────────┐
       2.1 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⠀⡧⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⡀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⠀⠀⠀⢀⠐⠂⠀⠀⠀⠀⠀⡇⠁⠐⠀⠀⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⠀│
           │⠠⠀⠀⠀⠀⠀⠐⠀⠀⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⠀⠄⡇⢀⠀⠀⠀⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠒⠖⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠚⠒⠒⡗⠒⠒⠒⢒⠒⠒⠒⠒⠒⠒⠒⠖⠒⠒⠒⠒⠒⠢│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠀⠀⠀⠀⠠⠀⠁⡇⠀⠀⠀⣀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠠⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⡇⠀⠀⠀⠀⢀⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠄⠀⠀⠀⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      -2.2 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
           └────────────────────────────────────────┘
           -2.1                                   1.9

See also
=========

`Plot`, `lineplot`, `stairs`, `BrailleCanvas`, `BlockCanvas`, `AsciiCanvas`
"""
function scatterplot(
        x::AbstractVector{<:Number}, y::AbstractVector{<:Number};
        color::Symbol = :auto,
        name::AbstractString = "",
        canvas::Type = BrailleCanvas,
        kw...)
    X = convert(Vector{Float64}, collect(x))
    Y = convert(Vector{Float64}, collect(y))
    new_plot = Plot(X, Y, canvas; kw...)
    color = (color == :auto) ? next_color!(new_plot) : color
    name == "" || annotate!(new_plot, :r, name, color)
    points!(new_plot, X, Y, color)
end

function scatterplot!(
        plot::Plot{<:Canvas}, x::AbstractVector{<:Number}, y::AbstractVector{<:Number};
        color::Symbol = :auto,
        name::AbstractString = "",
        kw...)
    X = convert(Vector{Float64}, collect(x))
    Y = convert(Vector{Float64}, collect(y))
    color = (color == :auto) ? next_color!(plot) : color
    name == "" || annotate!(plot, :r, name, color)
    points!(plot, X, Y, color)
end

function scatterplot(X::AbstractVector; kw...)
    scatterplot(1:length(X), X; kw...)
end

function scatterplot!(plot::Plot{<:Canvas}, X::AbstractVector; kw...)
    scatterplot!(plot, 1:length(X), X; kw...)
end


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

name = "Ripser"
version = v"0.3.2"

# Collection of sources required to build Ripser
sources = [
    "https://github.com/scikit-tda/ripser.py.git" =>
    "b749d0b170f708b9c9b125bffe1c302c870961a2"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
$CXX -std=c++11 -Ofast -fPIC -shared -L. $LDFLAGS \
     -D USE_COEFFICIENTS -D ASSEMBLE_REDUCTION_MATRIX -D LIBRIPSER \
     ripser.py/ripser/ripser.cpp -o libripser.$dlext
if [[ $dlext == dll ]]; then
    LIB="bin"
else
    LIB="lib"
fi
mkdir $prefix/$LIB
mv libripser.$dlext $prefix/$LIB
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(),
    FreeBSD(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libripser", :libripser)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

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

max_stars_repo_name: mtsch/RipserBuilder

max_issues_repo_name: mtsch/RipserBuilder

max_forks_repo_name: mtsch/RipserBuilder


text:
module TestConstraintRSOC

using Test

using VecMathOptInterface
const MOI = VecMathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_RSOC()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.RSOC{Float64}(mock)
    MOI.Test.runtests(
        bridged_mock,
        config,
        include = [
            "test_basic_VectorOfVariables_RotatedSecondOrderCone",
            "test_basic_VectorAffineFunction_RotatedSecondOrderCone",
            "test_basic_VectorQuadraticFunction_RotatedSecondOrderCone",
        ],
    )
    MOI.empty!(bridged_mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [0.5, 1.0, 1 / √2, 1 / √2],
            (MOI.VariableIndex, MOI.EqualTo{Float64}) => [-√2, -1 / √2],
            (MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone) =>
                [[3 / 2, 1 / 2, -1.0, -1.0]],
        )
    MOI.Test.test_conic_RotatedSecondOrderCone_VectorOfVariables(
        bridged_mock,
        config,
    )
    MOI.empty!(bridged_mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1 / √2, 1 / √2],
            (MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone) =>
                [[3 / 2, 1 / 2, -1.0, -1.0]],
        )
    MOI.Test.test_conic_RotatedSecondOrderCone_VectorAffineFunction(
        bridged_mock,
        config,
    )
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.RotatedSecondOrderCone,
            }(),
        ),
    )
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        value = [√2, 1 / √2, -1.0, -1.0]
        MOI.set(bridged_mock, attr, ci, value)
        @test MOI.get(bridged_mock, attr, ci) ≈ value
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        2,
        ((MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone, 0),),
    )
    return
end

function test_SOCR()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.SOCR{Float64}(mock)
    MOI.Test.runtests(
        bridged_mock,
        config,
        include = [
            "test_basic_VectorOfVariables_SecondOrderCone",
            "test_basic_VectorAffineFunction_SecondOrderCone",
            "test_basic_VectorQuadraticFunction_SecondOrderCone",
        ],
    )
    MOI.empty!(bridged_mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1.0, 1 / √2, 1 / √2],
            (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone) => [[1 - 1 / √2, 1 + 1 / √2, -1]],
            (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[-√2]],
        )
    MOI.Test.test_conic_SecondOrderCone_VectorOfVariables(bridged_mock, config)
    MOI.empty!(bridged_mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1.0, 1 / √2, 1 / √2],
            (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone) => [[1 - 1 / √2, 1 + 1 / √2, -1]],
            (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[-√2]],
        )
    MOI.Test.test_conic_SecondOrderCone_VectorAffineFunction(
        bridged_mock,
        config,
    )
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.SecondOrderCone,
            }(),
        ),
    )
    @test !MOI.Utilities.is_canonical(
        MOI.get(bridged_mock, MOI.ConstraintFunction(), ci),
    )
    @test MOI.Utilities.is_canonical(
        MOI.get(bridged_mock, MOI.CanonicalConstraintFunction(), ci),
    )
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        value = [√2, 1 / √2, -1.0, -1.0]
        MOI.set(bridged_mock, attr, ci, value)
        @test MOI.get(bridged_mock, attr, ci) ≈ value
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        3,
        ((MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),),
    )
    return
end

end  # module

TestConstraintRSOC.runtests()


text:
function link_table(file::String)

    id = Symbol(ncgetatt(file, "dim_space", "id"))

    df_links = CSV.File(joinpath(dirname(pathof(IcanProj)), "..", "data", "df_links.csv")) |> DataFrame

    eval(:(df_nc = DataFrame($(id) = convert(Array{Int64}, ncread($(file), "id")),
                             nc = 1:length(ncread($(file), "id")))))

    df_final = join(df_links, df_nc, on = id)

    return df_final

end

function link_table(path::String, variable::String, space_res::String)

    file = get_filename(path, variable, space_res, 1)

    id = Symbol(ncgetatt(file, "id", "id"))

    df_links = CSV.File(joinpath(dirname(pathof(IcanProj)), "..", "data", "df_links.csv")) |> DataFrame

    eval(:(df_nc = DataFrame($(id) = convert(Array{Int64}, ncread($(file), "id")),
                             nc = 1:length(ncread($(file), "id")))))

    df_final = join(df_links, df_nc, on = id)

    return df_final

end

function link_results(file_fine::String, file_coarse::String) #, id_fine, id_coarse)

    # Metadata table

    df_meta = CSV.File(joinpath(dirname(pathof(IcanProj)), "..", "data", "df_links.csv")) |> DataFrame

    # Get attributes

    id_fine = ncgetatt(file_fine, "id", "id")

    if typeof(id_fine) == Nothing
        id_fine = ncgetatt(file_fine, "dim_space", "id")
    end

    id_coarse = ncgetatt(file_coarse, "id", "id")

    if typeof(id_coarse) == Nothing
        id_coarse = ncgetatt(file_coarse, "dim_space", "id")
    end

    # Fine resolution

    df_fine = DataFrame(id_fine = convert(Array{Int64}, ncread(file_fine, "id")),
                        nc_fine = 1:length(ncread(file_fine, "id")))

    # Coarse resolution

    df_coarse = DataFrame(id_coarse = convert(Array{Int64}, ncread(file_coarse, "id")),
                          nc_coarse = 1:length(ncread(file_coarse, "id")))

    # Link resolutions

    df_links = DataFrame(id_fine = df_meta[Symbol(id_fine)],
                         id_coarse = df_meta[Symbol(id_coarse)],
                         ind_julia = df_meta[:ind_julia],
                         elev = df_meta[:elev])

    df_links = join(df_links, df_fine, on = :id_fine)

    df_links = join(df_links, df_coarse, on = :id_coarse)

    # Compute number of fine scale grid cells in coarse grid

    tmp = DataFrame(nc_coarse = df_links[:nc_coarse],
                    n = 1)

    tmp = aggregate(tmp, :nc_coarse, sum)

    df_links = join(df_links, tmp, on = :nc_coarse)

    # Sort results

    sort!(df_links, [:id_fine])

    return df_links

end

"""
Aggregate fine scale results and return together with coarse scale results.
"""
function unify_results(file_fine, file_coarse, df_links, variable)

    var_fine = ncread(file_fine, variable)

    dim_space = ncread(file_fine, "id")  # convert to int?

    if size(var_fine, 2) != length(dim_space)
        var_fine = var_fine'
    end

    var_coarse = ncread(file_coarse, variable)

    dim_space = ncread(file_coarse, "id")  # convert to int?

    if size(var_coarse, 2) != length(dim_space)
        var_coarse = var_coarse'
    end

    var_agg = fill(0.0, size(var_coarse))

    ngrids = fill(0.0, size(var_coarse, 2))

    for row in eachrow(df_links)

        # This works even though some of the gridcells lack values

        tmp = var_fine[:, row[:nc_fine]]

        if !any(isnan.(tmp))

            var_agg[:, row[:nc_coarse]] += var_fine[:, row[:nc_fine]]

            ngrids[row[:nc_coarse]] += 1

        end

        # This should work if all grid cells have results

        # var_agg[:, row[:nc_coarse]] += var_fine[:, row[:nc_fine]] / row[:n_sum]

        # ngrids[row[:nc_coarse]] += 1

    end

    var_agg = var_agg ./ repeat(permutedims(ngrids), size(var_agg, 1))

    ncclose()

    return var_coarse, var_agg, ngrids

end

"""
Project results to a map using senorge extents.
"""
function project_results(values, df_links, on = :nc_coarse)

    map_values = fill(NaN, 1550, 1195)

    map_values[df_links[:ind_julia]] = values[df_links[on]]

    return map_values

end


text:
#=
surface_test:
- Julia version: 1.1.0
- Author: Jirik
- Date: 2019-03-19
=#
# # using Revise
using Test
using Logging
using SparseArrays
using Distributed
if nprocs() == 1
    addprocs(3)
end

using LarSurf
# using Plasm
# using LinearAlgebraicRepresentation
# Lar = LinearAlgebraicRepresentation

# TODO repair this test
@testset "Get Fchar test" begin

    println("generate data")
    segmentation = LarSurf.data234()
    block_size = [2,2,2]
    # Flin2 = __grid_get_surface_Fchar_per_block_parallel_pmap(segmentation, block_size)

end


text:
struct Liouvillian{T<:Real,M<:AbstractMatrix,N,F,D} <: QuObject
    L₀::M
    Lₙ::NTuple{N,M}
    fₙ::F
    pₙ::NTuple{N,Vector{T}}
    dims::Dims{D}
end

Liouvillian(dims::Dims{D},L₀::M,Lₙ=(),fₙ=(),pₙ=()) where {D,M<:AbstractMatrix} = Liouvillian{real(eltype(M)),M,length(Lₙ),typeof(fₙ),D}(L₀, Lₙ, fₙ, pₙ, dims)

data(L::Liouvillian) = (L.L₀, L.Lₙ, L.fₙ, L.pₙ)

# Liouvillian ODE Type, compatible with AbstractODEFunction interface
using OrdinaryDiffEq.DiffEqBase: AbstractODEFunction
struct LiouvillianODE{LT<:Liouvillian,TD,TJ,TP} <: AbstractODEFunction{true}
    L::LT
    tgrad::TD
    jac::TJ
    jac_prototype::TP
    sparsity::TP
    mass_matrix::UniformScaling{Bool}
end

function LiouvillianODE(L::TL) where {TL<:Liouvillian}
    tgrad = (dT,ψ,p,t) -> tgrad_liouvillian!(L, dT, ψ, t)
    jac = (J,ψ,p,t) -> jac_liouvillian!(L, J, t)
    jac_prot = similar(L.L₀); jac(jac_prot,nothing,nothing,zero(real(eltype(jac_prot))))
    return LiouvillianODE{TL,typeof(tgrad),typeof(jac),typeof(jac_prot)}(L, tgrad, jac, jac_prot, jac_prot, I)
end

(f::LiouvillianODE)(dψ,ψ,p,t) = mul_liouvillian!(f.L, dψ, ψ, t)

# ODE
mul_liouvillian!(L::Liouvillian,dψ,ψ,t) = mul_liouvillian!(dψ, ψ, t, L.L₀, L.Lₙ, L.fₙ, L.pₙ)

function mul_liouvillian!(dψ::AbstractArray,ψ::AbstractArray,t::Real,L₀::AbstractMatrix,Lₙ::Tuple,fₙ::Tuple,pₙ::Tuple)
    mul!(dψ, L₀, ψ, 1, 0)
    applyfun!(dψ, ψ, t, Lₙ, fₙ, pₙ)
end
@inline function applyfun!(dψ::AbstractArray,ψ::AbstractArray,t::Real,Lₙ::Tuple,fₙ::Tuple,pₙ::Tuple)
    mul!(dψ, first(Lₙ), ψ, first(fₙ)(t,first(pₙ)), 1)
    applyfun!(dψ, ψ, t, tail(Lₙ), tail(fₙ), tail(pₙ))
end
@inline applyfun!(dψ::AbstractArray,ψ::AbstractArray,t::Real,Lₙ::Tuple{},fₙ::Tuple{},pₙ::Tuple{}) = nothing

# Jacobian
jac_liouvillian!(L::Liouvillian,J::AbstractMatrix,t::Real) = jac_liouvillian!(J, t, L.L₀, L.Lₙ, L.fₙ, L.pₙ)

function jac_liouvillian!(J::AbstractMatrix,t::Real,L₀::AbstractMatrix,Lₙ::Tuple,fₙ::Tuple,pₙ::Tuple)
    copy!(J, L₀)
    applyjac!(J, t, Lₙ, fₙ, pₙ)
end
@inline function applyjac!(J::AbstractMatrix,t::Real,Lₙ::Tuple,fₙ::Tuple,pₙ::Tuple)
    J += first(fₙ)(t,first(pₙ)) * first(Lₙ)
    applyjac!(J, t, tail(Lₙ), tail(fₙ), tail(pₙ))
end
@inline applyjac!(J::AbstractMatrix,t::Real,Lₙ::Tuple{},fₙ::Tuple{},pₙ::Tuple{}) = nothing

# Time-derivative
tgrad_liouvillian!(L::Liouvillian,dT,ψ,t) = tgrad_liouvillian!(dT, ψ, t, L.L₀, L.Lₙ, L.fₙ, L.pₙ)

function tgrad_liouvillian!(dT::AbstractArray,ψ::AbstractArray,t::Real,L₀::AbstractMatrix,Lₙ::Tuple,fₙ::Tuple,pₙ::Tuple)
    fill!(dT,0)
    apply_tgrad!(dT, ψ, t, Lₙ, fₙ, pₙ)
end
@inline function apply_tgrad!(dT::AbstractArray,ψ::AbstractArray,t::Real,Lₙ::Tuple,fₙ::Tuple,pₙ::Tuple)
    mul!(dT, first(Lₙ), ψ, deriv(t->first(fₙ)(t,first(pₙ)), t), 1)
    apply_tgrad!(dT, ψ, t, tail(Lₙ), tail(fₙ), tail(pₙ))
end
apply_tgrad!(dT::AbstractArray,ψ::AbstractArray,t::Real,Lₙ::Tuple{},fₙ::Tuple{},pₙ::Tuple{}) = nothing

# ForwardDiff doesn't support complex-valued functions so use finite difference in the meantime
function deriv(f::F, x::T) where {F,T<:Real}
    r = sqrt(eps(T))
    h = max(r*abs(x), r)
    return (f(x+h)-f(x-h))/2h
end


text:
ENV["GKSwstype"] = "nul"  # For StatsPlots

import Pkg; Pkg.activate("../../../")
include("../../info.jl")
include("scenarios.jl")

using CytofDiffDensity; const cdd = CytofDiffDensity
using Util
import Random
using DrWatson
using MCMC
using BSON
using Distributions  # required for `BSON.load`
using StatsPlots
using StatsFuns

# Plot settings
plotsize = (450, 450)
Plots.scalefontsizes()
Plots.scalefontsizes(1.5)

simname = "v1"
make_resultsdir(sim) = mkpath("$(Info.resultsdir_simstudy)/$(simname)/$(savename(sim))")

function run(sim)
  sim = (; sim...)

  # Simulate data.
  simdata = scenarios(sim.snum, Ni=10000, seed=1)

  # Make results dir if needed.
  resultsdir = make_resultsdir(sim)

  # Create Priors.
  y = [simdata.yC; simdata.yT]

  # Plot and save data
  imgdir = mkpath(joinpath(resultsdir, "img"))
  plot(size=plotsize)
  histogram!(simdata.yC, normalize=true, label=nothing, la=0, alpha=.3, color=:blue)
  histogram!(simdata.yT, normalize=true, label=nothing, la=0, alpha=.3, color=:red)
  ygrid = cdd.make_ygrid(y, 100)
  plot!(ygrid, pdf.(simdata.mmC, ygrid), label=nothing, color=:blue)
  plot!(ygrid, pdf.(simdata.mmT, ygrid), label=nothing, color=:red)
  savefig(joinpath(imgdir, "data.pdf"))
  closeall()

  # TODO: Does this help?
  mu_prior = cdd.make_ordered_prior(y, sim.K,
                                    s=(maximum(y) - minimum(y)) / (2 * sim.K))

  Random.seed!(0)
  model = Gtilde(simdata.yC, simdata.yT, sim.K, sim.beta, mu=mu_prior)
  exclude = [:vC, :vT, :zetaC, :zetaT, :lambdaC, :lambdaT]

  init = MCMC.make_init_state(model)
  spl = make_sampler(model, init=init)

  nburn, nsamps, thin = (10000, 4000, 2)
  callback = make_callback(model, nburn=nburn, nsamps=nsamps, thin=thin)

  Util.redirect_stdout_to_file(joinpath(resultsdir, "log.txt")) do
    # Run MCMC.
    chain, metrics = mcmc(spl, nsamps, init=init, nburn=nburn, thin=thin,
                          callback=callback, exclude=exclude);

    # Save chain and metrics.
    BSON.bson(joinpath(resultsdir, "results.bson"),
              Dict(:chain => chain, :metrics => metrics, :simdata => simdata,
                   :model => model))
  end
end

function _postprocess(sim, resultsdir)
  sim = (; sim...)

  # Make image dir if needed.
  imgdir = mkpath(joinpath(resultsdir, "img"))

  # Load results.
  # chain, metrics, simdata, model.
  r = (; BSON.load(joinpath(resultsdir, "results.bson"))...)

  # Print summary stats
  cdd.printsummary(r.chain, r.metrics); println()

  # Print model info
  cdd.print_model_info(r.model)

  # Plot loglike
  plot(r.metrics[:loglike], label=nothing)
  savefig(joinpath(imgdir, "loglike.pdf"))
  closeall()

  # Plot posterior for each multivariate parameter.
  for sym in [:mu, :sigma, :nu, :phi, :etaC, :etaT]
    p = hcat(getindex.(r.chain, sym)...)'

    plot(p, label=nothing)
    savefig(joinpath(imgdir, "$(sym)-trace.pdf"))
    closeall()

    boxplot(p, label=nothing)
    savefig(joinpath(imgdir, "$(sym).pdf"))
    closeall()
  end

  # Plot posterior for each univariate parameter.
  for sym in [:tau]
    p = vec(hcat(getindex.(r.chain, sym)...))

    plot(p, label=nothing)
    savefig(joinpath(imgdir, "$(sym)-trace.pdf"))
    closeall()

    histogram(p, label=nothing, normalize=true)
    savefig(joinpath(imgdir, "$(sym).pdf"))
    closeall()
  end

  # TODO: Plot posterior density for C and T.
  y = [r.simdata.yC; r.simdata.yT]
  ygrid = cdd.make_ygrid(y, 100)
  plot(size=plotsize)
  cdd.plot_post_density!(r.chain, ygrid)
  cdd.plot_simtruth(r.simdata.mmC, r.simdata.mmT, ygrid, label=nothing)
  savefig(joinpath(imgdir, "post-density.pdf"))
  closeall()
end

function postprocess(sim)
  # Make results dir if needed.
  resultsdir = make_resultsdir(sim)

  Util.redirect_stdout_to_file(joinpath(resultsdir, "info.txt")) do
    _postprocess(sim, resultsdir)
  end
end

function load_results(snum, K, beta)
  d = Dict(:snum => snum, :K => K, :beta => beta)
  resdir = make_resultsdir(d)
  return (; BSON.load(joinpath(resdir, "results.bson"))...)
end

function combine_results(snum::Integer, K::Integer; p::Real=0.1)
  (0 < p < 1) || error("0 < p < 1 required!")

  # Results dir
  resdir = make_resultsdir(Dict(:snum => snum, :K => K))
  imgdir = mkpath(joinpath(resdir, "img"))

  # Load results
  r0 = load_results(snum, K, 0)
  r1 = load_results(snum, K, 1)
  length(r0.chain) == length(r1.chain) || error("Chain lengths need to be equal!")

  # Loglikelihoods
  ll0 = r0.metrics[:loglike]
  ll1 = r1.metrics[:loglike]
  length(ll0) == length(ll1) || error("Loglikelihood lengths need to be equal!")

  # Number of MCMC samples
  B = length(ll0)

  # Log BF for Gtilde
  log_bf_gtilde = MCMC.log_bayes_factor(ll0, ll1)

  # Log bf for PZero
  pzero = Pzero(NC=r0.simdata.NC, NT=r0.simdata.NT,
                QC=r0.simdata.QC, QT=r0.simdata.QT)
  log_bf_pzero = cdd.compute_log_bf(pzero, B)

  # Compute P(β=1 | data)
  pm1 = logistic(log_bf_gtilde + log_bf_pzero + logit(p))

  # Combine chains and loglikelihoods.
  beta = pm1 .> rand(B)
  chain = [beta[b] ? r1.chain[b] : r0.chain[b] for b in 1:B]
  loglike = [beta[b] ? ll1[b] : ll0[b] for b in 1:B]
  m0_gamma = cdd.infer_Pzero0(pzero, B)
  m1_gamma = cdd.infer_Pzero1(pzero, B)
  gammaC = [beta[b] ? m1_gamma.gammaC_samples[b] : m0_gamma.gammaC_samples[b] for b in 1:B]
  gammaT = [beta[b] ? m1_gamma.gammaT_samples[b] : m0_gamma.gammaT_samples[b] for b in 1:B]

  # Plot gamma
  plot(size=plotsize)
  histogram!(gammaC, normalize=true, alpha=0.3, la=0, label=nothing, color=:blue)
  histogram!(gammaT, normalize=true, alpha=0.3, la=0, label=nothing, color=:red)
  savefig(joinpath(imgdir, "gamma.pdf"))
  closeall()

  # Plot Gtilde
  y = [r0.simdata.yC; r0.simdata.yT]
  ygrid = cdd.make_ygrid(y)
  plot(size=plotsize)
  cdd.plot_post_density!(chain, ygrid)
  cdd.plot_simtruth(r0.simdata.mmC, r0.simdata.mmT, ygrid, label=nothing)
  savefig(joinpath(imgdir, "post-density.pdf"))
  closeall()

  # DIC
  merged_dic = MCMC.dic(loglike)

  # Write results
  open(joinpath(resdir, "metrics.txt"), "w") do io
    write(io, "log BF(Gtilde) in favor of model 1: $(log_bf_gtilde) \n")
    write(io, "log BF(PZero) in favor of model 1: $(log_bf_pzero) \n")
    write(io, "p: $(p) \n")
    write(io, "pm1: $(pm1) \n")
    write(io, "DIC: $(merged_dic) \n")
  end

  return (chain=chain, loglike=loglike, dic=merged_dic)
end


text:
import CameraModel
import SLAMData
# import Surfels

import DenseMapPosterior

data_seq_name_w_ext = "livingroom2.klg"
data_seq_name = splitext(data_seq_name_w_ext)[1]

fused_model_fpath = "/tmp/fused_model.ply"

data_dir = joinpath(dirname(dirname(pathof(DenseMapPosterior))), "data", "livingroom2")

data_file = joinpath(data_dir, data_seq_name_w_ext)
cam = CameraModel.Aug_ICL_NUIM_CamParam()
rgbd_reader = SLAMData.KlgReader(data_file, cam=cam);

# traj_path = joinpath(data_dir, "m1.log")   # -14.06029388
traj_path = joinpath(data_dir, "m2.log")   # -37.356114
# traj_path = joinpath(data_dir, "m3.log")   # -48.314537
# traj_path = joinpath(data_dir, "m4.log")   # -982.4293978
# traj_path = joinpath(data_dir, "m5.log")   # -2252.5785928

est_name = splitext(basename(traj_path))[1]
traj = SLAMData.load_redwood_traj(traj_path)

# ------------------- fuse model -------------------------------------------
# create opengl context

# viewer_ = Surfels.viewer()

DenseMapPosterior.fuse_RGBD_to_3D_model(
    rgbd_reader,
    traj.poses, #id2Twc
    rgbd_reader.cam,
    fused_model_fpath
)

# -------------------- eval DMP ---------------------------------------------

dmp = DenseMapPosterior.DMP_RGBD(
    rgbd_reader,
    fused_model_fpath,
    traj.poses,
    cam
)

println("DMP of $est_name on $(data_seq_name) is: $dmp")


text:
@everywhere module PCBM  # defined on all workers
    type Model
        rxns
        mets
        genes
        S
        lb
        ub
        c
        b
        csense
        rxn_gene_mat
        rxn_name
        rxn_rules
        rxn_subsystem
        rxn_extra
        met_formula
        met_name
        met_extra
        gene_name
        gene_extra
        description
    end

    solvers = map(x -> in(x, readdir(Pkg.dir())) ? x : "" , ["GLPK", "CPLEX"])

    function set_solver(solver)
        Solver = ""
        if solver == "glpk"
            Solver = "GLPK"
        elseif solver == "cplex"
            Solver = "CPLEX"
        else
            warn("choose one of the following solvers:")
            warn("glpk")
            warn("cplex")
            return
        end

        global settings_default_lp = solver
    end

    if in("GLPK", solvers)
        using GLPK
        include(Pkg.dir() * "/CBM/src/core/solvers/GLPK.jl")
        global settings_default_lp = "glpk"
    end

    if in("CPLEX", solvers)
        using CPLEX
        include(Pkg.dir() * "/CBM/src/core/solvers/CPLEX.jl")
    end

    if !isdefined(:settings_default_lp)
        if in("GLPK", solvers)
            global settings_default_lp = "glpk"
        elseif in("CPLEX", solvers)
            global settings_default_lp = "cplex"
        else
            warn("Neither GLPK or CPLEX is installed")
        end
    end

    function global_setup_lp{T <: Number}(c::Array{T}, lb::Array{T}, ub::Array{T}, b_lower::Array{T}, b_upper::Array{T}, S::SparseMatrixCSC, objective::AbstractString = "max")
        global lp = PCBM.setup_lp(c,lb,ub,b_lower, b_upper, S, objective)
    end

    function setup_lp{T <: Number}(c::Array{T}, lb::Array{T}, ub::Array{T}, b_lower::Array{T}, b_upper::Array{T}, S::SparseMatrixCSC, objective::AbstractString = "max")
        if settings_default_lp == "glpk"
            return PCBM.setup_glpk(c,lb,ub,b_lower, b_upper, S, objective)
        elseif settings_default_lp == "cplex"
            return PCBM.setup_cplex(c,lb,ub,b_lower, b_upper, S, objective)
        else
            warn("No solver")
        end
    end

    ###############################################################################
    # stuff for FVA in Parallel
    ###############################################################################
    # THINK setup_workers is perhaps uneccesary check when you have time
    function setup_workers(model)
        @sync for pid in workers()
            @async remotecall_fetch(PCBM.global_setup_lp, pid,
                zeros(Float64, length(model.rxns)),
                    model.lb,
                    model.ub,
                    model.b,
                    model.b,
                    model.S)
        end
    end

    function fast_fva_objective(x)
        res = Array(Float64, length(x))

        for (i,v) in enumerate(x)
            PCBM.set_objective(lp, v, 1.0)

            PCBM.solve(lp)
            res[i] = PCBM.get_objective_value(lp)

            PCBM.set_objective(lp, v, 0.0)
        end

        return res
    end

    function fast_fva_solution(x)
        res = Array(Float64, length(get_variable_bounds(lp)[1]), length(x))

        for (i,v) in enumerate(x)
            PCBM.set_objective(lp, v, 1.0)

            PCBM.solve(lp)
            res[:,i] = PCBM.get_solution(lp)

            PCBM.set_objective(lp, v, 0.0)
        end

        return res
    end

    function fast_fva(model, optPercentage::Number = 1, flux_matrix::Bool = false)
        nrxns = length(model.rxns)
        x = collect(1:nrxns)

        lpmin = setup_lp(model.c, model.lb, model.ub, model.b,model.b, model.S)
        solve(lpmin)
        wt_min = get_objective_value(lpmin)

        lb = deepcopy(model.lb)
        lb[find(model.c)] *= wt_min

        @sync for pid in workers()
            @async remotecall_fetch(PCBM.global_setup_lp, pid,
                zeros(Float64, length(model.rxns)),
                    lb,
                    model.ub,
                    model.b,
                    model.b,
                    model.S)
        end

        stepsize = div(nrxns, 3*length(workers()))
        checklist = map(x -> collect(x:x+stepsize-1), 1:stepsize:nrxns)
        checklist[end] = checklist[end][checklist[end] .<= nrxns]

        maxF = flux_matrix ?  pmap(PCBM.fast_fva_solution, checklist) : pmap(PCBM.fast_fva_objective, checklist)

        map(x -> remotecall_fetch(() -> PCBM.set_direction(lp, "min"), x), workers())

        minF = flux_matrix ?  pmap(PCBM.fast_fva_solution, checklist) : pmap(PCBM.fast_fva_objective, checklist)

        if flux_matrix
            maxV = hcat(maxF...)
            minV = hcat(minF...)
            maxF = diag(maxV)
            minF = diag(minV)
            return minF, maxF, minV, maxV
        end

        return vcat(minF...), vcat(maxF...)
    end
    ###############################################################################
    ###############################################################################

    ###############################################################################
    # stuff for gene deletion in Parallel
    ###############################################################################
    function setup_workers_for_gene_deletion(model)
        @sync for pid in workers()
            @async remotecall_fetch(PCBM.global_setup_lp, pid,
                model.c,
                model.lb,
                model.ub,
                model.b,
                model.b,
                model.S)
        end
    end

    function pgd_analysis(rxns)
        #progress_meter = ProgressMeter.Progress(length(rxns), 0.5, "\n\u1b[1FFVA:", 50)
        combos = sort(rxns, lt=lexless)
        num_combos = length(combos)

        solve(lp)
        wt = get_objective_value(lp)

        lb, ub = get_variable_bounds(lp)

        flow = Array(Float64, num_combos)

        ####
        for d in combos[1]
            set_col_bounds(lp, d, 0.0, 0.0)
        end

        solve(lp)
        flow[1] = get_objective_value(lp)

        disable = []
        for i in 2:num_combos
            enable = setdiff(combos[i-1], combos[i])
            disable = setdiff(combos[i], combos[i-1])

            for e in enable
                set_col_bounds(lp, e, lb[e], ub[e])
            end

            for d in disable
                set_col_bounds(lp, d, 0.0, 0.0)
            end
            solve(lp)
            flow[i] = get_objective_value(lp)
            #next!(progress_meter)
        end

        for e in disable
            set_col_bounds(lp, e, lb[e], ub[e])
        end

        flow /= wt
        flow = round(flow, 16)

        result = Dict(zip(combos, flow))

        map(x -> set_col_bounds(lp, x, lb[x], ub[x]), 1:length(lb))

        return result
    end

    function pgd(model, rxns)
        listlen = length(rxns)
        setup_workers_for_gene_deletion(model)

        stepsize = div(listlen, 3*length(workers()))
        checklist = map(x -> collect(x:x+stepsize-1), 1:stepsize:listlen)
        checklist[end] = checklist[end][checklist[end] .<= listlen]

        rxnlist = map(x -> rxns[x], checklist)

        res = pmap(PCBM.pgd_analysis, rxnlist)
    end
    ###############################################################################
    ###############################################################################

end

max_stars_repo_name: isebarn/CBM

max_issues_repo_name: isebarn/CBM

max_forks_repo_name: isebarn/CBM


text:
struct WrappedLearner <: Learner
    learner
    stop
    askhook
    tellhook
end
WrappedLearner(learner,stop) = WrappedLearner(learner, stop, (x,y)->nothing, (x,y)->nothing)

function ask!(learner::WrappedLearner,input)
    if learner.stop(learner.learner)
        val = nothing
    else
        val = ask!(learner.learner,input)
    end
    learner.askhook(learner.learner,val)
    return val
end

function tell!(learner::WrappedLearner,message)
    tell!(learner.learner,message)
    learner.tellhook(learner.learner,message)
end

### The interface probabily needs to be a little different
mutable struct IgnorantLearner <: Learner
    iter
    state
end
IgnorantLearner(iter) = IgnorantLearner(iter,nothing)

function ask!(learner::IgnorantLearner,input)

    if input==nothing
        return nothing
    end

    temp = learner.state==nothing ? iterate(learner.iter) : iterate(learner.iter,learner.state)
    if temp==nothing
        return nothing
    else
        val,learner.state = temp
        return val
    end
end

tell!(learner::IgnorantLearner,m) = nothing

### Something which would help to reproduce the state of the Learner.
mutable struct HistoryMaster <: Master
    tasks
    results
    np
end

nparalel(master::HistoryMaster) = master.np

function HistoryMaster(output,np)
    N = length(output)
    tasks = Channel{Any}(N)
    results = Channel{Tuple{Any,Any}}(N)

    for oi in output
        put!(results,oi)
    end

    HistoryMaster(tasks,results,np)
end


text:
#functions for calling the matrix vector products
#i.e. function A(y_in,Nx,Ny,h) computes the product A*u and stores it in y

#Right now only_out[for  the second-order accurate SBP operators for constant coefficients from
#Mattsson and Nordstrom, 2004.

#Written for the 2D domain (x,y) \in (a, b) \times (c, d),
#stacking the grid function in the vertical direction,
#so that, e.g., \partial u / \partial x \approx D \kron I and
#\partial u / \partial y_out[\approx I \kron D
#the faces of the 2D domain are label bottom to top, left to right, i.e.
#side 1 is y_out[= c, side 2 is y_out[= d
#side 3 is x = a, side 4 is x = b

#D2x and D2y_out[compute approximations to
#\partial^2 u / \partial x^2 and
#\partial^2 u / \partial y^2, respectively

#FACEtoVOL(face,y_in_face,Nx,Ny) maps the face value y_in_face to a full length vector at the
#nodes corresponding to face

#BxSx \approx the traction on faces 3 and 4
#BySy_out[\approx the traction on faces 1 and 2

#BxSx_tran and BySy_tran are the transposes of BxSx and BySy

function D2x(y_in, Nx, Ny, h)
	N = Nx*Ny
	y_out = zeros(N)
	idx = 1:Ny
	y_out[idx] = (y_in[idx] - 2 .* y_in[Ny.+ idx] + y_in[2*Ny .+ idx]) ./ h^2

	idx1 = Ny+1:N-Ny
	y_out[idx1] = (y_in[idx1 .- Ny] - 2 .* y_in[idx1] + y_in[idx1 .+ Ny]) ./ h^2

	idx2 = N-Ny+1:N
	y_out[idx2] = (y_in[idx2 .- 2*Ny] -2 .* y_in[idx2 .- Ny] + y_in[idx2]) ./ h^2

	return y_out
end

function D2y(y_in, Nx, Ny, h)
	N = Nx*Ny
	y_out = zeros(N)
	idx = 1:Ny:N-Ny+1
	y_out[idx] = (y_in[idx] - 2 .* y_in[idx .+ 1] + y_in[idx .+ 2]) ./ h^2

	idx1 = Ny:Ny:N
	y_out[idx1] = (y_in[idx1 .- 2] - 2 .* y_in[idx1 .- 1] + y_in[idx1]) ./ h^2

	for j = 1:Nx
		idx = 2+(j-1)*Ny:j*Ny-1
		y_out[idx] = (y_in[idx .- 1] - 2 .* y_in[idx] + y_in[idx .+ 1]) ./ h^2
	end

	return y_out

end

function Dx(y_in, Nx, Ny, h)
	N = Nx*Ny
	y_out = zeros(N)

	idx = 1:Ny
	y_out[idx] = (y_in[idx .+ Ny] - y_in[idx]) ./ h

	idx1 = Ny+1:N-Ny
	y_out[idx1] = (y_in[idx1 .+ Ny]-y_in[idx1 .- Ny]) ./ (2*h)

	idx2 = N-Ny+1:N
	y_out[idx2] = (y_in[idx2]-y_in[idx2 .- Ny]) ./ h

	return y_out
end

function Dy(y_in, Nx, Ny, h)
	N = Nx*Ny
	y_out = zeros(N)

	idx = 1:Ny:N-Ny+1
	y_out[idx] = (y_in[idx .+ 1] - y_in[idx]) ./ h

	idx = Ny:Ny:N
	y_out[idx] = (y_in[idx] - y_in[idx .- 1]) ./h

	for j = 1:Nx
		idx = 2+(j-1)*Ny:j*Ny-1
		y_out[idx] = (y_in[idx .+ 1] - y_in[idx .- 1]) ./ (2*h)
	end

	return y_out
end

function Hxinv(y_in, Nx, Ny, h)
	N = Nx*Ny
	y_out = zeros(N)

	idx = 1:Ny
	y_out[idx] = (2*y_in[idx]) .* (1/h)

	idx = Ny+1:N-Ny
	y_out[idx] = (1*y_in[idx]) .* (1/h)

	idx = N-Ny+1:N
	y_out[idx] = (2*y_in[idx]) .* (1/h)

	return y_out
end

function Hyinv(y_in, Nx, Ny, h)
	N = Nx*Ny
	y_out = zeros(N)

	idx = 1:Ny:N-Ny+1
	y_out[idx] = (2*y_in[idx]) .* (1/h)

	idx = Ny:Ny:N
	y_out[idx] = (2*y_in[idx]) .* (1/h)

	for i = 1:Nx
		idx = 2+(i-1).*Ny:i*Ny-1
		y_out[idx] = y_in[idx] .* (1/h)
	end

	return y_out

end

function Hx(y_in, Nx, Ny, h)
	N = Nx*Ny
    y_out = zeros(N)

    idx = 1:Ny
	y_out[idx] = h .* (1/2)*y_in[idx]

    idx = Ny+1:N-Ny
    y_out[idx] = h .* 1*y_in[idx]

    idx = N-Ny+1:N
	y_out[idx] = h .* (1/2)*y_in[idx]

    return y_out

end

function Hy(y_in, Nx, Ny, h)
	N = Nx*Ny
    y_out = zeros(N)

    idx = 1:Ny:N-Ny+1
	y_out[idx] = h .* (1/2)*y_in[idx]

    idx = Ny:Ny:N
	y_out[idx] = h .* (1/2)*y_in[idx]

    for i = 1:Nx
        idx = 2+(i-1).*Ny:i*Ny-1
        y_out[idx] = h .* y_in[idx]
    end

    return y_out

end

function FACEtoVOL(y_in_face, face, Nx, Ny)
	N = Nx*Ny
	y_out = zeros(N)

	if face == 1
		idx = 1:Ny:N-Ny+1
	elseif face == 2
		idx = Ny:Ny:N
	elseif face == 3
		idx = 1:Ny
	elseif face == 4
		idx = N-Ny+1:N
	else
	end

	y_out[idx] = y_in_face

	return y_out

end

function VOLtoFACE(y_in, face, Nx, Ny)
	N = Nx*Ny
    y_out = zeros(N)

    if face == 1
        idx = 1:Ny:N-Ny+1
    elseif face == 2
        idx = Ny:Ny:N
    elseif face == 3
        idx = 1:Ny
    elseif face == 4
        idx = N-Ny+1:N
    else
    end

	y_out[idx] = y_in[idx]
    return y_out
end

function Bx(Nx,Ny)
	N = Nx*Ny
	y_out = zeros(N)

	idx = 1:Ny
	y_out[idx] = -1 .* ones(Ny)

	idx = N-Ny+1:N
	y_out[idx] = 1 .* ones(Ny)
	return y_out
end

function By(Nx,Ny)
	N = Nx*Ny
	y_out = zeros(N)

	idx = 1:Ny:N-Ny+1
	y_out[idx] = -1 .* ones(Ny)

	idx = Ny:Ny:N
	y_out[idx] = 1 .* ones(Ny)
	return y_out
end

function BxSx(y_in, Nx, Ny, h)
	N = Nx*Ny
	y_out = zeros(N)

	idx = 1:Ny
	y_out[idx] = (1/h) .* (1.5 .* y_in[idx] - 2 .* y_in[idx .+ Ny] + 0.5 .* y_in[idx .+ 2*Ny])
	y_out[N-Ny.+ idx] = (1/h) .* (0.5 .* y_in[N-3*Ny .+ idx] - 2 .* y_in[N-2*Ny .+ idx] + 1.5 .* y_in[N-Ny .+ idx])

	return y_out

end

function BySy(y_in, Nx, Ny, h)
	N = Nx*Ny
	y_out = zeros(N)

	idx = 1:Ny:N-Ny+1
	y_out[idx] = (1/h) .* (1.5 .* y_in[idx] - 2 .* y_in[idx .+ 1] + 0.5 .* y_in[idx .+ 2])

	idx = Ny:Ny:N
	y_out[idx] = (1/h) .* (0.5 .* y_in[idx .- 2] - 2 .* y_in[idx .- 1] + 1.5 .* y_in[idx])

	return y_out
end

function BxSx_tran(y_in, Nx, Ny, h)
	N = Nx*Ny
	y_out = zeros(N)

	idx1 = 1:Ny
	y_out[idx1] += (1.5 .* y_in[idx1]) .* (1/h)
	idx = Ny+1:2*Ny
	y_out[idx] += (-2 .* y_in[idx1]) .* (1/h)
	idx  = 2*Ny+1:3*Ny
	y_out[idx] += (0.5 .* y_in[idx1]) .* (1/h)

	idxN = N-Ny+1:N
	y_out[idxN] += (1.5 .* y_in[idxN]) .* (1/h)
	idx = N-2*Ny+1:N-Ny
	y_out[idx] += (-2 .* y_in[idxN]) .* (1/h)
	idx = N-3*Ny+1:N-2*Ny
	y_out[idx] += (0.5 .* y_in[idxN]) .* (1/h)

	return y_out
end

function BySy_tran(y_in, Nx, Ny, h)
	N = Nx*Ny
	y_out = zeros(N)

	idx1 = 1:Ny:N-Ny+1
	y_out[idx1] += (1.5 .* y_in[idx1]) .* (1/h)
	idx = 2:Ny:N-Ny+2
	y_out[idx] += (-2 .* y_in[idx1]) .* (1/h)
	idx = 3:Ny:N-Ny+3
	y_out[idx] += (0.5 .* y_in[idx1]) .* (1/h)

	idxN = Ny:Ny:N
	y_out[idxN] += (1.5 .* y_in[idxN]) .* (1/h)
	idx = Ny-1:Ny:N-1
	y_out[idx] += (-2 .* y_in[idxN]) .* (1/h)
	idx = Ny-2:Ny:N-2
	y_out[idx] += (0.5 .* y_in[idxN]) .* (1/h)

	return y_out
end


text:
import StatsFuns: log2π

@node MvNormalMeanPrecision Stochastic [ out, (μ, aliases = [ mean ]), (Λ, aliases = [ invcov, precision ]) ]

conjugate_type(::Type{ <: MvNormalMeanPrecision }, ::Type{ Val{ :out } }) = MvNormalMeanPrecision
conjugate_type(::Type{ <: MvNormalMeanPrecision }, ::Type{ Val{ :μ } })   = MvNormalMeanPrecision
conjugate_type(::Type{ <: MvNormalMeanPrecision }, ::Type{ Val{ :Λ } })   = Wishart

@average_energy MvNormalMeanPrecision (q_out::Any, q_μ::Any, q_Λ::Any) = begin
    (m_mean, v_mean) = mean(q_μ), cov(q_μ)
    (m_out, v_out)   = mean(q_out), cov(q_out)
    0.5 * (ndims(q_out) * log2π + logdet(cholinv(mean(q_Λ))) + tr(mean(q_Λ)*(v_out + v_mean + (m_out - m_mean)*(m_out - m_mean)')))
end

@average_energy MvNormalMeanPrecision (q_out_μ::Any, q_Λ::Any) = begin
    (m, V) = mean(q_out_μ), cov(q_out_μ)
    d = Int64(ndims(q_out_μ)/2)
    @views 0.5*(d*log2π + logdet(cholinv(mean(q_Λ))) + tr(mean(q_Λ)*( V[1:d,1:d] - V[1:d,d+1:end] - V[d+1:end,1:d] + V[d+1:end,d+1:end] + (m[1:d] - m[d+1:end])*(m[1:d] - m[d+1:end])' ) ))
end


text:
import Base.zero

# trace type for pairwise alignment
const Trace = Int8

# trace values
const TRACE_NONE         = Trace(0)
const TRACE_MATCH        = Trace(1)
const TRACE_INSERT       = Trace(2)
const TRACE_DELETE       = Trace(3)
const TRACE_CODON_INSERT = Trace(4)
const TRACE_CODON_DELETE = Trace(5)

const OFFSETS = ([1, 1],  # sub
                 [1, 0],  # insertion
                 [0, 1],  # deletion
                 [3, 0],  # codon insertion
                 [0, 3])  # codon deletion

function offset_forward(move::Trace, i::Int, j::Int)
    (a, b) = OFFSETS[move]
    return (i + a, j + b)
end

function offset_backward(move::Trace, i::Int, j::Int)
    (a, b) = OFFSETS[move]
    return (i - a, j - b)
end

function update_helper(final_score::Score, final_move::Trace,
                       move_score::Score, move::Trace,
                       newcols::Array{Score,2},
                       A::BandedArray{Score},
                       i::Int, j::Int, acol::Int)
    prev_i, prev_j = offset_backward(move, i, j)
    rangecol = min(prev_j, size(A)[2])
    if inband(A, prev_i, rangecol)
        score = if acol < 1 || prev_j <= acol
            A[prev_i, prev_j] + move_score
        else
            newcols[prev_i, prev_j - acol] + move_score
        end
        if score > final_score
            return score, move
        end
    end
    return final_score, final_move
end

function update(A::BandedArray{Score},
                i::Int, j::Int,
                s_base::DNA, t_base::DNA,
                pseq::RifrafSequence;
                newcols::Array{Score,2}=Array{Score}((0, 0)),
                doreverse::Bool=false,
                acol::Int=-1, trim::Bool=false,
                skew_matches::Bool=false)
    final_score = Score(-Inf)
    final_move = TRACE_NONE

    nrows, ncols = size(A)
    seqlen = length(pseq)
    # TODO: this cannot make mismatches preferable to codon indels
    seq_i = doreverse ? min(seqlen, seqlen - (i - 1) + 1) : max(i - 1, 1)
    del_i = doreverse ? nrows - i + 1 : i
    match_score = (s_base == t_base) ? pseq.match_scores[seq_i] : pseq.mismatch_scores[seq_i]
    ins_score = pseq.ins_scores[seq_i]
    del_score = pseq.del_scores[del_i]

    if skew_matches && s_base != t_base
        match_score *= 0.99
    end
    # allow terminal insertions for free
    if trim && ((j == 1) || (j == ncols))
        ins_score = 0.0
    end
    final_score, final_move = update_helper(final_score, final_move,
                                            match_score, TRACE_MATCH,
                                            newcols, A, i, j, acol)
    final_score, final_move = update_helper(final_score, final_move,
                                            ins_score, TRACE_INSERT,
                                            newcols, A, i, j, acol)
    final_score, final_move = update_helper(final_score, final_move,
                                            del_score, TRACE_DELETE,
                                            newcols, A, i, j, acol)

    if do_codon_moves(pseq)
        if do_codon_ins(pseq) && i > CODON_LENGTH
            codon_i = i - CODON_LENGTH
            if doreverse
                codon_i = length(pseq.codon_ins_scores) - codon_i + 1
            end
            codon_ins_score = pseq.codon_ins_scores[codon_i]
            final_score, final_move = update_helper(final_score, final_move,
                                                    codon_ins_score, TRACE_CODON_INSERT,
                                                    newcols, A, i, j, acol)
        end
        if do_codon_del(pseq) && j > CODON_LENGTH
            codon_del_score = pseq.codon_del_scores[del_i]
            final_score, final_move = update_helper(final_score, final_move,
                                                    codon_del_score, TRACE_CODON_DELETE,
                                                    newcols, A, i, j, acol,)
        end
    end
    if final_score == Score(-Inf)
        error("new score is invalid")
    end
    if final_move == TRACE_NONE
        error("failed to find a move")
    end
    return final_score, final_move
end

function forward_moves!(t::DNASeq, s::RifrafSequence,
                        result::BandedArray{Score},
                        moves::BandedArray{Trace};
                        trim::Bool=false,
                        skew_matches::Bool=false)
    new_shape = (length(s) + 1, length(t) + 1)
    newbandwidth!(result, s.bandwidth)
    newbandwidth!(moves, s.bandwidth)
    resize!(result, new_shape)
    resize!(moves, new_shape)
    result[1, 1] = Score(0.0)
    moves[1, 1] = TRACE_NONE
    nrows, ncols = new_shape
    for j = 1:ncols
        start, stop = row_range(result, j)
        for i = start:stop
            if i == 1 && j == 1
                continue
            end
            sbase = i > 1 ? s.seq[i - 1] : DNA_Gap
            tbase = j > 1 ? t[j - 1] : DNA_Gap
            result[i, j], moves[i, j] = update(result, i, j, sbase, tbase, s;
                                               doreverse=false,
                                               trim=trim,
                                               skew_matches=skew_matches)
        end
    end
end

"""Does backtracing to find best alignment."""
function forward_moves(t::DNASeq, s::RifrafSequence;
                       trim::Bool=false,
                       skew_matches::Bool=false)
    result = BandedArray(Score, (length(s) + 1, length(t) + 1), s.bandwidth;
                         default=-Inf, initialize=false)
    moves = BandedArray(Trace, size(result), s.bandwidth, initialize=false)
    forward_moves!(t, s, result, moves,
                   trim=trim, skew_matches=skew_matches)
    return result, moves
end

function forward!(t::DNASeq, s::RifrafSequence,
                  result::BandedArray{Score};
                  doreverse=false,
                  trim::Bool=false,
                  skew_matches::Bool=false)
    new_shape = (length(s) + 1, length(t) + 1)
    newbandwidth!(result, s.bandwidth)
    resize!(result, new_shape)
    result[1, 1] = Score(0.0)
    nrows, ncols = new_shape
    for j = 1:ncols
        start, stop = row_range(result, j)
        for i = start:stop
            if i == 1 && j == 1
                continue
            end
            sbase = i > 1 ? s.seq[doreverse ? length(s) - (i - 1) + 1 : i - 1] : DNA_Gap
            tbase = j > 1 ? t[doreverse ? length(t) - (j - 1) + 1 : j - 1] : DNA_Gap
            result[i, j], _ = update(result, i, j, sbase, tbase, s;
                                     doreverse=doreverse,
                                     trim=trim,
                                     skew_matches=skew_matches)
        end
    end
end

"""
F[i, j] is the log probability of aligning s[1:i-1] to t[1:j-1].

"""
function forward(t::DNASeq, s::RifrafSequence;
                 doreverse::Bool=false,
                 trim::Bool=false,
                 skew_matches::Bool=false)
    result = BandedArray(Score, (length(s) + 1, length(t) + 1), s.bandwidth;
                         initialize=false, default=-Inf)
    forward!(t, s, result; doreverse=doreverse, trim=trim,
             skew_matches=skew_matches)
    return result
end

function backward!(t::DNASeq, s::RifrafSequence,
                   result::BandedArray{Score})
    # this actually seems faster than a dedicated backwards alignment
    # algorithm. possibly because of cache effects.
    forward!(t, s, result, doreverse=true)
    flip!(result)
end

"""
B[i, j] is the log probability of aligning s[i:end] to t[j:end].

"""
function backward(t::DNASeq, s::RifrafSequence)
    result = forward(t, s, doreverse=true)
    flip!(result)
    return result
end

function backtrace_indices(moves::BandedArray{Trace};
                           start::Tuple{Int,Int}=(0, 0))
    result = Tuple{Int,Int}[]
    i, j = start
    if i == 0 || j == 0
        i, j = size(moves)
    end
    while i > 1 || j > 1
        m = moves[i, j]
        i, j = offset_backward(m, i, j)
        push!(result, (i, j))
    end
    return reverse(result)
end

function backtrace(moves::BandedArray{Trace})
    taken_moves = Trace[]
    i, j = size(moves)
    while i > 1 || j > 1
        m = moves[i, j]
        push!(taken_moves, m)
        i, j = offset_backward(m, i, j)
    end
    return reverse(taken_moves)
end

function count_errors(arr::BandedArray{Trace},
                      t::DNASeq, s::DNASeq)
    moves = backtrace(arr)
    a, b = moves_to_aligned_seqs(moves, t, s)
    return sum([1 for (x, y) in zip(a, b) if x != y])
end

function count_errors(t::DNASeq, s::RifrafSequence)
    _, Amoves = forward_moves(t, s, skew_matches=true)
    return count_errors(Amoves, t, s.seq)
end

"""should be close to Levenshtein"""
function edit_distance(t::DNASeq, s::DNASeq)
    log_ps = fill(-1.0, length(s))
    bandwidth = Int(ceil(min(length(t), length(s)) * 0.5))
    scores = Scores(ErrorModel(1.0, 1.0, 1.0))
    seq = RifrafSequence(s, log_ps, bandwidth, scores)
    _, Amoves = forward_moves(t, seq, skew_matches=true)
    return count_errors(Amoves, t, s)
end

function band_tolerance(Amoves::BandedArray{Trace})
    nrows, ncols = size(Amoves)
    dist = nrows
    i, j = nrows, ncols
    while i > 1 || j > 1
        start, stop = row_range(Amoves, j)
        if start > 1
            dist = min(dist, abs(i - start))
        end
        if stop < nrows
            dist = min(dist, abs(i - stop))
        end
        i, j = offset_backward(Amoves[i, j], i, j)
    end
    start, stop = row_range(Amoves, j)
    if start > 1
        dist = min(dist, abs(i - start))
    end
    if stop < nrows
        dist = min(dist, abs(i - stop))
    end
    return dist
end

function moves_to_aligned_seqs(moves::Vector{Trace},
                               t::DNASeq, s::DNASeq)
    aligned_t = DNASequence()
    aligned_s = DNASequence()
    i, j = (0, 0)
    for move in moves
        i, j = offset_forward(move, i, j)
        if move == TRACE_MATCH
            push!(aligned_t, t[j])
            push!(aligned_s, s[i])
        elseif move == TRACE_INSERT
            push!(aligned_t, DNA_Gap)
            push!(aligned_s, s[i])
        elseif move == TRACE_DELETE
            push!(aligned_t, t[j])
            push!(aligned_s, DNA_Gap)
        elseif move == TRACE_CODON_INSERT
            append!(aligned_t, DNASequence([DNA_Gap, DNA_Gap, DNA_Gap]))
            append!(aligned_s, DNASequence([s[i - 2], s[i - 1], s[i]]))
        elseif move == TRACE_CODON_DELETE
            append!(aligned_t, DNASequence([t[j - 2], t[j - 1], t[j]]))
            append!(aligned_s, DNASequence([DNA_Gap, DNA_Gap, DNA_Gap]))
        end
    end
    return aligned_t, aligned_s
end

"""
    moves_to_indices(moves, tlen, slen)

Compute index vector mapping positions in `t` to positions in `s`.

If a position in `t` is an insertion relative to `s`, it maps to the
last mapped base in `s` (or to `0`, if it is the first base in `t`).

"""
function moves_to_indices(moves::Vector{Trace},
                          tlen::Int, slen::Int)
    result = Int[]
    i, j = (0, 0)
    last_j = 0
    for move in moves
        i, j = offset_forward(move, i, j)
        if j > last_j
            push!(result, i)
            last_j = j
        end
    end
    return result
end

function align_moves(t::DNASeq, s::RifrafSequence;
                     trim::Bool=false,
                     skew_matches::Bool=false)
    A, Amoves = forward_moves(t, s;
                              trim=trim,
                              skew_matches=skew_matches)
    return backtrace(Amoves)
end

function align(t::DNASeq, s::RifrafSequence;
               trim::Bool=false,
               skew_matches::Bool=false)
    moves = align_moves(t, s;
                        trim=trim,
                        skew_matches=skew_matches)
    return moves_to_aligned_seqs(moves, t, s.seq)
end

function align(t::DNASeq, s::DNASeq, phreds::Vector{Phred},
               scores::Scores,
               bandwidth::Int;
               trim::Bool=false,
               skew_matches::Bool=false)
    moves = align_moves(t, RifrafSequence(s, phreds, bandwidth, scores),
                        trim=trim, skew_matches=skew_matches)
    return moves_to_aligned_seqs(moves, t, s)
end

function align_score(t::DNASeq, s::DNASeq, phreds::Vector{Phred},
                     scores::Scores,
                     bandwidth::Int;
                     trim::Bool=false,
                     skew_matches::Bool=false)
    A, Amoves = forward_moves(t, RifrafSequence(s, phreds, bandwidth, scores);
                              trim=trim,
                              skew_matches=skew_matches)
    moves = backtrace(Amoves)
    return moves_to_aligned_seqs(moves, t, s), A[end, end]
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

const => = Pair

"""
    Pair(x, y)
    x => y

Construct a `Pair` object with type `Pair{typeof(x), typeof(y)}`. The elements
are stored in the fields `first` and `second`. They can also be accessed via
iteration (but a `Pair` is treated as a single "scalar" for broadcasting operations).

See also [`Dict`](@ref).

# Examples
```jldoctest
julia> p = "foo" => 7
"foo" => 7

julia> typeof(p)
Pair{String, Int64}

julia> p.first
"foo"

julia> for x in p
           println(x)
       end
foo
7
```
"""
Pair, =>

eltype(p::Type{Pair{A, B}}) where {A, B} = Union{A, B}
iterate(p::Pair, i=1) = i > 2 ? nothing : (getfield(p, i), i + 1)
indexed_iterate(p::Pair, i::Int, state=1) = (getfield(p, i), i + 1)

hash(p::Pair, h::UInt) = hash(p.second, hash(p.first, h))

==(p::Pair, q::Pair) = (p.first==q.first) & (p.second==q.second)
isequal(p::Pair, q::Pair) = isequal(p.first,q.first) & isequal(p.second,q.second)

isless(p::Pair, q::Pair) = ifelse(!isequal(p.first,q.first), isless(p.first,q.first),
                                                             isless(p.second,q.second))
getindex(p::Pair,i::Int) = getfield(p,i)
getindex(p::Pair,i::Real) = getfield(p, convert(Int, i))
reverse(p::Pair{A,B}) where {A,B} = Pair{B,A}(p.second, p.first)

firstindex(p::Pair) = 1
lastindex(p::Pair) = 2
length(p::Pair) = 2
first(p::Pair) = p.first
last(p::Pair) = p.second

convert(::Type{Pair{A,B}}, x::Pair{A,B}) where {A,B} = x
function convert(::Type{Pair{A,B}}, x::Pair) where {A,B}
    Pair{A,B}(convert(A, x[1]), convert(B, x[2]))
end

promote_rule(::Type{Pair{A1,B1}}, ::Type{Pair{A2,B2}}) where {A1,B1,A2,B2} =
    Pair{promote_type(A1, A2), promote_type(B1, B2)}

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
# these require decoding to work correctly

isinteger(a::Float65) = isinteger(reflect(a.fp))

function (-){T<:Float65}(a::T)
    value = a.fp
    value = reflect(value)
    value = -value
    if getstate(a)
        value = setstate(value)
    end
    value
end

for fn in (:abs,)
   @eval begin

      function ($fn){T<:Float65}(a::T)
          value = (T)( ($fn)(reflect(a.fp)) )
          if getstate(a)
              value = setstate(value)
          end
          value
      end
   end
end

const sqrt65 = Dict( Tiny => Tiny, Huge => Huge )
function sqrt{T<:Float65}(a::T)
    get(sqrt65, a, sqrt65fn(a))
end
function sqrt65fn{T<:Float65}(a::T)
    value = (T)( ($fn)(reflect(a.fp)) )
    if getstate(a)
        value = setstate(value)
    end
    (T)( value )
end

for op in (:(+), :(-), :(*), :(/), :(\), :(%), :(^),
           :div, :fld, :rem, :mod, :mod1, :rem1, :fld1)

    @eval begin

        function ($op){T<:Float65}(a::T, b::T)
            value = (T)( ($op)(reflect(a.fp), reflect(b.fp)) )
            state = getstate(a) | getstate(b)
            if state
               value = setstate(value)
            end
            value
         end
    end
end

for op in (:(+), :(-), :(*), :(/), :(\), :(%), :(^),
           :div, :fld, :rem, :mod, :mod1, :rem1, :fld1)
    @eval begin

        function ($op){T<:Float65}(a::T, b::Float64)
            value = (T)( ($op)(reflect(a.fp), b) )
            if getstate(a)
               value = setstate(value)
            end
            value
         end

        function ($op){T<:Float65}(a::Float64, b::T)
            value = (T)( ($op)(a, reflect(b.fp)) )
            if getstate(b)
               value = setstate(value)
            end
            value
         end

        ($op){T<:Float65,I<:Union{Int64,Int32}}(a::T, b::I) =
            getstate(a) ? setstate(Float65( $op(reflect(a.fp), b) )) : Float65( $op(reflect(a.fp), b) )
        ($op){T<:Float65,I<:Union{Int64,Int32}}(a::I, b::T) =
            getstate(b) ? setstate(Float65( $op(a, reflect(b.fp)) )) : Float65( $op(a, reflect(b.fp)) )
        ($op){T<:Float65}(a::T, b::Bool) =
            getstate(a) ? setstate(Float65( $op(reflect(a.fp), convert(Float64,b)) )) : Float65( $op(reflect(a.fp), convert(Float64,b)) )
        ($op){T<:Float65}(a::Bool, b::T) =
            getstate(b) ? setstate(Float65( $op(convert(Float64,a), reflect(b.fp)) )) : Float65( $op(convert(Float64,a), reflect(b.fp)) )
        ($op){T<:Float65,R<:Real}(a::T, b::R) =
            getstate(a) ? setstate(Float65( $op(reflect(a.fp), convert(Float64,b)) )) : Float65( $op(reflect(a.fp), convert(Float64,b)) )
        ($op){T<:Float65,R<:Real}(a::R, b::T) =
            getstate(b) ? setstate(Float65( $op(convert(Float64,a), reflect(b.fp)) )) : Float65( $op(convert(Float64,a), reflect(b.fp)) )
    end
end

function (eps){T<:Float64}(a::T)
    value = (T)( (eps)(reflect(a.fp)) )
    if getstate(a)
       value = setstate(value)
    end
    value
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 59c9b880-c5a5-11eb-3fbb-8f6c9389cbbd
md"""
### Test Code for Module KNN
"""

# ╔═╡ 44bde682-a0d6-4e06-8123-ff380bec0266
module script include("./script.jl") end

# ╔═╡ 4ba5bce5-61f5-4e34-9627-2971112ce3fc
KNN = script.KNN

# ╔═╡ 0a12f2ac-459b-4ead-ae71-ded368bcbb5a
import Random

# ╔═╡ 32bb3911-1d61-4ca4-bba6-af98b3d6e545
import Plots

# ╔═╡ 6febcd4a-facd-4301-9183-10853be03633
module tools include("../tools.jl") end

# ╔═╡ 4aa17a82-5872-42d1-9315-620678eb1634
JuTools = tools.JuTools

# ╔═╡ 541ebf6c-0059-4305-ae2b-a16331fa5232
begin
# prepare data, no need to scale
linear_func = m -> 2 * m - 20
X_data, Y_data = JuTools.data_generate_linear_2d(linear_func=linear_func, data_size=500,
    range_min=0.0, range_max=100.0, random_scale=20.0)
size(X_data), size(Y_data)
end

# ╔═╡ 9b8a6f72-ac55-4567-bc4c-387602df923e
let
# plot the data
plot_X1 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_Y1 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_X2 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_Y2 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_X3 = [i for i in 0:100]
plot_Y3 = [linear_func(m) for m in plot_X3]
Plots.gr()
Plots.scatter(plot_X1, plot_Y1, leg=false, c="red", background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot_X2, plot_Y2, leg=false, c="blue")
Plots.plot!(plot_X3, plot_Y3, leg=false, c="orange", linewidth=5)
Plots.xlims!((0.0, 100.0))
Plots.ylims!((0.0, 100.0))
end

# ╔═╡ 62b49cde-f3fc-42b4-bbfb-ce2d01248925
md"""
### Naive Approach
"""

# ╔═╡ 0a4cc82c-4e14-4224-a172-d0fd3643db27
let
plot_X1 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_Y1 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_X2 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_Y2 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_X3 = [i for i in 0:100]
plot_Y3 = [linear_func(m) for m in plot_X3]
# plot contour
X_test_1 = 0:1:100
X_test_2 = 0:1:100
f_test(X1, X2) = begin
   KNN.predict_naive([X1, X2], X_data, Y_data, K=5)[1]
end
Plots.gr()
Plots.contour(X_test_1, X_test_2, f_test, leg=false, fill=true, c=:matter,
    linewidth=0, background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot_X1, plot_Y1, leg=false, c="red")
Plots.scatter!(plot_X2, plot_Y2, leg=false, c="blue")
Plots.plot!(plot_X3, plot_Y3, leg=false, c="orange", linewidth=5)
Plots.xlims!((0.0, 100.0))
Plots.ylims!((0.0, 100.0))
end

# ╔═╡ ecfeb65a-b7bd-4941-8b83-179860c3bf23
let
plot_X1 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_Y1 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_X2 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_Y2 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_X3 = [i for i in 0:100]
plot_Y3 = [linear_func(m) for m in plot_X3]
# test with different Ks
X_test_1 = 0:1:100
X_test_2 = 0:1:100
f_test(X1, X2) = begin
   KNN.predict_naive([X1, X2], X_data, Y_data, K=2)[1]
end

Plots.gr()
Plots.contour(X_test_1, X_test_2, f_test, leg=false, fill=true, c=:matter,
    linewidth=0, background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot_X1, plot_Y1, leg=false, c="red")
Plots.scatter!(plot_X2, plot_Y2, leg=false, c="blue")
Plots.plot!(plot_X3, plot_Y3, leg=false, c="orange", linewidth=5)
Plots.xlims!((0.0, 100.0))
Plots.ylims!((0.0, 100.0))
end

# ╔═╡ cdc0071f-a7c6-4d17-b559-ddbf777b9969
let
plot_X1 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_Y1 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_X2 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_Y2 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_X3 = [i for i in 0:100]
plot_Y3 = [linear_func(m) for m in plot_X3]
X_test_1 = 0:1:100
X_test_2 = 0:1:100
f_test(X1, X2) = begin
   KNN.predict_naive([X1, X2], X_data, Y_data, K=10)[1]
end

Plots.gr()
Plots.contour(X_test_1, X_test_2, f_test, leg=false, fill=true, c=:matter,
    linewidth=0, background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot_X1, plot_Y1, leg=false, c="red")
Plots.scatter!(plot_X2, plot_Y2, leg=false, c="blue")
Plots.plot!(plot_X3, plot_Y3, leg=false, c="orange", linewidth=5)
Plots.xlims!((0.0, 100.0))
Plots.ylims!((0.0, 100.0))
end

# ╔═╡ 86bcd18d-4562-4ba6-bdb8-c45b2a935540
md"""
With larger `K`, the algorithm can find the boundary more easily
"""

# ╔═╡ bcaaf451-1430-43d5-b737-87129dd0579b
md"""
Now compare with euclidean distance metric
"""

# ╔═╡ 22719609-1a57-4fe2-a1d5-0d7494c54aaf
let
plot_X1 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_Y1 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_X2 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_Y2 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_X3 = [i for i in 0:100]
plot_Y3 = [linear_func(m) for m in plot_X3]
X_test_1 = 0:1:100
X_test_2 = 0:1:100
f_test(X1, X2) = begin
   KNN.predict_naive_fun([X1, X2], X_data, Y_data, K=10)[1]
end

Plots.gr()
Plots.contour(X_test_1, X_test_2, f_test, leg=false, fill=true, c=:matter,
    linewidth=0, background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot_X1, plot_Y1, leg=false, c="red")
Plots.scatter!(plot_X2, plot_Y2, leg=false, c="blue")
Plots.plot!(plot_X3, plot_Y3, leg=false, c="orange", linewidth=5)
Plots.xlims!((0.0, 100.0))
Plots.ylims!((0.0, 100.0))
end

# ╔═╡ 8bba7dab-bbf6-439f-a091-b5f4cfe41dbf
let
plot_X1 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_Y1 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_X2 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_Y2 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_X3 = [i for i in 0:100]
plot_Y3 = [linear_func(m) for m in plot_X3]
X_test_1 = 0:1:100
X_test_2 = 0:1:100
f_test(X1, X2) = begin
   KNN.predict_naive_fun([X1, X2], X_data, Y_data, K=2)[1]
end

Plots.gr()
Plots.contour(X_test_1, X_test_2, f_test, leg=false, fill=true, c=:matter,
    linewidth=0, background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot_X1, plot_Y1, leg=false, c="red")
Plots.scatter!(plot_X2, plot_Y2, leg=false, c="blue")
Plots.plot!(plot_X3, plot_Y3, leg=false, c="orange", linewidth=5)
Plots.xlims!((0.0, 100.0))
Plots.ylims!((0.0, 100.0))
end

# ╔═╡ 8f6f1ae4-9b17-4972-9fa0-d9684e79771b
let
plot_X1 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_Y1 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_X2 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_Y2 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_X3 = [i for i in 0:100]
plot_Y3 = [linear_func(m) for m in plot_X3]
X_test_1 = 0:1:100
X_test_2 = 0:1:100
f_test(X1, X2) = begin
   KNN.predict_naive_fun([X1, X2], X_data, Y_data, K=30)[1]
end

Plots.gr()
Plots.contour(X_test_1, X_test_2, f_test, leg=false, fill=true, c=:matter,
    linewidth=0, background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot_X1, plot_Y1, leg=false, c="red")
Plots.scatter!(plot_X2, plot_Y2, leg=false, c="blue")
Plots.plot!(plot_X3, plot_Y3, leg=false, c="orange", linewidth=5)
Plots.xlims!((0.0, 100.0))
Plots.ylims!((0.0, 100.0))
end

# ╔═╡ 64a60947-e048-4d53-8b9c-14dc76e5cd0b
md"""
#### Note

From the comparison, we see the difference between euclidean distance and vector angle in terms of measuring data point similarity

The following `KdTree` and `BallTree` are all measured based on euclidean distance
So we will see similar boundary shape
"""

# ╔═╡ 5ee1cef1-376c-4401-b197-29e242403bb0
md"""
### K-d Tree Approach
"""

# ╔═╡ 96296a72-8f2c-4ea0-a476-31f929dda576
kdtree = KNN.create_kdtree(X_data, Y_data)

# ╔═╡ 18bbb469-96b2-4f19-8c7a-4b863a73871a
let
plot_X1 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_Y1 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_X2 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_Y2 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_X3 = [i for i in 0:100]
plot_Y3 = [linear_func(m) for m in plot_X3]
X_test_1 = 0:1:100
X_test_2 = 0:1:100
f_test(X1, X2) = begin
   KNN.predict_kdtree([X1, X2], kdtree, K=10)[1]
end

Plots.gr()
Plots.contour(X_test_1, X_test_2, f_test, leg=false, fill=true, c=:matter,
    linewidth=0, background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot_X1, plot_Y1, leg=false, c="red")
Plots.scatter!(plot_X2, plot_Y2, leg=false, c="blue")
Plots.plot!(plot_X3, plot_Y3, leg=false, c="orange", linewidth=5)
Plots.xlims!((0.0, 100.0))
Plots.ylims!((0.0, 100.0))
end

# ╔═╡ 8cce776f-5853-4246-a039-47eedeaea266
let
plot_X1 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_Y1 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_X2 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_Y2 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_X3 = [i for i in 0:100]
plot_Y3 = [linear_func(m) for m in plot_X3]
X_test_1 = 0:1:100
X_test_2 = 0:1:100
f_test(X1, X2) = begin
   KNN.predict_kdtree([X1, X2], kdtree, K=2)[1]
end

Plots.gr()
Plots.contour(X_test_1, X_test_2, f_test, leg=false, fill=true, c=:matter,
    linewidth=0, background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot_X1, plot_Y1, leg=false, c="red")
Plots.scatter!(plot_X2, plot_Y2, leg=false, c="blue")
Plots.plot!(plot_X3, plot_Y3, leg=false, c="orange", linewidth=5)
Plots.xlims!((0.0, 100.0))
Plots.ylims!((0.0, 100.0))
end

# ╔═╡ d5a9ebda-e772-4888-bc8f-3faff0413aad
let
plot_X1 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_Y1 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_X2 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_Y2 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_X3 = [i for i in 0:100]
plot_Y3 = [linear_func(m) for m in plot_X3]
X_test_1 = 0:1:100
X_test_2 = 0:1:100
f_test(X1, X2) = begin
   KNN.predict_kdtree([X1, X2], kdtree, K=30)[1]
end

Plots.gr()
Plots.contour(X_test_1, X_test_2, f_test, leg=false, fill=true, c=:matter,
    linewidth=0, background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot_X1, plot_Y1, leg=false, c="red")
Plots.scatter!(plot_X2, plot_Y2, leg=false, c="blue")
Plots.plot!(plot_X3, plot_Y3, leg=false, c="orange", linewidth=5)
Plots.xlims!((0.0, 100.0))
Plots.ylims!((0.0, 100.0))
end

# ╔═╡ d311739d-2455-448f-b343-06ecdf9d725b
md"""
With larger `K`, the boundary is more linear
"""

# ╔═╡ ee593cd1-21bc-4b33-bf87-fd5254676e9d
md"""
### Ball Tree Approach
"""

# ╔═╡ dc3f28e4-7acc-45bf-b720-fe33d8d6bb3f
balltree = KNN.create_balltree(X_data, Y_data)

# ╔═╡ 6642011a-15ee-4b5b-afb6-fd7ac990b503
let
plot_X1 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_Y1 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_X2 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_Y2 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_X3 = [i for i in 0:100]
plot_Y3 = [linear_func(m) for m in plot_X3]
X_test_1 = 0:1:100
X_test_2 = 0:1:100
f_test(X1, X2) = begin
   KNN.predict_balltree([X1, X2], balltree, K=10)[1]
end

Plots.gr()
Plots.contour(X_test_1, X_test_2, f_test, leg=false, fill=true, c=:matter,
    linewidth=0, background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot_X1, plot_Y1, leg=false, c="red")
Plots.scatter!(plot_X2, plot_Y2, leg=false, c="blue")
Plots.plot!(plot_X3, plot_Y3, leg=false, c="orange", linewidth=5)
Plots.xlims!((0.0, 100.0))
Plots.ylims!((0.0, 100.0))
end

# ╔═╡ 9e702dff-6ab7-4e9b-8dcb-de2afd85ef80
let
plot_X1 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_Y1 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_X2 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_Y2 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_X3 = [i for i in 0:100]
plot_Y3 = [linear_func(m) for m in plot_X3]
X_test_1 = 0:1:100
X_test_2 = 0:1:100
f_test(X1, X2) = begin
   KNN.predict_balltree([X1, X2], balltree, K=2)[1]
end

Plots.gr()
Plots.contour(X_test_1, X_test_2, f_test, leg=false, fill=true, c=:matter,
    linewidth=0, background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot_X1, plot_Y1, leg=false, c="red")
Plots.scatter!(plot_X2, plot_Y2, leg=false, c="blue")
Plots.plot!(plot_X3, plot_Y3, leg=false, c="orange", linewidth=5)
Plots.xlims!((0.0, 100.0))
Plots.ylims!((0.0, 100.0))
end

# ╔═╡ 1720df3b-7290-4b74-b4b7-6eb7b68e551d
let
plot_X1 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_Y1 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_X2 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_Y2 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_X3 = [i for i in 0:100]
plot_Y3 = [linear_func(m) for m in plot_X3]
X_test_1 = 0:1:100
X_test_2 = 0:1:100
f_test(X1, X2) = begin
   KNN.predict_balltree([X1, X2], balltree, K=30)[1]
end

Plots.gr()
Plots.contour(X_test_1, X_test_2, f_test, leg=false, fill=true, c=:matter,
    linewidth=0, background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot_X1, plot_Y1, leg=false, c="red")
Plots.scatter!(plot_X2, plot_Y2, leg=false, c="blue")
Plots.plot!(plot_X3, plot_Y3, leg=false, c="orange", linewidth=5)
Plots.xlims!((0.0, 100.0))
Plots.ylims!((0.0, 100.0))
end

# ╔═╡ 655441e7-ff51-4f7c-839d-fa60339edadf
# compare time
@time KNN.predict_naive(X_data, X_data, Y_data, K=10)

# ╔═╡ 2af2c142-8cd3-49f2-a646-450cb49bf080
@time KNN.predict_kdtree(X_data, kdtree, K=10)

# ╔═╡ 8d4d270f-b569-42bf-ad4b-215975a98bc6
@time KNN.predict_balltree(X_data, balltree, K=10)

# ╔═╡ ff118ee9-12cc-40d7-9219-c654256abdfe
md"""
We see that K-d Tree approach is faster than naive approach
Ball Tree is similar to K-d Tree, but should be faster on larger dataset
"""

# ╔═╡ Cell order:
# ╟─59c9b880-c5a5-11eb-3fbb-8f6c9389cbbd
# ╠═44bde682-a0d6-4e06-8123-ff380bec0266
# ╠═4ba5bce5-61f5-4e34-9627-2971112ce3fc
# ╠═0a12f2ac-459b-4ead-ae71-ded368bcbb5a
# ╠═32bb3911-1d61-4ca4-bba6-af98b3d6e545
# ╠═6febcd4a-facd-4301-9183-10853be03633
# ╠═4aa17a82-5872-42d1-9315-620678eb1634
# ╠═541ebf6c-0059-4305-ae2b-a16331fa5232
# ╠═9b8a6f72-ac55-4567-bc4c-387602df923e
# ╟─62b49cde-f3fc-42b4-bbfb-ce2d01248925
# ╠═0a4cc82c-4e14-4224-a172-d0fd3643db27
# ╠═ecfeb65a-b7bd-4941-8b83-179860c3bf23
# ╠═cdc0071f-a7c6-4d17-b559-ddbf777b9969
# ╟─86bcd18d-4562-4ba6-bdb8-c45b2a935540
# ╟─bcaaf451-1430-43d5-b737-87129dd0579b
# ╠═22719609-1a57-4fe2-a1d5-0d7494c54aaf
# ╠═8bba7dab-bbf6-439f-a091-b5f4cfe41dbf
# ╠═8f6f1ae4-9b17-4972-9fa0-d9684e79771b
# ╟─64a60947-e048-4d53-8b9c-14dc76e5cd0b
# ╟─5ee1cef1-376c-4401-b197-29e242403bb0
# ╠═96296a72-8f2c-4ea0-a476-31f929dda576
# ╠═18bbb469-96b2-4f19-8c7a-4b863a73871a
# ╠═8cce776f-5853-4246-a039-47eedeaea266
# ╠═d5a9ebda-e772-4888-bc8f-3faff0413aad
# ╟─d311739d-2455-448f-b343-06ecdf9d725b
# ╟─ee593cd1-21bc-4b33-bf87-fd5254676e9d
# ╠═dc3f28e4-7acc-45bf-b720-fe33d8d6bb3f
# ╠═6642011a-15ee-4b5b-afb6-fd7ac990b503
# ╠═9e702dff-6ab7-4e9b-8dcb-de2afd85ef80
# ╠═1720df3b-7290-4b74-b4b7-6eb7b68e551d
# ╠═655441e7-ff51-4f7c-839d-fa60339edadf
# ╠═2af2c142-8cd3-49f2-a646-450cb49bf080
# ╠═8d4d270f-b569-42bf-ad4b-215975a98bc6
# ╟─ff118ee9-12cc-40d7-9219-c654256abdfe

max_stars_repo_name: teamclouday/JuliaLearn

max_issues_repo_name: teamclouday/JuliaLearn

max_forks_repo_name: teamclouday/JuliaLearn


text:
module Islander

export Enmap, enfft, enifft, enfft!, enifft!
include("Pixell.jl")

using DelimitedFiles
using Distributions
using FFTW
"""
Draws an inverse gamma distribution given the sky map spectrum.
"""
function inverse_gamma(sigma_ell::Array{T,1}) where {T}
    lmax = size(sigma_ell, 1)
    gauss = Normal()
    Cl_realization = zeros(T, lmax)
    for ell in 1:lmax
        rho = 0.0
        for i in 1:2*ell-1
            rho += convert(T,rand(gauss)^2)
        end
        Cl_realization[ell] = (2 * ell + 1) * sigma_ell[ell] / rho
    end
    return Cl_realization
end

"""
Utility function to load ``C_{\\ell}^{TT}`` when working in the notebooks.
"""
function util_get_cltt(lmax=10000)
    camb = readdlm("../data/camb_theory.dat", skipstart=1)
    camb_ells = convert.(Int64, camb[1:lmax-1,1])
    camb_cltt = camb[1:lmax-1,2]
    theory_ells = collect(1:lmax)
    theory_cltt = zeros(lmax)
    theory_cltt[ camb_ells ] = (camb_cltt
        .* 2π ./ camb_ells ./ (camb_ells .+ 1))
    return theory_ells, theory_cltt
end

"""
Utility function which prepends some zeros to an array. It makes a copy instead
of modifying the input.
"""
function ellpad(arr::Array{T,N}; nzeros=1) where {T,N}
    result = arr[:]
    pushfirst!(result, zeros(T, nzeros)...)
    return result
end

"""
Quick and dirty binned, flat-sky power spectra.
"""
function ell_binpow( power, ell, Δell, ell_max=2500; use_Dell=false )
    nbin = convert(Int, ceil(ell_max / Δell)) + 1
    bin_centers = zeros(nbin)
    bin_means = zeros(nbin)
    for i in 1:nbin
        bin_left = (i-1) * Δell
        bin_right = bin_left + Δell
        bin_centers[i] = (bin_left + bin_right) / 2.0
        filt = (ell .> bin_left) .& (ell .<= bin_right)
        bin_data = power[filt]
        if use_Dell
            bin_data .*= ell[filt]
        end
        bin_means[i] = sum(bin_data) / size(bin_data,1)
    end

    bin_centers, bin_means
end

end # module


text:
using MinimalRides
using Plots
using DataFrames
using Dates
using Underscores

import MinimalRides: make_subset, results, DistanceSmoothing

RAW_DATA_URL_UPDATED = ("https://dl.dropboxusercontent.com" *
    "/sh/hibzl6fkzukltk9/AABSMicBJlwMnlmA3ljt1uY5a" *
    "/data_samples-json2.zip")
CACHE_DATA = "/media/win/Data/data_samples-json2.zip"

rides = MinimalRides.load(RAW_DATA_URL_UPDATED, CACHE_DATA);

DAYS = Date.(["2020-05-21"])
TYPES = [:freight]

subset = make_subset(rides, DAYS, TYPES)

df = results(subset; smoothalg = DistanceSmoothing(0.5), search_radius = 2.5 * 1.2) |> DataFrame

pairs_df = @_ sort(df, :op, rev = true) |> filter((_.op > 0.75) & (_.len > 400), __)

println("Best pairs: ")
println(pairs_df)

for i in axes(pairs_df, 1)
    plot(map(x -> (x.coord.lat, x.coord.lon), subset[pairs_df[i, :i]].route),
        title = "Rides ($(pairs_df[i, :i]), $(pairs_df[i, :j]))",
        label = "Ride $(pairs_df[i, :i])")
    d = plot!(map(x -> (x.coord.lat, x.coord.lon), subset[pairs_df[i, :j]].route),
        label = "Ride $(pairs_df[i, :j])")
    display(d)
end

using BenchmarkTools

# 37 routes, 666 (37*36/2) pairs
println("Benchmark results: ")
show(stdout, MIME"text/plain"(), @benchmark(results($subset)))