text:
module SamplingTools

    using LinearAlgebra
    using Statistics
    using Base64
    using Interpolations

    function hj_info(h::Dict, J::Dict)
        n = length(h)
        M = zeros(Float64, (n, n))
        degrees = zeros(Int64, n)
        j_full = Dict()
        for ((x, y), val) in J
            M[sort([x, y])...] = val
            degrees[x] += 1
            degrees[y] += 1
        end
        M += M'
        for (key, val) in h
            M[key, key] = val
            j_full[(key, key)] = val
        end
        return Array{Float64, 2}(UpperTriangular(M)), merge(J, j_full), degrees
    end

    function hj_info!(M::Array{Float64, 2}, j_full::Dict,
                      degrees::Array{Int64, 1}, h::Dict, J::Dict)
        n = length(h)
        fill!(degrees, 0.0)
        for ((x, y), val) in J
            M[x, y] = val
            degrees[x] += 1
            degrees[y] += 1
        end
        M += M'
        for (key, val) in h
            M[key, key] = val
            j_full[(key, key)] = val
        end
        return M, j_full, degrees
    end

    function generate_nbr(nspins::Int,
                          J::Dict,
                          max_degree::Int)
        # define vars
        nbs_vals = zeros(Float64, (max_degree, nspins))
        nbs_idx = zeros(Int64, (max_degree, nspins))
        ispin = 1
        ipair = 1
        x, y = 1, 1
        val = 0.0
        for ispin = 1:nspins
            ipair = 1
            for ((x, y), val) in J
                if x == ispin
                    nbs_idx[ipair, ispin] = y
                    nbs_vals[ipair, ispin] = val
                    ipair += 1
                elseif y == ispin
                    nbs_idx[ipair, ispin] = x
                    nbs_vals[ipair, ispin] = val
                    ipair += 1
                end
            end
        end
        return nbs_idx, nbs_vals
    end

    function generate_nbr!(nbs_idx::Array{Int64, 2}, nbs_vals::Array{Float64, 2},
                           nspins::Int,
                           J::Dict,
                           max_degree::Int)
        # define vars
        ispin = 1
        ipair = 1
        x, y = 1, 1
        val = 0.0
        for ispin = 1:nspins
            ipair = 1
            for ((x, y), val) in J
                if x == ispin
                    nbs_idx[ipair, ispin] = y
                    nbs_vals[ipair, ispin] = val
                    ipair += 1
                elseif y == ispin
                    nbs_idx[ipair, ispin] = x
                    nbs_vals[ipair, ispin] = val
                    ipair += 1
                end
            end
        end
        return nbs_idx, nbs_vals
    end

    function generate_nbr_dict(h::Dict{Int64, Float64},
                               J::Dict{Tuple{Int64, Int64}, Float64})
        # define vars
        nbs = Dict{Int64, Dict{Int64, Float64}}()
        nbs_edges = Dict()
        nspins = length(h)

        for ispin = 1:nspins
            local_nbr = Dict{Int64, Float64}(ispin => h[ispin])
            local_edges = []
            for ((x, y), val) in J
                if x == ispin
                    push!(local_nbr, y => val)
                    push!(local_edges, (x, y))
                elseif y == ispin
                    push!(local_nbr, x => val)
                    push!(local_edges, (x, y))
                end
            end
            push!(nbs, ispin => local_nbr)
            push!(nbs_edges, ispin => local_edges)
        end
        return nbs, nbs_edges
    end

    function shuffle(sorted_array::Array{Int64, 1}, n::Int64)
        t = 0
        j = 0
        @fastmath @inbounds for i = n:-1:1
            j = (rand(sorted_array) % i) + 1
            t = sorted_array[i]
            sorted_array[i] = sorted_array[j]
            sorted_array[j] = t
        end
    end

    function default_temp_range(M::Array{Float64, 2})
        JM = map(abs, M)
        min_de = minimum(diag(JM))
        max_de = maximum(sum(JM, dims=2))

        hot_temp = max_de / log(2)
        cold_temp = min_de / log(1000)
        return hot_temp, cold_temp
    end

    function ising_energy(spin_vector::Array{Float64, 1},
                          JM::Array{Float64, 2}; b_field=1.0
        )::Float64
        b = diag(JM)
        return b_field * ((spin_vector'* ((JM - diagm(b)) * spin_vector)) + (b' * spin_vector))
    end

    function discrete_aggregation(spin_vectors::Array{Array{Float64, 1}, 1},
                                  num_reads, JM::Array{Float64, 2})
        # spin_vectors = [spin_vectors[:, i] for i =1:num_reads]
        unq_vectors = unique(spin_vectors)
        num_unique = size(unq_vectors, 1)
        energies_unq = zeros(Float64, num_unique)
        occurences_unq = zeros(Int, num_unique)
        for i = 1:num_unique
            uvec = unq_vectors[i]
            energies_unq[i] = ising_energy(uvec, JM)
            occurences_unq[i] = count(x -> x == uvec, spin_vectors)
        end
        return unq_vectors, energies_unq, occurences_unq
    end

    function rotor_aggregation(rotor_vectors::Array{Array{Float64, 1}, 1}, num_reads,
                               JM::Array{Float64, 2}; cluster=false)
        if cluster
            spin_vectors = map(x-> sign.(x), map(x-> sin.(x), rotor_vectors))
        else
            spin_vectors = map(x-> sign.(x), map(x-> cos.(x), rotor_vectors))
        end
        unq_vectors = unique(spin_vectors)
        num_unique = size(unq_vectors, 1)
        energies_unq = zeros(Float64, num_unique)
        occurences_unq = zeros(Int, num_unique)
        for i = 1:num_unique
            uvec = unq_vectors[i]
            energies_unq[i] = ising_energy(uvec, JM)
            occurences_unq[i] = count(x -> x == uvec, spin_vectors)
        end
        return unq_vectors, energies_unq, occurences_unq
    end

    function pimc_aggregation(spin_lattice::Array{Array{Float64, 2}, 1},
                              num_reads, JM::Array{Float64, 2})
        # spin_vectors = [spin_vectors[:, i] for i =1:num_reads]
        nslices = size(spin_lattice[1], 2)
        spin_vectors = Array{Array{Float64, 1}, 1}()
        for i = 1:num_reads
            v, e, o = discrete_aggregation([spin_lattice[i][:, j] for j =1:nslices], nslices, JM)
            push!(spin_vectors, collect(v[findmin(e)[2]]))
        end
        return discrete_aggregation(spin_vectors, num_reads, JM)
    end

    function bootstrap(sample_set::Array{Float64, 1},
                       bootstraps=1000, return_bounds=true)
        # for finding a bootstrap median
        n = length(sample_set)
        sample_set = repeat(sample_set, outer=(1,bootstraps))
        @inbounds bootstrap_samples = sample_set[rand(1:n, n, bootstraps)]
        medians = vcat(median(bootstrap_samples, dims=1)...)
        m = median(medians)
        if return_bounds
            bounds = quantile(medians, [0.05, 0.95])
            return m, bounds
        end
        return m
    end

    function add_cross_talk(h::Dict, j::Dict, susceptibility=-0.05)
        nbs, nbs_edges = generate_nbr_dict(h, j)
        j_p = merge(j, Dict([((y, x), val) for ((x, y), val) in j]))
        h_c = copy(h)
        j_c = copy(j)
        j_p_keys = collect(keys(j_p))
        j_c_keys = collect(keys(j_c))
        for (node, n_edges) in nbs_edges
            h_off = 0.0
            ne = length(n_edges)
            for i1 = 1:ne
                (x, y) = n_edges[i1]
                if node == y
                    k = x
                else
                    k = y
                end
                h_off += j[(x, y)] * h[k]
                for i2 = 1:ne
                    (u, v) = n_edges[i2]
                    if ((x, y) != (u, v)) & (i1 < i2)
                        if v == node
                            a = u
                        else
                            a = v
                        end
                        if k < a
                            jij = (k, a)
                        else
                            jij = (a, k)
                        end
                        j_off = susceptibility * get(j_p, (node, a), 0.0) * get(j_p, (node, k), 0.0)
                        if (jij in j_c_keys) == false
                            push!(j_c_keys, jij)
                            push!(j_c, jij => 0.0)
                        end
                        j_c[jij] += j_off
                    end
                end
            end
            h_c[node] += susceptibility * h_off
        end
        return h_c, j_c
    end

    function add_noise(h, j, h_std=0.023, j_std=0.03, h_mean=0.0, j_mean=0.0)
        h_noise = (h_std .* randn(length(h))) .+ h_mean
        j_noise = (j_std .* randn(length(j))) .+ j_mean
        h_new = Dict(Tuple(zip(keys(h), values(h) .+ h_noise)))
        j_new = Dict(Tuple(zip(keys(j), values(j) .+ j_noise)))
        return h_new, j_new
    end

    function PFC_dict(d::Float64; N=2, R=1.0)
        h = Dict{Int64, Float64}([(i, 0.0) for i = 1:(2*N)])
        j = Dict{Tuple, Float64}()
        @assert R > d
        for i = 1:N
            h[i] = R - d
            h[i+N] = -1.0 * R
            j[(i, i+N)] = -1.0 * R
            if i < N
                j[(i, i+1)] = -1.0 * R
            end
        end
        return h, j
    end

    function energy_ranges(h::Dict{Int64, Float64}, j::Dict{Tuple, Float64},
                           max_field=1.0)
        n = length(h)
        @assert n >= 1
        abs_h = []
        abs_dict = Dict{Int64, Float64}()
        abs_j = []
        # bias contrib
        for (key, hh) in h
            abs_h_dict[key] = abs(hh)
            if hh != 0.0
                push!(abs_h, abs(hh))
            end
        end
        # j contrib
        for ((i1, i2), jj) in j
            ajj = abs(jj)
            abs_dict[i1] += ajj
            abs_dict[i2] += ajj
            if jj != 0.0
                push!(abs_j, ajj)
            end
        end

        abs_vals = vcat(abs_h, abs_j)
        if length(abs_vals) == 0
            return [0.0, 0.0]
        end

        min_de = minimum(abs_vals)
        max_de = maximum(values(abs_dict))
        return [min_de, max_de]
    end

    function DQA_schedule(k::Int, nspins::Int, a_schedule::Array{Float64, 1},
                          b_schedule::Array{Float64, 1}; sx=0.2, cx=0.0, c1=0.0)
        num_steps = length(a_schedule)
        new_steps = convert(Int64, ceil(num_steps / (1.0 - sx)))
        s_old = LinRange(sx, 1.0, num_steps)
        s = collect(LinRange(0.0, 1.0, new_steps))
        ak_sch = collect(LinRange(cx, c1, num_steps))
        a_dqa_dict = Dict()
        b_dqa_dict = Dict()
        for i = 1:nspins
            a_sch = deepcopy(a_schedule)
            b_sch = deepcopy(b_schedule)
            if i == k
                afn = monotonic_interpolation(s_old, ak_sch)
                a_dqa_dict[i] = afn.(s)
                b_dqa_dict[i] = LinRange(-sx/(1.0-sx), 1.0, new_steps)
            else
                afn = monotonic_interpolation(s_old, a_sch)
                bfn = monotonic_interpolation(s_old, b_sch)
                a_dqa_dict[i] = afn.(s)
                b_dqa_dict[i] = bfn.(s)
            end
        end
        return a_dqa_dict, b_dqa_dict
    end

    function monotonic_interpolation(x, y)
        return extrapolate(interpolate(x, y, SteffenMonotonicInterpolation()), Flat())
    end

    export rotor_aggregation, discrete_aggregation, ising_energy,
            generate_nbr, default_temp_range, hj_info, shuffle,
            generate_nbr_dict, bootstrap, add_cross_talk, add_noise,
            hj_info!, generate_nbr!, pimc_aggregation, PFC_dict,
            energy_ranges, DQA_schedule, monotonic_interpolation
end

max_stars_repo_name: dtoconnor/Annealing-Samplers

max_issues_repo_name: dtoconnor/Annealing-Samplers

max_forks_repo_name: dtoconnor/Annealing-Samplers


text:
export TabularQ, TabularV, AggregationV, PolynomialV, FourierV, Tiling, TilingsV, TilingsQ, LinearV, NeuralNetworkQ

include("tabular_Q.jl")
include("tabular_V.jl")
include("aggregation_V.jl")
include("polynomial_V.jl")
include("fourier_V.jl")
include("tilings.jl")
include("linear_V.jl")
include("neural_networks.jl")

ext: jl


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Foundations of machine learning")

# # Modeling data 2

#-

# ## Building a model
#
# Recall that in notebook 3, we saw that we could use a mathematical function to classify an image as an apple or a banana, based on the average amount of green in an image:
#
#
#
#
#

#-

# A common function for performing this kind of **classification** is the sigmoid that we saw in the last notebook, and that we will now extend by adding two **parameters**, $w$ and $b$:
#
# $$\sigma(x; w, b) := \frac{1}{1 + \exp(-wx + b)}$$
#
# $$ x = \mathrm{data} $$
#
# \begin{align}
# \sigma(x;w,b) &\approx 0 \implies \mathrm{apple} \\
# \sigma(x;w,b) &\approx 1 \implies \mathrm{banana}
# \end{align}

#-

# In our mathematical notation above, the `;` in the function differentiates between the **data** and the **parameters**. `x` is the data and is determined from the image. The parameters, `w` and `b`, are numbers which we choose to make our function match the results it should be modeling.
#
# Note that in the code below, we don't distinguish between data and parameters - both are just inputs to our function, σ!

using Images, Statistics

apple = load(datapath("data/10_100.jpg"))
banana = load(datapath("data/104_100.jpg"))

apple_green_amount = mean(Float64.(green.(apple)))
banana_green_amount = mean(Float64.(green.(banana)))

println("Average green for apple = $apple_green_amount")
println("Average green for banana = $banana_green_amount")

#-

σ(x, w, b) = 1 / (1 + exp(-w * x + b))

# What we want is that when we give σ as input the average green for the apple, roughly `x = 0.3385`, it should return as output something close to 0, meaning "apple". And when we give σ the input `x = 0.8808`, it should output something close to 1, meaning "banana".
#
# By changing the parameters of the function, we can change the shape of the function, and hence make it represent, or **fit**, the data better!

#-

# ## Data fitting by varying parameters

#-

# We can understand how our choice of `w` and `b` affects our model by seeing how our values for `w` and `b` change the plot of the $\sigma$ function.

using Plots; gr()   # GR works better for interactive manipulations

