text:
include("./spp_dm1.jl")

# ---------------------------------------------------------------------------
# GRASP

function greedyRandomizedConstruction(C, A, a)
    x = zeros(Int64, length(C))
    S = 1:length(C)
    M = 1:size(A, 1)
    E = Array{Array{Int64,1},1}(undef, length(C))
    Einit = Array{Array{Int64,1},1}(undef, length(C))
    F = Array{Array{Int64,1},1}(undef, length(C))
    U = Array{Float64,1}(undef, length(C))     # profit
    firstloop = true
    while (length(S) != 0)
        fill!(E, Int64[])
        fill!(F, Int64[])
        fill!(U, -Inf)
        if length(S) == 1
            iSel = S[1]
        else
            #mettre à jour E et F
            for j in M
                current = intersect(findall(isequal(1), A[j, :]), S)
                for i in current
                    E[i] = union(E[i], setdiff(current, i))
                    F[i] = union(F[i], j)
                end
            end
            #calculer uMax, uMin et Ui
            uMax = -Inf
            uMin = Inf
            for i in S
                if length(E[i]) != 0
                    U[i] = C[i] - sum(C[v] for v in E[i])
                else
                    U[i] = C[i]
                end
                if U[i] < uMin
                    uMin = U[i]
                end
                if U[i] > uMax
                    uMax = U[i]
                end
            end
            #calculer limit, RCL et choisir iSel aléatoire
            limit = uMin + a*(uMax - uMin)
            RCL = findall(>=(limit),U)
            iSel = RCL[rand(1:length(RCL))]
        end
        if firstloop
            Einit = copy(E)
            firstloop = false
        end
        #mettre a jour le problème
        x[iSel] = 1
        S = setdiff(S, union(E[iSel], iSel))
        M = setdiff(M, F[iSel])
    end
    z = sum(C[i]x[i] for i in 1:length(C))
    return x, z, Einit
end

function GRASP(C, A, a, finish)
    zMax = 0
    xMax = zeros(Int64, length(C))
    t = 0
    while (t <= finish)
        t = t + @elapsed begin
            xInit, zInit, Einit = greedyRandomizedConstruction(C, A, a)
            x, z = GreedyImprovement(C, A, xInit, zInit, Einit)
        end
        if z > zMax
            zMax = z
            xMax = x
        end
    end
    return xMax, zMax
end

# ---------------------------------------------------------------------------
# ReactiveGRASP

# Fonction tirer aléatoire un élément suivant son poids
function draw(a,p)
    k = findfirst(cumsum(p) .> rand())
    return k
end

function reactiveGRASP(C, A, a, finish, N)
    zMin = 0
    m = length(a)
    zTotal = Array{Pair{Int64, Int64}}(undef, m)        #pair la sommes de z[k] avec nb fois utilise a[k]
    fill!(zTotal, 0 => 0)
    zMax = 0
    xMax = zeros(Int64, length(C))
    p = Array{Float64,1}(undef, m)
    fill!(p, 1/m)
    q = Array{Float64,1}(undef, m)
    nbIter = 1
    t = 0
    while (t <= finish)
        #modifier p
        if (nbIter - 1)%N == 0 && nbIter != 1
            if zMax == zMin             #optimal found
                break
            else
                sumQ = 0
                for k in 1:m
                    q[k] = ((zTotal[k].first / zTotal[k].second) - zMin)/(zMax - zMin)
                    sumQ = sumQ + q[k]
                end
                for k in 1:m
                    p[k] = q[k]/sumQ
                end
            end
        end
        #tirer a[k]
        if nbIter <= m  #assurer que chaque a[k] a utilisé au moins une fois
            k = nbIter
        else
            k = draw(a,p)
        end
        #GRASP
        t = t + @elapsed begin
            xInit, zInit, Einit = greedyRandomizedConstruction(C, A, a[k])
            x, z = GreedyImprovement(C, A, xInit, zInit, Einit)
        end
        #mettre à jours les valeurs de x, z et nbIter
        if nbIter == 1
            zMax = z
            xMax = x
            zMin = z
        elseif z > zMax
            zMax = z
            xMax = x
        elseif z < zMin
            zMin = z
        end
        zTotal[k] = zTotal[k].first + z => zTotal[k].second + 1
        nbIter += 1
    end
    return xMax, zMax
end

# -------------------------------------------------------------------------------------------
# Version GRASP et ReactiveGRASP avec la condition d'arrêt est nombre des itérations
# L'ultiliser pour le graphique

function graspSPP(fname, a, nbIterationGrasp)
    target = "../../../Data"
    C, A = loadSPP(string(target,"/",fname))

    zconstruction = zeros(Int64,nbIterationGrasp)
    zamelioration = zeros(Int64,nbIterationGrasp)
    zbest = zeros(Int64,nbIterationGrasp)
    zbetter = 0

    for i=1:nbIterationGrasp
        xInit, zInit, Einit = greedyRandomizedConstruction(C, A, a)
        x, z = GreedyImprovement(C, A, xInit, zInit, Einit)
        zconstruction[i] = zInit
        zamelioration[i] = z
        zbetter = max(zbetter, zamelioration[i])
        zbest[i] = zbetter
    end
    return zconstruction, zamelioration, zbest
end

function regraspSPP(fname, a, nbIterationGrasp, N)
    target = "../../../Data"
    C, A = loadSPP(string(target,"/",fname))

    zconstruction = zeros(Int64,nbIterationGrasp)
    zamelioration = zeros(Int64,nbIterationGrasp)
    zbest = zeros(Int64,nbIterationGrasp)
    zbetter = 0

    zMin = 0
    m = length(a)
    zTotal = Array{Pair{Int64, Int64}}(undef, m)        #pair la sommes de z[k] avec nb fois utilise a[k]
    fill!(zTotal, 0 => 0)
    zMax = 0
    xMax = zeros(Int64, length(C))
    p = Array{Float64,1}(undef, m)
    fill!(p, 1/m)
    q = Array{Float64,1}(undef, m)
    nbIter = 1
    while (nbIter <= nbIterationGrasp)
        #modifier p
        if (nbIter - 1)%N == 0 && nbIter != 1
            if zMax == zMin             #optimal found
                break
            else
                sumQ = 0
                for k in 1:m
                    q[k] = ((zTotal[k].first / zTotal[k].second) - zMin)/(zMax - zMin)
                    sumQ = sumQ + q[k]
                end
                for k in 1:m
                    p[k] = q[k]/sumQ
                end
            end
        end
        #tirer a[k]
        if nbIter <= m  #assurer que chaque a[k] a utilisé au moins une fois
            k = nbIter
        else
            k = draw(a,p)
        end
        #GRASP
        xInit, zInit, Einit = greedyRandomizedConstruction(C, A, a[k])
        x, z = GreedyImprovement(C, A, xInit, zInit, Einit)
        zconstruction[nbIter] = zInit
        zamelioration[nbIter] = z
        zbetter = max(zbetter, zamelioration[nbIter])
        zbest[nbIter] = zbetter
        #mettre à jours les valeurs de x, z et nbIter
        if nbIter == 1
            zMax = z
            xMax = x
            zMin = z
        elseif z > zMax
            zMax = z
            xMax = x
        elseif z < zMin
            zMin = z
        end
        zTotal[k] = zTotal[k].first + z => zTotal[k].second + 1
        nbIter += 1
    end
    return zconstruction, zamelioration, zbest, p
end

max_stars_repo_name: ducanh-le/metaSPPstu

max_issues_repo_name: ducanh-le/metaSPPstu

max_forks_repo_name: ducanh-le/metaSPPstu


text:
using CLArrays, GLVisualize, GPUArrays, GLAbstraction, GeometryTypes

# source: https://github.com/johnfgibson/julia-pde-benchmark/blob/master/1-Kuramoto-Sivashinksy-benchmark.ipynb
function inner_ks(n, IFFT!, FFT!, Nt, Nn, Nn1, u, G, A_inv, B, dt2, dt32, uslice, U)
    Nn1 .= Nn       # shift nonlinear term in time
    Nn .= u         # put u into Nn in prep for comp of nonlinear term

    IFFT! * Nn

    # plotting
    uslice .= real.(Nn) ./ 10f0
    U[1:Nt, n] = reshape(Array(uslice), (Nt, 1)) # copy from gpu to opengl gpu not implemented for now

    # transform Nn to gridpt values, in place
    Nn .= Nn .* Nn   # collocation calculation of u^2
    FFT!*Nn        # transform Nn back to spectral coeffs, in place

    Nn .= G .* Nn    # compute Nn == -1/2 d/dx (u^2) = -u u_x

    # loop fusion! Julia translates the folling line of code to a single for loop.
    u .= A_inv .* (B .* u .+ dt32 .* Nn .- dt2 .* Nn1)
end

T = Float32; AT = CLArray
N = 1000
Lx = T(64*pi)
Nx = T(N)
dt = T(1/16)

x = Lx*(0:Nx-1)/Nx
u = T.(cos.(x) + 0.1*sin.(x/8) + 0.01*cos.((2*pi/Lx)*x))

u = AT((T(1)+T(0)im)*u)             # force u to be complex
Nx = length(u)                      # number of gridpoints
kx = T.(vcat(0:Nx/2-1, 0:0, -Nx/2+1:-1))# integer wavenumbers: exp(2*pi*kx*x/L)
alpha = T(2)*pi*kx/Lx                  # real wavenumbers:    exp(alpha*x)

D = T(1)im*alpha                       # spectral D = d/dx operator

L = alpha.^2 .- alpha.^4            # spectral L = -D^2 - D^4 operator

G = AT(T(-0.5) .* D)               # spectral -1/2 D operator, to eval -u u_x = 1/2 d/dx u^2

# convenience variables
dt2  = T(dt/2)
dt32 = T(3*dt/2)
A_inv = AT((ones(T, Nx) - dt2*L).^(-1))
B = AT(ones(T, Nx) + dt2*L)

# compute in-place FFTW plans
FFT! = plan_fft!(u)
IFFT! = plan_ifft!(u)

# compute nonlinear term Nn == -u u_x
powed = u .* u
Nn = G .* fft(powed);    # Nn == -1/2 d/dx (u^2) = -u u_x
Nn1 = copy(Nn);        # Nn1 = Nn at first time step
FFT! * u;

uslice = real(u)
U = zeros(Float32, N, N)

w = glscreen(resolution = (1920, 1080))

robj = visualize(
    U, :surface, color_norm = Vec2f0(-0.1, 0.1),
    ranges = ((-3f0, 3f0), (-3f0, 3f0))
).children[]
Ugpu = robj[:position_z]

# setup camera and view object
_view(robj, position = Vec3f0(-4.33, 3.8, -2.7), lookat = Vec3f0(-0.23, 0.5, 0.7))
cam = w.cameras[:perspective]
push!(cam.up, Vec3f0(0.3, -0.33, -0.9))

io, buffer = GLVisualize.create_video_stream(homedir()*"/Desktop/pd2.mkv", w)
for n in 1:N
    isopen(w) || break
    inner_ks(n, (IFFT!), (FFT!), N, Nn, Nn1, u, G, A_inv, B, dt2, dt32, uslice, Ugpu)
    GLWindow.poll_glfw()
    GLWindow.reactive_run_till_now()
    GLWindow.render_frame(w)
    GLWindow.swapbuffers(w)
    (n % 4 == 0) && GLVisualize.add_frame!(io, w, buffer)
end
close(io)
GLWindow.destroy!(w)

ext: jl


text:
# Functions for simulating data from GBMs
#
# Copyright (c) 2020: Jeffrey W. Miller.
# This file is released under the MIT "Expat" License.

module Sim

export simulate_parameters_and_data, permute_to_match_latent_dimensions!

using Distributions
using Statistics: mean
using LinearAlgebra: pinv, qr

# NegativeBinomial distribution with mean=mu and dispersion=1/r.
#   r = inverse dispersion
#   p = r/(mu+r)  # success probability
# Note that Julia uses the Wolfram definition which is that r is the number of successes until stopping.
NegBin(mu,r) = NegativeBinomial(r, r/(mu+r))

# Standardize x to zero mean, unit variance
center_scale(x,dim=1) = (z = x .- mean(x,dims=dim); z./sqrt.(mean(z.^2,dims=dim)))

# Enforce bound on absolute values of x
bound(x,maxabsval) = max(min(x,maxabsval),-maxabsval)