# Run the code in the next cell. You should see two "sliders" appear, one for `w` and one for `b`.
#
# **Game**:
# Change w and b around until the blue curve, labeled "model", which is the graph of the `\sigma` function, passes through *both* of the data points at the same time.

w = 10.0 # try manipulating w between -10 and 30
b = 10.0 # try manipulating b between 0 and 20

plot(x -> σ(x, w, b), xlim=(-0,1), ylim=(-0.1,1.1), label="model", legend=:topleft, lw=3)

scatter!([apple_green_amount],  [0.0], label="apple", ms=5)   # marker size = 5
scatter!([banana_green_amount], [1.0], label="banana", ms=5)

# Notice that the two parameters do two very different things. The **weight**, `w`, determines *how fast* the transition between 0 and 1 occurs. It encodes how trustworthy we think our data  actually is, and in what range we should be putting points between 0 and 1 and thus calling them "unsure". The **bias**, `b`, encodes *where* on the $x$-axis the switch should take place. It can be seen as shifting the function left-right. We'll come to understand these *parameters* more in notebook 6.
#
# Here are some parameter choices that work well:

w = 25.58; b = 15.6

plot(x -> σ(x, w, b), xlim=(0,1), ylim=(-0.1,1.1), label="model", legend=:topleft, lw=3)

scatter!([apple_green_amount], [0.0], label="apple")
scatter!([banana_green_amount],[1.0], label="banana")

# (Note that in this problem there are many combinations of `w` and `b` that fit the data well.)

#-

# Once we have a model, we have a computational representation for how to choose between "apple" and "banana". So let's pull in some new images and see what our model says about them!

apple2 = load(datapath("data/107_100.jpg"))

#-

green_amount = mean(Float64.(green.(apple2)))
@show green_amount

scatter!([green_amount], [0.0], label="new apple")

# Our model successfully says that our new image is an apple! Pat yourself on the back: you've actually just trained your first neural network!

#-

# #### Exercise 1
#
# Load the image of a banana in `data/8_100.jpg` as `mybanana`. Edit the code below to calculate the amount of green in `mybanana` and to overlay data for this image with the existing model and data points.

#-

# #### Solution

#-

# To get the desired overlay, the code we need is

mybanana = load(datapath("data/8_100.jpg"))
mybanana_green_amount = mean(Float64.(green.(banana)))
scatter!([mybanana_green_amount], [1.0], label="my banana")

# ## Closing remarks: bigger models, more data, more accuracy
#
# That last apple should start making you think: not all apples are red; some are yellow. "Redness" is one attribute of being an apple, but isn't the whole thing. What we need to do is incorporate more ideas into our model by allowing more inputs. However, more inputs would mean more parameters to play with. Also, we would like to have the computer start "learning" on its own, instead of modifying the parameters ourselves until we think it "looks right". How do we take the next step?
#
# The first thing to think about is, if you wanted to incorporate more data into the model, how would you change the sigmoid function? Play around with some ideas. But also, start thinking about how you chose parameters. What process did you do to finally end up at good parameters? These two problems (working with models with more data and automatically choosing parameters) are the last remaining step to understanding deep learning.

max_stars_repo_path: Courses/Foundations of machine learning/0500.Building-models.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
# # Counting residues
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__notebooks/10_Information.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__notebooks/10_Information.ipynb)
#
# MIToS was designed to perform fast counting of residues. To achieve that,
# each `Residue` is encoded as an integer that can be used to index.
#

using MIToS.MSA
res = Residue('C')

#-

Int(res)

# Let's count residues in the *UBA domain* Pfam family.

using MIToS.Pfam
const msa_file = downloadpfam("PF16577")

#-

const msa = read(msa_file, Stockholm)

# We can access the MSA as a matrix to get a particular column or sequence:
# ```julia
# msa[:, column]
# msa[sequence, :]
# ```

col_i = msa[:, 10]

# MIToS Information module has functions to get residue frequencies (counts or
# fractions) and to return contingency tables:

using MIToS.Information

#-

Ni = count(col_i)

#-

Pi = probabilities(col_i)

# `Counts` and `Probabilities` objects wrap a `ContingencyTable`:

Pi.table

# They can be quickly indexed by using `Residue`s:

Pi[Residue('P')]

# You can create a bidimensional `ContingencyTable` by using two MSA columns or
# sequences (there is no limit to the number of dimensions):

col_j = msa[:, 15]

Nij = count(col_i, col_j)

#-

Nij[Residue('P'), Residue('R')]

# You can also get probabilities by normalizing counts:

Pij = normalize(Nij.table)

#-

Pij[Residue('P'), Residue('R')]

# ## Plotting counts
#
# We use a `heatmap` to plot a bidimensional table:
#

using Plots
gr()

# We can list the available color libraries:

clibraries()

# Visualize the color maps defined in a library:

showlibrary(:Plots)

# And select the color library we want to use:

clibrary(:cmocean)

# For this plot, we use the interactive *Plotly* backend

plotlyjs()

# We need a vector of strings for the ticks labels...

x = string.(Residue.(UngappedAlphabet()))

# ...because `xticks` and `xticks` keyword arguments take a tuple:
# `(positions, labels)`

heatmap(Nij,
		color=:tempo,
		ratio=:equal,
		xticks=(1:length(x), x),
		yticks=(1:length(x), x))

# ## Functions taking counts and probabilities
#
# Information defines different functions that take contingency tables as an
# input, for example:

entropy(Ni)

#-

mutual_information(Nij)

# ### Conservation
#
# You can use the MIToS' Information module to calculate the Shannon entropy of
# each MSA column:

using MIToS.Information

const alphabet = UngappedAlphabet()

# You also have `GappedAlphabet` and `ReducedAlphabet`:

ReducedAlphabet("(AILMV)(RHK)(NQST)(DE)(FWY)CGP")

# Empty contingency tables can be defined by using
# `ContingencyTable(element_type, Val{dimensions}, alphabet)`:

const table = ContingencyTable(Int, Val{1}, alphabet)

# The `mapcolfreq!`, `mapseqfreq!`, `mapcolpairfreq!` and `mapseqpairfreq!`
# functions from the Information module apply a function on the table filled
# using each column, sequence, column pair or sequence pair of the MSA,
# respectively:

Hx = mapcolfreq!(entropy, msa, Counts(table))

# We use **Plots.jl** to visualize the entropy of each column to find variable
# and conserved regions

gr()

#-

plot_entropy = plot(vec(Hx),
					color=:orange,
					fill=0,
					fillalpha=0.5,
					legend=false,
					ylab="Entropy")

#-

plot_msa = plot(msa, legend=false)

#-

plot(plot_entropy,
	 plot_msa,
	 layout=grid(2, 1, heights=[0.2, 0.8]),
	 link=:x)

# ## Multivariate mutual information
#
# We are going to measure mutual information in MSA column triplets. We use
# Julia parallelism because the number of combinations is high.
#

using Distributed

#-

nprocs()

#-

addprocs(3)

#-

nprocs()

# We are going to use [Combinatorics](https://github.com/JuliaMath/Combinatorics.jl)
# everywhere (in all the process) to get a generator of the combinations and
# [ProgressMeter](https://github.com/timholy/ProgressMeter.jl) to see the
# progress of the parallel loop.

@everywhere using Combinatorics
@everywhere using ProgressMeter
@everywhere using MIToS.MSA
@everywhere using MIToS.Information

# #### Exercise 1
#
# Fill the missing parts of the function to calculate mutual information
# between MSA columns i, j and k.

function parallel_mmi(msa)
	ncols = ncolumns(msa)
	triplets = combinations(1:ncols, 3)
	@showprogress pmap(triplets) do (i, j, k)
##           ...to complete...
##           (i, j, k) => ...to complete...
#src		table = probabilities(msa[:, i], msa[:, j], msa[:, k])
#src		(i, j, k) => mutual_information(table)
	end
end

#-

mi_values = parallel_mmi(msa)

max_stars_repo_name: diegozea/JuliaForBioinformatics

max_issues_repo_name: diegozea/JuliaForBioinformatics

max_forks_repo_name: diegozea/JuliaForBioinformatics


text:
## regression loss

export mae
export maeLoss, maeCost
export L1Loss, L1Cost

export mse
export mseLoss, mseCost
export L2Loss, L2Cost

export Lp, LpLoss, LpCost

"""
    mae(x::Variable{T}, label::Variable{T}) -> y::Variable{T}

mean absolute error (mae) between each element in the input `x` and target `label`. Also called L1Loss. i.e. ⤦\n
    y = |xᵢ - lᵢ|
"""
function mae(x::Variable{T}, label::Variable{T}) where T
    @assert (x.shape == label.shape)
    backprop = (x.backprop || label.backprop)
    y = Variable{T}(abs.(ᵛ(x) - ᵛ(label)), backprop)
    if backprop
        function maeBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* sign.(ᵛ(y))
            end
            ifNotKeepδThenFreeδ!(y)
        end
        push!(graph.backward, maeBackward)
    end
    return y
end

maeLoss(x::Variable{T}, label::Variable{T}; reduction::String="sum") where T = loss( mae(x, label), reduction=reduction )
L1Loss(x::Variable{T},  label::Variable{T}; reduction::String="sum") where T = loss( mae(x, label), reduction=reduction )
maeCost(x::Variable{T}, label::Variable{T}) where T = cost( mae(x, label) )
L1Cost(x::Variable{T},  label::Variable{T}) where T = cost( mae(x, label) )

"""
    mse(x::Variable{T}, label::Variable{T}) -> y::Variable{T}

mean sqrt error (mse) between each element in the input `x` and target `label`. Also called L2Loss. i.e. ⤦\n
    y = (xᵢ - lᵢ)²
"""
function mse(x::Variable{T}, label::Variable{T}) where T
    @assert (x.shape == label.shape)
    backprop = (x.backprop || label.backprop)
    𝟚 = eltype(x)(2.0f0)
    y = Variable{T}((ᵛ(x) - ᵛ(label)).^𝟚, backprop)
    if backprop
        function mseBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* 𝟚 .* (ᵛ(x) - ᵛ(label))
            end
            ifNotKeepδThenFreeδ!(y)
        end
        push!(graph.backward, mseBackward)
    end
    return y
end

mseLoss(x::Variable{T}, label::Variable{T}; reduction::String="sum") where T = loss( mse(x, label), reduction=reduction )
L2Loss(x::Variable{T},  label::Variable{T}; reduction::String="sum") where T = loss( mse(x, label), reduction=reduction )
mseCost(x::Variable{T}, label::Variable{T}) where T = cost( mse(x, label) )
L2Cost(x::Variable{T},  label::Variable{T}) where T = cost( mse(x, label) )

"""
    Lp(x::Variable{T}, label::Variable{T}; p=3) -> y::Variable{T}

absolute error's `p`-th power between each element in the input `x` and target `label`. Also called LpLoss. i.e. ⤦\n
    y = |xᵢ - lᵢ|ᵖ
"""
function Lp(x::Variable{T}, label::Variable{T}; p=3) where T
    @assert (x.shape == label.shape)
    backprop = (x.backprop || label.backprop)
    Δ = ᵛ(x) - ᵛ(label)
    y = Variable{T}(Δ .^ p, backprop)
    if backprop
        function LpBackward()
            if need2computeδ!(x)
                i = (Δ .!= eltype(T)(0.0))
                x.delta[i] .+= y.delta[i] .* y.value[i] ./ Δ[i] .* p
                # δ(x) .+= δ(y) .* ᵛ(y) ./ Δ .* p
            end
            ifNotKeepδThenFreeδ!(y)
        end
        push!(graph.backward, LpBackward)
    end
    return y
end

LpLoss(x::Variable{T}, label::Variable{T}; p=3, reduction::String="sum") where T = loss( Lp(x, label; p=p), reduction=reduction )
LpCost(x::Variable{T}, label::Variable{T}; p=3) where T = cost( Lp(x, label; p=p) )

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
module InferGMM
using Flux, Flux.Tracker
using Flux: ADAM
using StatsFuns: logsumexp
using Distributions
using Formatting
using Random # randperm, MersenneTwister
using LinearAlgebra  #: cholesky, logdet, diag, inv, triu
using NNlib: softmax
using AxUtil # Math, array, Flux extensions
using BSON, Logging, ProgressMeter
using Pkg

export GMM, importance_sample, responsibilities, importance_sample, gmm_fit

include("misc.jl")
using .Misc

include("llh.jl")
include("gmm.jl")
using .gmm

include("variational.jl")
using .variational

include("scm.jl")
using .scm

include("amis.jl")

end


text:
"""
Primitive Gaussian basis function.
"""
struct PGBF
    center::Tuple3F64
    power::Tuple3I64
    expn::Float64
    NORM::Float64
end

"""
Construct and instace of PGBF with exponent `expn`,
centered at `center` and angular momentum `power`.
"""
function PGBF( expn::Float64, center::Tuple3F64, power::Tuple3I64 )
    NORM = 1.0 / sqrt(overlap(expn,center,power))
    p = PGBF( center, power, expn, NORM )
    return p
end

"""
Construct an instance of PGBF with exponent `expn` and
center at `center`, with
angular momentums: (0, 0, 0) as the default value.
"""
function PGBF( expn::Float64, center::Tuple3F64 )
    power = (0,0,0)
    NORM = 1.0 / sqrt(overlap(expn,center,power))
    p = PGBF( center, power, expn, NORM )
    return p
end

"""
Construct an instance of PGBF with exponent `expn` with the following
default values:
- centered at (0.0, 0.0, 0.0)
- with angular momentums: (0, 0, 0)
"""
function PGBF( expn::Float64 )
    center = (0.0,0.0,0.0)
    power = (0,0,0)
    NORM = 1.0 / sqrt(overlap(expn,center,power))
    p = PGBF( center, power, expn, NORM )
    return p
end

"""
Evaluate `bf::PGBF` at point `(x,y,z)`
"""
function evaluate( bf::PGBF, x::Float64, y::Float64, z::Float64 )
    xo = bf.center[1]
    yo = bf.center[2]
    zo = bf.center[3]
    I = bf.power[1]
    J = bf.power[2]
    K = bf.power[3]
    expn = bf.expn
    NORM = bf.NORM

    dx, dy, dz = x-xo, y-yo, z-zo
    r2 = dist2(dx, dy, dz)
    return NORM*(dx^I)*(dy^J)*(dz^K)*exp(-expn*r2)
end

import Base: println
"""
Display basic information about PGBF
"""
function println( bf::PGBF )
    @printf("\n")
    @printf("Primitive Gaussian:\n")
    @printf("Exponent: %18.10f\n", bf.expn)
    @printf("Center  : (%18.10f,%18.10f,%18.10f)\n", bf.center[1], bf.center[2], bf.center[3])
    @printf("AM      : (%2d,%2d,%2d)\n", bf.power[1], bf.power[2], bf.power[3])
    @printf("Norm    : %18.10f\n", bf.NORM)
    @printf("\n")
end


text:
using GridPotentialFlow
using Plots
Lx = 10.0
Ly = 4.0
xlim = (-Lx/2,Lx/2)
ylim = (-Ly/2,Ly/2)
g = PhysicalGrid(xlim, ylim, 0.04)

Rc = 1
circle = Circle(Rc,2*cellsize(g))
pfb = PotentialFlowBody(circle)

Δs = dlength(circle);

y∞ = Rc/2
v1 = Vortex(-Lx/2+10*cellsize(g),y∞,1.0);
v2 = Vortex(-Lx/2+10*cellsize(g),-y∞,-1.0);

x(r) = sqrt(r^2-(y∞^2*(r^2-1)^2)/((r^2-1)^2-4*y∞^2))
y(r) = y∞*(r^2-1)/sqrt((r^2-1)^2-4*y∞^2)
s = 0:0.001:π/2
rmin = 1.4375649
c1 = rmin/(-Lx/2+rmin)
c2 = -Lx/2*c1
r = c2 ./ (cos.(s).-c1)
x_trajectory = -reverse(x.(r)); append!(x_trajectory,x.(r))
y_trajectory_upper = reverse(y.(r)); append!(y_trajectory_upper,y.(r))
y_trajectory_lower = -y_trajectory_upper;

Px_func(x,y) = y - y/(x^2+y^2) + y - y/(x^2+y^2);
Py_func(x,y) = -x + x/(x^2+y^2) + x - x/(x^2+y^2);
Px_exact = Px_func.(x_trajectory,y_trajectory_upper);
Py_exact = Py_func.(x_trajectory,y_trajectory_upper);

model = VortexModel(g,bodies=[pfb],vortices=[v1,v2]);

function rhs(X,model,t)
    setvortexpositions!(model,X)
    Ẋ = vortexvelocities!(model)
    return Ẋ
end

import OrdinaryDiffEq
X = getvortexpositions(model)
prob = OrdinaryDiffEq.ODEProblem(rhs,X,(0.0,64.0),model);
sol = OrdinaryDiffEq.solve(prob,dt=0.1,OrdinaryDiffEq.RK4(),dense=false,adaptive=false);

Px_numerical_hist = []
Py_numerical_hist = []
for u in sol.u
    setvortexpositions!(model,u)
    Px, Py = impulse(model)
    push!(Px_numerical_hist,Px)
    push!(Py_numerical_hist,Py)
end

plot(circle,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2,xlabel="x",ylabel="y")
scatter!(model.vortices.x,model.vortices.y,color=:red)
plot!(x_trajectory,y_trajectory_upper,linecolor=:red,label="exact")
plot!(x_trajectory,y_trajectory_lower,linecolor=:red,label="")
plot!(map(s->s.u[1],sol.u),map(s->s.v[1],sol.u),color=:blue,linestyle=:dash,label="simulated")
plot!(map(s->s.u[2],sol.u),map(s->s.v[2],sol.u),color=:blue,linestyle=:dash,label="")

plot(x_trajectory,Px_exact,color=:red,label="exact",xlabel="x",ylabel="Px")
plot!(map(s->s.u[1],sol.u),Px_numerical_hist,color=:blue,linestyle=:dash,label="simulated")

c = 1.0;   # chord length

ċ = -1.0;  # translational velocity
α = -π/3;   # angle of attack

σLE = 0.0;
σTE = 0.0;

Δt = 5e-2;
tf = 1.0;
T = 0.0:Δt:tf;

using PotentialFlow
function compute_ẋ!(ẋ, x, t)
    plate, ambient_sys = x
    motion = ẋ[1]
    motion.ċ, motion.c̈, motion.α̇ = motion.kin(t)
    Plates.enforce_no_flow_through!(plate, motion, ambient_sys, t)
    reset_velocity!(ẋ, x)
    self_induce_velocity!(ẋ, x, t)
end
function shed_new_vorticity!(blobs, plate, motion, t, lesp = 0.0, tesp = 0.0)
    z₊ = (blobs[end-1].z + 2plate.zs[end])/3
    z₋ = (blobs[end].z + 2plate.zs[1])/3
    blob₊ = PotentialFlow.Vortex.Blob(z₊, 1.0, δ)
    blob₋ = PotentialFlow.Vortex.Blob(z₋, 1.0, δ)
    Plates.enforce_no_flow_through!(plate, motion, blobs, t)
    Γ₊, Γ₋, _, _ = Plates.vorticity_flux!(plate, blob₊, blob₋, t, lesp, tesp);
    push!(blobs, PotentialFlow.Vortex.Blob(z₊, Γ₊, blobs[1].δ), PotentialFlow.Vortex.Blob(z₋, Γ₋, blobs[1].δ))
end
N = 128 # number of plate control points (distributed along a extrema Chebyshev grid)
δ = 0.01
plate = PotentialFlow.Plate(N, c, zero(ComplexF64), α)
motion = Plates.RigidBodyMotion(ċ, 0.0);

Δz₀ = im*3Δt*exp(im*plate.α) # vectors perpendicular to the plate
z₋, z₊ = plate.zs[[1,N]] # LE and TE
blobs = PotentialFlow.Vortex.Blob.(Δz₀ .+ [z₊, z₋], 1.0, δ) # First two point vortices are placed close to the LE and TE with unit strength
Plates.enforce_no_flow_through!(plate, motion, (), 0)
Γ₊, Γ₋, _, _ = Plates.vorticity_flux!(plate, blobs[1], blobs[2], 0.0, σLE, σTE); # Determine strength of first two vortices
blobs = PotentialFlow.Vortex.Blob.(Δz₀ .+ [z₊, z₋], [Γ₊, Γ₋], δ) # Create first two point vortices now with calculated strengths
sys₀ = (plate, blobs)
sys = deepcopy(sys₀)
sys₊ = deepcopy(sys₀) # Used for storage during time-marching
ẋs = (motion, allocate_velocity(blobs))
imp = ComplexF64[];
t_hist = Float64[]
global t = 0
push!(t_hist,t)

for tloc in T[2:end]
    global t += Δt
    global sys
    global sys₊
    push!(imp,Elements.impulse(sys))
    push!(t_hist,t)
    local plate, ambient_ω = sys
    local motion, ambient_u = ẋs
    resize!(sys₊[2], length(sys[2]))
    forward_euler!(sys₊, sys, t, Δt, compute_ẋ!, advect!, ẋs)
    sys, sys₊ = sys₊, sys
    shed_new_vorticity!(sys[2], sys[1], ẋs[1], t, σLE, σTE)
end
force = -diff(imp)/Δt;

Δx = 0.01
xlim = (-0.5,2)
ylim = (-1,1)
g = PhysicalGrid(xlim,ylim,Δx);

dsdx = 2
plate = RigidBodyTools.Plate(c,dsdx*Δx)
pfb = PotentialFlowBody(plate,edges=[1,length(plate)],σ=[SuctionParameter(σLE),SuctionParameter(σTE)])
transform = RigidTransform((0.0,0.0),-π/3)
transform(plate);
Δs = dlength(plate)
maximum(Δs/Δx);

firstvLE = GridPotentialFlow.Vortex(plate.x[1]+3Δt*(-ċ)*cos(plate.α+π/2),plate.y[1]+3Δt*(-ċ)*sin(plate.α+π/2),0.0);
firstvTE = GridPotentialFlow.Vortex(plate.x[end]+3Δt*(-ċ)*cos(plate.α+π/2),plate.y[end]+3Δt*(-ċ)*sin(plate.α+π/2),0.0);

function createsheddedvortices(plate,oldvortices)

    vLE = GridPotentialFlow.Vortex(2/3*plate.x[1]+1/3*oldvortices[end-1].x,2/3*plate.y[1]+1/3*oldvortices[end-1].y,0.0)
    vTE = GridPotentialFlow.Vortex(2/3*plate.x[end]+1/3*oldvortices[end].x,2/3*plate.y[end]+1/3*oldvortices[end].y,0.0)

    return vLE,vTE
end

model = VortexModel(g,bodies=[pfb],vortices=[firstvLE,firstvTE],U∞=(-ċ,0.0));

Px_hist = Float64[];
Py_hist = Float64[];
sol = solve(model);
for tloc in T[2:end]
    X = getvortexpositions(model) # gets bigger every time step because we add vortices
    Ẋ = deepcopy(X)

    solve!(sol, model)
    setvortexstrengths!(model, sol.δΓ_vec, length(X.u)-1:length(X.u))
    subtractcirculation!(model.bodies, sol.δΓ_vec)
    Px, Py = impulse(model)

    vortexvelocities!(Ẋ, model, sol.ψ)
    X .= X .+ Ẋ*Δt
    setvortexpositions!(model, X)

    vLE, vTE = createsheddedvortices(plate,model.vortices)
    pushvortices!(model,vLE,vTE)
    push!(Px_hist,Px)
    push!(Py_hist,Py)
end

Fx_hist = -diff(Px_hist)/Δt;
Fy_hist = -diff(Py_hist)/Δt;

plot(plate,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2,xlim=xlim,ylim=ylim,xlabel="x",ylabel="y")
scatter!(real.((v->v.z).(sys[2])).-tf*ċ,imag.((v->v.z).(sys[2])),color=:red,markersize=4,label="PotentialFlow.jl")
scatter!(model.vortices.x,model.vortices.y,color=:blue,markersize=2,label="GridPotentialFlow.jl")

plot(xlabel="t",ylabel="Py")
plot!(T[1:end-1],imag.(imp),color=:blue,label="PotentialFlow.jl")
plot!(T[1:end-1],Py_hist,color=:red,label="GridPotentialFlow.jl")

plot(xlabel="t",ylabel="Fy")
plot!(T[2:end-1],imag.(force),color=:blue,label="PotentialFlow.jl")
plot!(T[2:end-1],Fy_hist,color=:red,label="GridPotentialFlow.jl")

Lx = 2.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,0.01)
Δx = cellsize(g);

a = 0.5
b = 0.25
ellipse = Ellipse(a,b,Δx)
pfb = PotentialFlowBody(ellipse)
model = VortexModel(g,bodies=[pfb])
M = GridPotentialFlow.addedmass(model)

using Test
@test isapprox(M[1,1], π*b^2, rtol=1e-1)

@test isapprox(M[2,2], π*a^2, atol=1e-1)

GRratios = [0.1,0.2,0.3,0.4,0.5,0.6,0.8,1.0,1.2,1.4,1.6];

function rectangulararray(rows,columns,spacing)
     centers = VectorData(rows*columns)
     for j in 1:columns
        for i in 1:rows
            xc = (i-1)*spacing
            yc = (j-1)*spacing
            centers.u[(j-1)*rows+i] = xc
            centers.v[(j-1)*rows+i] = yc
        end
    end
    centers.u .-= mean(centers.u)
    centers.v .-= mean(centers.v)
    return centers
 end

n = 100
rows = 3
columns = 3
N = rows*columns
R = 0.20
𝒱 = π*R^2
bodies = fill(Circle(R,2*cellsize(g)),N);
Δs = minimum(dlength(bodies[1]));

using Statistics: mean
using LinearAlgebra: eigen

λoverMratios = zeros(length(GRratios))
for idx in 1:length(GRratios)
    global bodies
    gap = GRratios[idx]*R
    spacing = 2*R+gap
    bodycenters = rectangulararray(rows,columns,spacing)
    for i in 1:N
        Tf = RigidTransform((bodycenters.u[i],bodycenters.v[i]),0.0)
        global bodies[i] = Tf(deepcopy(bodies[i]))
    end
    body_list = PotentialFlowBody.(bodies)
    vm = VortexModel(g,bodies=body_list);
    M = GridPotentialFlow.addedmass(vm)/𝒱
    eigM = eigen(M);
    max_eig_value_coef = maximum(real(eigM.values))
    max_self_added_mass_coef = maximum(M)
    λoverMratios[idx] = max_eig_value_coef/max_self_added_mass_coef
end

nineRodArrayChen = [2.3637,2.2092,2.1007,2.0120,1.9350,1.8665,1.7494,1.6531,1.5732,1.5066,1.4508];
plot(xlabel="GR",ylabel="λ/max(Mij)")
plot!(GRratios,nineRodArrayChen,label="Chen1975")
plot!(GRratios,λoverMratios,label="GridPotentialFlow.jl")

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

# This is the builder for the library Dierckx in Netlib
# (http://www.netlib.org/dierckx/), which doesn't really have a versin number.
name = "Dierckx"
version = v"0.0.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/scipy/scipy/releases/download/v1.5.1/scipy-1.5.1.tar.xz",
                  "0728bd66a5251cfeff17a72280ae5a40ec14add217f94868d1415b3c469b610a"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/scipy-*/scipy/interpolate/fitpack
mkdir -p "${libdir}"
gfortran -o "${libdir}/libddierckx.${dlext}" -O3 -shared -fPIC *.f
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_gfortran_versions(supported_platforms())

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
]

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

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
"""
    AbstractGlobalGeometry

Determines the conversion from local coordinates and vector bases to a Cartesian basis.
"""
abstract type AbstractGlobalGeometry end

Cartesian123Point(pt::AbstractPoint, global_geometry::AbstractGlobalGeometry) =
    Cartesian123Point(CartesianPoint(pt, global_geometry))

(::Type{<:CartesianVector{<:Any, I}})(
    u::AxisVector,
    global_geometry::AbstractGlobalGeometry,
    local_geometry::LocalGeometry,
) where {I} = transform(
    CartesianAxis{I}(),
    CartesianVector(u, global_geometry, local_geometry),
)

"""
    CartesianGlobalGeometry()

Specifies that the local coordinates align with the Cartesian coordinates, e.g.
`XYZPoint` aligns with `Cartesian123Point`, and `UVWVector` aligns with
`Cartesian123Vector`.
"""
struct CartesianGlobalGeometry <: AbstractGlobalGeometry end

# coordinates
CartesianPoint(pt::XPoint{FT}, ::CartesianGlobalGeometry) where {FT} =
    Cartesian1Point{FT}(pt.x)