function simulate_covariate_matrix(I,K; maxabsval=100, cdist="Normal")
    if cdist=="Normal"
        cdistn = Normal()
    elseif cdist=="Bernoulli"
        cdistn = Bernoulli()
    elseif cdist=="Gamma"
        cdistn = Gamma(2,1/sqrt(2))
    else
        @error("Unknown value of cdist for simulating data.")
    end

    if cdist=="Normal"
        X = randn(I,K)*randn(K,K)
    else
        # The following line is mathematically equivalent to the approach used, and is algorithmically identical to the Normal case due to the use of center_scale:
        # W = randn(I,K); Q = randn(K,K); C = Q'*Q; sigma = sqrt.(diag(C)); X = (W*Q)./sigma'

        Q = randn(K,K); C = Q'*Q; sigma = sqrt.(diag(C)); C = (C./sigma)./sigma'; C = (C+C')/2  # generate a random correlation matrix
        X = rand(MvNormal(C),I)'  # feature covariates
        X = quantile.(cdistn,cdf.(Normal(),X))  # copula transformation
    end
    X = bound.(X,maxabsval)  # enforce bound on magnitude
    X[:,1] .= 1.0  # constant for intercepts
    X[:,2:end] = center_scale(X[:,2:end],1)  # standardize the covariates
    return X
end

# Simulate true parameters and data
function simulate_parameters_and_data(I,J,K,L,M; maxabsval=100, omega0=-2.3, odist="NegativeBinomial", cdist="Normal", pdist="Normal")
    if pdist=="Normal"
        pdistn = Normal()
    elseif pdist=="Gamma"
        pdistn = Gamma(2,1/sqrt(2))
    else
        @error("Unknown value of pdist for simulating data.")
    end

    # Simulate covariates
    X = simulate_covariate_matrix(I,K; maxabsval=100, cdist=cdist)  # feature covariates
    Z = simulate_covariate_matrix(J,L; maxabsval=100, cdist=cdist)  # sample covariates

    # Simulate parameters
    A0 = rand(pdistn,J,K)/(2*sqrt(K))  # true feature coefficients
    B0 = rand(pdistn,I,L)/(2*sqrt(L))  # true sample coefficients
    C0 = rand(pdistn,K,L)/sqrt(K*L); if K>0; C0[1,1] = C0[1,1] + 3.0; end  # true intercepts and interaction coefficients
    B0 = B0 - X*(pinv(X)*B0)  # enforce constraint that X'*B = 0 by projecting onto nullspace of X'.
    A0 = A0 - Z*(pinv(Z)*A0)  # enforce constraint that Z'*A = 0 by projecting onto nullspace of Z'.

    U0 = Matrix(qr(randn(I,M)).Q)  # generate U0 uniformly from the Stiefel manifold
    V0 = Matrix(qr(randn(J,M)).Q)  # generate V0 uniformly from the Stiefel manifold
    if M > 0
        U0 = U0 - X*(pinv(X)*U0)  # enforce constraint that X'*U = 0 by projecting onto nullspace of X'.
        V0 = V0 - Z*(pinv(Z)*V0)  # enforce constraint that Z'*V = 0 by projecting onto nullspace of Z'.
    end
    D0 = (sqrt(I)+sqrt(J))*(collect(0:M-1)/max(M-1,1) .+ 1)  # use evenly spaced D values following the Marchenko-Pastur scaling

    # omega0 = true overall log-dispersion
    S0 = randn(I,1); S0 = S0 .- log(mean(exp.(S0)))  # true feature-specific log-dispersion offsets
    T0 = randn(J,1); T0 = T0 .- log(mean(exp.(T0)))  # true sample-specific log-dispersion offsets

    # Simulate data
    Mu0 = exp.(X*A0' + B0*Z' + X*C0*Z' + U0*(D0.*V0'))
    r0 = exp.(-S0.-T0'.-omega0)
    if odist=="NB"
        Y = rand.(NegBin.(Mu0,r0))
    elseif odist=="Poisson"
        Y = rand.(Poisson.(Mu0))
    elseif odist=="Geometric"
        Y = rand.(Geometric.(1.0./(Mu0 .+ 1)))
    elseif odist=="LNP"
        Sigma2 = log.(1.0./r0 .+ 1.0)
        epsilon = rand.(Normal.(-0.5*Sigma2, sqrt.(Sigma2)))
        Y = rand.(Poisson.(Mu0.*exp.(epsilon)))
    else
        @error("Unknown value of odist for simulating data.")
    end

    return Y,X,Z,A0,B0,C0,D0,U0,V0,S0,T0,omega0
end

# Match latent dimensions and latent factor signs
function permute_to_match_latent_dimensions!(D,U,V,D0,U0,V0)
    M = length(D)
    D1,U1,V1 = copy(D),copy(U),copy(V)
    mlist0 = collect(1:M)
    mlist1 = collect(1:M)
    for m = 1:M
        Cor = cor(U0[:,mlist0],U1[:,mlist1])
        index = findmax(abs.(Cor))[2]
        i0,i1 = index[1],index[2]
        sgn = sign(Cor[i0,i1])
        m0,m1 = mlist0[i0],mlist1[i1]
        U[:,m0] = U1[:,m1]*sgn
        V[:,m0] = V1[:,m1]*sgn
        D[m0] = D1[m1]
        deleteat!(mlist0,i0)
        deleteat!(mlist1,i1)
    end
end

end # module Sim


text:
include("Shape.jl")

"""
    morph_to(to_func::Function; object=:stroke)

A closure for the [`_morph_to`](@ref) function.
This makes it easier to write the function inside an `Object`.

Currently morphing is quite simple and only works for basic shapes.
It especially does not work with functions which produce more than one polygon
or which produce filled polygons.
Blending between fills of polygons is definitely coming at a later stage.

**Important:** The functions itself should not draw the polygon
i.e. use `circle(Point(100,100), 50)` instead of `circle(Point(100,100), 50, :stroke)`

# Arguments
- `to_func::Function`: Same as `from_func` but it defines the "result" polygon,
                       which will be displayed at the end of the Object

# Keywords
- `do_action::Symbol` defines whether the object has a fill or just a stroke. Defaults to `:stroke`.

# Example

This creates a star that morphs into a circle and back.

```
astar(args...; do_action=:stroke) = star(O, 50, 5, 0.5, 0, do_action)
acirc(args...; do_action=:stroke) = circle(Point(100,100), 50, do_action)

video = Video(500, 500)
back = Background(1:20, ground)
star_obj = Object(1:10, astar)
act!(star_obj, Action(linear(), morph_to(acirc)))
circle_obj = Object(11:20, acirc)
act!(circle_obj, Action(:same, morph_to(astar)))
```
"""
function morph_to(to_func::Function; do_action = :stroke)
    return (video, object, action, frame) ->
        _morph_to(video, object, action, frame, to_func; do_action = do_action)
end

"""
    _morph_to(video::Video, object::Object, action::Action, frame, to_func::Function; do_action=:stroke)

Internal version of [`morph_to`](@ref) but described there.
"""
function _morph_to(
    video::Video,
    object::Object,
    action::Action,
    frame,
    to_func::Function;
    do_action = :stroke,
)
    newpath()
    object.func(video, object, frame; do_action = :none)
    closepath()
    from_polys = pathtopoly()

    newpath()
    to_func(video, object, frame; do_action = :none)
    closepath()
    to_polys = pathtopoly()

    return morph_between(video, action, frame, from_polys, to_polys; do_action = do_action)
end

"""
    morph_between(video::Video, action::Action, frame,
        from_polys::Vector{Vector{Point}}, to_polys::Vector{Vector{Point}};
        do_action=:stroke)

Internal version of [`morph_to`](@ref) after the from poly is defined.
"""
function morph_between(
    video::Video,
    action::Action,
    frame,
    from_polys::Vector{Vector{Point}},
    to_polys::Vector{Vector{Point}};
    do_action = :stroke,
)
    cs = get_current_setting()
    cs.show_object = false

    # computation of the polygons and the best way to morph in the first frame
    if frame == first(get_frames(action))
        save_morph_polygons!(action, from_polys, to_polys)
    end

    # obtain the computed polygons. These polygons have the same number of points.
    number_of_shapes = length(action.defs[:from_shape])
    bool_move = ones(Bool, number_of_shapes)
    bool_appear = zeros(Bool, number_of_shapes)

    for si in 1:number_of_shapes
        from_shape = action.defs[:from_shape][si]
        to_shape = action.defs[:to_shape][si]
        inter_shape = action.defs[:inter_shape][si]

        if isempty(from_shape) || isempty(to_shape)
            if isempty(to_shape)
                action.defs[:inter_shape][si] = from_shape
            else
                action.defs[:inter_shape][si] = to_shape
                bool_appear[si] = true
            end
            bool_move[si] = false
            continue
        end

        # compute the interpolation variable `t` for the current frame
        t = get_interpolation(action, frame)
        interpolate_shape!(inter_shape, from_shape, to_shape, t)
        draw_shape(inter_shape, do_action)
    end

    # let new paths appear
    t = get_interpolation(action, frame)
    setopacity(t)

    shape_ids = [si for si in 1:number_of_shapes if (!bool_move[si] && bool_appear[si])]
    shapes = action.defs[:inter_shape][shape_ids]
    draw_shape.(shapes, do_action)

    # let old paths disappear
    t = get_interpolation(action, frame)
    setopacity(1 - t)

    shape_ids = [si for si in 1:number_of_shapes if (!bool_move[si] && !bool_appear[si])]
    shapes = action.defs[:inter_shape][shape_ids]
    draw_shape.(shapes, do_action)

    setopacity(1)
end

"""
    match_num_points(poly_1::Vector{Point}, poly_2::Vector{Point})

This is a helper function for [`morph_to`](@ref).
Given two polygons `poly_1` and `poly_2` points are added to the polygon with less points
until both polygons have the same number of points.
The polygon with less points gets mutated during this process.

# Arguments
- `poly_1::Vector{Point}`: The points which define the first polygon
- `poly_2::Vector{Point}`: The points which define the second polygon
"""
function match_num_points(poly_1::Vector{Point}, poly_2::Vector{Point})
    l1 = length(poly_1)
    l2 = length(poly_2)
    # if both have the same number of points => we are done
    l1 == l2 && return poly_1, poly_2

    new_poly_1 = simplify(poly_1)
    new_poly_2 = simplify(poly_2)
    l1 = length(new_poly_1)
    l2 = length(new_poly_2)

    # poly_1 should have less points than poly_2 so we flip if this is not the case
    flipped = false
    if l1 > l2
        new_poly_1, new_poly_2 = new_poly_2, new_poly_1
        l1, l2 = l2, l1
        flipped = true
    end
    # the difference of the length of points
    missing_nodes = l2 - l1

    add_points!(new_poly_1, missing_nodes)

    if flipped
        new_poly_1, new_poly_2 = new_poly_2, new_poly_1
    end
    @assert length(new_poly_1) == length(new_poly_2)
    return new_poly_1, new_poly_2
end

"""
    add_points!(poly, missing_nodes)

Add #`missing_nodes` to poly.
"""
function add_points!(poly, missing_nodes)
    pdists = polydistances(poly)
    every_dist = pdists[end] / missing_nodes
    pdiffs = diff(pdists)

    ct = 0.0
    poly_orig = copy(poly)
    npi = 1
    for pi in 1:length(poly_orig)
        add_nodes = convert(Int, round(pdiffs[pi] / every_dist))
        t = pdiffs[pi] / (add_nodes + 1)
        ct_local = ct
        for i in 1:add_nodes
            ct_local += t
            new_point = get_polypoint_at(poly_orig, ct_local / pdists[end]; pdist = pdists)
            npi += 1
            insert!(poly, npi, new_point)
        end
        npi += 1
        missing_nodes -= add_nodes
        ct += pdiffs[pi]
        every_dist = (pdists[end] - ct) / missing_nodes
    end
end

"""
    save_morph_polygons!(action::Action, from_func::Vector{Vector{Point}},
                                         to_func::Vector{Vector{Point}})

Calls the functions to polygons and calls [`match_num_points`](@ref)
such that both polygons have the same number of points.
This is done once inside [`_morph_to`](@ref).
Saves the two polygons inside `action.defs[:from_poly]` and `action.defs[:to_poly]`.

**Assumption:** Both functions create the same number of polygons.
"""
function save_morph_polygons!(
    action::Action,
    from_polys::Vector{Vector{Point}},
    to_polys::Vector{Vector{Point}},
)
    # delete polygons with less than 2 points
    for i in length(from_polys):-1:1
        length(from_polys[i]) <= 1 && splice!(from_polys, i)
    end
    for i in length(to_polys):-1:1
        length(to_polys[i]) <= 1 && splice!(to_polys, i)
    end

    if length(from_polys) != length(to_polys)
        try_merge_polygons(from_polys)
        try_merge_polygons(to_polys)
    end

    from_shapes = create_shapes(from_polys)
    to_shapes = create_shapes(to_polys)

    if length(from_shapes) > 1
        from_shapes, to_shapes = reorder_match(from_shapes, to_shapes)
    end

    action.defs[:from_shape] = Vector{Shape}()
    action.defs[:to_shape] = Vector{Shape}()
    action.defs[:inter_shape] = Vector{Shape}()

    counter = 0
    for (from_shape, to_shape) in zip(from_shapes, to_shapes)
        counter += 1
        if isempty(from_shape) || isempty(to_shape)
            if isempty(from_shape)
                push!(action.defs[:from_shape], EmptyShape())
                push!(action.defs[:to_shape], to_shape)
                push!(
                    action.defs[:inter_shape],
                    Shape(
                        length(to_shape.points),
                        [length(subpath) for subpath in to_shape.subpaths],
                    ),
                )
            else
                push!(action.defs[:from_shape], from_shape)
                push!(action.defs[:to_shape], EmptyShape())
                push!(
                    action.defs[:inter_shape],
                    Shape(
                        length(from_shape.points),
                        [length(subpath) for subpath in from_shape.subpaths],
                    ),
                )
            end
        else
            from_shape, to_shape = prepare_to_interpolate(from_shape, to_shape)

            push!(action.defs[:from_shape], from_shape)
            push!(action.defs[:to_shape], to_shape)
            push!(
                action.defs[:inter_shape],
                Shape(
                    length(from_shape.points),
                    [length(subpath) for subpath in from_shape.subpaths],
                ),
            )
        end
    end
end

"""
    try_merge_polygons(polys)

Try to merge polygons together to match the number of polygons that get morphed.
The only example I encountered is that the `[` of a 3xY matrix consists of 3 parts which
are merged together.
"""
function try_merge_polygons(polys)
    for pi in 1:length(polys)
        !ispolyclockwise(polys[pi]) && continue
        for pj in 1:length(polys)
            pi >= pj && continue
            !ispolyclockwise(polys[pj]) && continue
            smallest_dist = Inf
            found_Ap = O
            found_Bp = O
            start_with_A = true
            for (pA, pointA1, pointA2) in
                zip(1:length(polys[pi]), polys[pi], circshift(polys[pi], -1))
                for (pB, pointB1, pointB2) in
                    zip(1:length(polys[pj]), polys[pj], circshift(polys[pj], -1))
                    dist1 = distance(pointA1, pointB2)
                    dist2 = distance(pointA2, pointB1)
                    if dist1 <= 3 && dist2 <= 3
                        if isinside((pointB1 + pointB2) / 2, polys[pi]; allowonedge = true)
                            pB1 = pB + 1 <= length(polys[pj]) ? pB + 1 : 1
                            polys[pj][pB1] = pointA1 - (pointB2 - pointA1)
                            polys[pj][pB] = pointA2 - (pointB1 - pointA2)
                            pointB2 = polys[pj][pB1]
                            pointB1 = polys[pj][pB]
                        end
                        if ispolyclockwise([pointA1, pointA2, pointB1, pointB2])
                            start_with_A = false
                            smallest_dist = dist1
                            found_Ap = pA
                            found_Bp = pB + 1 <= length(polys[pj]) ? pB + 1 : 1
                        else
                            start_with_A = true
                            smallest_dist = dist2
                            found_Ap = pA + 1 <= length(polys[pi]) ? pA + 1 : 1
                            found_Bp = pB + 1
                        end
                    end
                end
            end
            if smallest_dist <= 3
                if start_with_A
                    tmp_polyA = circshift(polys[pi], length(polys[pi]) - found_Ap + 1)
                    tmp_polyB = circshift(polys[pj], length(polys[pj]) - found_Bp)
                    polys[pi] = vcat(tmp_polyA, tmp_polyB[1:(end - 1)])
                    polys[pj] = Point[]
                else
                    tmp_polyB = circshift(polys[pj], length(polys[pj]) - found_Bp + 1)
                    tmp_polyA = circshift(polys[pi], length(polys[pi]) - found_Ap)
                    polys[pi] = vcat(tmp_polyB, tmp_polyA)
                    polys[pj] = Point[]
                end
            end
        end
    end
    for pi in length(polys):-1:1
        if isempty(polys[pi])
            splice!(polys, pi)
        end
    end
end

"""
    compute_shortest_morphing_dist(from_poly::Vector{Point}, to_poly::Vector{Point})

Rotates `from_poly` internally to check which mapping produces the smallest morphing distance.
It returns the start index of the rotation of `from_poly` as well as the smallest distance.
"""
function compute_shortest_morphing_dist(from_poly::Vector{Point}, to_poly::Vector{Point})
    # find the smallest morphing distance to match the points in a more natural way
    # smallest_i holds the best starting point of from_path
    smallest_i = 1
    smallest_distance = typemax(Float64)

    for i in 1:length(from_poly)
        overall_distance = 0.0
        for j in 1:length(from_poly)
            p1 = from_poly[mod1(j + i - 1, length(from_poly))]
            p2 = to_poly[j]
            overall_distance += distance(p1, p2)
        end
        if overall_distance < smallest_distance
            smallest_distance = overall_distance
            smallest_i = i
        end
    end
    return smallest_i, smallest_distance
end

"""
    reorder_match(from_shapes::Vector{Shape}, to_shapes::Vector{Shape})

Computes the similiarty of the shapes and finds the best mapping such that the sum of similiarty
is maximized.

Additionally it creates empty shapes when needed such that
`reordered_from` and `reordered_to` contain the same number of shapes.

# Returns
- `reordered_from::Vector{Shape}`
- `reordered_to::Vector{Shape}`
"""
function reorder_match(from_shapes::Vector{Shape}, to_shapes::Vector{Shape})
    num_shapes = max(length(from_shapes), length(to_shapes))

    similiarity_matrix = zeros(length(from_shapes), length(to_shapes))
    for fi in 1:length(from_shapes)
        from_shape = from_shapes[fi]
        for ti in 1:length(to_shapes)
            to_shape = to_shapes[ti]
            similiarity_matrix[fi, ti] = -get_similarity(from_shape, to_shape)
        end
    end
    assignment, cost = hungarian(similiarity_matrix)

    new_from_shapes = Vector{Shape}(undef, num_shapes)
    for i in 1:num_shapes
        new_from_shapes[i] = EmptyShape()
    end

    new_to_shapes = Vector{Shape}(undef, num_shapes)
    for i in 1:num_shapes
        if i <= length(to_shapes)
            new_to_shapes[i] = to_shapes[i]
        else
            new_to_shapes[i] = EmptyShape()
        end
    end

    ptr = length(to_shapes) + 1
    for i in 1:length(from_shapes)
        if assignment[i] == 0
            new_from_shapes[ptr] = from_shapes[i]
            ptr += 1
        else
            new_from_shapes[assignment[i]] = from_shapes[i]
        end
    end
    return new_from_shapes, new_to_shapes
end


text:
import YAML

# This would be better in a general library!
atomicmass = Dict{AbstractString,Float64}(
"H"=>1.00794, "C" => 12.01, "N" => 14.01, "S" => 32.07, "Zn" => 65.38,
"I" => 126.9, "Br" => 79.904, "Cl" => 35.45,
"Pb" => 207.2)
#Zn is actually Sn; stupid work-around for Pymol seeing 'Sn' as 'S'
#The initial test cases which I am covering, are phonons of CH3NH3.PbI3 and SnS

mesh = YAML.load(open("mesh.yaml"))     #Phonopy mesh.yaml file; with phonons

# Native VASP POSCAR reader
P=readdlm(open("POSCAR","r"))

lattice=[ P[l,f]::Float64 for l=3:5,f=1:3 ]
println(STDERR,lattice)
species=[ P[6,f] for f=1:length(P[7,:]) ]
speciescount=[ P[7,f]::Int for f=1:length(P[7,:]) ]
NATOMS=sum(speciescount)
println(STDERR,species)
positions=[ P[l,f]::Float64 for l=9:9+NATOMS,f=1:3 ]
# The following is probably overkill, but reads the VASP atom formats + expands
# into a one dimensional string vector
#     species:    C    N    H    Pb   I
#     speciescount:   1     1     6     1     3
atomnames=AbstractString[]
for (count,specie) in zip(speciescount,species)
    for i=1:count push!(atomnames,specie) end
end
println(STDERR,atomnames)

# SUPERCELL definition
supercellexpansions=[ a*lattice[1,:] + b*lattice[2,:] + c*lattice[3,:] for a=0:1,b=0:1,c=0:1 ] #generates set of lattice vectors to apply for supercell expansions
println(STDERR,"supercellexpansions ==>",supercellexpansions)

function output_animated_xyz(eigenmode,eigenvector,freq,steps=32)
    filename= @sprintf("anim_%02d.xyz",eigenmode)
    anim=open(filename,"w")

    for phi=0:2*pi/steps:2*pi-1e-6 #slightly offset from 2pi so we don't repeat 0=2pi frame
#        projection= lattice[1][1]*positions + 2*realeigenvector*sin(phi) # this does all the maths
#        println("Lattice: ",lattice,"\n Eigenvec: ",realeigenvector)
        #projection=positions*lattice + 2*sin(phi)*realeigenvector
#        println("Projection: ",projection)

        # output routines to .xyz format
        @printf(anim,"%d\n\n",NATOMS*length(supercellexpansions)) # header for .xyz multi part files
        for i=1:NATOMS
#            println("u ==> ",projection[i,:])
#            println("Positions[",i,"]: ",positions[i,:])
#            println("Realeigenvector[",i,"]: ",realeigenvector[i,:])

# NB: Currently uncertain as to whether to apply mass weighting by:-
#   diving the eigenvector by sqrt(amu) - converting from Energy --> displacement
#   multiplying the eigenvector by sqrt(amu) - weighting the Dynamic matrix with atomic mass ?
            projection=lattice * (positions[i,:]' + 0.2 / sqrt(atomicmass[atomnames[i]]) * eigenvector[i,:]'*sin(phi))
            for supercellexpansion in supercellexpansions
                supercellprojection=projection+supercellexpansion'
                @printf(anim,"%s %f %f %f\n",atomnames[i],supercellprojection[1],supercellprojection[2],supercellprojection[3])
            end
        end
    end
    close(anim)

end

# This decomposes the amount that the different atomtypes contribute to each phonon mode, in the unit cell
function decompose_eigenmode_atomtype(eigenmode,realeigenvector,freq)
    print("Eigenmode: ",eigenmode)
    @printf("\tFreq: %.2f THz %03.1f (cm-1)\t",freq[2],freq[2]*33.36)

    #atomiccontribution = Dict{AbstractString,Float64}("Pb"=>0.0, "Br" => 0.0, "N" => 0.0, "C" => 0.0, "H" => 0.0)
    atomiccontribution=[species[i]=>0.0 for i in 1:length(species)]
    for i=1:NATOMS
        atomiccontribution[atomnames[i]]+= norm(realeigenvector[i,:])
    end

    totallength=sum(values(atomiccontribution))
    # Now weight every object
    for contri in keys(atomiccontribution) # Surely a better way that iterating over?
        atomiccontribution[contri]/=totallength
    end

    for contri in sort(species, by=x->atomicmass[x],rev=true) #Heaviest first, by lookup
        @printf("%s %.4f\t",contri,atomiccontribution[contri])
    end
#    println(atomiccontribution)
#    println(sum(values(atomiccontribution)))
    println()
end

# This outputs, for each atom in the unit cell, the contribution in terms of displacement and energy, for the phonon
function decompose_eigenmode_atom_contributions(eigenmode,realeigenvector)
    normsum=0.0
    normsummassweighted=0.0
    for i=1:NATOMS
        normsum+=norm(realeigenvector[i,:])
        normsummassweighted+=norm(realeigenvector[i,:])/sqrt(atomicmass[atomnames[i]])
    end
    println("Norm sum: ",normsum, " Norm sum(mass weighted): ",normsummassweighted)
    for i=1:NATOMS
        println("Mode: ",eigenmode," Atom: ",i," ",atomnames[i],
#        "\n",
#        " Norm: ", norm(realeigenvector[i,:]),
        " Norm(weighted): ",norm(realeigenvector[i,:])/normsum,
        " Norm(mass weighted): ",(norm(realeigenvector[i,:])/sqrt(atomicmass[atomnames[i]]))/normsummassweighted)
    end
end

# Header line for GNUPLOT, plotting mode decompositions.
@printf("Mode 1 Freq THz THz cm1 cm1 ")
for contri in sort(species, by=x->atomicmass[x],rev=true) #Heaviest first, by lookup
    @printf("%s %s\t",contri,contri)
end
@printf("\n")

# Data structure looks like: mesh["phonon"][1]["band"][2]["eigenvector"][1][2][1]
for (eigenmode,(eigenvector,freq)) in enumerate(mesh["phonon"][1]["band"])
#    println("freq (THz) ==> ",freq[2], "\tWavenumbers (cm-1, 3sf) ==> ",freq[2]*33.36)

## These functions used when figuring out form / normalisation of eigenvectors, from Phonopy mesh.yaml
#    println("phonon[\"eigenvector\"] ==>",phonon["eigenvector"])
#    println("eigenvector ==> ",eigenvector[2])
#    for atom in eigenvector[2]
#        #println ("atom ==> ",atom)
#        @printf("atom %e %e %e \n",atom[1][1],atom[2][1],atom[3][1])
#        disp=[atom[1][1],atom[2][1],atom[3][1]] #pull data out of YAML into sensible Julia Vector
#        println(disp)
#    end

    realeigenvector=[ eigenvector[2][n][d][1]::Float64 for n=1:NATOMS, d=1:3 ]
    # Array comprehension to reform mesh.yaml format into [n][d] shap
    realeigenvector=reshape(realeigenvector,NATOMS,3) # doesn't do anything?

    output_animated_xyz(eigenmode,realeigenvector,freq)

    decompose_eigenmode_atomtype(eigenmode,realeigenvector,freq)
    #decompose_eigenmode_atom_contributions(eigenmode,realeigenvector)

#=
    for I=10:12 # iodine indexes, hard coded
        println(show(positions))
        println(show(positions[I,:]))
        PbI=positions[I,:]-positions[9,:] # Vector from Pb to I
        show(PbI)
        @printf("I: %d PbI: %f %f %f Phonon-u: %f %f %f PbI.Phonon-u: %f\n",I,
        PbI[1],PbI[2],PbI[3],
        realeigenvector[I,1],realeigenvector[I,2],realeigenvector[I,3],
        dot(PbI::Float64,realeigenvector[I]::Float64) )
    end
=#
end

max_stars_repo_name: jarvist/Julia-Phonons

max_issues_repo_name: lucydot/Julia-Phonons

max_forks_repo_name: lucydot/Julia-Phonons


text:
#####
# Binning error analysis
#####

"""
    binning_error(X[; binsize=0, warnings=false])

Calculates statistical error (eff. standard deviation) for correlated data.
How: Binning of data and assuming statistical independence of bins
(i.e. R plateau has been reached). (Eq. 3.18 basically)

The default `binsize=0` indicates automatic binning.
"""
function binning_error(X::AbstractVector{T}; binsize=0, warnings=false) where T<:Real
    # Data: real numbers
    if binsize == 0
        binsize = floor(Int, length(X)/32)
        binsize = binsize==0?1:binsize
    end

    isinteger(length(X) / binsize) || !warnings ||
        warn("Non-integer number of bins $(length(X) / binsize). " *
             "Last bin will be smaller than all others.")

    bin_means = map(mean, Iterators.partition(X, binsize))
    return sqrt(1/length(bin_means) * var(bin_means))
end
function binning_error(X::AbstractVector{T}; binsize=0, warnings=false) where T<:Complex
    # Data: complex numbers
    sqrt(binning_error(real(X), binsize=binsize, warnings=warnings)^2 +
        binning_error(imag(X), binsize=binsize, warnings=warnings)^2)
end

# Data: arrays
function binning_error(X::AbstractArray{T}; binsize=0, warnings=false) where T<:Number
    ndimsX = ndims(X)
    mapslices(y->binning_error(y; binsize=binsize, warnings=warnings), X, ndimsX)[(Colon() for _ in 1:ndimsX-1)...,1]
end
function binning_error(X::AbstractVector{T}; binsize=0, warnings=false) where T<:(AbstractArray{S} where S)
    binning_error(cat(ndims(X[1])+1, X...); binsize=binsize, warnings=warnings)
end

#####
# Useful plots
#####
"""
Plots the binning error coefficient function `R(binsize)` for multiple bin sizes.
Ideally, this plot shows a plateau (Fig. 3.3), indicating that the bin averages have become independent.
For correlated data `R>≈1` and `sqrt(R)` quantifies how much one would have underestimated
the one-sigma errorbar.

Returns bss and R from `R_function(X)`.
"""
function plot_binning_R(X::Vector{T}; min_nbins=32, figsize=(4,3)) where T<:Real

    try PyPlot
    catch
        warn("PyPlot not loaded. Loading it now...")
        eval(Expr(:toplevel, Expr(:using, Symbol("PyPlot"))))
    end

    bss, R = R_function(X, min_nbins=min_nbins)
    figure(figsize=figsize)
    plot(bss, R, "m.-")
    xlabel("bin size")
    ylabel("R");

    return bss, R
end

"""
Plot histogram of X with mean and statistical errorbars (including correlation effects)
and autocorrelation function C(t).

Optional keyword `error` can be "tau_binning", "tau_integrated", or "tau_fitting", "binning",
and sets method for error estimation.
"""
function plot_error(X::Vector{T}; binsize=0, histbins=50, error="binning", figsize=(10,4), digits=3) where T<:Real

    try PyPlot
    catch
        warn("PyPlot not loaded. Loading it now...")
        eval(Expr(:toplevel, Expr(:using, Symbol("PyPlot"))))
    end

    fig, ax = subplots(1, 2, figsize=figsize)

    Xmean = mean(X)
    Xstd = std(X)

    if error == "binning"
        err = binning_error(X, binsize=binsize, warnings=false)
    elseif error == "tau_binning"
        err = tau_binning_error(X, binsize=binsize)
    elseif error == "tau_integrated"
        err = tau_integrated_error(X)
    elseif error == "tau_fitting"
        err = tau_fitting_error(X)
    end

    ax[1][:hist](X, histbins, color="gray", alpha=.5, normed=1)
    ax[1][:set_ylabel]("\$ P \$")
    ax[1][:set_xlabel]("\$ X \$")
    ax[1][:set_yticks]([])
    ax[1][:axvline](Xmean, color="black", label="\$ $(round.(Xmean, digits)) \$", linewidth=2.0)

    ax[1][:axvline](Xmean+err, color="r", label="\$ \\pm $(round.(err, digits)) \$ (binning σ)", linewidth=2.0)
    ax[1][:axvline](Xmean-err, color="r", linewidth=2.0)

    ax[1][:axvline](Xmean+2*err, color="r", alpha=.3, label="\$ \\pm $(round.(2*err, digits)) \$ (binning 2σ)", linewidth=2.0)
    ax[1][:axvline](Xmean-2*err, color="r", alpha=.3, linewidth=2.0)

    ax[1][:axvline](Xmean+Xstd, color="g", label="\$ \\pm $(round.(Xstd, digits)) \$ (std)", linewidth=2.0)
    ax[1][:axvline](Xmean-Xstd, color="g", linewidth=2.0)

    ax[1][:legend](frameon=false, loc="best")

    ax[2][:plot](autocor(X), "-", color="k", linewidth=2.0)
    ax[2][:set_xlabel]("\$ t \$")
    ax[2][:set_ylabel]("\$ C(t) \$")

    tight_layout()
    return ax
end

#####
# Calculation of error coefficient (function) R. (Ch. 3.4 in QMC book)
#####
"""
Groups datapoints in bins of varying size `bs`.
Returns the used binsizes `bss` and the error coefficient function `R(bss)` (Eq. 3.20) which should feature
a plateau, i.e. `R(bs_p) ~ R(bs)` for `bs >= bs_p`. (Fig. 3.3)

Optional keyword `min_nbins`. Only bin sizes used that lead to at least `min_nbins` bins.
"""
function R_function(X::Vector{T}; min_nbins=10) where T<:Real
    bss = Int[]
    N = length(X)
    Xmean = mean(X)
    Xvar = var(X)

    for bs in 1:N
        N%bs == 0 ? push!(bss, bs) : nothing
    end

    n_bins = Int.(N./bss)
    bss = bss[n_bins .>= min_nbins] # at least 500 bins for every binsize

    R = zeros(length(bss))
    for (i, bs) in enumerate(bss)

        blockmeans = vec(mean(reshape(X, (bs,n_bins[i])), 1))

        blocksigma2 = 1/(n_bins[i]-1)*sum((blockmeans - Xmean).^2)

        R[i] = bs * blocksigma2 / Xvar
    end

    return bss, R
end

"""
Groups datapoints in bins of fixed binsize and returns error coefficient R. (Eq. 3.20)
"""
function R_value(X::Vector{T}, binsize::Int) where T<:Real
    N = length(X)
    n_bins = div(N,binsize)
    lastbs = rem(N,binsize)
    lastbs == 0 || (lastbs >= binsize/2 || warn("Binsize leads to last bin having less than binsize/2 elements."))

    blockmeans = vec(mean(reshape(X[1:n_bins*binsize], (binsize,n_bins)), 1))
    if lastbs != 0
        vcat(blockmeans, mean(X[n_bins*binsize+1:end]))
        n_bins += 1
    end

    blocksigma2 = 1/(n_bins-1)*sum((blockmeans - mean(X)).^2)
    return binsize * blocksigma2 / var(X)
end


text:
@info "Generating Validation Table"
function generate_validation_table(filepath::AbstractString)
    descriptor = InfrastructureSystems.read_validation_descriptor(joinpath(PSYPATH,"descriptors","power_system_structs.json"))
    open(filepath, "w") do io
        write(io, "# Data Requirements\n\n")
        write(io, "|  Struct Name  |  Field Name  |  DataType  |  Min  |  Max  |  Action  |\n")
        write(io, "|---------------|--------------|------------|-------|-------|----------|\n")
        for item in descriptor
            for field in item["fields"]
                write(io, "|$(item["struct_name"])|$(field["name"])|$(field["data_type"])|")
                if haskey(field, "valid_range")
                    if field["valid_range"] isa Dict
                        valid_min = field["valid_range"]["min"]
                        valid_max = field["valid_range"]["max"]
                        for value in (valid_min, valid_max)
                            write(io, isnothing(value) ? "null" : string(value))
                            write(io, "|")
                        end
                    else
                        write(io, "$(field["valid_range"])|$(field["valid_range"])|")
                    end
                else
                    write(io, "-|-")
                end
                if haskey(field, "validation_action")
                    write(io, "$(field["validation_action"])|\n")
                else
                    write(io, "|-\n")
                end
            end
        end
    end
end

generate_validation_table(joinpath(PSYPATH, "../docs/src/man/data_requirements_table.md"))


text:
using JFSM
using NetCDF
using TimeSeries
using IcanProj
using ProgressMeter
using ExcelReaders
using DataFrames

"""
Run fsm and save results to netcdf.
"""
function run_simulation(path, tstart, tstop, am, cm, dm, em, hm, dt, res, iexp, pot_melt, comp_prec)

    # Read metadata

    id = ncread(joinpath(path, "forcings", "rainf_$(res).nc"), "id")

    id_desc = ncgetatt(joinpath(path, "forcings", "rainf_$(res).nc"), "dim_space", "id")

    # Space dimension

    dim_s_in = length(id)

    # Read time for forcings

    time_str_in = ncread(joinpath(path, "forcings", "tair_$(res).nc"), "time_str")

    time_in = DateTime.(time_str_in, "yyyy-mm-dd HH:MM:SS")

    # Starting position of forcings and time dimension

    istart = find(time_in .== tstart)[1]

    dim_t_in = sum(tstart .<= time_in .<= tstop)

    # Crop forcing time

    time_in = time_in[istart:(istart+dim_t_in-1)]

    # Split forcing time

    year_in = Dates.year.(time_in)

    month_in = Dates.month.(time_in)

    day_in = Dates.day.(time_in)

    hour_in = Dates.hour.(time_in)

    # Time and dimensions for outputs

    ts_3h = TimeArray(time_in, fill(0.0, length(time_in), 3))

    ts_24h = collapse(ts_3h, day, first, mean)

    time_str_out = Dates.format.(ts_24h.timestamp, "yyyy-mm-dd HH:MM:SS")

    dim_t_out = length(time_str_out)

    dim_s_out = dim_s_in

    # Meteorological inputs

    rainf = fill(0.0, dim_t_in)
    snowf = fill(0.0, dim_t_in)
    tair  = fill(0.0, dim_t_in)
    iswr  = fill(0.0, dim_t_in)
    ilwr  = fill(0.0, dim_t_in)
    pres  = fill(0.0, dim_t_in)
    rhum  = fill(0.0, dim_t_in)
    wind  = fill(0.0, dim_t_in)

    # Create output files

    mkpath(joinpath(path, "results_$(iexp)"))

    create_netcdf(joinpath(path, "results_$(iexp)", "hs_$(res).nc"), "hs", Dict("units" => "m"), dim_t_out, dim_s_out, time_str_out, id, id_desc)
    create_netcdf(joinpath(path, "results_$(iexp)", "swe_$(res).nc"), "swe", Dict("units" => "mm"), dim_t_out, dim_s_out, time_str_out, id, id_desc)
    create_netcdf(joinpath(path, "results_$(iexp)", "rof_$(res).nc"), "rof", Dict("units" => "unknown"), dim_t_out, dim_s_out, time_str_out, id, id_desc)

    # Read input and run model

    @showprogress 1 "Running..." for i in 1:dim_s_in

        # Read forcings

        ncread!(joinpath(path, "forcings", "rainf_$(res).nc"), "rainf", rainf, start = [istart, i], count=[dim_t_in, 1])
        ncread!(joinpath(path, "forcings", "snowf_$(res).nc"), "snowf", snowf, start = [istart, i], count=[dim_t_in, 1])
        ncread!(joinpath(path, "forcings", "tair_$(res).nc"),  "tair",  tair,  start = [istart, i], count=[dim_t_in, 1])
        ncread!(joinpath(path, "forcings", "iswr_$(res).nc"),  "iswr",  iswr,  start = [istart, i], count=[dim_t_in, 1])
        ncread!(joinpath(path, "forcings", "ilwr_$(res).nc"),  "ilwr",  ilwr,  start = [istart, i], count=[dim_t_in, 1])
        ncread!(joinpath(path, "forcings", "pres_$(res).nc"),  "pres",  pres,  start = [istart, i], count=[dim_t_in, 1])
        ncread!(joinpath(path, "forcings", "rhum_$(res).nc"),  "rhum",  rhum,  start = [istart, i], count=[dim_t_in, 1])
        ncread!(joinpath(path, "forcings", "wind_$(res).nc"),  "wind",  wind,  start = [istart, i], count=[dim_t_in, 1])

        # Compute snowfall and rainfall from precipitation and temperature

        if comp_prec == 1

            prec = snowf + rainf

            fill!(snowf, 0.0)
            fill!(rainf, 0.0)

            snowf[tair .<= 0.0] .= prec[tair .<= 0.0]
            rainf[tair .> 0.0] .= prec[tair .> 0.0]

        end

        metdata = [year_in month_in day_in hour_in iswr ilwr snowf/dt rainf/dt tair+273.15 rhum wind pres*1000]

        # Initilize states

        md = FsmType(am, cm, dm, em, hm, dt)

        # Run model

        if pot_melt == 1
            results, desc = run_fsm_special(md, metdata)
        else
            results, desc = run_fsm(md, metdata)
        end

        # Aggregate to daily

        ts_3h.values .= results

        ts_24h_mean = collapse(ts_3h, day, first, mean)

        ts_24h_sum = collapse(ts_3h, day, first, sum)

        # Save to netcdf

        ncwrite(ts_24h_mean.values[:,1], joinpath(path, "results_$(iexp)", "hs_$(res).nc"), "hs", start=[1, i], count=[-1, 1])
        ncwrite(ts_24h_mean.values[:,2], joinpath(path, "results_$(iexp)", "swe_$(res).nc"), "swe", start=[1, i], count=[-1, 1])
        ncwrite(ts_24h_sum.values[:,3], joinpath(path, "results_$(iexp)", "rof_$(res).nc"), "rof", start=[1, i], count=[-1, 1])

    end

    ncclose()

end

function run_all_resultions(path, tstart, tstop, dt, df_settings)

    for row in eachrow(df_settings)

        for res in ["50km", "25km", "10km", "5km", "1km"]

            run_simulation(path,
                           tstart,
                           tstop,
                           Int64(row[:am]),
                           Int64(row[:cm]),
                           Int64(row[:dm]),
                           Int64(row[:em]),
                           Int64(row[:hm]),
                           dt,
                           res,
                           Int64(row[:experiment]),
                           Int64(row[:pot_melt]),
                           Int64(row[:comp_prec]))

            print("Finished for $(res)\n")

        end

    end

end

# Settings common for all simulations

path = "/data02/Ican/vic_sim/fsm_past_1km/netcdf"

dt = 3600*3

tstart = DateTime(2002, 9, 1, 0, 0, 0)

tstop  = DateTime(2003, 8, 31, 21, 0, 0)

# Settings specific for each run

file = joinpath(dirname(pathof(IcanProj)), "..", "data", "simulations.xlsx")

df_settings = readxlsheet(DataFrame, file, "simulations")

# Run simulations

run_all_resultions(path, tstart, tstop, dt, df_settings)


text:
module Inputs

using QuantumESPRESSOBase.Inputs:
    optionof,
    optionpool,
    groupname,
    required_namelists,
    optional_namelists,
    required_cards,
    optional_cards,
    allnamelists,
    allcards
using QuantumESPRESSOParser.Inputs: InvalidInput

export InvalidInput,
    optionof,
    optionpool,
    groupname,
    required_namelists,
    optional_namelists,
    required_cards,
    optional_cards,
    allnamelists,
    allcards

module PWscf

using QuantumESPRESSOBase.Inputs.PWscf:
    ControlNamelist,
    SystemNamelist,
    ElectronsNamelist,
    IonsNamelist,
    CellNamelist,
    DosNamelist,
    BandsNamelist,
    AtomicSpecies,
    AtomicSpeciesCard,
    AtomicPosition,
    AtomicPositionsCard,
    CellParametersCard,
    AtomicForce,
    AtomicForcesCard,
    KPointsCard,
    KMeshCard,
    GammaPointCard,
    SpecialPointsCard,
    PWInput,
    VerbositySetter,
    ElectronicTemperatureSetter,
    ElecTempSetter,
    VolumeSetter,
    PressureSetter,
    CellParametersCardSetter,
    AtomicPositionsCardSetter,
    optconvert,
    getxmldir,
    wfcfiles,
    getpseudodir,
    getpotentials,
    exitfile,
    mkexitfile,
    optionof,
    optionpool,
    allnamelists,
    allcards,
    required_namelists,
    optional_namelists,
    required_cards,
    optional_cards
using QuantumESPRESSOParser.Inputs.PWscf
using QuantumESPRESSOFormatter.Inputs.PWscf: format_file, format_text

export ControlNamelist,
    SystemNamelist,
    ElectronsNamelist,
    IonsNamelist,
    CellNamelist,
    DosNamelist,
    BandsNamelist,
    AtomicSpecies,
    AtomicSpeciesCard,
    AtomicPosition,
    AtomicPositionsCard,
    CellParametersCard,
    AtomicForce,
    AtomicForcesCard,
    KPointsCard,
    KMeshCard,
    GammaPointCard,
    SpecialPointsCard,
    PWInput,
    VerbositySetter,
    ElectronicTemperatureSetter,
    ElecTempSetter,
    VolumeSetter,
    PressureSetter,
    CellParametersCardSetter,
    AtomicPositionsCardSetter,
    optconvert,
    getxmldir,
    wfcfiles,
    getpseudodir,
    getpotentials,
    optionof,
    optionpool,
    allnamelists,
    allcards,
    required_namelists,
    optional_namelists,
    required_cards,
    optional_cards,
    format_file,
    format_text

end # module PWscf

module PHonon

using QuantumESPRESSOBase.Inputs.PHonon:
    PhNamelist,
    Q2rNamelist,
    MatdynNamelist,
    DynmatNamelist,
    Q2rInput,
    DynmatInput,
    PhInput,
    MatdynInput,
    VerbositySetter,
    relayinfo
using QuantumESPRESSOParser.Inputs.PHonon

export PhNamelist,
    Q2rNamelist,
    MatdynNamelist,
    DynmatNamelist,
    Q2rInput,
    DynmatInput,
    PhInput,
    MatdynInput,
    VerbositySetter,
    relayinfo

end # module PHonon

end # module Inputs


text:
using RuntimeGeneratedFunctions, BenchmarkTools
using Test

function f(_du,_u,_p,_t)
    @inbounds _du[1] = _u[1]
    @inbounds _du[2] = _u[2]
    nothing
end

ex1 = :((_du,_u,_p,_t) -> begin
    @inbounds _du[1] = _u[1]
    @inbounds _du[2] = _u[2]
    nothing
end)

ex2 = :(function f(_du,_u,_p,_t)
    @inbounds _du[1] = _u[1]
    @inbounds _du[2] = _u[2]
    nothing
end)

ex3 = :(function (_du::T,_u::Vector{E},_p::P,_t::Any) where {T<:Vector,E,P}
    @inbounds _du[1] = _u[1]
    @inbounds _du[2] = _u[2]
    nothing
end)

f1 = RuntimeGeneratedFunction(ex1)
f2 = RuntimeGeneratedFunction(ex2)
f3 = RuntimeGeneratedFunction(ex3)

du = rand(2)
u = rand(2)
p = nothing
t = nothing

@test f1(du,u,p,t) === nothing
du == u
du = rand(2)
f2(du,u,p,t)
@test du == u
du = rand(2)
@test f3(du,u,p,t) === nothing
du == u

t1 = @belapsed $f($du,$u,$p,$t)
t2 = @belapsed $f1($du,$u,$p,$t)
t3 = @belapsed $f2($du,$u,$p,$t)
t4 = @belapsed $f3($du,$u,$p,$t)

@test t1 ≈ t2 atol = 3e-9
@test t1 ≈ t3 atol = 3e-9
@test t1 ≈ t4 atol = 3e-9

function no_worldage()
    ex = :(function f(_du,_u,_p,_t)
        @inbounds _du[1] = _u[1]
        @inbounds _du[2] = _u[2]
        nothing
    end)
    f1 = RuntimeGeneratedFunction(ex)
    du = rand(2)
    u = rand(2)
    p = nothing
    t = nothing
    f1(du,u,p,t)
end
@test no_worldage() === nothing


text:
# Weyl points, comparing w/ http://dx.doi.org/10.1103/PhysRevB.85.155118
cases = [214 => "P", 213 => "R", 212 => "R", 199 => "P", 198 => "R", 98 => "P", 96 => "A",
          92 => "A", 24 => "W", 19 => "R"]
cases = [182 => "K", 181 => "K", 180 => "K", 179 => "K", 178 => "K", 177 => "K", 154 => "K", 152 => "K", 150 => "K"]
for (sgnum, klab) in cases
    #klab="H"
    lgirsd = lgirreps(sgnum, Val(3))
    lgirs = lgirsd[klab]
    println(sgnum, " (", klab, "):")
    for lgir in lgirs
        println("   ", Crystalline.formatirreplabel(label(lgir)), ": ", KdotP.isweyl(lgir))
    end
end

sgnum, klab = 199, "PA"
lgirsd = lgirreps(sgnum, Val(3))
lgirs = lgirsd[klab]
lgir = lgirs[3]
kdotp(lgir)

for sgnum in 1:230
    timereversal=true
    lgirsd = lgirreps(sgnum, Val(3))
    printed_sg = false
    for (klab, lgirs) in lgirsd
        printed_k = false
        timereversal && (lgirs = realify(lgirs))
        for lgir in lgirs
            #println(num(lgir), ": ", label(lgir))
            if KdotP.isweyl(lgir; timereversal)
                printed_sg || (print("SG ", sgnum, ":\t"); printed_sg=true)
                if !printed_k
                    print(" | ")
                    printed_k=true
                else
                    print(", ")
                end
                print(Crystalline.formatirreplabel(label(lgir)))
            end
        end

    end
    printed_sg && println()
end

# Current discrepancise/"issues" (with time-reversal):
#   - SG 19, R₁: pseudo-real 4D irrep; bugs out on kdotp
#   - SG 80, P₁P₂: pseudo-real 2D irrep; not in Manes' tables.
#   - SG 92 & 96, A₁ & A₂: complex 4D irreps - kdotp bugs out on U
#   - SG 198, R₁₂₃: one complex and one pseudoreal 4D irreps - dunno about kdotp
#   - SG 212 & 213, R₁R₂: 4D complex irrep - dunno about kdotp
#   - SG 212 & 213, R₃: 4D irrep, not checked by `isweyl` cf. 4D nature
# All of Table 2 is good.

#----
# high-order degeneracies
sgnums_and_klabs = [(218, "R"), (220, "H"), (222, "R"), (223, "R"), (230, "H")]
for (sgnum, klab) in sgnums_and_klabs
    lgirsd = lgirreps(sgnum, Val(3))
    lgirs = realify(lgirsd[klab])
    lgir = lgirs[end]

    printstyled("─"^75, "\n\n", color=:light_black)
    println("━━━ SG ", sgnum, " ━━━\n")
    display(kdotp(lgir))
end


text:
using Performance
using Test

@testset "type abuse" begin
    ##
    ## Sometimes when people start getting excited about multiple dispatch, they start trying to
    ## to solve every problem as a dispatch problem. This works well in some circumstances and badly
    ## in others. This problem will illustrate that dispatch is not "magic": when types can't
    ## be predicted in advance, you're better off doing more at runtime and less via the type system.
    ##

    numbers = [15, 3, 6, 11, 4]
    ops = Operator[Plus(), Divide(), Times(), Minus()]
    opsstr = "+/*-"
    @test calculate_typed(numbers, ops) == 29
    @test calculate_branch(numbers, opsstr) == 29
    btyped = @benchmark(calculate_typed($numbers, $ops), samples=100, evals=10)
    bbranch = @benchmark(calculate_branch($numbers, $opsstr), samples=100, evals=10)
    println("type_abuse.jl: actual performance ratio: ", minimum(btyped).time / minimum(bbranch).time)

    # The other thing about excessive use of the type system is that it can increase compile-time latency,
    # since code has to be compiled for each different type. In some cases, it has been observed that one
    # can go from tens of minutes to a few seconds of compile time:
    # https://discourse.julialang.org/t/understanding-and-optimizing-compiler-time-just-a-bit/51759/13?u=tim.holy
    # Bottom line: types are great when they are predictable, but they can get in your way when not.
    # Ask yourself, "will the compiler be able to predict the next type in this computation?"
    # If the answer is "no," consider using branches or more conventional programming logic.
end


text:
normalise(d::Truncated{Normal, Continuous}) = d.untruncated.σ*(d.ucdf-d.lcdf)

lognormalise(d::Truncated{Normal, Continuous}) = log(normalise(d))


text:
# ---
# title: 1249. Minimum Remove to Make Valid Parentheses
# id: problem1249
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String, Stack
# link: <https://leetcode.com/problems/minimum-remove-to-make-valid-parentheses/description/>
# hidden: true
# ---
#
# Given a string s of `'('` , `')'` and lowercase English characters.
#
# Your task is to remove the minimum number of parentheses ( `'('` or `')'`, in
# any positions ) so that the resulting _parentheses string_ is valid and return
# **any** valid string.
#
# Formally, a _parentheses string_ is valid if and only if:
#
#   * It is the empty string, contains only lowercase characters, or
#   * It can be written as `AB` (`A` concatenated with `B`), where `A` and `B` are valid strings, or
#   * It can be written as `(A)`, where `A` is a valid string.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "lee(t(c)o)de)"
#     Output: "lee(t(c)o)de"
#     Explanation: "lee(t(co)de)" , "lee(t(c)ode)" would also be accepted.
#
#
# **Example 2:**
#
#
#
#     Input: s = "a)b(c)d"
#     Output: "ab(c)d"
#
#
# **Example 3:**
#
#
#
#     Input: s = "))(("
#     Output: ""
#     Explanation: An empty string is also valid.
#
#
# **Example 4:**
#
#
#
#     Input: s = "(a(b(c)d)"
#     Output: "a(b(c)d)"
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 10^5`
#   * `s[i]` is one of  `'('` , `')'` and lowercase English letters`.`
#
#
## @lc code=start
using LeetCode

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


text:
using Statistics
using LinearAlgebra
using ProgressMeter
using FileIO
using ArgParse

mutable struct SortedDat{F,I}
    num_categories::I
    num_bits::I
    features::Array{F,2}
    labels::Array{I,1}
    category_ranges::Array{UnitRange{I}}
end

mutable struct DecisionTreeModel{F,I}
    num_categories::I
    num_bits::I
    category_to_leaf::Array{I,1}
    leaf_to_category::Array{I,1}
    weights::Array{F,2}
    biases::Array{F,1}
    avg_lls::Array{F,1}
end

function SortedDat(dir_path; dimensions = nothing, subset = "train")
    if dimensions === nothing
        dimensions = "full"
    else
        dimensions = "k$dimensions"
    end

    features = load_features(joinpath(dir_path, "$subset-features-$dimensions.np"))
    labels, num_categories = load_labels(joinpath(dir_path, "$subset-labels-first.np"))
    @assert size(features, 2) == size(labels, 1)

    num_bits = Int32(ceil(log2(num_categories)))

    # Sort data set by label.
    permutation = sortperm(labels)
    labels = labels[permutation]
    features = features[:, permutation]

    # Remember where data points for each label begin and end.
    category_ranges = Array{UnitRange{Int32}}(undef, 2^num_bits)
    data_i = 1
    for cat = 1:2^num_bits
        start = data_i
        while data_i <= length(labels) && labels[data_i] == cat
            data_i += 1
        end
        category_ranges[cat] = start:(data_i - 1)
    end

    SortedDat{Float32,Int32}(num_categories, num_bits, features, labels, category_ranges)
end

function load_features(path)
    println("Loading features from \"$path\" ...")
    file = open(path, "r")
    n = read(file, Int32)
    k = read(file, Int32)
    println("=> dimensions: $k x $n")
    features = reinterpret(Float32, read(file, 4 * k * n))
    @assert eof(file)
    close(file)

    reshape(features, Int64(k), Int64(n))
end

function load_labels(path)
    println("Loading labels from \"$path\" ...")
    file = open(path, "r")
    num_dat = read(file, Int32)
    num_categories = read(file, Int32)
    println("=> $num_dat data points over $num_categories categories")
    labels = reinterpret(Int32, read(file, 4 * num_dat)) .+ 1
    @assert eof(file)
    @assert length(labels) == num_dat
    close(file)
    labels, num_categories
end

function leaves_for_node(node, num_bits)
    level = Int32(floor(log2(node)))
    start = (node << (num_bits - level)) - (1 << num_bits) + 1
    stop = start + (1 << (num_bits - level)) - 1
    start:stop
end

softplus(x) = log(1 + exp(x))

function expand_dims(a, i)
    s = size(a)
    reshape(a, s[1:i - 1]..., 1, s[i:end]...)
end

struct Yes end
struct No end

function grad_log_joint(weight_and_bias, features, label_index, is_left, regularizer; hess = No())
    k = length(weight_and_bias) - 1
    scores = features' * view(weight_and_bias, 1:k) .+ weight_and_bias[k + 1]
    for i = 1:length(scores)
        if is_left[label_index[i]]
            scores[i] = - scores[i]
        end
    end

    signed_sigmoid_neg_scores = 1 ./ (1 .+ exp.(scores))
    signed_sigmoid_pos_scores = if hess isa Yes
        1 .- signed_sigmoid_neg_scores
    else
        nothing
    end

    for i = 1:length(signed_sigmoid_neg_scores)
        if is_left[label_index[i]]
            signed_sigmoid_neg_scores[i] = - signed_sigmoid_neg_scores[i]
            if hess isa Yes
                signed_sigmoid_pos_scores[i] = - signed_sigmoid_pos_scores[i]
            end
        end
    end

    features_signed_sigmoid_neg_scores = features .* expand_dims(signed_sigmoid_neg_scores, 1)

    grad = -regularizer * weight_and_bias  # Gradient of log prior.
    grad[1:k] += reshape(sum(features_signed_sigmoid_neg_scores, dims = 2), k)
    grad[k + 1] += sum(signed_sigmoid_neg_scores)

    if hess isa Yes
        features_signed_sigmoid_pos_scores = features .* expand_dims(signed_sigmoid_pos_scores, 1)
        hess_mat = Matrix{Float32}(regularizer * LinearAlgebra.I, k + 1, k + 1)
        hess_mat[1:k, 1:k] += sum(expand_dims(features_signed_sigmoid_pos_scores, 1) .* expand_dims(features_signed_sigmoid_neg_scores, 2),
            dims = 3)
        hess_mat[end, 1:k] = sum(features_signed_sigmoid_pos_scores .* expand_dims(signed_sigmoid_neg_scores, 1), dims = 2)
        hess_mat[1:k, end] = view(hess_mat, k + 1, 1:k)'
        hess_mat[end, end] += sum(signed_sigmoid_pos_scores .* signed_sigmoid_neg_scores)
        grad, hess_mat
    else
        grad
    end
end

function split_in_half(dat, categories, cat_emb, max_half, regularizer)
    k = size(cat_emb, 1)

    if isempty(categories)
        return zeros(Float32, k), Float32(0), Int32[], Int32[]
    elseif length(categories) == 1
        # Trivial classifier where all except one leaf are padding. Put the one relevant leaf into right subtree and
        # return a zero weight vector and a very large bias so that the left subtree has almost zero probability.
        return zeros(Float32, k), Float32(1000), Int32[], Int32[1]
    end

    max_half = min(max_half, length(categories))
    min_half = length(categories) - max_half
    @assert 0 <= min_half <= max_half

    # Make initial split at center.
    left_count = (min_half + max_half) ÷ 2

    # Initialize weight vector with PCA of category embeddings
    cur_cat_emb = cat_emb[:, categories]
    centered_cat_emb = cur_cat_emb .- Statistics.mean(cur_cat_emb, dims = 2)
    cov = centered_cat_emb * centered_cat_emb'  # shape (emb_dim, emb_dim)
    centered_cat_emb = nothing

    weight = reshape(LinearAlgebra.eigen(LinearAlgebra.Symmetric(cov), k:k).vectors, size(cov, 1))

    # Sort scores and initialize partition.
    # We use the relation `log(sigmoid(-x)) = log(sigmoid(x)) - x`.
    cat_scores = cur_cat_emb' * weight
    permutation = sortperm(cat_scores)
    is_left = BitArray{1}(undef, length(categories))
    is_left[view(permutation, 1:left_count)] .= true
    is_left[view(permutation, left_count + 1:length(permutation))] .= false

    # Extract relevant points from data set.
    num_data_points = Int32(0)
    for cat = categories
        num_data_points += length(dat.category_ranges[cat])
    end
    features = Array{Float32}(undef, size(cat_emb, 1), num_data_points)
    label_index = Array{Int32}(undef, num_data_points)  # Index into `is_left`.
    i = 1
    for (idx, cat) = enumerate(categories)
        r = dat.category_ranges[cat]
        features[:,i:i + length(r) - 1] = dat.features[:,r]
        label_index[i:i + length(r) - 1] .= idx
        i += length(r)
    end
    @assert i == num_data_points + 1

    # Initialize bias such that scores have zero mean.
    # A roughly correct initial bias is important for Newton's descent to work because
    # the objective function has a very low curvature far away from the optimum.
    bias = -Statistics.mean(features' * weight)

    weight_and_bias = vcat(weight, Float32[bias])
    steps_since_change = 0
    steps_until_hess = 0
    hess = Array{Float32}(undef, 1, 1)
    total_step = 0

    while steps_since_change < 10
        # Continuous optimization: optimize until convergence but break if change is large enough to warrant discrete update
        path_length = Float32(0)
        while steps_since_change < 10 && (total_step < 10 || path_length < 0.01)
            if steps_until_hess == 0
                grad, hess = grad_log_joint(weight_and_bias, features, label_index, is_left, regularizer, hess = Yes())
                steps_until_hess = 8
            else
                grad = grad_log_joint(weight_and_bias, features, label_index, is_left, regularizer)
                steps_until_hess -= 1
            end

            step = hess \ grad
            step_norm = sqrt(sum(step.^2))
            path_length += step_norm
            if step_norm > 1.0
                steps_until_hess = 0
                step *= 1 / step_norm
            end
            weight_and_bias += step

            if step_norm > 1e-10
                steps_since_change = 0
            else
                steps_since_change += 1
            end
            total_step += 1
        end

        # Discrete optimization.
        # We use the relation `log(sigmoid(-x)) = log(sigmoid(x)) - x`.
        cat_scores = cur_cat_emb' * weight
        permutation = sortperm(cat_scores)
        sorted_scores = cat_scores[permutation]
        first_found_pos = findfirst(x->x >= 0, sorted_scores)
        last_negative = if first_found_pos === nothing
            length(sorted_scores)
        else
            first_found_pos - 1
        end
        left_count = max(min_half, min(max_half, last_negative))

        if (!all(is_left[view(permutation, 1:left_count)]) || any(is_left[view(permutation, left_count + 1:length(permutation))]))
            is_left[view(permutation, 1:left_count)] .= true
            is_left[view(permutation, left_count + 1:length(permutation))] .= false
            steps_since_change = 0
            steps_until_hess = 0
            if left_count == 0 || left_count == length(permutation)
                # Classifier became trivial as one child tree only contains paddings. Once we make this decision we go all in.
                return zeros(Float32, k), Float32(1000), Int32[], collect(Int32, 1:length(permutation))
            end
        end
    end

    view(weight_and_bias, 1:k), weight_and_bias[end], view(permutation, 1:left_count), view(permutation, left_count + 1:length(permutation))
end

function optimize_tree(dat, regularizer, scale_regularizer::Bool; normalize = false)
    cat_emb = zeros(Float64, size(dat.features, 1), 2^dat.num_bits)
    if normalize
        for (category, data_range) = enumerate(dat.category_ranges)
            if length(data_range) != 0
                cat_emb[:, category] = view(Statistics.mean(Array{Float64}(view(dat.features, :, data_range)), dims = 2), :, 1)
            end
        end
        cat_emb[:,1:dat.num_categories] .-= Statistics.mean(cat_emb[:,1:dat.num_categories], dims = 2)
    else
        for (category, data_range) = enumerate(dat.category_ranges)
            cat_emb[:,category] = view(sum(Array{Float64}(view(dat.features, :, data_range)), dims = 2), :, 1)
            # No need to center since feature matrix is already centered.
        end
    end

    @assert isapprox(Statistics.mean(cat_emb, dims = 2), zeros(size(dat.features, 1), 1), atol = 1e-8)

    leaf_to_category = Array(Int32(1):Int32(2^dat.num_bits))
    k = size(dat.features, 1)
    node_weights = Array{Float32}(undef, k, 2^dat.num_bits - 1)
    node_biases = Array{Float32}(undef, 2^dat.num_bits - 1)

    @showprogress for layer = 1:dat.num_bits
        for node = 2^(layer - 1):(2^layer - 1)
            cur_leaves = leaves_for_node(node, dat.num_bits)
            cur_categories_with_padding = leaf_to_category[cur_leaves]
            cur_categories = filter(x->x <= dat.num_categories, cur_categories_with_padding)
            cur_padding = filter(x->x > dat.num_categories, cur_categories_with_padding)
            reg = if scale_regularizer
                regularizer * length(cur_categories)
            else
                regularizer
            end
            node_weights[:, node], node_biases[node], indices_left, indices_right = split_in_half(
                dat, cur_categories, cat_emb, length(cur_leaves) ÷ 2, reg)
            leaf_to_category[cur_leaves] = vcat(cur_categories[indices_left], cur_padding, cur_categories[indices_right])
        end
    end

    # Check that `leaf_to_category` is a permutation.
    counts = zeros(Int32, (2^dat.num_bits,))
    for i = leaf_to_category
        counts[i] += 1
    end
    @assert all(counts .== 1)

    category_to_leaf = Array{Int32}(undef, dat.num_categories)
    for (leaf, category) = enumerate(leaf_to_category)
        if category <= dat.num_categories
            category_to_leaf[category] = leaf
        end
    end

    leaf_to_category, category_to_leaf, node_weights, node_biases
end

function average_category_log_likelihoods(dat::SortedDat, category_to_leaf, weights, biases)
    avg_lls = Array{Float32}(undef, dat.num_categories)
    bits = Int32[1<<i for i = (dat.num_bits - 1):-1:0]
    for (cat, (leaf, r)) in enumerate(zip(category_to_leaf, dat.category_ranges))
        decisions = ((leaf - 1) .& bits) .!= 0
        nodes = Int32[(leaf - 1 + 2^dat.num_bits) >> i for i = dat.num_bits:-1:1]
        cur_weights = weights[:, nodes]
        cur_biases = biases[nodes]
        for i in 1:dat.num_bits
            if decisions[i]
                cur_weights[:, i] .*= -1
                cur_biases[i] *= -1
            end
        end
        scores = view(dat.features, :, r)' * cur_weights .+ expand_dims(cur_biases, 1)  # shape (num_dat, num_bits)
        avg_lls[cat] = - sum(log.(1 .+ exp.(scores))) / length(r)
    end

    avg_lls
end

arg_parser = ArgParseSettings()
@add_arg_table arg_parser begin
    "input_dir"
        help = "Path to indput directory."
    "output_file"
        help = "Path to output file."
    "--reg"
        help = "Regularizer strength."
        arg_type = Float32
        default = Float32(1.0)
    "--frequency_reg"
        help = "Multiply regularizer strength with number of data points under node."
        action = :store_true
    "-d", "--embedding_dim"
        help = "Embedding dimension."
        arg_type = Int
        default = 16
end

args = parse_args(ARGS, arg_parser, as_symbols = true)

dat = SortedDat(args[:input_dir], dimensions = args[:embedding_dim])
println("Starting fit ...")
(leaf_to_category, category_to_leaf, weights, biases), elapsed, _ = @timed optimize_tree(dat, args[:reg], args[:frequency_reg])
println("done ($elapsed seconds).")

avg_lls = average_category_log_likelihoods(dat, category_to_leaf, weights, biases)
println("Average log likelihood per category on training set: $(Statistics.mean(avg_lls))")

save(args[:output_file], "model",
DecisionTreeModel(dat.num_categories, dat.num_bits, category_to_leaf, leaf_to_category, weights, biases, avg_lls))

max_stars_repo_name: mandt-lab/adversarial-negative-sampling

max_issues_repo_name: mandt-lab/adversarial-negative-sampling

max_forks_repo_name: mandt-lab/adversarial-negative-sampling


text:
using MultivariateFunctions
tol = 10*eps()

# pe functions
function test_pe(a::Float64,b::Float64,base::Float64,d::Int,x::Float64)
    return a * exp(b*(x-base)) * (x-base)^d
end
f1 = PE_Function(1.0, 2.0,3.0, 4)
f1_test_result = test_pe(1.0, 2.0,3.0, 4,5.0)
abs(evaluate(f1, 5.0) - f1_test_result) < tol
f2 = PE_Function(0.0, 2.0,3.0, 4)
f2_test_result = 0.0
abs(evaluate(f2, 5.0) - f2_test_result) < tol
f3 = PE_Function(1.0, 8.0,6.0, 8)
f3_test_result = test_pe(1.0, 8.0,6.0, 8,5.0)
abs(evaluate(f3, 5.0) - f3_test_result) < tol

# Sum of functions
sum0 =  Sum_Of_Functions([])
typeof(sum0) == MultivariateFunctions.Sum_Of_Functions
sum1 =  Sum_Of_Functions([f1])
typeof(sum1) == MultivariateFunctions.Sum_Of_Functions
sum3 =  Sum_Of_Functions([f2,f3])
typeof(sum3) == MultivariateFunctions.Sum_Of_Functions
sum4 =  Sum_Of_Functions([f1,f2,f3])
typeof(sum4) == MultivariateFunctions.Sum_Of_Functions
length(sum4.functions_) == 2
sum5 =  Sum_Of_Functions([sum3, sum4])
typeof(sum5) == MultivariateFunctions.Sum_Of_Functions
length(sum5.functions_) == 2
abs(evaluate(sum1, 5.0) - f1_test_result) < tol
abs(evaluate(sum5, 5.0) - 2*f3_test_result - f1_test_result) < 100* tol
abs(evaluate(sum0, 5.0)) < tol

coordinates = Dict(Symbol.(["x", "y", "z", "w"]) .=> [0.5, 2.0, 3.5, 4.0])

f = PE_Function(2.4, Dict{Symbol,PE_Unit}([:x,:y,:z,:w] .=> [ PE_Unit(2.0,3.0,2), PE_Unit(2.0,3.0,2) , PE_Unit(2.0,3.0,2), PE_Unit(2.0,3.0,2)   ]))
new_bases = Dict{Symbol,Float64}([:x,:y,:z, :q] .=> [ 2.8,3.0,1.0,3.0  ])
converted_f = change_base(f,new_bases)
abs(sum(evaluate.(converted_f, Ref(coordinates))) - evaluate(f, coordinates)) < tol

fl = 8.9
integ = 7

function test_result(func, expected_type, eval_to, len = 1)
    val_test = abs(evaluate(func, 5.0) - eval_to) < 1e-09
    if (!val_test)
        print("Failed Val Test")
    end
    type_test = typeof(func) == expected_type
    if (!type_test)
        print("Failed Type Test")
    end
    if typeof(func) == MultivariateFunctions.Sum_Of_Functions
        len_test = length(func.functions_) == len
        if (!len_test)
            print("Failed Length Test")
        end
        return all([val_test,type_test,len_test])
    else
        return all([val_test,type_test])
    end
end

# pe and Float
test_result(f1 + fl, MultivariateFunctions.Sum_Of_Functions, f1_test_result + fl,2 )
test_result(f1 - fl, MultivariateFunctions.Sum_Of_Functions, f1_test_result - fl,2 )
test_result(f1 * fl, MultivariateFunctions.PE_Function, (f1_test_result) * fl )
test_result(f1 / fl, MultivariateFunctions.PE_Function, (f1_test_result) / fl )
test_result(fl + f1, MultivariateFunctions.Sum_Of_Functions, f1_test_result + fl,2 )
test_result(fl - f1, MultivariateFunctions.Sum_Of_Functions, fl - f1_test_result,2 )
test_result(fl * f1, MultivariateFunctions.PE_Function, (f1_test_result) * fl )
# pe and Int
test_result(f1 + integ, MultivariateFunctions.Sum_Of_Functions, f1_test_result + integ,2 )
test_result(f1 - integ, MultivariateFunctions.Sum_Of_Functions, f1_test_result - integ,2 )
test_result(f1 * integ, MultivariateFunctions.PE_Function, (f1_test_result) * integ )
test_result(f1 / integ, MultivariateFunctions.PE_Function, (f1_test_result) / integ )
test_result(integ + f1, MultivariateFunctions.Sum_Of_Functions, f1_test_result + integ,2 )
test_result(integ - f1, MultivariateFunctions.Sum_Of_Functions, -f1_test_result + integ,2 )
test_result(integ * f1, MultivariateFunctions.PE_Function, (f1_test_result) * integ )

# Sum and Float
test_result(sum4 + fl, MultivariateFunctions.Sum_Of_Functions, f1_test_result + f3_test_result + fl, 3 )
test_result(sum4 - fl, MultivariateFunctions.Sum_Of_Functions, f1_test_result + f3_test_result - fl, 3 )
test_result(sum4 * fl, MultivariateFunctions.Sum_Of_Functions, (f1_test_result + f3_test_result) * fl, 2 )
test_result(sum4 / fl, MultivariateFunctions.Sum_Of_Functions, (f1_test_result + f3_test_result) / fl, 2 )
test_result(fl + sum4, MultivariateFunctions.Sum_Of_Functions, f1_test_result + f3_test_result + fl, 3 )
test_result(fl - sum4, MultivariateFunctions.Sum_Of_Functions, -f1_test_result -f3_test_result + fl, 3 )
test_result(fl * sum4, MultivariateFunctions.Sum_Of_Functions, (f1_test_result + f3_test_result) * fl, 2 )
# Sum and Int
test_result(sum4 + integ, MultivariateFunctions.Sum_Of_Functions, f1_test_result + f3_test_result + integ, 3 )
test_result(sum4 - integ, MultivariateFunctions.Sum_Of_Functions, f1_test_result + f3_test_result - integ, 3 )
test_result(sum4 * integ, MultivariateFunctions.Sum_Of_Functions, (f1_test_result + f3_test_result) * integ, 2 )
test_result(sum4 / integ, MultivariateFunctions.Sum_Of_Functions, (f1_test_result + f3_test_result) / integ, 2 )
test_result(integ + sum4, MultivariateFunctions.Sum_Of_Functions, f1_test_result + f3_test_result + integ, 3 )
test_result(integ - sum4, MultivariateFunctions.Sum_Of_Functions, -f1_test_result -f3_test_result + integ, 3 )
test_result(integ * sum4, MultivariateFunctions.Sum_Of_Functions, (f1_test_result + f3_test_result) * integ, 2 )

### Making sums with addition and subtraction.
evaluate(sum4, 5.0) == evaluate(f1 + f2 + f3, 5.0)
abs(evaluate(sum4, 5.0) - 2*evaluate(f3, 5.0) - evaluate(f1 + f2 - f3, 5.0)) < 0.001

### Changing of base
f1_unchanged = change_base(f1, Dict{Symbol,Float64}(:default => 3.0))
abs(f1.units_[:default].base_ - f1_unchanged[1].units_[:default].base_) < tol

f1_changed = change_base(f1, Dict{Symbol,Float64}(:default => 4.0))
abs(f1_changed[1].units_[:default].base_ - 4.0) < tol
abs(f1_changed[3].units_[:default].base_ - 4.0) < tol
abs(sum(evaluate.(f1_changed, 5.0)) - f1_test_result) < 100*tol # Changing bases should not change this.

f1_changed_again = change_base(f1, Dict{Symbol,Float64}(:default => -1.0))
abs(f1_changed_again[1].units_[:default].base_ + 1.0) < tol
abs(f1_changed_again[3].units_[:default].base_ + 1.0) < tol
abs(sum(evaluate.(f1_changed_again, 5.0)) - f1_test_result) < 1000000*tol # Changing bases should not change this.

f3_changed = change_base(f3, Dict{Symbol,Float64}(:default => 11.0))
abs(f3_changed[1].units_[:default].base_ - 11.0) < tol
abs(f3_changed[3].units_[:default].base_ - 11.0) < tol
abs(sum(evaluate.(f3_changed, 5.0)) - f3_test_result) < 1000000*tol # Changing bases should not change this.

### multiplication of functions
test_result( f1 * f3 , MultivariateFunctions.Sum_Of_Functions, (f1_test_result * f3_test_result), 9 )
test_result( f1 * PE_Function(1.0,2.0,3.0,4) , MultivariateFunctions.PE_Function, (f1_test_result * f1_test_result), 1 )
test_result( f1 * sum4 , MultivariateFunctions.Sum_Of_Functions, f1_test_result * (f1_test_result + f3_test_result), 10 )
test_result( sum4 * f1 , MultivariateFunctions.Sum_Of_Functions, f1_test_result * (f1_test_result + f3_test_result), 10 )
test_result( sum4 * sum4 , MultivariateFunctions.Sum_Of_Functions, (f1_test_result + f3_test_result)^2, 11 )

### Powers pe
test_result( f1 ^ 0 , MultivariateFunctions.PE_Function, 1.0 )
test_result( f1 ^ 1 , MultivariateFunctions.PE_Function, (f1_test_result) )
test_result( f1 ^ 2 , MultivariateFunctions.PE_Function, (f1_test_result * f1_test_result) )
abs(evaluate(f1 ^ 4 ,5.0) - (f1_test_result * f1_test_result * f1_test_result * f1_test_result)) < 1e-03 # This fails normal tol.

### Powers sums
test_result( sum4 ^ 0 , MultivariateFunctions.PE_Function, 1.0 )
test_result( sum4 ^ 1 , MultivariateFunctions.Sum_Of_Functions, (f1_test_result + f3_test_result), 2 )
test_result( sum4 ^ 2 , MultivariateFunctions.Sum_Of_Functions, (f1_test_result + f3_test_result)^2, 11 )
# Higher powers give an underflow problem with base changes.

## Testing Calculus.
pe_const    = PE_Function(2.0,0.0,5.0,0)
pe_lin      = PE_Function(2.0,0.0,5.0,1)
pe_quad     = PE_Function(2.0,0.0,2.0,2)
pe_exp      = PE_Function(2.0,2.0,1.0,0)
pe_exp_quad = PE_Function(7.0,2.0,2.0,2)

# Linear gradient constant
abs(evaluate(derivative(pe_lin),5.0) - evaluate(derivative(pe_lin),1.0) ) < tol
typeof(derivative(pe_lin)) == MultivariateFunctions.Sum_Of_Functions
# This is also linear
abs(evaluate(derivative(derivative(pe_quad)),5.0) -evaluate(derivative(derivative(pe_quad)),9.0) ) < tol
# derivative into a sum of functions
typeof(derivative(pe_exp_quad)) == MultivariateFunctions.Sum_Of_Functions
abs(evaluate(derivative(pe_exp_quad),5.0) - ( 7.0*exp(2.0*(5.0-2.0))*(5.0-2.0)*(2 + 2.0*(5.0-2.0)) ) ) < tol
# Derivaitve of sum
abs(evaluate(derivative(sum4),5.0) - evaluate(derivative(f1),5.0) - evaluate(derivative(f3),5.0) ) < 1e-09

# Integral of quadratic
abs(integral(pe_quad,1.0,2.0) - ((2/3)*2^3 - 4*2^2 + 8*2) + ((2/3)*1^3 - 4*1^2 + 8*1)) < tol
# Integral of exponential
abs(integral(pe_exp,0.2,3.8) - (exp(2.0*(3.8-1.0)) - exp(2.0*(0.2-1.0))  )) < tol
# Integral of combined.
combined_analytical_integral = PE_Function(3.5,2.0,2.0,2) - PE_Function(14.0/4.0,2.0,2.0,1) + PE_Function(14.0/8.0,2.0,2.0,0)
abs(integral(pe_exp_quad,0.2,3.8) - ( evaluate(combined_analytical_integral, 3.8) - evaluate(combined_analytical_integral, 0.2))) < 1000*tol

# left and right integrals
result = integral(pe_exp_quad,0.2,3.8)
left_integral = integral(pe_exp_quad, Dict{Symbol,Tuple{Float64,Symbol}}(:default => (0.2, :default_right)))
abs(evaluate(left_integral, Dict{Symbol,Float64}(:default_right => 3.8)) - result) < 100 * tol
right_integral = integral(pe_exp_quad, Dict{Symbol,Tuple{Symbol,Float64}}(:default => (:default_left, 3.8)))
abs(evaluate(right_integral, Dict{Symbol,Float64}(:default_left => 0.2)) - result) <  100 * tol
both_integral = integral(pe_exp_quad, Dict{Symbol,Tuple{Symbol,Symbol}}(:default => (:default_left, :default_right)))
abs(evaluate(both_integral, Dict{Symbol,Float64}([:default_left, :default_right] .=> [0.2, 3.8])) - result) <  100 * tol

## Testing linear rescaling
# For PE_Functions
test_func = PE_Function(0.5,1.2,0.9,2)
inp = 4.0
alpha = 0.435
beta = -1.52
rescaled_input = alpha*inp + beta
converted_test_func = convert_to_linearly_rescale_inputs(test_func, alpha, beta)
abs(evaluate(test_func, inp) - evaluate(converted_test_func, rescaled_input)) < 1e-10
# SumOfFunctions
sumFunc = test_func + pe_exp_quad + pe_exp
converted_sumFunc = convert_to_linearly_rescale_inputs(sumFunc, alpha, beta)
abs(evaluate(sumFunc, inp) - evaluate(converted_sumFunc, rescaled_input)) < 1e-10

# Multiplications with no units.
n1 = PE_Function(5.0)
n2 = PE_Function(12.0)
abs((n1 * n2).multiplier_ - 60.0) < tol
abs((n1 * pe_lin).multiplier_ - 10.0) < tol


text:
mutable struct SVAEvamp <: SVAE
	q
	g
	zdim
	hue
	μzfromhidden
	κzfromhidden
	pseudo_inputs

	"""
	SVAE_vamp(q, g, hdim, zdim, T) Constructor of the S-VAE where `zdim > 3` and T determines the floating point type (default Float32)
	"""
	SVAEvamp(q, g, hdim::Integer, zdim::Integer, num_pseudoinputs::Integer, T = Float32) = new(q, g, zdim, convert(T, huentropy(zdim)), Adapt.adapt(T, Chain(Dense(hdim, zdim), x -> normalizecolumns(x))), Adapt.adapt(T, Dense(hdim, 1, x -> σ.(x) .* 100)), Flux.param(Adapt.adapt(T, normalizecolumns(randn(size(q[1].W, 2), num_pseudoinputs)))))
end

function SVAEvamp(inputDim, hiddenDim, latentDim, numLayers, nonlinearity, layerType, num_pseudoinputs, T = Float32)
	encoder = Adapt.adapt(T, layerbuilder(inputDim, hiddenDim, hiddenDim, numLayers, nonlinearity, "", layerType))
    decoder = Adapt.adapt(T, layerbuilder(latentDim, hiddenDim, inputDim, numLayers + 1, nonlinearity, "linear", layerType))
	return SVAEvamp(encoder, decoder, hiddenDim, latentDim, num_pseudoinputs, T)
end

Flux.@treelike(SVAEvamp)

function sampleVamp(m::SVAEvamp, n::Integer)
	components = rand(1:size(m.pseudo_inputs, 2), n)
	(μz, κz) = zparams(m, m.pseudo_inputs[:, components])
	return samplez(m, μz, κz)
end

function wloss(m::SVAEvamp, x, β, d)
	(μz, κz) = zparams(m, x)
	z = samplez(m, μz, κz)
	zp = sampleVamp(m, size(z, 2))
	#prior = samplez(m, ones(size(μz)) .* normalizecolumns(m.priorμ), ones(size(κz)) .* m.priorκ)
	Ω = d(z, zp)
	xgivenz = m.g(z)
	return Flux.mse(x, xgivenz) + β * Ω
end

function pz(m::SVAEvamp, x)
	if size(x, 2) > 1
		xs = [x[:, i] for i in 1:size(x, 2)]
		return map(x -> pz_singleinstance(m, x), xs)
	else
		return pz_singleinstance(m, x)
	end
end

function pz_singleinstance(m::SVAEvamp, x)
	@assert size(x, 2) == 1
	μz = zparams(m, x)[1].data
	pseudoz = zparams(m, m.pseudo_inputs)[1].data
	p = 0
	for i in 1:size(m.pseudo_inputs, 2)
		p += exp.(log_vmf_wo_c(μz, pseudoz[:, i], 1))
	end
	p /= size(m.pseudo_inputs, 2)
end


text:
##############################################################
# Mandlebrot code
##############################################################

function calc(x, y, iterations)
    c0 = complex(x, y)
    c = 0
    for i = 1:iterations
        if abs(c) > 2
            return i
        end
        c = c * c + c0
    end
    return 0
end

function mandlebrot_fractal(min_x, min_y, max_x, max_y, canvas_width, canvas_height, iterations)
    pixels = fill(0, canvas_width, canvas_height)

    x_pixel_diff = (max_x - min_x) / canvas_width
    y_pixel_diff = (max_y - min_y) / canvas_height

    for x = 1:canvas_width
        for y = 1:canvas_height
            i = calc((x * x_pixel_diff) + min_x, (y * y_pixel_diff) + min_y, iterations)
            pixels[x,y] = (i << 21) + (i << 10) + i * 8
        end
    end

    return pixels
end

##############################################################
# Julia-set code
##############################################################

function julia_fractal(min_x, min_y, max_x, max_y, canvas_width, canvas_height, iterations)
    pixels = fill(0, canvas_width, canvas_height)

    x_pixel_diff = (max_x - min_x) / canvas_width
    y_pixel_diff = (max_y - min_y) / canvas_height

    cX, cY = -0.7, 0.27015

    for x = 1:canvas_width
        for y = 1:canvas_height
            zx = (x * x_pixel_diff) + min_x
            zy = (y * y_pixel_diff) + min_y
            i = iterations
            while ((((zx * zx) + (zy * zy)) < 4) && (i > 1))
                temp = (zx * zx) - (zy * zy) + cX
                zy = (2.0 * zx * zy) + cY
                zx = temp
                i -= 1
            end

            pixels[x,y] = (i << 21) + (i << 10) + i * 8
        end
    end

    return pixels
end

max_stars_repo_name: James-P-D/Fractals

max_issues_repo_name: James-P-D/Fractals

max_forks_repo_name: James-P-D/Fractals


text:
push!(LOAD_PATH, string(pwd(),"/src"))

import ModifiedStokesSolver
using Base.Test

mss = ModifiedStokesSolver

@testset "Split" begin

    srand(1)
    r = 1 .- 2.*(rand(), rand())
    f = 1 .- 2.*(rand(), rand())
    n = 1 .- 2.*(rand(), rand())
    rnorm = sqrt(r[1]^2+r[2]^2)

    cutoff = mss.POWERSERIES_ZMAX

    @testset "T1-T3" begin
        # Test split of stresslet kernel functions
        z = 0.987 # This value should work well for both
        (T1,T2,T3) = mss.Ti_direct(z)
        (T1S, T1L, T2S, T2L, T3S, T3L) =
            mss.Ti_split(z)
        E1 = (T1S + T1L*log(z) - T1)/T1
        E2 = (T2S + T2L*log(z) + 1/(8*pi*z^2) - 1/(pi*z^4)
              - T2)/T2
        E3 = (T3S + T3L*log(z) - T3)/T3

        @test abs(E1) < 1e-14
        @test abs(E2) < 1e-14
        @test abs(E3) < 1e-14

        if cutoff>0
            # Test with dispatcher for large and small arguments
            # 3.5 is the largest argument that works with explicit split
            for z=[cutoff/2 cutoff*2 1e-5 3.5]
                @show z
                (T1,T2,T3) = mss.Ti(z)
                (T1S, T1L, T2S, T2L, T3S, T3L) =
                    mss.Ti_split(z)
                E1 = (T1S + T1L*log(z) - T1)/T1
                E2 = (T2S + T2L*log(z) + 1/(8*pi*z^2) - 1/(pi*z^4)
                      - T2)/T2
                E3 = (T3S + T3L*log(z) - T3)/T3

                @test abs(E1) < 2e-14
                @test abs(E2) < 1e-14
                @test abs(E3) < 1e-14

            end
        end
    end

    @testset "Derivatives" begin
        @testset "Split" begin
            z = 0.9999*ModifiedStokesSolver.POWERSERIES_ZMAX_PRIME
            err = abs.( (mss.Ti_split_pow(z, true) .- mss.Ti_split_direct(z, true)) ./ mss.Ti_split_direct(z, true) )
            @test maximum(err) < 1e-13
        end
        @testset "Recombine" begin
            z = 0.99*ModifiedStokesSolver.POWERSERIES_ZMAX_PRIME
            err = abs.( (mss.Ti_prime(z) .- mss.Ti_direct_prime(z)) ./ mss.Ti_direct_prime(z) )
            @test maximum(err) < 3e-14
        end
    end

    @testset "Stresslet" begin
        # Test split stresslet
        r = rand(2)
        f = rand(2)
        n = rand(2)
        alpha = rand()
        u1 = mss.stresslet(r, f, n, alpha)
        (uS, uL, uC, uQ) = mss.stresslet_split(r, f, n, alpha)
        rnorm = norm(r)
        rdotf = dot(r,f)
        rdotn = dot(r,n)
        u2 = uS + uL*log(rnorm) + uC*rdotn/rnorm^2 +
            uQ*rdotf*r*rdotn/rnorm^4
        E = norm(u1-u2,Inf)/norm(u2,Inf)
        @test E < 1e-14
    end
end


text:
function string_code_unit(s::String, i::Integer)
    if i <= length(s)
        return codeunit(s,i)
    else
        return 0x00
    end
end

function stringsort!(vs::AbstractVector{T}) where T<:String
    l = length(vs)
    ts = similar(vs)

    # Init
    iters = maximum(length.(vs))

    bin = zeros(UInt32, 256, iters)
    # Histogram for each element, radix
    for i = 1:l
        vsi = vs[i]
        for j =1:iters
            idx = string_code_unit(vsi,j) + 1
            @inbounds bin[idx,j] += 1
        end
    end

    # Sort!
    swaps = 0
    for j = iters:-1:1
        # Unroll first data iteration, check for degenerate case
        idx = string_code_unit(vs[l],j) + 1

        # are all values the same at this radix?
        if bin[idx,j] == l;  continue;  end

        cbin = cumsum(bin[:,j])
        ci = cbin[idx]
        ts[ci] = vs[l]
        cbin[idx] -= 1

        # Finish the loop...
        @inbounds for i in l-1:-1:1
            idx = string_code_unit(vs[i],j) + 1
            ci = cbin[idx]
            ts[ci] = vs[i]
            cbin[idx] -= 1
        end
        vs,ts = ts,vs
        print(now())
        swaps += 1
    end

    if isodd(swaps)
        vs,ts = ts,vs
        for i = 1:l
            @inbounds vs[i] = ts[i]
        end
    end
    vs
end

const M=100_000_000; const K=100
srand(1)
@time svec1 = rand(["i"*dec(k,rand(1:7)) for k in 1:M÷K], M)
#@time svec1 = ["i"*dec(k,rand(1:7)) for k in 1:M÷K]
@time stringsort!(svec1)
issorted(svec1)

@code_warntype stringsort!(svec1)

sorttwo2!(load_bits.(svec1), svec1)

string_code_unit("id001",6)

length("id001")


text:
# Goal: Test the correctness of state initialization
using Test

a, b = 1.1, -0.5
f(u, p, t) = [a*u[1], b*u[2]]
u0 = [0.1, 1.0]
tspan = (0.0, 5.0)
t0, T = tspan
prob = ODEProblem(f, u0, tspan)
p = prob.p

# True Solutions and derivatives
u(t) = [a^0*u0[1] * exp(a*t), u0[2] * exp(b*t)]
du(t) = [a^1*u0[1] * exp(a*t), b * u0[2] * exp(b*t)]
ddu(t) = [a^2*u0[1] * exp(a*t), (b)^2 * u0[2] * exp(b*t)]
dddu(t) = [a^3*u0[1] * exp(a*t), (b)^3 * u0[2] * exp(b*t)]
ddddu(t) = [a^4*u0[1] * exp(a*t), (b)^4 * u0[2] * exp(b*t)]
dddddu(t) = [a^5*u0[1] * exp(a*t), (b)^5 * u0[2] * exp(b*t)]
ddddddu(t) = [a^6*u0[1] * exp(a*t), (b)^6 * u0[2] * exp(b*t)]
true_init_states = [u(t0); du(t0); ddu(t0); dddu(t0); ddddu(t0); dddddu(t0); ddddddu(t0)]

@testset "OOP state init" begin
    m0, P0 = ODEFilters.initialize_with_derivatives(prob.u0, prob.f, prob.p, prob.tspan[1], 6)
    @test m0 ≈ true_init_states
    @test all(P0 .== 0)
end

@testset "IIP state init" begin
    f!(du, u, p, t) = (du .= f(u, p, t))
    prob = ODEProblem(f!, u0, tspan)
    m0, P0 = ODEFilters.initialize_with_derivatives(prob.u0, prob.f, prob.p, prob.tspan[1], 6)
    @test m0 ≈ true_init_states
    @test all(P0 .== 0)
end


text:
export train_rvq, quantize_rvq, experiment_rvq

"""
    quantize_rvq(X, C, V=false) -> B, singletons

Given data and full-dimensional codebooks, quantize.

# Arguments
- `X::Matrix{T}`: `d`-by-`n` data to quantize
- `C::Vector{Matrix{T}}`: `m`-long vector with `d`-by-`h` matrix entries. Each matrix is a codebook.
- `V::Bool`: Whether to print progress

# Returns
- `B::Matrix{Int16}`: `m`-by-`n` codes that approximate `X`
- `singletons::Vector{Matrix{T}}`: `m` matrices with unused codebook entries
"""
function quantize_rvq(
  X::Matrix{T},         # d-by-n. Data to encode
  C::Vector{Matrix{T}}, # codebooks
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  d, n = size( X )
  m    = length( C )
  h    = size( C[1], 2 )
  B    = Vector{Vector{Int}}(undef, m) # codes
  for i = 1:m; B[i] = zeros(Int,n); end # Allocate codes

  singletons = Vector{Matrix{T}}(undef, m)

  # Residual after encoding the ith codebook
  Xr = copy(X)

  for i = 1:m
    if V print("RVQ encoding on codebook $i / $m... ") end

    # Find distances from X to codebook
    dmat = Distances.pairwise(Distances.SqEuclidean(), C[i], Xr)

    # Update the codes
    costs     = zeros(T, n)    # Out. The cost of assignment
    counts    = zeros(Int, h)  # Out. Number of points assigned to this centre
    to_update = ones(Bool, h)  # Out. Whether this centre needs update
    unused    = Vector{Int}()  # Out. Indices of centres that have no assignments

    Clustering.update_assignments!(dmat, true, B[i], costs, counts, to_update, unused)

    # Create new codebook entries that we are missing
    if !isempty(unused)
      temp_codebook = similar(C[i]) # Out. The new codebooks will be added here
      Clustering.repick_unused_centers(Xr, costs, temp_codebook, unused)
      singletons[i] = temp_codebook[:, unused]
    end

    # Update the residual
    Xr .-= C[i][:, B[i]]

    if V println("done"); end
  end

  B = hcat(B...)
  B = convert(Matrix{Int16}, B)
  B = collect(B')

  B, singletons
end

"""
    train_rvq(X, m, h, niter=25, V=false) -> C, B, error

Train a residual quantizer.

# Arguments
- `X::Matrix{T}`: `d`-by-`n` data to quantize
- `m::Integer`: Number of codebooks
- `h::Integer`: Number of entries in each codebook (typically 256)
- `niter::Integer`: Number of iterations to use
- `V::Bool`: Whether to print progress

# Returns
- `C::Vector{Matrix{T}}`: `m`-long vector with `d`-by-`h` matrix entries. Each matrix is a codebook of size approximately `d/m`-by-`h`.
- `B::Matrix{Int16}`: `m`-by-`n` matrix with the codes
- `error::T`: The quantization error after training
"""
function train_rvq(
  X::Matrix{T},       # d-by-n. Data to learn codebooks from
  m::Integer,         # number of codebooks
  h::Integer,         # number of entries per codebook
  niter::Integer=25,  # Number of k-means iterations for training
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  d, n = size( X )

  C = Vector{Matrix{T}}(undef, m) # codebooks
  B = zeros(Int16, n, m) # codes

  # Residual
  Xr = copy(X)

  for i = 1:m
    if V print("Working on codebook $i / $m... "); end
    # FAISS uses 25 iterations by default
    # See https://github.com/facebookresearch/faiss/blob/master/Clustering.cpp#L28
    cluster = kmeans(Xr, h, init=:kmpp, maxiter=niter)
    C[i], B[:,i] = cluster.centers, cluster.assignments

    # Update the residual
    Xr .-= C[i][:,B[:,i]]

    if V
      println("done.")
      println("  Ran for $(cluster.iterations) iterations")
      println("  Error after codebook $(i) is $(cluster.totalcost ./ n)")
      println("  Converged: $(cluster.converged)")
    end
  end

  B = collect(B')
  error = qerror(X, B, C)

  C, B, error
end

function experiment_rvq(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  Xb::Matrix{T}, # d-by-n. Base set
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of k-means iterations for training
  knn::Integer=1000,
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  # === RVQ train ===
  d, _ = size(Xt)
  C, B, train_error = Rayuela.train_rvq(Xt, m, h, niter, V)
	norms_B, norms_C = get_norms_codebook(B, C)

  # === Encode the base set ===
  B_base, _ = Rayuela.quantize_rvq(Xb, C, V)
  base_error = qerror(Xb, B_base, C)
  if V; @printf("Error in base is %e\n", base_error); end

  # Compute and quantize the database norms
  B_base_norms, _ = quantize_norms(B_base, C, norms_C)
  db_norms        = vec(norms_C[ B_base_norms ])

  # === Compute recall ===
  # B_base       = convert(Matrix{UInt8}, B_base-1)
  # B_base_norms = convert(Vector{UInt8}, B_base_norms-1)

  if V; print("Querying m=$m ... "); end
  @time dists, idx = linscan_lsq(B_base, Xq, C, db_norms, Matrix{Float32}(1.0I, d, d), knn)
  if V; println("done"); end

  recall = eval_recall(gt, idx, knn)
  return C, B, train_error, B_base, recall
end

function experiment_rvq_query_base(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of k-means iterations for training
  knn::Integer=1000,
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  # === RVQ train ===
  d, _ = size(Xt)
  C, B, train_error = Rayuela.train_rvq(Xt, m, h, niter, V)
	norms_B, norms_C = get_norms_codebook(B, C)
  db_norms     = vec(norms_C[ norms_B ])

  # === Compute recall ===
  # B_base       = convert(Matrix{UInt8}, B_base-1)
  # B_base_norms = convert(Vector{UInt8}, B_base_norms-1)
  if V; print("Querying m=$m ... "); end
  @time dists, idx = linscan_lsq(B, Xq, C, db_norms, Matrix{Float32}(1.0I, d, d), knn)
  if V; println("done"); end

  recall = eval_recall(gt, idx, knn)
  return C, B, train_error, recall
end