CartesianPoint(pt::ZPoint{FT}, ::CartesianGlobalGeometry) where {FT} =
    Cartesian3Point{FT}(pt.z)
CartesianPoint(pt::XYPoint{FT}, ::CartesianGlobalGeometry) where {FT} =
    Cartesian12Point{FT}(pt.x, pt.y)
CartesianPoint(pt::XZPoint{FT}, ::CartesianGlobalGeometry) where {FT} =
    Cartesian13Point{FT}(pt.x, pt.z)
CartesianPoint(pt::XYZPoint{FT}, ::CartesianGlobalGeometry) where {FT} =
    Cartesian123Point{FT}(pt.x, pt.y, pt.z)

# vectors
CartesianVector(
    u::CartesianVector,
    ::CartesianGlobalGeometry,
    ::LocalGeometry,
) = u
function CartesianVector(
    u::AxisVector,
    ::CartesianGlobalGeometry,
    local_geometry::LocalGeometry{I},
) where {I}
    u_local = LocalVector(u, local_geometry)
    AxisVector(CartesianAxis{I}(), components(u_local))
end

#=
LocalVector(u::CartesianVector{T,I}, ::CartesianGlobalGeometry) where {T,I} =
    AxisVector(LocalAxis{I}(), components(u))
=#

"""
    SphericalGlobalGeometry(radius)

Specifies that the local coordinates are specified in reference to a sphere of
radius `radius`. The `x1` axis is aligned with the zero longitude line.

The local vector basis has `u` in the zonal direction (with east being
positive), `v` in the meridonal (north positive), and `w` in the radial
direction (outward positive). For a point located at the pole, we take the limit
along the zero longitude line:
- at the north pole, this corresponds to `u` being aligned with the `x2`
  direction, `v` being aligned with the negative `x1` direction, and `w` being
  aligned with the `x3` direction.
- at the south pole, this corresponds to `u` being aligned with the `x2`
  direction, `v` being aligned with the `x1` direction, and `w` being aligned
  with the negative `x3` direction.
"""
struct SphericalGlobalGeometry{FT} <: AbstractGlobalGeometry
    radius::FT
end

# coordinates
function CartesianPoint(pt::LatLongPoint, global_geom::SphericalGlobalGeometry)
    r = global_geom.radius
    x1 = r * cosd(pt.long) * cosd(pt.lat)
    x2 = r * sind(pt.long) * cosd(pt.lat)
    x3 = r * sind(pt.lat)
    Cartesian123Point(x1, x2, x3)
end
function LatLongPoint(pt::Cartesian123Point, ::SphericalGlobalGeometry)
    ϕ = atand(pt.x3, hypot(pt.x2, pt.x1))
    # IEEE754 spec states that atand(±0.0, −0.0) == ±180, however to make the UV
    # orienation consistent, we define the longitude to be zero at the poles
    if abs(ϕ) == 90
        λ = zero(ϕ)
    else
        λ = atand(pt.x2, pt.x1)
    end
    LatLongPoint(ϕ, λ)
end

"""
    great_circle_distance(pt1::LatLongPoint, pt2::LatLongPoint, global_geometry::SphericalGlobalGeometry)

Compute the great circle (spherical geodesic) distance between `pt1` and `pt2`.
"""
function great_circle_distance(
    pt1::LatLongPoint,
    pt2::LatLongPoint,
    global_geom::SphericalGlobalGeometry,
)
    r = global_geom.radius
    ϕ1 = pt1.lat
    λ1 = pt1.long
    ϕ2 = pt2.lat
    λ2 = pt2.long
    Δλ = λ1 - λ2
    return r * atan(
        hypot(
            cosd(ϕ2) * sind(Δλ),
            cosd(ϕ1) * sind(ϕ2) - sind(ϕ1) * cosd(ϕ2) * cosd(Δλ),
        ),
        cosd(ϕ1) * cosd(ϕ2) * cosd(Δλ) + sind(ϕ1) * sind(ϕ2),
    )
end

# vectors
CartesianVector(
    u::CartesianVector,
    ::SphericalGlobalGeometry,
    ::LocalGeometry,
) = u
CartesianVector(
    u::AxisVector,
    global_geometry::SphericalGlobalGeometry,
    local_geometry::LocalGeometry,
) = CartesianVector(
    UVWVector(u, local_geometry),
    global_geometry,
    local_geometry.coordinates,
)
function CartesianVector(
    u::UVWVector,
    ::SphericalGlobalGeometry,
    coord::LatLongPoint,
)
    ϕ = coord.lat
    λ = coord.long
    G = @SMatrix [
        -sind(λ) -sind(ϕ)*cosd(λ) cosd(λ)*cosd(ϕ)
        cosd(λ) -sind(ϕ)*sind(λ) sind(λ)*cosd(ϕ)
        0 cosd(ϕ) sind(ϕ)
    ]
    AxisVector(Cartesian123Axis(), G * Geometry.components(u))
end


text:
"""
    eigsROT_Distance(P::Matrix{Float64}, W::SparseMatrixCSC{Float64, Int64}, X::Matrix{Float64}; α::Float64 = 1.0)

computes the sROT distance matrix of P's column vectors on a (unweighted) tree.

# Input Argument
- `P::Matrix{Float64}`: a matrix whose columns are vector measures with the same total mass.
- `W::SparseMatrixCSC{Float64, Int64}`: the weight matrix of the tree.
- `X::Matrix{Float64}`: the node positions (i-th row represents node `i`'s location)
- `α::Float64`: default is 1.0. ROT parameter.

# Output Argument
- `dist_sROT::Matrix{Float64}`: distance matrix, dist_sROT[i,j] = d_{sROT}(pᵢ, pⱼ; α).
- `Ws::SparseMatrixCSC{Float64, Int64}`: the weight matrix of the simplified tree.
- `Xs::Matrix{Float64}`: the node locations of the simplified tree
- `𝚯::Matrix{Float64}`: the shortened pmfs from `P`.

"""
function eigsROT_Distance(P::Matrix{Float64}, W::SparseMatrixCSC{Float64, Int64}, X::Matrix{Float64}; α::Float64 = 1.0)
    G = SimpleGraph(W)
    # check if the input weight matrix forms a tree
    if ne(G) + 1 != nv(G) || !is_connected(G)
        @error("input graph is not a tree.")
    end
    d = degree(G)
    # find index of junction nodes
    ijc = findall(d .> 2)
    # cut the graph into several disconnected subgraphs
    Wc = deepcopy(W)
    for i in ijc; Wc[i, :] .= 0; Wc[:, i] .= 0; end
    # find the indices for each subgraph
    Ind = find_subgraph_inds(Wc)
    # low dimensional pmfs
    𝚯 = Ind' * P
    # the centroids of subgraphs
    Xs = Diagonal(1 ./ sum(Ind, dims = 1)[:]) * Ind' * X
    # build Gs, i.e., the graph of subgraphs or the simplified tree
    Gs = Graph(size(Ind, 2))
    Ns = nv(Gs)
    # index of Gs's junction nodes
    ijcs = []
    for k = 1:Ns
        supportind = findall(Ind[:, k] .== 1)
        if issubset(supportind, ijc)
            push!(ijcs, k)
        end
    end
    # connect the branches to the junctions to form the simplified tree
    for k in setdiff(1:Ns, ijcs)
        supportind = findall(Ind[:,k] .== 1)
        for i = 1:length(ijc)
            if sum(W[supportind, ijc[i]]) > 0
                add_edge!(Gs, Edge(k, ijcs[i]))
            end
        end
    end
    Ws = 1.0 * adjacency_matrix(Gs)
    # compute the sROT distance matrix
    Qs = incidence_matrix(Gs; oriented = true)
    dist_sROT = eigROT_Distance(𝚯, Qs; edge_length = 1, α = α)
    return dist_sROT, Ws, Xs, 𝚯
end

"""
    find_subgraph_inds(Wc::SparseMatrixCSC{Float64, Int64})

find all subgraph indices of a tree. (subgraph includes: branches and junctions)

# Input Argument
- `Wc::SparseMatrixCSC{Float64, Int64}`: the weight matrix of the tree chopped by
    the junctions.

# Output Argument
- `Ind::Matrix{Float64}`: a matrix whose columns represent the subgraph node
    indices in binary form.

"""
function find_subgraph_inds(Wc::SparseMatrixCSC{Float64, Int64})
    N = size(Wc, 1)
    Dc = Diagonal(sum(Wc, dims=1)[:])
    Lc = Matrix(Dc - Wc)
    𝛌c, 𝚽c = eigen(Lc)
    p = findall(𝛌c .< 100 * eps())
    Uc = 𝚽c[:, p]
    # use heat diffusion at t = ∞ and the siginal `1:N` to figure out all branch indices
    a = round.(1e8 * Uc * Uc' * [k for k = 1:N])
    ind = unique(a)
    Ind = zeros(N, length(ind))
    for k = 1:length(ind); Ind[findall(a .== ind[k]), k] .= 1; end
    return Ind
end


text:
include("file_io/read_input.jl")
include("file_io/read_key.jl")
include("file_io/write_output.jl")
include("file_io/write_key.jl")

include("initialize_system.jl")
include("evolve_system.jl")
include("compress_system.jl")

include("generate_random_unit_quaternion.jl")
include("characteristic_matrix_ellipsoid.jl")
include("characteristic_matrix_ellipse.jl")
include("rotation_matrix.jl")

include("generate_proposal_position.jl")
include("generate_proposal_orientation.jl")
include("overlap_ellipsoid.jl")
include("overlap_ellipse.jl")

include("overlap_cuboid.jl")
include("overlap_cuboid_binary.jl")

include("signed_distance_mod.jl")
include("quaternion_mult.jl")

include("vtk/voxel_structure.jl")
include("vtk/write_vtk.jl")

function run_generation()
	# Inititalization of random number generation device.
	random_seed::Int64 = convert(Int64, time_ns())
	srand(random_seed)

	# Start time.
	t_start_ns::Int64 = convert(Int64, time_ns())

	# Change current folder to the folder where this script lies.
	(program_file_dir::String, program_file_name::String) = splitdir(PROGRAM_FILE)
	program_file_dir = abspath(program_file_dir)
	cd(program_file_dir)

	# Assert that input is file and store path.
	input_file_path::String = ""
	if isfile(ARGS[1])
		input_file_path = ARGS[1]
	else
		println("No input file specified or specified input file does not exist. Aborting.")
		return nothing
	end

	# Read input from file.
	println(join(("Reading input from file ", input_file_path, "...")))
	(	particle_type::String,
		R::Array{Float64, 2},
		Lx::Float64,
		Ly::Float64,
		Lz::Float64,
		phi::Float64,
		phi_target::Float64,
		position_constraint_axis::String,
		position_constraint_lower::Float64,
		position_constraint_upper::Float64,
		orientation_axis::Array{Float64, 1},
		orientation_constraint_axis::Array{Float64, 1},
		orientation_constraint_lower::Float64,
		orientation_constraint_upper::Float64,
		sigma_translation_max::Float64,
		sigma_rotation_max::Float64,
		sigma_ratio::Float64,
		number_of_equilibration_sweeps::Int64,
		delta_phi::Float64,
		output_file_path::String) = read_input(input_file_path)

	# Initialize system.
	(	Lx,
		Ly,
		Lz,
		phi,
		X::Array{Float64, 1},
		Y::Array{Float64, 1},
		Z::Array{Float64, 1},
		Q0::Array{Float64, 1},
		Q1::Array{Float64, 1},
		Q2::Array{Float64, 1},
		Q3::Array{Float64, 1},
		A11::Array{Float64, 1},
		A12::Array{Float64, 1},
		A13::Array{Float64, 1},
		A21::Array{Float64, 1},
		A22::Array{Float64, 1},
		A23::Array{Float64, 1},
		A31::Array{Float64, 1},
		A32::Array{Float64, 1},
		A33::Array{Float64, 1}) = initialize_system(	particle_type,
														R,
														Lx,
														Ly,
														Lz,
														phi,
														position_constraint_axis,
														position_constraint_lower,
														position_constraint_upper,
														orientation_axis,
														orientation_constraint_axis,
														orientation_constraint_lower,
														orientation_constraint_upper)

	# Translation and rotation speeds.
	sigma_translation::Float64 = sigma_translation_max
	sigma_rotation::Float64 = 0.0
	if sigma_ratio >= 0.0
		sigma_rotation = sigma_ratio * sigma_translation
	else
		sigma_rotation = sigma_rotation_max
	end

	# Relax and equlibrate system.
	number_of_relaxation_sweeps_max::Int64 = typemax(Int64)
	(	X,
		Y,
		Z,
		Q0,
		Q1,
		Q2,
		Q3,
		A11,
		A12,
		A13,
		A21,
		A22,
		A23,
		A31,
		A32,
		A33,
		sigma_translation,
		sigma_rotation,
		is_relaxed::Bool) = evolve_system(	particle_type,
										R,
										Lx,
										Ly,
										Lz,
										X,
										Y,
										Z,
										Q0,
										Q1,
										Q2,
										Q3,
										A11,
										A12,
										A13,
										A21,
										A22,
										A23,
										A31,
										A32,
										A33,
										position_constraint_axis,
										position_constraint_lower,
										position_constraint_upper,
										orientation_axis,
										orientation_constraint_axis,
										orientation_constraint_lower,
										orientation_constraint_upper,
										sigma_translation,
										sigma_translation_max,
										sigma_rotation,
										sigma_rotation_max,
										sigma_ratio,
										number_of_relaxation_sweeps_max,
										number_of_equilibration_sweeps)

	# Compress system.
	number_of_relaxation_sweeps_max = 1000
	if particle_type != "ellipse"
		if phi_target > phi
			(	Lx,
				Ly,
				Lz,
				phi,
				X,
				Y,
				Z,
				Q0,
				Q1,
				Q2,
				Q3,
				A11,
				A12,
				A13,
				A21,
				A22,
				A23,
				A31,
				A32,
				A33,
				sigma_translation,
				sigma_rotation) = compress_system(	particle_type,
												R,
												Lx,
												Ly,
												Lz,
												phi_target,
												X,
												Y,
												Z,
												Q0,
												Q1,
												Q2,
												Q3,
												A11,
												A12,
												A13,
												A21,
												A22,
												A23,
												A31,
												A32,
												A33,
												position_constraint_axis,
												position_constraint_lower,
												position_constraint_upper,
												orientation_axis,
												orientation_constraint_axis,
												orientation_constraint_lower,
												orientation_constraint_upper,
												sigma_translation,
												sigma_translation_max,
												sigma_rotation,
												sigma_rotation_max,
												sigma_ratio,
												delta_phi,
												number_of_relaxation_sweeps_max)
		else
			println("Target volume fraction not larger than initial volume fraction. Not compressing system.")
		end
	else
		if phi_target > phi
			println("System with particle type 'ellipse' not eligible for compression. Not compressing system.")
		end
	end

	# Write output.
	t_exec_ns::Int64 = convert(Int64, time_ns()) - t_start_ns
	t_exec_s::Float64 = convert(Float64, t_exec_ns) / 1e9
	write_output(	output_file_path,
				particle_type,
				R,
				Lx,
				Ly,
				Lz,
				phi,
				X,
				Y,
				Z,
				Q0,
				Q1,
				Q2,
				Q3,
				t_exec_s)

	nothing
end

run_generation()

max_stars_repo_name: roding/whitefish

max_issues_repo_name: roding/whitefish

max_forks_repo_name: roding/whitefish


text:
import Base: elsize, size, ndims, convert, showarray, vec, flipdim, vcat, hcat, cat, reshape, permutedims, circshift, repeat

export AFInfo, replace!, mergeArrays, getDevicePointer, getDataRefCount, deviceMemInfo

immutable Dim4
    dim1::Integer
    dim2::Integer
    dim3::Integer
    dim4::Integer
end

new_ptr() = Base.RefValue{Ptr{Void}}(C_NULL)

sizeof{T}(a::AFArray{T}) = elsize(a) * length(a)

function convert{T,N}(::Type{Array{T,N}}, x::AFArray{T,N})
    ret = Array(UInt8, sizeof(x))
    err = ccall((:af_get_data_ptr, af_lib),
                Cint, (Ptr{T}, Ptr{Ptr{Void}}),
                pointer(ret), x.ptr)
    err == 0 || throwAFerror(err)
    #af_get_data_ptr!(ret, x, T)
    ret = reinterpret(T, ret)
    ret = reshape(ret, size(x)...)
    ret
end

@compat (::Type{Array}){T,N}(a::AFArray{T,N}) = convert(Array{T,N}, a)

function size(a::AFArray)
    dim1 = Base.RefValue{Cuint}(0)
    dim2 = Base.RefValue{Cuint}(0)
    dim3 = Base.RefValue{Cuint}(0)
    dim4 =Base.RefValue{Cuint}(0)
    af_get_dims!(dim1, dim2, dim3, dim4, a)
    n = ndims(a)
    dim4_to_dims(Dim4(dim1[], dim2[], dim3[], dim4[]), n)
end

function dim4_to_dims(d::Dim4, n::Integer)
    if n == 1
        return (Int(d.dim1),)
    elseif n == 2
        return (Int(d.dim1), Int(d.dim2))
    elseif n == 3
        return (Int(d.dim1), Int(d.dim2), Int(d.dim3))
    elseif n == 4
        return (Int(d.dim1), Int(d.dim2), Int(d.dim3), Int(d.dim4))
    else
        throw(ArgumentError("Too many dimensions for an ArrayFire Array"))
    end
end

function get_all_dims(a::Array)
    n = ndims(a)
    dims = [size(a)...]
    for i = n+1:4
        push!(dims, 1)
    end
    dims
end

function ndims(a::AFArray)
    n = Base.RefValue{Cuint}(0)
    af_get_numdims!(n, a.ptr)
    Int(n[])
end

function ndims(ptr::Ptr{Void})
    n = Base.RefValue{Cuint}(0)
    af_get_numdims!(n, ptr)
    Int(n[])
end

@compat (::Type{AFArray{T}}){T}(ptr::Ptr{Void}) = AFArray{T, ndims(ptr)}(ptr)

@compat function Base.show(io::IO, m::MIME"text/plain", a::AFArray)
    print(io, summary(a))
    !isempty(a) && println(io,":")
    if VERSION > v"0.5.0-dev+4993"
        Base.showarray(io::IO, Array(a), false; header = false)
    else
        Base.showarray(io::IO, Array(a); header=false)
    end
end

#=function showarray{T}(io::IO, X::AFAbstractArray{T};
                   header::Bool=true, kwargs...)
    header && print(io, summary(X))
    !isempty(X) && println(io,":")
    showarray(io, convert(Array{T},X); header = false, kwargs...)
end=#

function backend_eltype(p::Ptr{Void})
    t = Base.Ref{Cuint}(0)
    af_get_type(t, p)
    jltype(t[])
end

backend_eltype(a::AFArray) = backend_eltype(a.ptr)

function AFInfo()
    af_info()
    nothing
end

# Array Utils

function vec{T}(a::AFArray{T})
    out = new_ptr()
    af_flat(out, a)
    AFArray{T}(out[])
end

function flipdim{T}(a::AFArray{T}, dim::Int)
    out = new_ptr()
    af_flip(out, a, Cuint(dim-1))
    AFArray{T}(out[])
end

function vcat{T,S}(a::AFArray{T}, b::AFArray{S})
    out = new_ptr()
    af_join(out, 0, a, b)
    AFArray{af_promote(T,S)}(out[])
end

function hcat{T,S}(a::AFArray{T}, b::AFArray{S})
    out = new_ptr()
    af_join(out, 1, a, b)
    AFArray{af_promote(T,S)}(out[])
end

function cat(dim::Int, a::AFArray...)
    out = new_ptr()
    n = length(a)
    b = [a[i].ptr for i = 1:n]
    af_join_many(out, dim - 1, Cuint(n), b)
    AFArray{backend_eltype(out[])}(out[])
end

function reshape{T}(a::AFArray{T}, dims::Int...)
    out = new_ptr()
    l = length(dims)
    dims = [dims...]
    af_moddims(out, a, Cuint(l), dims)
    AFArray{T}(out[])
end

function permutedims{T}(a::AFArray{T}, perm::Vector{Int})
    make4Dperm!(perm)
    out = new_ptr()
    perm = perm - 1
    af_reorder(out, a, Cuint(perm[1]), Cuint(perm[2]), Cuint(perm[3]), Cuint(perm[4]))
    AFArray{T}(out[])
end

function make4Dperm!(perm::Vector{Int})
    l = length(perm)
    for i = l+1:4
        push!(perm, i)
    end
end

function replace!(a::AFArray, b::AFArray, cond::AFArray{Bool})
    af_replace(a, cond, b)
    nothing
end

function replace!(a::AFArray, b::Real, cond::AFArray{Bool})
    af_replace_scalar(a, cond, b)
    nothing
end

function mergeArrays{T,S}(a::AFArray{T}, b::AFArray{S}, cond::AFArray{Bool})
    out = new_ptr()
    af_select(out, cond, a, b)
    AFArray{af_promote(T,S)}(out[])
end

function mergeArrays{T<:Number,S}(a::T, b::AFArray{S}, cond::AFArray{Bool})
    out = new_ptr()
    af_select_scalar_l(out, cond, a, b)
    AFArray{af_promote(T,S)}(out[])
end

function mergeArrays{T,S<:Number}(a::AFArray{T}, b::S, cond::AFArray{Bool})
    out = new_ptr()
    af_select_scalar_r(out, cond, a, b)
    AFArray{af_promote(T,S)}(out[])
end

function circshift{T}(a::AFArray{T}, shifts::Vector{Int})
    make4Dshift!(shifts)
    out = new_ptr()
    af_shift(out, a, shifts...)
    AFArray{T}(out[])
end

function make4Dshift!(s::Vector{Int})
    l = length(s)
    for i = l+1:4
        push!(s, 0)
    end
end

function repeat{T}(a::AFArray{T}; outer::Vector{Int} = [1,1,1,1])
    make4Dtile!(outer)
    out = new_ptr()
    outer = map(Cuint, outer)
    af_tile(out, a, outer...)
    AFArray{T}(out[])
end

function make4Dtile!(d::Vector{Int})
    l = length(d)
    for i = l+1:4
        push!(d, 1)
    end
end

function convert{T}(::Type{AFArray{T}}, a::AFArray)
    out =  new_ptr()
    af_cast(out, a, T)
    AFArray{T}(out[])
end

function getDevicePointer(a::AFArray)
    out = new_ptr()
    af_get_device_ptr(out, a)
    out[]
end

function getDataRefCount(a::AFArray)
    out = Base.Ref{Cint}(0)
    af_get_data_ref_count(out, a)
    out[]
end

function deviceMemInfo()
    alloc_bytes = Base.Ref{Csize_t}(0)
    alloc_buffers = Base.Ref{Csize_t}(0)
    lock_bytes = Base.Ref{Csize_t}(0)
    lock_buffers = Base.Ref{Csize_t}(0)
    af_device_mem_info(alloc_bytes, alloc_buffers, lock_bytes, lock_buffers)
    Int(alloc_bytes[]), Int(alloc_buffers[]), Int(lock_bytes[]), Int(lock_buffers[])
end


text:
struct dmrgIterable{TMPS <: AbstractMPS, TMPO <: AbstractMPO}
    initmps::TMPS
    mpo::TMPO
end

struct dmrgState{TMPS <: AbstractMPS, TE <: AbstractTensor}
    mps::TMPS
    envs::Vector{TE}
    energies::Vector{Float64}
end

function dmrgiterable(mps::TMPS, mpo::TMPO) where {TMPS <: AbstractMPS, TMPO <: AbstractMPO}
    cmps = normalize!(canonicalize(mps,1))
    return dmrgIterable{typeof(cmps), TMPO}(cmps , mpo)
end

function iterate(iter::dmrgIterable{TMPS}) where TMPS
    @unpack initmps, mpo = iter
    envs = initenvironments(initmps, mpo)
    TE = typeof(envs[1])
    state = dmrgState{TMPS,TE}(initmps, envs, [real(inner(initmps,mpo,initmps))])
    return state, state
end

function iterate(iter::dmrgIterable, state::dmrgState)
    @unpack mps, envs, energies = state
    @unpack mpo = iter
    E = dmrgsweep!(mps, mpo, envs)
    push!(energies, real(E))
    return state, state
end

function initenvironments(mps::CanonicalMPS{L,T,TE,TB}, mpo::MPO{L,T}) where {L,T,TE,TB}
    centerlink(mps) == 1 || throw(ArgumentError("centerlink needs to be set to 1"))
    envs = TB[]
    foldr((x,y) -> push!(envs, contractsites(x,y))[end],
        collect(zip(mps[2:L], mpo[2:L], mps[2:L]')), init = ())
    return reverse(envs)
end

function updateenvironments!(envs::Vector{TB},
     mps::CanonicalMPS{L,T,TE,TB}, mpo, i, up = true) where {L,T,TE,TB}
    centerlink(mps) == i || centerlink(mps) + 1 == i || throw(
        ArgumentError("illegal environment update"))
    if i == 1 #down -> up
        envs[1]   = contractsites((), (mps[1],mpo[1],mps[1]'))
    elseif i == L #up -> down
        envs[L-1] = contractsites((mps[L],mpo[L],mps[L]'),())
    elseif up
        envs[i]   = contractsites(envs[i-1], (mps[i], mpo[i], mps[i]'))
    else
        envs[i-1] = contractsites((mps[i], mpo[i], mps[i]'), envs[i])
    end
    return envs
end

"""
    dmrg(ansatz::AbstractMPS, mpo::AbstractMPO [; tol=1e-6, verbose=false,minit=2,period=1,maxit=10^6])
returns the ground-state as a `CanonicalMPS` of the MPO `mpo` as well as a list of the energies
calculated after each full sweep of the single-site dmrg algorithm.
The algorithm finishes either if `maxit` iterations are reached, or at least `minit` iterations
have been done _and_ the absolute change of energy between the two most recent steps is less
than `tol`.
If `tol` is negative, `maxit` iterations will be done.

If `verbose=true`, information about the execution will be printed every `period` iteration.
"""
function dmrg(ansatz::AbstractMPS{L,T}, mpo::AbstractMPO{L,T};
                tol::Float64 = 1e-6,
                verbose::Bool = false,
                minit::Int = 2,
                period::Int = 1,
                maxit::Int = 1_000_000,) where {L,T}

    stop(state) = length(state.energies) > minit &&
                    abs(state.energies[end] - state.energies[end-1]) < tol
    disp(state) = @printf("%5d \t| %.3e | %.3e\n", state[2][1]-1, state[1]/1e9, state[2][2].energies[end])

    iter = dmrgiterable(ansatz, mpo)
    tol > 0 && (iter = halt(iter, stop) )
    iter = take(iter, maxit)
    iter = enumerate(iter)

    if verbose
        @printf("\tn \t| time (s)\t| E \n")
        iter = sample(iter, period)
        iter = stopwatch(iter)
        iter = tee(iter, disp)
        (_, (n, state)) = loop(iter)
    else
        (n, state) = loop(iter)
    end
    return (state.mps, state.energies[2:end])
end

function dmrgsweep!(mps::AbstractMPS{L}, mpo, envs) where L
    E = 0
    for i in 1:L
        E = minimisesite!(mps, mpo, i, envs)
        updateenvironments!(envs, mps, mpo, i, true)
    end
    for i in L-1:-1:2
        E = minimisesite!(mps, mpo, i, envs)
        canonicalize!(mps,i-1)
        updateenvironments!(envs, mps, mpo, i, false)
    end
    return E
end

function minimisesite!(mps::CanonicalMPS{L}, mpo::AbstractMPO, i, envs) where {L}
    canonicalize!(mps, min(i, L-1))
    #absorb zero site
    if centerlink(mps) >= i
        mps[i] = contractsite(mps[i], zerosite(mps))
    else
        mps[i] = contractsite(zerosite(mps), mps[i])
    end

    #smallest realpart eigenvector
    mps[i], E, info = minevec(mps, mpo, i, envs)

    #canonical form again
    if i == L
        U, S, mps[L] = onesitesvd2left(mps[L])
        mps.zerosite = contractsite(U, S)
    else
        mps[i], S, V = onesitesvd2right(mps[i])
        mps.zerosite = contractsite(S, V)
    end
    return E
end

function minevec(mps::CanonicalMPS{L}, mpo::AbstractMPO{L}, i, envs) where L
    heffmapfun = heffmapconst(mpo, envs, i)
    v0 = mps[i]
    res = eigsolve(heffmapfun, v0, 1, :SR, ishermitian=true)
    return (res[2][1], res[1][1], res[3]) #evec, eval, info
end


text:
@testset "competing Poisson" begin
    r = [0.2, 0.7, 3.9]
    cp = CompetingPoisson(r; events = 2:4)

    # printing
    @test repr(cp) == "Competing Poisson with total rate 4.8; P(2)=0.0417 P(3)=0.146 P(4)=0.812"

    # properties
    @test cp.rates ≈ r
    @test cp.total_rate ≈ sum(r)
    @test cp.probabilities ≈ normalize(r, 1)
    @test cp.events == 2:4
    @test Set(propertynames(cp)) == Set((:total_rate, :rates, :probabilities, :events))
    @test Set(propertynames(cp, true)) == Set((propertynames(cp)..., :index_distribution))

    # empirical count
    c = zeros(Int, 4)           # event counts
    d = zeros(10000)            # durations
    for i in axes(d, 1)
        x = @inferred rand(cp)
        d[i] = x.duration
        c[x.event] += 1
    end
    e = ecdf(d)
    p = (1:9)./10
    q = quantile.(Exponential(1 / sum(r)), p)
    @test mean(d) ≈ (1 ./ sum(r)) atol = 0.01
    @test all(isapprox.(e.(q), p; atol = 0.02))
    @test normalize(c, 1) ≈ vcat([0], normalize(r, 1)) atol = 0.02
end


text:
"""
    Atmos

Atmospheric parameters.
"""
module Atmos

module SubgridScale

export C_smag,
       C_drag,
       inv_Pr_turb,
       Prandtl_air,
       c_a_KASM,
       c_e1_KASM,
       c_e2_KASM,
       c_1_KASM,
       c_2_KASM,
       c_3_KASM

""" Smagorinsky Coefficient (dimensionless) """
function C_smag end

""" Drag coefficient (dimensionless) """
function C_drag end

""" Turbulent Prandtl Number (dimensionless) """
function inv_Pr_turb end

""" Molecular Prandtl Number, dry air (dimensionless) """
function Prandtl_air end

""" ``c_{a}`` KASM, 2006 (dimensionless) """
function c_a_KASM end

""" ``c_{e1}`` KASM, 2006 (dimensionless) """
function c_e1_KASM end

""" ``c_{e2}`` KASM, 2006 (dimensionless) """
function c_e2_KASM end

""" ``c_{1}``  KASM, 2006 (dimensionless) """
function c_1_KASM end

""" ``c_{2}`` KASM, 2006 (dimensionless) """
function c_2_KASM end

""" ``c_{3}`` KASM, 2006 (dimensionless) """
function c_3_KASM end

end # module SubgridScale

module EDMF

export c_λ,
       c_ε,
       c_δ,
       c_γ,
       β,
       μ_0,
       χ,
       w_min,
       lim_ϵ,
       lim_amp,
       a_min,
       a_surf,
       κ_star²,
       ψϕ_stab,
       α_d,
       α_a,
       α_b,
       H_up_min,
       c_d,
       c_m,
       c_b,
       a1,
       a2,
       ω_pr,
       Pr_n,
       Ri_c,
       smin_ub,
       smin_rm

""" Entrainment TKE scale (dimensionless) """
function c_λ end

""" Entrainment factor (dimensionless) """
function c_ε end

""" Detrainment factor (dimensionless) """
function c_δ end

""" Turbulent Entrainment factor (dimensionless) """
function c_γ end

""" Detrainment RH power (dimensionless) """
function β end

""" Logistic function scale (s``^{-1}``) """
function μ_0 end

""" Updraft mixing fraction (dimensionless) """
function χ end

""" Minimum updraft velocity (m/s) """
function w_min end

""" Exponential area limiter scale (dimensionless) """
function lim_ϵ end

""" Exponential area limiter amplitude(dimensionless) """
function lim_amp end

""" Minimum area fraction for any subdomain (dimensionless) """
function a_min end

""" Updraft surface area fraction (dimensionless) """
function a_surf end

""" Square ratio of rms turbulent velocity to friction velocity (dimensionless) """
function κ_star² end

""" Surface covariance stability coefficient (dimensionless) """
function ψϕ_stab end

""" Pressure drag factor (dimensionless) """
function α_d end

""" Advective pressure factor (dimensionless """
function α_a end

""" Buoyancy pressure factor (dimensionless) """
function α_b end

""" Minimum updraft height for closures (m) """
function H_up_min end

""" Environmental dissipation factor (dimensionless) """
function c_d end

""" Eddy viscosity factor (dimensionless) """
function c_m end

""" Static stability factor (dimensionless) """
function c_b end

""" Empirical stability function coefficient (dimensionless) """
function a1 end

""" Empirical stability function coefficient (dimensionless) """
function a2 end

""" Prandtl number empirical coefficient (dimensionless) """
function ω_pr end

""" Neutral Prandtl number (dimensionless) """
function Pr_n end

""" Critical Richardson number (dimensionless) """
function Ri_c end

""" Smooth minimum's fractional upper bound (dimensionless) """
function smin_ub end

""" Smooth minimum's regularization length minimum (m) """
function smin_rm end

end # module EDMF

module Microphysics_0M

export τ_precip,
       qc_0,
       S_0

""" precipitation removal timescale (s) """
function τ_precip end

""" precipitation removal threshold expressed in condensate specific humidity (kg/kg) """
function qc_0 end

""" precipitation removal threshold expressed in supersaturation (dimensionless) """
function S_0 end

end # module Microphysics_0M

module Microphysics

export C_drag,
       K_therm,
       D_vapor,
       ν_air,
       N_Sc,
       τ_cond_evap,
       τ_sub_dep,
       r_ice_snow,
       n0_ice,
       r0_ice,
       me_ice,
       m0_ice,
       χm_ice,
       Δm_ice,
       q_liq_threshold,
       τ_acnv_rai,
       a_vent_rai,
       b_vent_rai,
       n0_rai,
       r0_rai,
       me_rai,
       ae_rai,
       ve_rai,
       m0_rai,
       a0_rai,
       χm_rai,
       Δm_rai,
       χa_rai,
       Δa_rai,
       χv_rai,
       Δv_rai,
       q_ice_threshold,
       τ_acnv_sno,
       a_vent_sno,
       b_vent_sno,
       μ_sno,
       ν_sno,
       r0_sno,
       me_sno,
       ae_sno,
       ve_sno,
       m0_sno,
       a0_sno,
       v0_sno,
       χm_sno,
       Δm_sno,
       χa_sno,
       Δa_sno,
       χv_sno,
       Δv_sno,
       microph_scaling,
       microph_scaling_dep_sub,
       microph_scaling_melt,
       E_liq_rai,
       E_liq_sno,
       E_ice_rai,
       E_ice_sno,
       E_rai_sno

""" Marshall-Palmer distribution `n_0` coeff for rain or snow (1/m``^4``) """
function n0_ice end
""" Marshall-Palmer distribution `n_0` coeff for rain or snow (1/m``^4``) """
function n0_rai end

""" coefficient to compute Marshall-Palmer distribution coefficient `n_0(ρq_snow/ρ0)` for snow (1/m``^4``) """
function μ_sno end

""" coefficient to compute Marshall-Palmer distribution coefficient `n_0(ρq_snow/ρ0)` for snow (-) """
function ν_sno end

""" assumed length scale for water drop or ice crystal (m) """
function r0_rai end
""" assumed length scale for water drop or ice crystal (m) """
function r0_ice end
""" assumed length scale for water drop or ice crystal (m) """
function r0_sno end

""" coefficient in mass, `radius/r0`, (kg) """
function m0_ice end
""" coefficient in mass, `radius/r0`, (kg) """
function m0_rai end
""" coefficient in mass, `radius/r0`, (kg) """
function m0_sno end

""" coefficient in mass, `radius/r0`, (-) """
function χm_ice end
""" coefficient in mass, `radius/r0`, (-) """
function χm_rai end
""" coefficient in mass, `radius/r0`, (-) """
function χm_sno end

""" exponent in mass, `radius/r0`, (-) """
function me_sno end
""" exponent in mass, `radius/r0`, (-) """
function me_ice end
""" exponent in mass, `radius/r0`, (-) """
function me_rai end

""" exponent in mass, `radius/r0`, (-) """
function Δm_ice end
""" exponent in mass, `radius/r0`, (-) """
function Δm_rai end
""" exponent in mass, `radius/r0`, (-) """
function Δm_sno end

""" coefficient in area, `radius/r0`, relation (m``^2``) """
function a0_rai end
""" coefficient in area, `radius/r0`, relation (m``^2``) """
function a0_sno end

""" coefficient in area, `radius/r0`, (-) """
function χa_rai end
""" coefficient in area, `radius/r0`, (-) """
function χa_sno end

""" exponent in area, `radius/r0`, relation """
function ae_rai end
""" exponent in area, `radius/r0`, relation """
function ae_sno end

""" exponent in area, `radius/r0`, (-) """
function Δa_rai end
""" exponent in area, `radius/r0`, (-) """
function Δa_sno end

""" coefficient in velocity, `radius/r0`, (m/s) """
function v0_sno end

""" coefficient in velocity, `radius/r0`, (-) """
function χv_rai end
""" coefficient in velocity, `radius/r0`, (-) """
function χv_sno end

""" exponent in velocity, `radius/r0`, """
function ve_rai end
""" exponent in velocity, `radius/r0`, """
function ve_sno end

""" exponent in velocity, `radius/r0`, (-) """
function Δv_rai end
""" exponent in velocity, `radius/r0`, (-) """
function Δv_sno end

""" drag coefficient for rain drops (-) """
function C_drag end

""" condensation/evaporation timescale (s) """
function τ_cond_evap end

""" sublimation/deposition timescale (s) """
function τ_sub_dep end

""" rain autoconversion threshold `∈(0.5, 1) * 1e-3` (-) """
function q_liq_threshold end

""" rain autoconversion timescale `∈(1e3, 1e4)` (s) """
function τ_acnv_rai end

""" snow autoconversion threshold (-) """
function q_ice_threshold end

""" snow autoconversion timescale (s) """
function τ_acnv_sno end

""" threshold between ice and snow (m) """
function r_ice_snow end

""" collision efficiency (-) """
function E_liq_rai end
""" collision efficiency (-) """
function E_liq_sno end
""" collision efficiency (-) """
function E_ice_rai end
""" collision efficiency (-) """
function E_ice_sno end
""" collision efficiency (-) """
function E_rai_sno end

""" ventilation factor coefficient for rain (-) """
function a_vent_rai end
""" ventilation factor coefficient for snow (-) """
function a_vent_sno end

""" ventilation factor coefficient for rain (-) """
function b_vent_rai end
""" ventilation factor coefficient for snow (-) """
function b_vent_sno end

""" rain evaporation scaling factor for coupling with EDMF (-) """
function microph_scaling end
""" snow deposition and sublimation scaling factor for coupling with EDMF (-) """
function microph_scaling_dep_sub end
""" snow melting scaling factor for coupling with EDMF (-) """
function microph_scaling_melt end

""" thermal conductivity of air (J/m/s/K) """
function K_therm end

""" diffusivity of water vapor (m``^2``/s) """
function D_vapor end

""" kinematic viscosity of air (m``^2``/s) """
function ν_air end

""" Schmidt number (-) """
function N_Sc end

end

end # module Atmos


text:
using Plots
x = [1,2,5]; y=x.^2

plot(x,y, label = "My First Graph")

max_stars_repo_name: dshehtanian/julia-learn

max_issues_repo_name: dshehtanian/julia-learn

max_forks_repo_name: dshehtanian/julia-learn


text:
# # Tensor Operations
# #-------------------
# scalar(t::Tensor)=iscnumber(space(t)) ? t.data[1] : throw(SpaceError("Not a scalar"))
#
# function tensorcopy!(t1::Tensor,labels1,t2::Tensor,labels2)
#     # Replaces tensor t2 with t1
#     N1=numind(t1)
#     perm=indexin(labels2,labels1)
#
#     length(perm) == N1 || throw(TensorOperations.LabelError("invalid label specification"))
#     isperm(perm) || throw(TensorOperations.LabelError("invalid label specification"))
#     for i = 1:N1
#         space(t1,i) == space(t2,perm[i]) || throw(SpaceError())
#     end
#     N1==0 && (t2.data[1]=t1.data[1]; return t2)
#     perm==[1:N1] && return copy!(t2,t1)
#     TensorOperations.tensorcopy_native!(t1.data,t2.data,perm)
#     return t2
# end
# function tensoradd!(alpha::Number,t1::Tensor,labels1,beta::Number,t2::Tensor,labels2)
#     # Replaces tensor t2 with beta*t2+alpha*t1
#     N1=numind(t1)
#     perm=indexin(labels2,labels1)
#
#     length(perm) == N1 || throw(TensorOperations.LabelError("invalid label specification"))
#     isperm(perm) || throw(TensorOperations.LabelError("invalid label specification"))
#     for i = 1:N1
#         space(t1,perm[i]) == space(t2,i) || throw(SpaceError("incompatible index spaces of tensors"))
#     end
#     N1==0 && (t2.data[1]=beta*t2.data[1]+alpha*t1.data[1]; return t2)
#     perm==[1:N1] && return (beta==0 ? scale!(copy!(t2,t1),alpha) : Base.LinAlg.axpy!(alpha,t1,scale!(t2,beta)))
#     beta==0 && (TensorOperations.tensorcopy_native!(t1.data,t2.data,perm);return scale!(t2,alpha))
#     TensorOperations.tensoradd_native!(alpha,t1.data,beta,t2.data,perm)
#     return t2
# end
# function tensortrace!(alpha::Number,A::Tensor,labelsA,beta::Number,C::Tensor,labelsC)
#     NA=numind(A)
#     NC=numind(C)
#     (length(labelsA)==NA && length(labelsC)==NC) || throw(LabelError("invalid label specification"))
#     NA==NC && return tensoradd!(alpha,A,labelsA,beta,C,labelsC) # nothing to trace
#
#     oindA=indexin(labelsC,labelsA)
#     clabels=unique(setdiff(labelsA,labelsC))
#     NA==NC+2*length(clabels) || throw(LabelError("invalid label specification"))
#
#     cindA1=Array(Int,length(clabels))
#     cindA2=Array(Int,length(clabels))
#     for i=1:length(clabels)
#         cindA1[i]=findfirst(labelsA,clabels[i])
#         cindA2[i]=findnext(labelsA,clabels[i],cindA1[i]+1)
#     end
#     isperm(vcat(oindA,cindA1,cindA2)) || throw(LabelError("invalid label specification"))
#
#     for i = 1:NC
#         space(A,oindA[i]) == space(C,i) || throw(SpaceError("space mismatch"))
#     end
#     for i = 1:div(NA-NC,2)
#         space(A,cindA1[i]) == dual(space(A,cindA2[i])) || throw(SpaceError("space mismatch"))
#     end
#
#     TensorOperations.tensortrace_native!(alpha,A.data,beta,C.data,oindA,cindA1,cindA2)
#     return C
# end
# function tensorcontract!(alpha::Number,A::Tensor,labelsA,conjA::Char,B::Tensor,labelsB,conjB::Char,beta::Number,C::Tensor,labelsC;method=:BLAS,buffer::TCBuffer=defaultcontractbuffer)
#     # Get properties of input arrays
#     NA=numind(A)
#     NB=numind(B)
#     NC=numind(C)
#
#     # Process labels, do some error checking and analyse problem structure
#     if NA!=length(labelsA) || NB!=length(labelsB) || NC!=length(labelsC)
#         throw(TensorOperations.LabelError("invalid label specification"))
#     end
#     ulabelsA=unique(labelsA)
#     ulabelsB=unique(labelsB)
#     ulabelsC=unique(labelsC)
#     if NA!=length(ulabelsA) || NB!=length(ulabelsB) || NC!=length(ulabelsC)
#         throw(TensorOperations.LabelError("tensorcontract requires unique label for every index of the tensor, handle inner contraction first with tensortrace"))
#     end
#
#     clabels=intersect(ulabelsA,ulabelsB)
#     numcontract=length(clabels)
#     olabelsA=intersect(ulabelsA,ulabelsC)
#     numopenA=length(olabelsA)
#     olabelsB=intersect(ulabelsB,ulabelsC)
#     numopenB=length(olabelsB)
#
#     if numcontract+numopenA!=NA || numcontract+numopenB!=NB || numopenA+numopenB!=NC
#         throw(LabelError("invalid contraction pattern"))
#     end
#
#     # Compute and contraction indices and check size compatibility
#     cindA=indexin(clabels,ulabelsA)
#     oindA=indexin(olabelsA,ulabelsA)
#     oindCA=indexin(olabelsA,ulabelsC)
#     cindB=indexin(clabels,ulabelsB)
#     oindB=indexin(olabelsB,ulabelsB)
#     oindCB=indexin(olabelsB,ulabelsC)
#
#     # check size compatibility
#     spaceA=space(A)
#     spaceB=space(B)
#     spaceC=space(C)
#
#     ospaceA=spaceA[oindA]
#     ospaceB=spaceB[oindB]
#
#     conjA=='C' || conjA=='N' || throw(ArgumentError("conjA should be 'C' or 'N'."))
#     conjB=='C' || conjB=='N' || throw(ArgumentError("conjB should be 'C' or 'N'."))
#
#     if conjA == conjB
#         for (i,j) in zip(cindA,cindB)
#             spaceA[i] == dual(spaceB[j]) || throw(SpaceError("incompatible index space for label $(ulabelsA[i])"))
#         end
#     else
#         for (i,j) in zip(cindA,cindB)
#             spaceA[i] == dual(conj(spaceB[j])) || throw(SpaceError("incompatible index space for label $(ulabelsA[i])"))
#         end
#     end
#     for (i,j) in zip(oindA,oindCA)
#         spaceC[j] == (conjA=='C' ? conj(spaceA[i]) : spaceA[i]) || throw(SpaceError("incompatible index space for label $(ulabelsA[i])"))
#     end
#     for (i,j) in zip(oindB,oindCB)
#         spaceC[j] == (conjB=='C' ? conj(spaceB[i]) : spaceB[i]) || throw(SpaceError("incompatible index space for label $(ulabelsB[i])"))
#     end
#
#     if method==:BLAS
#         TensorOperations.tensorcontract_blas!(alpha,A.data,conjA,B.data,conjB,beta,C.data,buffer,oindA,cindA,oindB,cindB,oindCA,oindCB)
#     elseif method==:native
#         if NA>=NB
#             TensorOperations.tensorcontract_native!(alpha,A.data,conjA,B.data,conjB,beta,C.data,oindA,cindA,oindB,cindB,oindCA,oindCB)
#         else
#             TensorOperations.tensorcontract_native!(alpha,B.data,conjB,A.data,conjA,beta,C.data,oindB,cindB,oindA,cindA,oindCB,oindCA)
#         end
#     else
#         throw(ArgumentError("method should be :BLAS or :native"))
#     end
#     return C
# end
#
# function tensorproduct!(alpha::Number,A::Tensor,labelsA,B::Tensor,labelsB,beta::Number,C::Tensor,labelsC)
#     # Get properties of input arrays
#     NA=numind(A)
#     NB=numind(B)
#     NC=numind(C)
#
#     # Process labels, do some error checking and analyse problem structure
#     if NA!=length(labelsA) || NB!=length(labelsB) || NC!=length(labelsC)
#         throw(TensorOperations.LabelError("invalid label specification"))
#     end
#     NC==NA+NB || throw(TensorOperations.LabelError("invalid label specification for tensor product"))
#
#     tensorcontract!(alpha,A,labelsA,'N',B,labelsB,'N',beta,C,labelsC;method=:native)
# end
#


text:
pa=FBD.simple_problem()

function adjtest()
	x=randn(size(F,2))
	y=randn(size(F,1))
	a=LinearAlgebra.dot(y,F*x)
	b=LinearAlgebra.dot(x,adjoint(F)*y)
	c=LinearAlgebra.dot(x, transpose(F)*F*x)
	println("adjoint test: ", a, "\t", b)
	@test isapprox(a,b,rtol=1e-6)
	println("must be positive: ", c)
	@test c>0.0
end

@testset "IBD" begin
	global p=pa.pfibd
	for attrib in [FBD.S(), FBD.G()]
		global F=FBD.operator(p, attrib);
		adjtest()
	end
end

@testset "BD" begin
	global p=pa.plsbd
	for attrib in [FBD.S(), FBD.G()]
		global F=FBD.operator(p, attrib);
		adjtest()
	end
end

# adj tests of filter
pa=FBD.P_bandpass(Float64, fmin=0.1, fmax=0.4, nt=101)
global F=FBD.create_operator(pa)
@testset "bandpass" begin
	adjtest()
end


text:
using ProximalAlgorithms
using ProximalOperators
using LinearAlgebra
using Random
using Test

@testset "Nonconvex QP (tiny, $T)" for T in [Float64]
    Q = Matrix(Diagonal(T[-0.5, 1.0]))
    q = T[0.3, 0.5]
    low = T(-1.0)
    upp = T(+1.0)

    f = Quadratic(Q, q)
    g = IndBox(low, upp)

    n = 2

    Lip = maximum(diag(Q))
    gamma = T(0.95) / Lip

    @testset "PANOC" begin
        x0 = zeros(T, n)
        solver = ProximalAlgorithms.PANOC{T}()
        x, it = solver(x0, f=f, g=g)
        z = min.(upp, max.(low, x .- gamma .* (Q*x + q)))
        # println(it, " ", 0.5*dot(x, Q*x) + dot(q, x))
        @test norm(x - z, Inf)/gamma <= solver.tol
    end

    @testset "ZeroFPR" begin
        x0 = zeros(T, n)
        solver = ProximalAlgorithms.ZeroFPR{T}()
        x, it = solver(x0, f=f, g=g)
        z = min.(upp, max.(low, x .- gamma .* (Q*x + q)))
        # println(it, " ", 0.5*dot(x, Q*x) + dot(q, x))
        @test norm(x - z, Inf)/gamma <= solver.tol
    end

    @testset "LiLin" begin
        x0 = zeros(T, n)
        solver = ProximalAlgorithms.LiLin{T}(gamma=gamma)
        x, it = solver(x0, f=f, g=g)
        z = min.(upp, max.(low, x .- gamma .* (Q*x + q)))
        # println(it, " ", 0.5*dot(x, Q*x) + dot(q, x))
        @test norm(x - z, Inf)/gamma <= solver.tol
    end
end

@testset "Nonconvex QP (small, $T)" for T in [Float64]
    @testset "Random problem $k" for k in 1:5
        Random.seed!(k)

        n = 100
        A = randn(T, n, n)
        U, R = qr(A)
        eigenvalues = T(2) .* rand(T, n) .- T(1)
        Q = U*Diagonal(eigenvalues)*U'
        Q = 0.5*(Q + Q')
        q = randn(T, n)

        low = T(-1.0)
        upp = T(+1.0)

        f = Quadratic(Q, q)
        g = IndBox(low, upp)

        Lip = maximum(abs.(eigenvalues))
        gamma = T(0.95) / Lip

        TOL = 1e-4

        @testset "PANOC" begin
            x0 = zeros(T, n)
            solver = ProximalAlgorithms.PANOC{T}(tol=TOL)
            x, it = solver(x0, f=f, g=g)
            z = min.(upp, max.(low, x .- gamma .* (Q*x + q)))
            # println(it, " ", 0.5*dot(x, Q*x) + dot(q, x))
            @test norm(x - z, Inf)/gamma <= solver.tol
        end

        @testset "ZeroFPR" begin
            x0 = zeros(T, n)
            solver = ProximalAlgorithms.ZeroFPR{T}(tol=TOL)
            x, it = solver(x0, f=f, g=g)
            z = min.(upp, max.(low, x .- gamma .* (Q*x + q)))
            # println(it, " ", 0.5*dot(x, Q*x) + dot(q, x))
            @test norm(x - z, Inf)/gamma <= solver.tol
        end

        @testset "LiLin" begin
            x0 = zeros(T, n)
            solver = ProximalAlgorithms.LiLin(gamma=gamma, tol=TOL)
            x, it = solver(x0, f=f, g=g)
            z = min.(upp, max.(low, x .- gamma .* (Q*x + q)))
            # println(it, " ", 0.5*dot(x, Q*x) + dot(q, x))
            @test norm(x - z, Inf)/gamma <= solver.tol
        end
    end
end


text:
using FANN
using Base.Test

# test simple linear problem for which we know the solution
net = ANN([20, 1], [:linear], b=0.1)

X = randn(20, 5000)
y = mean(X, 1)

dset = DataSet(X, y)

# we just test these work
# ~~~~~~~~~~~~~~~~~~~~~~
#  ~~~ QUICKPROP ~~~~~
setup_qprop!(net; mu=1.7, decay=-0.00001)

# train
train!(net, dset,
	   max_epochs=1000,
	   epochs_between_reports=0,
	   desired_error=1e-8)

# ~~~~~~~~~~~~~~~~~~~~~~
#  ~~~~~ RPROP ~~~~~~~~~
setup_rprop!(net; increase_factor=1.5, decrease_factor=0.9, delta_min=0.0, delta_zero=0.001, delta_max=10.0)

# train
train!(net, dset,
	   max_epochs=1000,
	   epochs_between_reports=0,
	   desired_error=1e-8)

# ~~~~~~~~~~~~~~~~~~~~~~
#  ~~~~~ BATCH ~~~~~~~~~
setup_bprop!(net; batch=true, learning_rate=0.5, learning_momentum=0.1)

# train
train!(net, dset,
	   max_epochs=1000,
	   epochs_between_reports=0,
	   desired_error=1e-8)

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#  ~~~~~ INCREMENTAL ~~~~~~~~~
setup_bprop!(net; batch=false, learning_rate=0.5, learning_momentum=0.1)

# train
train!(net, dset,
	   max_epochs=1000,
	   epochs_between_reports=0,
	   desired_error=1e-8)


text:
srcdir = "/Users/anray/Desktop/TDGP/src"
any(srcdir .== LOAD_PATH) || push!(LOAD_PATH, srcdir)
any(pwd() .== LOAD_PATH) || push!(LOAD_PATH, pwd())
using PyPlot, Random, Statistics, LinearAlgebra, Test, Revise,
MCMC_Driver, GP
##
Random.seed!(11)
xtest = collect(0:0.005:1)'
λ² = [0.1].^2
δ = 0.01
δtry = 10δ
p = 2
demean = true
K = GP.Mat32()

f = zeros(size(xtest,2))
change1 = round(Int, 1/2*length(f))
f[change1:end] .= 1
f .-= 0.5
ynoisy = f + δ*randn(length(f))
dec = 8
ytrain = ynoisy[1:dec:end]'
xtrain = xtest[1:dec:end]'

# f1 = figure()
# plot(xtest[:][:],f, "-k",label="true")
# plot(xtrain[:],ytrain, "+m", label="train",markersize=10)
#
ytest, σ2, σ_prior  = GP.GPfit(K, ytrain, xtrain, xtest, λ², δtry, p=p, demean=demean)
# plot(xtest[:],ytest,label="testing")
# plot(xtest[:],ytest+2sqrt.(diag(σ2)),"--r",alpha=0.2)
# plot(xtest[:],ytest-2sqrt.(diag(σ2)),"--b",alpha=0.2)

## non stationary
λ²test = λ²[1]*ones(size(xtest))
λ²train = λ²[1]*ones(Float64, size(xtrain))
ytest_ns, σ2_ns, σ_prior_ns  = GP.GPfit(K, ytrain, xtrain, xtest, λ²test, λ²train, δtry, p=p, demean=demean)
@test (norm(mean(ytest - ytest_ns)) < 1e-12)
f2, ax2 = plt.subplots(2,1, sharex=true, figsize=(7,7))
ax2[1].plot(xtest[:],f, "-k",label="true", linewidth=2)
ax2[1].plot(xtrain[:],ytrain, "xm", label="train",markersize=10, alpha=0.8)
ax2[1].plot(xtest[:],ytest_ns,label="fixed λ", linewidth=3)
# plot(xtest[:],ytest_ns+2sqrt.(diag(σ2_ns)),"--r",alpha=0.2)
# plot(xtest[:],ytest_ns-2sqrt.(diag(σ2_ns)),"--b",alpha=0.2)
##
lscale_max² = λ²[1]
dropoff = 0.00025
λ²test = lscale_max²*(1 .+eps() .- exp.(-1/dropoff*abs.(xtest.-xtest[change1+1]).^2)).^2
λ²train[:] = λ²test[1:dec:end]
ytest_ns, σ2_ns, σ_prior_ns  = GP.GPfit(K, ytrain, xtrain,
    xtest, λ²test, λ²train, δtry, p=2, demean=demean, nogetvars=false)

ax2[1].plot(xtest[:],ytest_ns,label="variable λ", linewidth=3, alpha=0.7, "-r")
ax2[2].plot(xtest[:],sqrt.(λ²test[:]), linewidth=3, color="r", alpha=0.7)
ax2[1].set_ylabel(L"f(x)")
ax2[2].set_xlabel(L"x")
ax2[2].set_ylabel(L"\lambda")
ax2[1].grid()
ax2[2].grid()
MCMC_Driver.nicenup(gcf(), fsize=14)
savefig("variable.png", dpi=300)
## timing
@info "meshkernel with broadcast"
@time for i = 1:1000
    GP.meshkernel(K, xtrain, xtest, λ²test, λ²train, p)
end
@info "mapkernel with pairwise"
@time for i = 1:1000
    GP.mapkernel(K, xtrain, xtest, λ²test, λ²train)
end
@info "loops and full expression"
@time for i = 1:1000
    GP.makekernel(K, xtrain, xtest, λ²test, λ²train, p)
end
## unit tests
@testset "Non stationary kernel tests" begin
    @testset "broadcast and map" begin
        @test norm(mean(GP.meshkernel(K, xtrain, xtest, λ²test, λ²train, p) -
                    GP.mapkernel(K, xtrain, xtest, λ²test, λ²train))) < 1e-12
        end
    @testset "map and full kernel" begin
        @test norm(mean(GP.makekernel(K, xtrain, xtest, λ²test, λ²train, p) -
                GP.mapkernel(K, xtrain, xtest, λ²test, λ²train)')) < 1e-12
        end
end
## now to 2D
Random.seed!(11)
x, y = xtest, xtest[1:round(Int, length(xtest)/2)]
x2D_test = zeros(2,length(x)*length(y))
λ²2D_test = zeros(2,length(x)*length(y))
for i in 1:size(x2D_test,2)
    xid, yid = Tuple(CartesianIndices((length(x),length(y)))[i])
    x2D_test[:,i] = [x[xid]; y[yid]]
    λ²xid, λ²yid = Tuple(CartesianIndices((length(x),length(y)))[i])
    λ²2D_test[1,i] = λ²test[λ²xid]
end
λ²2D_test[2,:] .= abs(y[end]-y[1])^2
y2D       = [f[i] for i in 1:length(x), j in 1:length(y)]
y2D_noisy = y2D + δ*randn(size(y2D))
y2D_train = y2D_noisy[1:dec:end, 1:dec:end][:]'
lidx = LinearIndices((1:length(x),1:length(y)))[1:dec:end, 1:dec:end][:]
λ²2D_train = λ²2D_test[:,lidx]
x2D_train = x2D_test[:,lidx]
## fit GPs
y2D_test_ns, = GP.GPfit(K, y2D_train, x2D_train,
    x2D_test, λ²2D_test, λ²2D_train, δtry, p=2, demean=demean, nogetvars=true)
y2D_test, = GP.GPfit(K, y2D_train, x2D_train,x2D_test,
            [λ²[1], λ²2D_test[2,1]],
            δtry, p=2, demean=demean, nogetvars=true)
## plot
fig, ax = plt.subplots(1,3, sharex=true, sharey=true)
im1 = ax[1].imshow(y2D_noisy)
ax[1].set_title("data");
colorbar(im1, ax=ax[1])
im2 = ax[2].imshow(reshape(y2D_test, size(y2D, 1), size(y2D, 2)))
ax[2].set_title("stationary GP")
colorbar(im2, ax=ax[2])
im3 = ax[3].imshow(reshape(y2D_test_ns, size(y2D, 1), size(y2D, 2)))
ax[3].set_title("non stationary GP")
colorbar(im3, ax=ax[3])
## see how fast you can change a row
Kstar = GP.mapkernel(K, x2D_train, x2D_test, λ²2D_test, λ²2D_train)
K_y  = GP.mapkernel(K, x2D_train, x2D_train, λ²2D_train, λ²2D_train) +
    Matrix(δtry^2*I, (size(x2D_train,2)), (size(x2D_train,2)))
A = copy(Kstar)
B = copy(K_y)
# change a training point
n = rand(1:length(y2D_train))
x2D_train[:,n] = x2D_test[:,n]
λ²2D_train[:,n] = λ²2D_test[:,n]
y2D_train[n]   = y2D_noisy[CartesianIndices((length(x),length(y)))[n]]
D, = GP.GPfit(K, y2D_train, x2D_train,
    x2D_test, λ²2D_test, λ²2D_train, δtry, p=2, demean=demean, nogetvars=true)
@time begin
Kstarv = @view Kstar[:,n]
map!(x²->x²,Kstarv,GP.colwise(K, x2D_train[:,n], x2D_test, λ²2D_train[:,n], λ²2D_test))
K_yv = @view K_y[n,:]
map!(x²->x²,K_yv,GP.colwise(K, x2D_train[:,n], x2D_train, λ²2D_train[:,n], λ²2D_train))
K_y[:,n] = K_y[n,:]
K_y[n,n] = K_y[n,n] + δtry^2
my = mean(y2D_train, dims=2)
y2D_train = y2D_train .- my
U = cholesky(K_y).U
C = my' .+ Kstar*(U\(U'\y2D_train'))
end
@test norm(mean(C - D)) < 1e-12


text:
using BinaryBuilder, Pkg

name = "Pathfinder"
version = v"0.5.0"

sources = [
    GitSource("https://github.com/servo/pathfinder.git",
              "7281a607a80dacc5470d7d3ec208341b6f0b0727"),
    DirectorySource("./bundled")
]

script = raw"""
cd ${WORKSPACE}/srcdir/pathfinder*/c

# by default, pathfinder will compile to a static library,
# so we change "staticlib" to "cdylib" to compile to a dynamic library
atomic_patch -d .. -p1 ../../patches/Compile-to-dynamic-lib.patch

# generate the header file
cbindgen \
    --lang c \
    --crate pathfinder_c \
    --config cbindgen.toml \
    --output "${includedir}/pathfinder.h"

# build pathfinder
cargo build --release --lib --target=${rust_target}

# install the library
install -Dvm 755 ../target/${rust_target}/release/*pathfinder.${dlext} "${libdir}/libpathfinder.${dlext}"

# install the licenses
install_license ../LICENSE-APACHE ../LICENSE-MIT
"""

platforms = supported_platforms()

# font-kit fails to compile on arm and i686 linux
filter!(p -> arch(p) ∉ ("armv6l", "armv7l"), platforms)
filter!(p -> !Sys.islinux(p) || arch(p) != "i686", platforms)
# cdylib apparently doesn't support musl
filter!(p -> libc(p) != "musl", platforms)
# Rust toolchain for i686 Windows is unusable
filter!(p -> !Sys.iswindows(p) || arch(p) != "i686", platforms)

products = [
    LibraryProduct("libpathfinder", :libpathfinder),
]

not_windows = filter(!Sys.iswindows, platforms)

dependencies = [
    Dependency("Fontconfig_jll", platforms=not_windows),
    Dependency("FreeType2_jll",  platforms=not_windows),
    Dependency("HarfBuzz_jll",   platforms=not_windows),
    HostBuildDependency(PackageSpec(; name="cbindgen_jll", uuid="a52b955f-5256-5bb0-8795-313e28591558")),
]

compilers = [:rust, :c]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; compilers, julia_compat="1.6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
export PDToken
export tokenize_space

mutable struct PDToken
    data
    children::Vector{PDToken}
    page::Int
    fcoord::Rectangle
    gcoord::Rectangle
end

function PDToken(chars::Vector{PDChar})
    str = join(map(c -> c.str, chars))
    fcoord = merge(map(c -> c.fcoord, chars))
    gcoord = merge(map(c -> c.gcoord, chars))
    PDToken(chars, str, fcoord, gcoord)
end

function tokenize_space(chars::Vector{PDChar})
    tokens = PDToken[]
    buffer = [chars[1]]
    average = chars[1].fcoord.w
    for i = 2:length(chars)
        prev, curr = chars[i-1], chars[i]
        expected = prev.fcoord.x + prev.fcoord.w + 0.3*average
        if prev.page != curr.page || prev.fcoord.x > curr.fcoord.x || expected < curr.fcoord.x
            push!(tokens, PDToken(buffer))
            buffer = [curr]
            average = curr.fcoord.w
        else
            push!(buffer, curr)
            average = (average + curr.fcoord.w) / 2
        end
    end
    isempty(buffer) || push!(tokens,PDToken(buffer))
    tokens
end

function writepdftxt(filename::String, tokens::Vector{PDToken})
    open(filename, "w") do io
        for t in tokens
            println(io, string(t))
        end
    end
end

export maketable
function maketable(chars::Vector{PDChar})
    tokens = tokenize_space(chars)
    maketable(tokens)
end

function maketable(tokens::Vector{PDToken})
    sorted = sort(tokens, by=t->t.gcoord.x)
    groups = [tokens[1]]
    for i = 2:length(tokens)
        t = sorted[i]
        prev = groups[end]
        if t.gcoord.x <= prev.gcoord.x+prev.gcoord.w
            push!(prev, t)
        else
            push!(groups, t)
        end
    end
end


text:
workspace()

 function postOrder(par,child)
        N = length(par)
        order = zeros(Int64,N)
        root = find(x->x==0,par)[1]
        stack = Array{Int64,1}(0)
        iii = N
        push!(stack,root)
        while !isempty(stack)
            v = pop!(stack)
            order[v] = iii
            iii-=1
            push!(stack,child[v]...)
        end

        return order
    end

N = 6
par = [2;4;4;5;0;5]
child = [Array{Int64}(0) for i=1:N]
for iii=1:N
    par[iii] != 0 && (push!(child[par[iii]],iii))
end

postOrder(par,child)


text:
struct RecursiveVec{T<:Union{Tuple,AbstractVector}}
    vecs::T

end
function RecursiveVec(arg1::AbstractVector{T}) where T
    if isbitstype(T)
        return RecursiveVec((arg1,))
    else
        return RecursiveVec{typeof(arg1)}(arg1)
    end
end
RecursiveVec(arg1, args...) = RecursiveVec((arg1, args...))

Base.getindex(v::RecursiveVec, i) = v.vecs[i]

Base.iterate(v::RecursiveVec) = iterate(v.vecs)
Base.iterate(v::RecursiveVec, s) = iterate(v.vecs, s)
Base.IteratorEltype(::Type{<:RecursiveVec}) = Base.EltypeUnknown() # since `eltype` is not the eltype of the iterator
Base.IteratorSize(::Type{<:RecursiveVec}) = Base.HasLength()
Base.length(v::RecursiveVec) = length(v.vecs)

Base.first(v::RecursiveVec) = first(v.vecs)
Base.last(v::RecursiveVec) = last(v.vecs)

Base.eltype(v::RecursiveVec) = eltype(typeof(v))
Base.eltype(::Type{RecursiveVec{T}}) where {T<:Tuple} = _eltype(T)
Base.eltype(::Type{RecursiveVec{T}}) where {T<:AbstractVector} = eltype(eltype(T))

_eltype(::Type{Tuple{T}}) where {T} = eltype(T)
function _eltype(::Type{TT}) where {TT<:Tuple}
    T = eltype(Base.tuple_type_head(TT))
    T2 = _eltype(Base.tuple_type_tail(TT))
    T == T2 ? T : error("all elements of a `RecursiveVec` should have same `eltype`")
end

Base.:-(v::RecursiveVec) = RecursiveVec(map(-,v.vecs))
Base.:+(v::RecursiveVec, w::RecursiveVec) = RecursiveVec(map(+,v.vecs, w.vecs))
Base.:-(v::RecursiveVec, w::RecursiveVec) = RecursiveVec(map(-,v.vecs, w.vecs))
Base.:*(v::RecursiveVec, a) = RecursiveVec(map(x->x*a, v.vecs))
Base.:*(a, v::RecursiveVec) = RecursiveVec(map(x->a*x, v.vecs))
Base.:/(v::RecursiveVec, a) = RecursiveVec(map(x->x/a, v.vecs))
Base.:\(a, v::RecursiveVec) = RecursiveVec(map(x->a\x, v.vecs))

function Base.similar(v::RecursiveVec, ::Type{T} = eltype(v)) where {T}
    RecursiveVec(map(x->similar(x,T), v.vecs))
end

function Base.fill!(v::RecursiveVec, a)
    for x in v.vecs
        fill!(x, a)
    end
    return v
end

function Base.copyto!(w::RecursiveVec, v::RecursiveVec)
    @assert length(w.vecs) == length(v.vecs)
    @inbounds for i = 1:length(w.vecs)
        copyto!(w.vecs[i], v.vecs[i])
    end
    return w
end

function LinearAlgebra.mul!(w::RecursiveVec, a, v::RecursiveVec)
    @assert length(w.vecs) == length(v.vecs)
    @inbounds for i = 1:length(w.vecs)
        mul!(w.vecs[i], a, v.vecs[i])
    end
    return w
end

function LinearAlgebra.mul!(w::RecursiveVec, v::RecursiveVec, a)
    @assert length(w.vecs) == length(v.vecs)
    @inbounds for i = 1:length(w.vecs)
        mul!(w.vecs[i], v.vecs[i], a)
    end
    return w
end

function LinearAlgebra.rmul!(v::RecursiveVec, a)
    for x in v.vecs
        rmul!(x, a)
    end
    return v
end

function LinearAlgebra.axpy!(a, v::RecursiveVec, w::RecursiveVec)
    @assert length(w.vecs) == length(v.vecs)
    @inbounds for i = 1:length(w.vecs)
        axpy!(a, v.vecs[i], w.vecs[i])
    end
    return w
end
function LinearAlgebra.axpby!(a, v::RecursiveVec, b, w::RecursiveVec)
    @assert length(w.vecs) == length(v.vecs)
    @inbounds for i = 1:length(w.vecs)
        axpby!(a, v.vecs[i], b, w.vecs[i])
    end
    return w
end

LinearAlgebra.dot(v::RecursiveVec{T}, w::RecursiveVec{T}) where {T} = sum(dot.(v.vecs, w.vecs))
LinearAlgebra.norm(v::RecursiveVec) = norm(norm.(v.vecs))


text:
module DiseaseOutbreak

include("distancing.jl")
export makeCitizens, infectCitizens!, followCitizens!
export uninfected, newInfections, cumulativeInfections, infectionSpread

include("dynamics.jl")
export Dynamics, Population

include("evolution.jl")
export evolve

include("sirx.jl")
include("sei3r.jl")
export SEI3R
export SIRX
export getParams, initialize
export caseModel

include("casemodelfit.jl")
export CaseModelFitResult
export fitCaseModel, summary, fitted, estimatedStates

include("sirplot.jl")
export plotEvolution, plotFit
export pyplotEvolution, pyplotFit
end # module