text:
###############################################
# Version of the objective function to be used
# by BlackBoxOptim and the simulated annealing
# algorithm
##############################################
function objective_function_2(parameters, Weighting_matrix::Matrix, random_numbers::Vector, b0)
    ##################
    # Function that returns the inverse of lognormal distribution
    function logninv(p,mu,sigma)
        #output: scalar x
        logx0 = -sqrt(2)*erfcinv(2*p);
        x = exp(sigma*logx0 + mu);
        return x
    end
    #################
    ##########################
    #define match productivity:
    function yi_m(y_index,x_index, ygrid, xgrid)
        match = ygrid[y_index,1]*xgrid[x_index,1]
        return match
    end
    ##########################
    ##########################
    # define opportunity cost of employment:
    function zi_m(y_index,x_index, ygrid, xgrid, z0, alpha)
        op_cost = z0 + alpha*(yi_m(y_index,x_index, ygrid, xgrid) - z0) #call the function defined above
        return op_cost
    end
    ##########################
    # Exit from unemployment rate: ft
    function ft_f(St_m,ut_m,ut,lm,lambda0,M)
    #St_m a vector of dimension M times 1
        ft = 0;
        for i=1:M
            #calculate the sum
            if St_m[1,i]> 0
                ft = ft + ut_m[1,i]*lm[i,1];
            end
        end
        #multiply by lamba0 and divide by the overall level of unemployment ut:
        ft = lambda0*ft/ut;
        return ft
    end

    #function that calculates the job destruction rate
    function st_f(St_m,ut_m,ut,lm, delta, M)
        st = 0;
        for i=1:M
        #calculate the integral
            if St_m[1,i] <= 0
                st = st + (1 - ut_m[1,i])*lm[i,1];
            end
        end
        #finish the calculation:
        st = delta + (1 - delta)*st/(1 - ut);
        return st
    end

    function ut_plus1_f(St_m, ut_m, delta, lambda0, M)
    #Calculates the next period's unemployment rate
        #initialization:
        ut_plus1_m = zeros(1,M);

        for i=1:M
            #if the surplus is negative, unemployment is 1:
            if St_m[1,i] <= 0
                ut_plus1_m[1,i]= 1;
            else
                ut_plus1_m[1,i]= ut_m[1,i] + delta*(1 - ut_m[1,i]) - lambda0*ut_m[1,i];
            end
        end
        return ut_plus1_m
    end

    ############
    # Parameters
    ############

    #parameters to estimate:
    """
    z0 = parameters[1];
    sigma = parameters[2];
    pho = parameters[3];
    lambda0 = parameters[4];
    eta =  parameters[5];
    mu =  parameters[6];
    x_lower_bound = parameters[7];
    delta = parameters[8];

    k = 0.12;
    lambda1  = k*lambda0;
    #x_lower_bound = 0.73
    alpha = 0.64
    tau = 0.5
    #delta = 0.046 #exogenous layoff rate
    r = 0.05/4 # interest rate
    discount = (1 - delta)/(1 + r)
    epsilon = 0.002
    """

    """
    z0 = parameters[1];
    sigma = parameters[2];
    lambda0 = parameters[3];
    eta =  parameters[4];
    mu =  parameters[5];
    x_lower_bound = parameters[6] #lower bound of the ability level
    alpha = parameters[7]
    delta = parameters[8]

    pho = 0.94
    k = 0.12;
    lambda1  = k*lambda0;
    #x_lower_bound = 0.73
    #alpha = 0.64
    tau = 0.5
    #delta = 0.042 # exogenous layoff rate

    r = 0.05/4 # interest rate
    discount = (1 - delta)/(1 + r)
    epsilon = 0.002
    """

    z0 = parameters[1];
    sigma = parameters[2];
    lambda0 = parameters[3];
    eta =  parameters[4];
    mu =  parameters[5];
    x_lower_bound = parameters[6] #lower bound of the ability level

    pho = 0.94
    k = 0.12;
    lambda1  = k*lambda0;
    #x_lower_bound = 0.73
    alpha = 0.64
    tau = 0.5
    delta = 0.042 # exogenous layoff rate

    r = 0.05/4 # interest rate
    discount = (1 - delta)/(1 + r)
    epsilon = 0.002

    # grid parameters
    N = 100 #number of states
    M = 500 #number of ability levels

    # Define the grid for ability x
    xgrid_space = linspace(x_lower_bound,x_lower_bound+1,M); # column vector with M rows

    xgrid = Array{Float64}(M,1);
    for i=1:M
   		xgrid[i,1] = xgrid_space[i];
    end

    # Define the grid for match ability yi:
        # 1. a:
        agrid = linspace((0+epsilon),(1-epsilon),N);

        # 2. yi:
        #intialization:
        ygrid = Array{Float64}(N,1);

        for i=1:N
            ygrid[i] = logninv(agrid[i],0,sigma); #calculate the inverse of the lognormal
        end

    #Markov transition matrix:
    Markov = zeros(N,N);

    # number of dimension
    # AR1 parameter:

    cop = Copulas.NormalCopula(2,pho)

    #Gaussian copula pdf:
    for j=1:N #move along the column:
        for i=1:N #move along the lines:
            Markov[i,j] = Copulas.dnormCopula([agrid[i] agrid[j]], cop)[];
        end
    end

    #normalize so that each row sum to 1:
    for i = 1:N
        Markov[i,:] = Markov[i,:]./sum(Markov[i,:]);
    end

    # distribution of workers on the grid x
    lm = Array{Float64}(M,1);

    # load the beta pdf
    b = Beta(eta, mu)

    for i=1:M
        lm[i,1] = pdf(b, xgrid[i] - x_lower_bound); #evaluation of the beta pdf at a given point
    end

    #normalization so that mass of workers sum up to 1:
    lm[:,1] = lm[:,1]./sum(lm[:,1]);

    # Pre calculate the opportunity cost along the grid of ability:
    # N times M matrix
    # rows = aggregate staes
    # columns = ability level
    zi_m_grid = zeros(N,M);

    for i=1:N
        for m=1:M
            zi_m_grid[i,m] = zi_m(i,m, ygrid, xgrid, z0, alpha);
        end
    end

    ##########################################################
    # Calculate the match surplus by value function iteration:
    # N times M matrices
    # N = aggregate level
    # M = ability level
    G = zeros(N,M); #intialization
    p = zeros(N,M);

    # Pre calculate values:
    for i = 1:N
        for m = 1:M
            G[i,m] = yi_m(i,m, ygrid, xgrid) - zi_m_grid[i,m];
            p[i,m] = yi_m(i,m, ygrid, xgrid);
        end
    end

    # N times M matrix
    # column index indicates the workers' type
    # row index indicates the shock index

        Si_m = zeros(N,M);

        #tol = 0.0001;
        tol = 0.01;
        maxits = 300;
        dif = tol+tol;
        its = 1;

        #initialization:
        up = ones(N,M);
        up_plus1 = zeros(N,M);
        compare = zeros(N,M);

        while dif>tol

            up_plus1 =  G + discount*Markov*max(up,compare)
            dif = norm(up_plus1 - up);
            up = up_plus1;
            its = its + 1;

            if its > maxits
                break
            end
        end

        Si_m = up;

###########################
# Simulation of the economy
###########################
n_years = 2500; # years to simulate
number_periods = n_years*4; # one period = one quarter

discard = round(Int,floor(number_periods/10)); #get rid of the first 10th observations

y_index_r = Array{Int16}(number_periods,1); #store the indexes of the shock
y_r = Array{Float64}(number_periods,1) #store value of the shock

#initial shock:
y_index_r[1,1] = 50; #initial shock index
y_r[1,1] = ygrid[y_index_r[1,1], 1]; #initial shock value

#Initialization:
ut_m_r = ones(number_periods,M); #row for the period, column for ability
St_m_r = Array{Float64}(number_periods,M); #row for the period, column for ability

ut_r = Array{Float64}(number_periods,1);

ft_r = Array{Float64}(number_periods,1);
qt_r = Array{Float64}(number_periods,1);
st_r = Array{Float64}(number_periods,1);

######################
# Loop over the economy:
for t=1:(number_periods-1)

    # Calculate the aggregate unemployment
    ut_r[t,1] = dot(ut_m_r[t,:],lm[:,1]); #dot product.

    # Calculate the surplus given the actual value of the shock
    St_m_r[t,:] = Si_m[y_index_r[t,1],:];

    #Remark: everytime I select a row, Julia treats it as a column vector
    #That is why I have to transpose my inputs in the next function
    # Exit rate from unemployment:
    ft_r[t,1] = ft_f(St_m_r[t,:]', ut_m_r[t,:]', ut_r[t,1], lm, lambda0, M);

    #Job destruction rate:
    st_r[t,1] = st_f(St_m_r[t,:]', ut_m_r[t,:]', ut_r[t,1], lm, delta, M);

    # Law of motion of unemployment:
    ut_m_r[t+1,:] = ut_plus1_f(St_m_r[t,:]', ut_m_r[t,:]', delta, lambda0, M);

    # New shock from the markov transition matrix:
    #r = rand(); #draw a number from the uniform distribution on 0, 1
    r = random_numbers[t]

    # I use the Markov transition matrix previously defined:
    prob = Markov[y_index_r[t,1],:];

    #I stock the index in the Markov matrix, as well as the value of the
    #shock y
    y_index_r[t+1,1] = sum(r .>= cumsum([0; prob]));

    y_r[t+1,1] = ygrid[y_index_r[t+1,1], 1];
end

#Check for odd values. These entities cannot be negative or zero:
for i = discard:(number_periods-1)
    if ut_r[i,1]<=0
        return Inf #+infinity of type Float64
    elseif ft_r[i,1]<=0
        return Inf
    elseif st_r[i,1]<=0
        return Inf
    end
end

#Moments to match:
mean_prod = mean(y_r[discard:(number_periods-1),1]) # Productivity yt
std_prod = std(log(y_r[discard:(number_periods-1),1]))

mean_unemployment = mean(ut_r[discard:(number_periods-1),1]) #Unemployment ut
std_unemployment = std(log(ut_r[discard:(number_periods-1),1]))

mean_exit_rate = mean(ft_r[discard:(number_periods-1),1]) #Exit rate from unemployment ft
std_exit_rate = std(log(ft_r[discard:(number_periods-1),1]))

mean_job_destruction_rate = mean(st_r[discard:(number_periods-1),1]) #Job separation rate st
std_job_destruction_rate = std(log(st_r[discard:(number_periods-1),1]))

#Check for odd values. Don't want an economy with constant values:

if std_prod == 0
    return Inf
elseif std_unemployment == 0
    return Inf
elseif std_exit_rate == 0
    return Inf
elseif std_job_destruction_rate == 0
    return Inf
end

#kurtosis_unemployment = kurtosis(log(ut_r[discard:(number_periods-1),1])) + 3 #returns excess kurtosis, so I add 3

#b1 = [mean_prod; std_prod; mean_unemployment; std_unemployment; kurtosis_unemployment; mean_exit_rate]; #simulated moments
#b0 = [ 1; 0.0226;  0.058;  0.22; 2.65; 0.78]; #moments to match from the original paper
#b0 = [ 1; 0.023;  0.059;  0.210; 2.374; 0.761] #with my empirical momemnts

#b1 = [mean_prod; std_prod; mean_unemployment; std_unemployment; mean_exit_rate; std_exit_rate ]; #simulated moments
#b0 = [ 1; 0.023;  0.059;  0.218; 0.762; 0.2358] #without kurtosis, with the standard deviation of the exit rate

#b1 = [mean_prod; std_prod; mean_unemployment; std_unemployment; kurtosis_unemployment; mean_exit_rate; std_exit_rate]; #simulated moments
#b0 = [ 1.0; 0.023;  0.059;  0.210; 2.374; 0.761; 0.085]; #moments to match

# mean productivity, std productivity, mean unemployment rate, std unemployment rate, mean job finding rate, std job finding rate, mean job separation rate, std separation rate
#b1 = [mean_prod; std_prod; mean_unemployment; std_unemployment; mean_exit_rate; std_exit_rate; mean_job_destruction_rate ; std_job_destruction_rate] #simulated moments

# mean productivity, std productivity, mean unemployment rate, std unemployment rate, mean job finding rate, std job finding rate, mean job separation rate
b1 = [mean_prod; std_prod; mean_unemployment; std_unemployment; mean_exit_rate; std_exit_rate; mean_job_destruction_rate] #simulated moments

db = (b1-b0)
distance = sum(transpose(db)*Weighting_matrix*db) #use the function "sum" because the output has to be a scalar, not a 1 times 1 matrix.

return distance
end

"""
using PyPlot, Distributions, StatsFuns, JLD, HDF5, DataFrames, Copulas, NLopt, BlackBoxOptim
path_main = "/home/julien/Final-Project-Julien-Pascal"
path_table = string(path_main,"/tables/")  #the path tables
path_surpluses = string(path_main,"/surpluses/")#path to precalculated objects, to speed up the process
file_random_numbers =  jldopen(string(path_surpluses,"random_numbers.jld"), "r")
random_numbers = read(file_random_numbers, "random_numbers")

#Weighting_matrix = eye(8,8)

#params_opt = [0.77, 0.023, 0.94, 0.99, 2.00, 5.56, 0.73, 0.046]

#["z0", "sigma", "lambda0", "eta", "mu", "x_lower_bar", "alpha"]
params_opt = [0.77, 0.023, 0.99, 2.00, 5.56, 0.73, 0.64]

#params_opt = [0.850196, 0.897381, 0.738404, 0.890517, 0.892101,  0.822158, 0.738211 ]

b0 = [1.0; 0.023;  0.059;  0.210; 0.761; 0.085; 0.046; 0.156]

params_opt = [0.77, 0.023, 0.99, 2.00, 5.56, 0.73, 0.64]
Weighting_matrix = eye(8,8)./(b0.*b0)

objective_function_2(params_opt, Weighting_matrix,random_numbers,b0)
"""


text:
using Test
using Flux
using Flux: DataLoader
using ParametrisedConvexApproximators
const PCA = ParametrisedConvexApproximators
using UnPack
using Transducers
using Convex
using Random
using Plots, StatsPlots
ENV["GKSwstype"]="nul"  # deactivate X server needs
using DataFrames
using Statistics
using LinearAlgebra
using SCS

__dir_save = "main/basic"

function target_function(x, u)
    n_x = length(x)
    n_u = length(u)
    [0.5 * ( -(1/n_x)*x'*x + (1/n_u)*u'*u )]
end

function sample(min, max)
    sampled = min + (max - min) .* rand(size(min)...)
end

function approximator_type(approximator)
    approx_type = nothing
    if typeof(approximator) == FNN
        approx_type = "FNN"
    elseif typeof(approximator) == MA
        approx_type = "MA"
    elseif typeof(approximator) == LSE
        approx_type = "LSE"
    elseif typeof(approximator) == PMA
        approx_type = "PMA"
    elseif typeof(approximator) == PLSE
        approx_type = "PLSE"
    elseif typeof(approximator) == PICNN
        approx_type = "PICNN"
    else
        error("Specify save directory")
    end
    approx_type
end

function infer_test(approximator)
    @unpack n, m = approximator
    # single data inference
    x = rand(n)
    u = rand(m)
    @test approximator(x, u) |> size  == (1,)
    # multiple data inference
    d = 10
    _xs = rand(n, d)
    _us = rand(m, d)
    @test approximator(_xs, _us) |> size  == (1, d)
end

function optimise_test(approximator, data, min_nt, max_nt)
    @unpack n, m = approximator
    d = data.x |> length
    println("optimise test; with $(d) test data")
    if typeof(approximator) <: ParametrisedConvexApproximator
        x = rand(n)
        u = Convex.Variable(m)
        @test approximator(x, u) |> size == (1,)  # inference with Convex.jl
    end
    _xs = hcat(data.x...)
    # begin  # TODO: remove it; only for test
    #     u = Convex.Variable(m)
    #     problem = minimize(approximator(_xs[:, 1], u)[1])
    #     problem.constraints += [u >= min_nt.u]
    #     problem.constraints += [u <= max_nt.u]
    #     @time solve!(problem, opt; verbose=true, silent_solver=false)
    # end
    println("Optimise $(d) points")
    @time res_timed = @timed optimise(approximator, _xs;
                                      u_min=min_nt.u, u_max=max_nt.u,
                                      multithreading=false,
                                     )  # with box constraints
    res = res_timed.value
    bchmkr = res_timed.time / d  # average time
    @test res.minimiser |> size == (m, d)  # optimise; minimiser
    @test res.optval |> size == (1, d)  # optimise; optval
    # compare true and estimated minimisers and optvals
    minimisers_true = 1:d |> Map(i -> zeros(m)) |> collect
    optvals_true = 1:d |> Map(i -> target_function(data.x[i], minimisers_true[i])) |> collect
    minimisers_estimated = 1:d |> Map(i -> res.minimiser[:, i]) |> collect
    optvals_estimated = 1:d |> Map(i -> res.optval[:, i]) |> collect
    minimisers_diff_norm = skipmissing(1:d |> Map(i -> minimisers_estimated[i] == [nothing] ? missing : norm(minimisers_estimated[i] - minimisers_true[i])) |> collect)  # make each element 'Number'
    optvals_diff_abs = skipmissing(1:d |> Map(i -> optvals_estimated[i] == [-Inf] || optvals_estimated[i] == [Inf] ? missing : abs((optvals_estimated[i] - optvals_true[i])[1])) |> collect)  # make each element 'Number'
    println("norm(estimated minimiser - true minimiser)'s mean (only for success cases): $(mean(minimisers_diff_norm))")
    println("norm(estimated optval - true optval)'s mean (only for success cases): $(mean(optvals_diff_abs))")
    _result = (;
               benchmark = bchmkr,
               minimisers_diff_norm=minimisers_diff_norm,
               optvals_diff_abs=optvals_diff_abs,
              )
    return _result
end

function training_test(approximator, data_train, data_test, epochs)
    _data_train = (;
                   x = hcat(data_train.x...),
                   u = hcat(data_train.u...),
                   f = hcat(data_train.f...),
                  )
    _data_test = (;
                   x = hcat(data_test.x...),
                   u = hcat(data_test.u...),
                   f = hcat(data_test.f...),
                  )
    loss(d) = Flux.Losses.mse(approximator(d.x, d.u), d.f)
    opt = ADAM(1e-3)
    ps = Flux.params(approximator)
    dataloader = DataLoader(_data_train; batchsize=16, shuffle=true, partial=false)
    println("Training $(epochs) epoch...")
    for epoch in 0:epochs
        if epoch != 0
            for d in dataloader
                train_loss, back = Flux.Zygote.pullback(() -> loss(d), ps)
                gs = back(one(train_loss))
                Flux.update!(opt, ps, gs)
                # projection for PICNN
                if typeof(approximator) == PICNN
                    for layer in approximator.NN.layers
                        if isdefined(layer, :Wz)
                            project_nonnegative!(layer.Wz)
                            if !all(layer.Wz .>= 0.0)
                                error("Projection seems not work")
                            end
                        end
                    end
                end
            end
        end
    end
end

function test_all(approximator, data, epochs, min_nt, max_nt)
    # split data
    xs_us_fs = zip(data.x, data.u, data.f) |> collect
    xs_us_fs_train, xs_us_fs_test = partitionTrainTest(xs_us_fs, 0.9)  # 90:10
    data_train = (;
                  x=xs_us_fs_train |> Map(xuf -> xuf[1]) |> collect,
                  u=xs_us_fs_train |> Map(xuf -> xuf[2]) |> collect,
                  f=xs_us_fs_train |> Map(xuf -> xuf[3]) |> collect,
                 )
    data_test = (;
                  x=xs_us_fs_test |> Map(xuf -> xuf[1]) |> collect,
                  u=xs_us_fs_test |> Map(xuf -> xuf[2]) |> collect,
                  f=xs_us_fs_test |> Map(xuf -> xuf[3]) |> collect,
                )
    @show typeof(approximator)
    @show number_of_parameters(approximator)
    # tests
    infer_test(approximator)
    @time training_test(approximator, data_train, data_test, epochs)
    _result = optimise_test(approximator, data_test, min_nt, max_nt)
    # figures
    fig_surface = nothing
    @unpack n, m = approximator
    approx_type = approximator_type(approximator)
    if n == 1 && m == 1
        fig_surface = plot_approximator(approximator, min_nt, max_nt)
    else
        println("plotting surface is ignored for high-dimensional cases")
        fig_surface = plot()
    end
    result = (; _result..., fig_surface=fig_surface, approx_type=approximator_type(approximator))
end

function plot_approximator(approximator, min_nt, max_nt; kwargs...)
    l = 20
    _xs = range(min_nt.x[1], stop=max_nt.x[1], length=l)
    _us = range(min_nt.u[1], stop=max_nt.u[1], length=l)
    _f(x, u) = approximator([x], [u])[1]
    fig = plot(_xs, _us, _f;
               st=:wireframe,
               xlim=(min_nt.x[1], max_nt.x[1]),
               ylim=(min_nt.u[1], max_nt.u[1]),
               zlim=(-1.0, 1.0),
               xticks=-0.5:0.5:1.0,
               aspect_ratio=:equal,
               camera=(60, 45),
               xlabel="x",
               ylabel="u",
               colorbar=false,
               xtickfontsize=13,
               ytickfontsize=13,
               ztickfontsize=13,
               xguidefontsize=13,
               yguidefontsize=13,
               zguidefontsize=13,
               legendfontsize=13,
               kwargs...
              )
    fig
end

"""
# References
[1] T. Bian and Z.-P. Jiang, “Value Iteration, Adaptive Dynamic Programming, and Optimal Control of Nonlinear Systems,” in 2016 IEEE 55th Conference on Decision and Control (CDC), Las Vegas, NV, USA, Dec. 2016, pp. 3375–3380. doi: 10.1109/CDC.2016.7798777.
[2] B. L. Stevens, F. L. Lewis, and E. N. Johnson, Aircraft Control and Simulation: Dynamics, Controls Design, and Autonomous Systems, Third edition. Hoboken, N.J: John Wiley & Sons, 2016.
[3] gym, "Humanoid-v2", https://gym.openai.com/envs/Humanoid-v2/
# dimensions
# n/m = 1/1: for visualisation, 13/4: [2, p.185, Table 3.5-2, F-16 Model Test Case], 376/17: [3]
"""

# @testset "basic" begin
# end

function main(n, m; epochs=100, seed=2021)
    @warn("Note: if you change the target function, you may have to change the true minimisers manually (in the function `optimise_test`).")
    df = DataFrame()
    # tests
    Random.seed!(seed)
    # training data
    d = 5_000
    println("No. of data points: $(d)")
    min_nt = (; x = -1*ones(n), u = -1*ones(m))
    max_nt = (; x = 1*ones(n), u = 1*ones(m))
    xs = 1:d |> Map(i -> sample(min_nt.x, max_nt.x)) |> collect
    us = 1:d |> Map(i -> sample(min_nt.u, max_nt.u)) |> collect
    fs = zip(xs, us) |> MapSplat((x, u) -> target_function(x, u)) |> collect
    data = (; x=xs, u=us, f=fs)
    println("#"^10 * " " * "n = $(n), m = $(m), epochs = $(epochs)" * " " * "#"^10)
    i_max = 30
    T = 1e-1
    h_array = [64, 64]  # fix it
    z_array = [64, 64]  # for PICNN
    h_array_fnn = [64, 64]
    multiplication_factor = 1
    u_array_0 = 64
    u_array = vcat(u_array_0, z_array...)  # for PICNN; length(u_array) != length(z_array) + 1
    act = Flux.leakyrelu
    α_is = 1:i_max*multiplication_factor |> Map(i -> Flux.glorot_uniform(n+m)) |> collect
    β_is = 1:i_max*multiplication_factor |> Map(i -> Flux.glorot_uniform(1)) |> collect
    # generate approximators
    fnn = FNN(n, m, h_array_fnn, act)
    ma = MA(α_is, β_is; n=n, m=m)
    lse = LSE(α_is, β_is, T; n=n, m=m)
    picnn = PICNN(n, m, u_array, z_array, act, act)
    pma = PMA(n, m, i_max, h_array, act)
    plse = PLSE(n, m, i_max, T, h_array, act)
    approximators = (;
                     fnn=fnn,
                     ma=ma,
                     lse=lse,
                     picnn=picnn,
                     pma=pma,
                     plse=plse,
                    )
    _dir_save = joinpath(__dir_save, "n=$(n)_m=$(m)_epochs=$(epochs)")
    mkpath(_dir_save)
    results = approximators |> Map(approximator -> test_all(approximator, data, epochs, min_nt, max_nt)) |> collect
    for (approximator, result) in zip(approximators, results)
        @unpack benchmark, minimisers_diff_norm, optvals_diff_abs, = result
        push!(df, (;
                   n=n, m=m, epochs=epochs,
                   approximator=approximator_type(approximator),
                   optimise_time_mean=benchmark |> mean,
                   minimisers_diff_norm_mean=minimisers_diff_norm |> mean,
                   optvals_diff_abs_mean=optvals_diff_abs |> mean,
                   no_of_minimiser_success=minimisers_diff_norm |> collect |> length,
                   no_of_optval_success=optvals_diff_abs |> collect |> length,
                   number_of_parameters=approximator |> number_of_parameters,
                  );
              cols=:union,
              promote=true,
             )
    end
    # violin plots
    minimisers_diff_norm = results |> Map(result -> result.minimisers_diff_norm) |> collect  # skipmissing's
    optvals_diff_abs = results |> Map(result -> result.optvals_diff_abs) |> collect  # skipmissing's
    approx_types = approximators |> Map(approximator_type) |> collect
    box_minimiser_diff_norm = plot()
    for (approx_type, minimiser_diff_norm) in zip(approx_types, minimisers_diff_norm)
        violin!(box_minimiser_diff_norm,
                [approx_type], minimiser_diff_norm |> collect;  # remove missing's
                legend=nothing,
                xtickfontsize=13,
                ytickfontsize=13,
                xguidefontsize=13,
                yguidefontsize=13,
                legendfontsize=13,
               )
    end
    savefig(box_minimiser_diff_norm, joinpath(_dir_save, "minimiser_diff_norm.png"))
    savefig(box_minimiser_diff_norm, joinpath(_dir_save, "minimiser_diff_norm.pdf"))
    box_optval_diff_abs = plot()
    for (approx_type, optval_diff_abs) in zip(approx_types, optvals_diff_abs)
        violin!(box_optval_diff_abs,
                [approx_type], optval_diff_abs |> collect;  # remove missing's
                legend=nothing,
                xtickfontsize=13,
                ytickfontsize=13,
                xguidefontsize=13,
                yguidefontsize=13,
                legendfontsize=13,
               )
    end
    savefig(box_optval_diff_abs, joinpath(_dir_save, "optval_diff_abs.png"))
    savefig(box_optval_diff_abs, joinpath(_dir_save, "optval_diff_abs.pdf"))
    # plotting
    if n == 1 && m == 1
        for result in results
            @unpack fig_surface, approx_type = result
            savefig(fig_surface, joinpath(_dir_save, "surface_" * approx_type * ".png"))
            savefig(fig_surface, joinpath(_dir_save, "surface_" * approx_type * ".pdf"))
        end
        fig_surface_true = plot_approximator(target_function, min_nt, max_nt)
        savefig(fig_surface_true, joinpath(_dir_save, "surface_true.png"))
        savefig(fig_surface_true, joinpath(_dir_save, "surface_true.pdf"))
    end
    title_minimiser_diff_norm = plot(title="2-norm of minimiser errors",
                                     framestyle=nothing,showaxis=false,xticks=false,yticks=false,margin=0Plots.px,
                                    )
    @show df
end

ext: jl


text:
function example2(n)
    h = function hello()
        println(n)
    end
    h
end
ourfunc = example2(5)
ourfunc()


text:
function kinshipcoef(father, mother, start, finish)
#
# Author: Kenneth Lange
#
# This function computes the kinship matrix of a pedigree.
# Children must appear after parents. The pedigree occurs
# in a block and is defined implicitly by giving fathers
# and mothers. Person i is a pedigree founder if and only if
# father[i] = 0.
#
   ped_size = finish - start + 1
   kinship = zeros(ped_size, ped_size)
   q = start - 1
   for i = 1 : ped_size
      if father[i + q] == 0
         kinship[i, i] = 0.5
      else
         j = father[i + q] - q
         k = mother[i + q] - q
         kinship[i, i] = 0.5 + 0.5 * kinship[j, k]
         for m = 1 : i - 1
            kinship[i, m] = 0.5 * (kinship[j, m] + kinship[k, m])
            kinship[m, i] = kinship[i, m]
         end
      end
   end
   return kinship
end


text:
function couple_terms(t1::Term, t2::Term)
    # It is assumed that t1 and t2 originate from non-equivalent
    # electrons, since the vector model does not predict correct term
    # couplings for equivalent electrons; some of the generated terms
    # would violate the Pauli principle; cf. Cowan p. 108–109.
    L1 = t1.L
    L2 = t2.L
    S1 = t1.S
    S2 = t2.S
    p = t1.parity * t2.parity
    sort(vcat([[Term(L, S, p) for S in abs(S1-S2):(S1+S2)]
               for L in abs(L1-L2):(L1+L2)]...))
end

function couple_terms(t1s::Vector{<:Term}, t2s::Vector{<:Term})
    ts = map(t1s) do t1
        map(t2s) do t2
            couple_terms(t1, t2)
        end
    end
    sort(unique(vcat(vcat(ts...)...)))
end

final_terms(ts::Vector{<:Vector{<:T}}) where {T<:Union{Term,Real}} =
    foldl(couple_terms, ts)

couple_terms(J1::T, J2::T) where {T <: Union{Integer,HalfInteger}} =
    collect(abs(J1-J2):(J1+J2))
couple_terms(J1::Real, J2::Real) =
    couple_terms(convert(HalfInteger, J1), convert(HalfInteger, J2))

term_type(::Type{IntermediateTerm}) = Term
term_type(::Type{T}) where {T<:Real} = T

function intermediate_couplings(its::Vector{T}, t₀::T=zero(T)) where {T<:Union{Term,Integer,HalfInteger}}
    ts = Vector{Vector{T}}()
    for t in couple_terms(t₀, its[1])
        if length(its) == 1
            push!(ts, [t₀, t])
        else
            for ts′ in intermediate_couplings(its[2:end], t)
                push!(ts, vcat(t₀, ts′...))
            end
        end
    end
    ts
end

intermediate_couplings(its::Vector{IntermediateTerm}, t₀::Term=zero(Term)) =
    intermediate_couplings(map(t -> t.term, its), t₀)

intermediate_couplings(J::Vector{IT}, j₀::T=zero(T)) where {IT <: Real, T <: Real} =
    intermediate_couplings(convert.(HalfInteger, J), convert(HalfInteger, j₀))

export couple_terms, intermediate_couplings


text:
"""
    readSRF(fname::AbstractString)

Import surface data in the BESA loc format

fname   = location and name of file to import
verbose = should the function output user feedback

surface = readSRF(verbose=false)
"""
function readSRF(fname::AbstractString = Pkg.dir("Private", "test", "data", "Default50Brain.srf"); verbose::Bool=false)

    debug("Surface file: $(fname)")

    # Open file
    fid = open(fname, "r")

    # Read header information from file
    version         = read(fid, Float32, 1)[1]
    reserved        = read(fid, Int32, 1)[1]
    numVertices     = read(fid, Int32, 1)[1]
    numTriangles    = read(fid, Int32, 1)[1]
    MeshCenterX     = read(fid, Float32, 1)[1]
    MeshCenterY     = read(fid, Float32, 1)[1]
    MeshCenterZ     = read(fid, Float32, 1)[1]
    Xcords          = read(fid, Float32, numVertices)
    Ycords          = read(fid, Float32, numVertices)
    Zcords          = read(fid, Float32, numVertices)
    Xcomps          = read(fid, Float32, numVertices)
    Ycomps          = read(fid, Float32, numVertices)
    Zcomps          = read(fid, Float32, numVertices)
    rConvColor      = read(fid, Float32, 1)[1]
    gConvColor      = read(fid, Float32, 1)[1]
    bConvColor      = read(fid, Float32, 1)[1]
    aConvColor      = read(fid, Float32, 1)[1]
    rConcColor      = read(fid, Float32, 1)[1]
    gConcColor      = read(fid, Float32, 1)[1]
    bConcColor      = read(fid, Float32, 1)[1]
    aConcColor      = read(fid, Float32, 1)[1]
    meshColor       = read(fid, Int32, numVertices)

    if verbose
        println("Version:            $(version)")
        println("Reserved:           $(reserved)")
        println("Vertices #:         $(numVertices)")
        println("Triangles #:        $(numTriangles)")

    end

    i = 1
    while i <= numVertices
        number = read(fid, Int32, 1)[1]
        # Just overwriting. Not saving data. TODO: Fix this
        neighbour = read(fid, Int32, number)
        i = i + 1
    end

    Triangles = read(fid, Int32, 3*numTriangles)
    numTriStripElements = read(fid, Int32, 1)[1]
    seqStripElements = read(fid, Int32, numTriStripElements)

    close(fid)

    return Xcords, Ycords, Zcords

end


text:
# Used in finite difference derivatives
function compute_timestep(params, t0, points, h, ib, it)
    rstar = params[end]
    xc = components(points[ib,it,:,1]./rstar, h)
    yc = components(points[ib,it,:,2]./rstar, h)
    trans = transit_init(params[ib - 1], zero(typeof(t0)), params[end - 2:end-1], false)
    return compute_flux(t0, t0, xc, yc, trans)
end

# Used in finite difference derivatives
function compute_timestep_nbody(coords, s_copy, tt, trans, ic, intr, it, ib)
    s = deepcopy(s_copy)
    coords = reshape(coords, 7, s.n)
    for (i, col) in enumerate(eachcol(coords))
        s.x[:,i] .= col[1:3]
        s.v[:,i] .= col[4:6]
        s.m[i]    = col[end]
    end

    pd = compute_pd(s, ic, tt, intr, grad=false);
    rstar = big(0.00465047 * 0.1192) # Trappist-1 (Rstar/AU)
    normalize_points!(pd.points, rstar);

    t0 = pd.times[it]
    xc = components(pd.points[ib,it,:,1], pd.h)
    yc = components(pd.points[ib,it,:,2], pd.h)

    return compute_flux(t0, t0, xc, yc, trans)
end

function test_compute_timestep_derivatives(n)
    # Setup initial conditions and compute points
    BJD = 7250.0; t0_ic = 7.0
    tmax = 2.0
    k = get_radius_ratios_trappist(n);
    u_n = get_limdark_coeffs_trappist();
    ic = setup_ICs(n, BJD, t0_ic);
    intr = setup_integrator(ic, tmax);
    tt = compute_transit_times(ic, intr);
    pd = compute_pd(ic, tt, intr, grad=true);

    # Normalize points to stellar radius
    rstar = 0.00465047 * 0.1192 # Trappist-1 (Rstar/AU)
    normalize_points!(pd.points, rstar);
    normalize_points!(pd.dpoints, rstar);

    for it in 1:length(pd.times)
        ib = pd.bodies[it]
        t0 = pd.times[it]

        # Get a transit structure
        trans = transit_init(k[ib - 1], 0.0, u_n, false)

        # Compute x and y components for impact parameter calculation
        xc = components(pd.points[ib,it,:,1], pd.h)
        yc = components(pd.points[ib,it,:,2], pd.h)

        # Compute the flux at t0
        tol = 1e-15
        maxdepth = 40
        flux = compute_flux(t0, t0, xc, yc, trans)

        # Make bigfloat variants
        ic_big = setup_ICs(n, big(BJD), big(t0_ic));
        intr_big = setup_integrator(ic_big, big(tmax));
        tt_big = compute_transit_times(ic_big, intr_big);
        pd_big = compute_pd(ic_big, tt_big, intr_big, grad=false);

        # Make sure we get the same answer
        params = [k...,u_n...,rstar]
        flux_big = compute_timestep(big.(params), big(t0), copy(pd_big.points), pd_big.h, ib, it)
        @test flux_big ≈ flux

        # Compute finite diff derivatives with respect to the transit parameters
        grad_num = grad(central_fdm(5, 1), p -> compute_timestep(p, big(t0), copy(pd_big.points), pd_big.h, ib, it), big.(params))[1]

        # Compute analytic derivatives
        inds = [1,2,3,4,5,6,7]
        dxc = [components(pd.dpoints[ib,it,:,1,k,i], pd.h) for i in inds, k in 1:ic.nbody][:]
        dyc = [components(pd.dpoints[ib,it,:,2,k,i], pd.h) for i in inds, k in 1:ic.nbody][:]
        dbdq0 = zeros(7 * ic.nbody)
        ki = ib - 1
        trans_grad = transit_init(k[ki], 0.0, u_n, true)
        n_params = length(dbdq0) + length(k) + length(u_n) + 1 # flux
        lc = Lightcurve(0.0, [0.0], [0.0], [0.0], u_n, k, rstar, 7*(length(k)+1))
        compute_flux!(t0, t0, xc, yc, dxc, dyc, lc, trans_grad, 1, ki, inv(rstar))

        # Check flux again
        @test lc.flux[1] ≈ flux

        # Check derivative wrt radius ratios
        @test isapprox(lc.dfdk[1,:], Float64.(grad_num[1:n-1]), norm=x -> maximum(abs.(x)))

        # Now the derivatives wrt limbdark coefficients
        @test isapprox(lc.dfdu[1,:], Float64.(grad_num[n:end-1]), norm=x -> maximum(abs.(x)))

        # Now the derivatives wrt the stellar radius
        @test isapprox(lc.dfdr[1], Float64.(grad_num[end]), norm=x -> maximum(abs(x)), atol=1e-8)

        # Finally, derivatives wrt the Nbody initial conditions
        s_copy = deepcopy(State(ic_big))
        params_nbody = vcat([[x...,v...,m] for (x, v, m) in zip(eachcol(s_copy.x), eachcol(s_copy.v), s_copy.m)]...)

        # Again, check flux
        trans = transit_init(k[ib - 1], 0.0, u_n, false)
        flux_nbody = compute_timestep_nbody(params_nbody, State(ic), tt, trans, deepcopy(ic), intr, it, ib)
        @test flux_nbody ≈ flux

        # Compute finite FiniteDifferences
        trans_big = transit_init(big(k[ki]), big(0.0), big.(u_n), false)
        grad_nbody_num = grad(central_fdm(12, 1), p -> compute_timestep_nbody(p, s_copy, tt_big, trans_big, ic_big, intr_big, it, ib), big.(params_nbody))[1]

        # Compare with analytic
        @test isapprox(lc.dfdq0[1,:], Float64.(grad_nbody_num), norm=x -> maximum(abs.(x)))
    end
end

@testset "Non-Integrated Flux" begin
    test_compute_timestep_derivatives(3) # This runs for a while...
end


text:
module GuiData

using Gtk

import ..GtkUtilities.Link: AbstractState, disconnect, set!, set_quietly!

export guidata, trigger

struct GUIData{T}
    data::Dict{T,Dict{Symbol,Any}}
end

GUIData() = GUIData(Dict{Any,Dict{Symbol,Any}}())

const empty_guidata = Dict{Symbol,Any}()

# guidata[w, :name] strips the State wrapper, unless requested by keyword
# Likewise, guidata[w, :name] = val set!s the State rather than replacing it,
# unless requested by keyword.

Base.getindex(wd::GUIData, w) = wd.data[w]
function Base.getindex(wd::GUIData, w, s::Symbol; raw::Bool=false)
    ret = wd.data[w][s]
    if !raw && isa(ret, AbstractState)
        return get(ret)
    end
    ret
end
Base.getindex(wd::GUIData{T}, ws::Tuple{T,Symbol}) where T = wd[ws[1], ws[2]]

function Base.setindex!(wd::GUIData, val, w, s::Symbol; raw::Bool=false)
    d = get(wd.data, w, empty_guidata)
    if d == empty_guidata
        d = wd.data[w] = Dict{Symbol,Any}()
        if isa(w, Gtk.GtkWidget)
            signal_connect(destroy_callback, w, "destroy", Nothing, (), false, (w,wd))
        end
    end
    if !raw
        obj = get(d, s, nothing)
        if isa(obj, AbstractState)
            set!(obj, val)
            return val
        end
    end
    d[s] = val
end

@guarded function destroy_callback(widgetptr::Ptr, user_data)
    widget, wd = user_data
    dct = wd.data[widget]
    for (k,v) in dct
        if isa(v, AbstractState)
            disconnect(v)
        end
    end
    delete!(wd.data, widget)
    nothing
end

Base.setindex!(wd::GUIData, val, ws::Tuple{T,Symbol}) where T = wd.data[ws[1],ws[2]] = val

Base.haskey(wd::GUIData, key) = haskey(wd.data, key)

function Base.get(wd::GUIData, ws::Tuple{T,Symbol}, default; raw::Bool=false) where T
    d = get(wd.data, ws[1], empty_guidata)
    val = d == empty_guidata ? default : get(d, ws[2], default)
    if !raw && isa(val, AbstractState)
        return get(val)
    end
    val
end

function Base.delete!(wd::GUIData, ws::Tuple{T,Symbol}) where T
    d = get(wd.data, ws[1], empty_guidata)
    d == empty_guidata ? wd : (delete!(d, ws[2]); wd)
end

Base.delete!(wd::GUIData, w) = delete!(wd.data, w)
Base.show(io::IO, wd::GUIData) = print(io, "GUIdata")

const guidata = GUIData()

"""
`trigger(widgets, symbols)` is used when you want to set several
state variables simultaneously, but don't want to refresh the screen
more frequently than necessary. You can set the `.value` parameter of
the state variables directly, then call `trigger` to synchronize the
GUI.
"""
function trigger(widgets, syms)
    dct = Dict{Gtk.GtkWidget,Any}()
    canvases = Set{Gtk.GtkCanvas}()
    for widget in widgets, sym in syms
        state = get(guidata, (widget, sym), nothing; raw=true)
        if state != nothing
            for w in state.widgets
                dct[w] = get(state)
            end
            for c in state.canvases
                push!(canvases, c)
            end
        end
    end
    for (w,val) in dct
        set_quietly!(w, val)
    end
    for c in canvases
        draw(c)
    end
end

trigger(widget::Gtk.GtkWidget, sym::Symbol) = trigger((widget,), (sym,))
trigger(widget::Gtk.GtkWidget, syms)        = trigger((widget,), syms)
trigger(widgets,               sym::Symbol) = trigger(widgets, (sym,))

"""
Given a widget (Button, Canvas, Window, etc.) or other graphical object
`w`, a value `val` can be associated with ("stored in") `w` using
```
guidata[w, :name] = val
```
where `:name` is the name (a Symbol) you've assigned to `val` for the
purposes of storage.

The value can be retrieved with
```
val = guidata[w, :name]
```
Here are some other things you can do with `guidata`:
```
alldata = guidata[w]           # fetch all data associated with w
val = get(guidata, (w,:name), default)   # returns default if :name not defined
delete!(guidata, (w,:name))    # deletes the value associated with :name
delete!(guidata, w)            # deletes all data associated with w

Example:

    c = Canvas()
    panzoom(c, (0,1), (0,1))
    xview = guidata[c, :xview]
```

Note that if `:name` corresponds to a `State` object, `guidata` will
get/set the **value** of the State object, transparently stripping the
wrapper itself.  This means that you change the value without
destroying any links you have set up.  If you need the state object
itself, use the `raw` keyword:
```jl
    state = getindex(guidata, c, :name; raw=true)     # retrieves a state object
    setindex!(guidata, newstate, c, :name; raw=true)  # replaces the old state object
```
"""
guidata

end


text:
module MatMul

using Random
using CachedArrays

collection(f, num, sz) = [f(randn(Float32, sz, sz)) for _ in 1:num]

go(num, sz; kw...) = go(identity, num, sz; kw...)

function go(f, num, sz; iters = 1000)
    Random.seed!(123)
    matrices = collection(f, num, sz)

    # Show the total size of all the matrices.
    @show sum(sizeof, matrices)

    @time for _ in 1:iters
        x = rand(matrices) * rand(matrices)
        i = rand(1:length(matrices))
        matrices[i] = x
    end
    return matrices
end

function go_inplace(f, num, sz; iters = 1000)
    Random.seed!(123)
    matrices = collection(f, num, sz)

    # Show the total size of all the matrices.
    @show sum(sizeof, matrices)
    @time for _ in 1:iters
        i = rand(1:length(matrices))
        matrices[i] .= rand(matrices) * rand(matrices)
    end
    return matrices
end

#####
##### Record tests here because of the ever increasing number of steps required.
#####

function largish_test()
    GC.gc()
    manager = CachedArrays.GlobalManager[]
    resize!(manager, 20_000_000_000)
    CachedArrays.resize_remote!(manager, 20_000_000_000)
    manager.flushpercent = 0.8

    return (
        () -> go_inplace(identity, 10, 10000; iters = 30),
        () -> go_inplace(LockedCachedArray, 10, 10000; iters = 30),
    )
end

end


text:
export checkpoint!, restoreModel!, continueSim!, find_t_for_tr

function checkpoint!(M)
    for s in [:pars, :init, :spec]
        M.originalState[s] = deepcopy(M.(s))
    end
    return(M)
end

function restoreModel!(M)
    for s in [:pars, :init, :spec]
        M.(s) = deepcopy(M.originalState[s])
    end
    return(M)
end

#function to find t corresponding to desired tr
function find_t_for_tr(M, simname, query_time)
    T = M.sims[simname].D["t"]
    TR = M.sims[simname].D["tr"]
    total = T[findfirst([tr >= query_time for tr in TR])]
    return(total)
end

function continueSim!(M, t1, t2, name, pars2; inits2 = false)
    # M.spec["total"] = t1
    n1 = name * "_1"
    n2 = name * "_2"
    runSimulation!(M, n1)
    ix = findfirst([t >= t1 for t in M.sims[n1].D["t"]])
    M.init = Dict([v => M.sims[n1].D[v][ix] for v in keys(M.odes)])
    if inits2 != false
        for k in keys(inits2)
            M.init[k] = inits2[k]
        end
    end
    M.pars = pars2
    M.spec["total"] = t2
    runSimulation!(M, n2)
    M.sims[n2].D["t"] = M.sims[n2].D["t"] .+ t1
    M.sims[name] = SimulationData(M, name)
    for v in keys(M.sims[name].D)
        append!(M.sims[name].D[v], M.sims[n1].D[v][1:ix])
        append!(M.sims[name].D[v], M.sims[n2].D[v])
    end
    restoreModel!(M)
    return(M)
end

max_stars_repo_name: novakgroupoxford/XPPjl

max_issues_repo_name: novakgroupoxford/XPPjl

max_forks_repo_name: novakgroupoxford/XPPjl


text:
using FEMAssembler

## scalar parameters
E  = 2e11  # Young modulus in Pa
A  = 5e-3  # Cross-section area in m^2
d  = 1.0   # horizontal distance in m
h  = 1.0   # height in m
Fx = 0     # horizontal load in N
Fy = -1e3  # vertical   load in N

## set structs
E1 = E
E2 = E

steel1  = materialsData( E1, 0.3 );
steel2  = materialsData( E2, 0.3 );
section = elementsData("truss", A);

# The coordinate matrix is given by
nodalCoords = [  0.  0.  0. ;
                 d   0.  h  ;
                2d   0.  0. ]

# matrix with MEBI indexes of each element (on each row)
# MEBI = Materials / Elements / BoundaryConditions / InitialConditions
MEBIValsMat = [ 0 1 1 0  ; # no material / first element / first BC / no IC
0 1 2 0  ;
1 2 0 0  ;
2 2 0 0 ] ;

# the connectivity matrix is given by a vector of vectors.
elemNodalConnec = [ [ 1 ],  [ 2 ], [ 3 ], [ 1, 2],  [ 2, 3] ]
MEBIVec         = [     1,      2,     1,       3,       4  ]

BCs = [ boundaryCondsData([ 1, 3, 5],[0, 0, 0],[     ], [ ]       ) ,
        boundaryCondsData([ 3    ]  ,[0      ],[1, 5 ], [Fx, Fy ] ) ]

# assemble the sitiffness equation
KG, FG, Kmatrices, neumDofs = assembler( nodalCoords, elemNodalConnec, MEBIValsMat, MEBIVec, [steel1, steel2], [[],section], BCs, [] ) ;

print("KG:\n")
display(KG)
display(FG)
print("neumDofs", neumDofs)
print("\n")

KGred = copy( KG )
FGred = copy( FG )

KGred = KGred[ neumDofs, neumDofs ]
FGred = FGred[ neumDofs           ]

print("KGred", KGred,"\n")
# the system is solved.
UGred = KGred \ FGred

UG = zeros( size( FG ) )
UG[  neumDofs ] = UGred


text:
"""
	magnitude_to_moment(m::T) where T<:Real

Converts moment magnitude to seismic moment (in dyne-cm).

# Examples
```julia-repl
	m = 6.0
	M0 = magnitude_to_moment(m)
```
"""
function magnitude_to_moment(m::T) where T<:Real
	# equivalent to: return 10.0^( 1.5*( m + 10.7 ) )
    return 10.0^( 1.5m + 16.05 )
end

"""
	corner_frequency_brune(m::S, Δσ::T, β::Float64=3.5) where {S<:Real, T<:Real}

Computes the corner frequency using the Brune model.
	- `m` is the moment magnitude
	- `Δσ` is the stress drop in units of bars
	- `β` is the shear-wave velocity at the source in units of km/s

# Examples
```julia-repl
	m = 6.0
	Δσ = 100.0
	β = 3.5
	fc = corner_frequency_brune(m, Δσ, β)
```
"""
function corner_frequency_brune(m::S, Δσ::T, β::Float64=3.5) where {S<:Real, T<:Real}
    Mo = magnitude_to_moment(m)
    return 4.9058e6 * β * ( Δσ / Mo )^(1/3)
end

"""
    corner_frequency_atkinson_silva_2000(m::T) where T<:Real

Computes the corner frequencies, `fa`, `fb`, and the mixing parameter `ε` from the Atkinson & Silva (2000) double corner frequency model.
This is the default source corner frequency model used by Boore & Thompson (2014) to define their source duration. But note that they just use fa. This function returns fb and ε also

# Examples
```julia-repl
    m = 6.0
    fa, fb, ε = corner_frequency_atkinson_silva_2000(m)
```

See also: [`corner_frequency`](@ref)
"""
function corner_frequency_atkinson_silva_2000(m::T) where T<:Real
	fa = 10.0^( 2.181 - 0.496*m )
	fb = 10.0^( 2.410 - 0.408*m )
	ε = 10.0^( 0.605 - 0.255*m )
    return fa, fb, ε
end

"""
    corner_frequency(m::U, src::SourceParameters{S,T}) where {S<:Float64, T<:Real, U<:Real}

Computes a 3-tuple of corner frequency components, depending upon source spectrum type.
By default the single-corner Brune spectrum is considered, but if `src.model` equals `:Atkinson_Silva_2000` then the components of the double-corner spectrum are returned.
If some other symbol is passed then the Brune model is returned.

# Examples
```julia-repl
    m = 6.0
    Δσ = 100.0
    κ0 = 0.035
    fas = FASParams(Δσ, κ0)
    # compute single corner frequency
	src.model = :Brune
    fc, tmp1, tmp2 = corner_frequency(m, src)
    # compute double corner frequencies
	src.model = :Atkinson_Silva_2000
    fa, fb, ε = corner_frequency(m, src)
```
"""
function corner_frequency(m::U, src::SourceParameters{S,T}) where {S<:Float64, T<:Real, U<:Real}
    if src.model == :Atkinson_Silva_2000
        fa, fb, ε = corner_frequency_atkinson_silva_2000(m)
		if T <: Dual
			return T(fa), T(fb), T(ε)
		else
			return fa, fb, ε
		end
    else
		# default to Brune
		fc = corner_frequency_brune(m, src.Δσ, src.β)
		V = typeof(fc)
		return fc, V(NaN), V(NaN)
    end
end

corner_frequency(m, fas::FourierParameters) = corner_frequency(m, fas.source)


text:
using Modia
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@testset "RunSimulations" begin

    Modia.ModiaLogging.resetTestStatus()

    include("models/TestVariableTypes.jl")
    include("models/TestEquations.jl")

    include("models/TestElectrical.jl")
    include("models/TestFilter.jl")
    include("models/TestArrayOfComponents.jl")
    include("models/TestConditionalComponents.jl")
    include("models/TestConditionalEquations.jl")

    include("models/TestSpatialDiscretization.jl")

    #include("models/TestCoupledInertias.jl")
    #include("models/TestPendulum.jl")
    #include("models/TestStateSelection.jl")
    #include("models/TestRetranslationOmega.jl")

    Modia.ModiaLogging.printTestStatus()

end
nothing


text:
using Oceananigans.Utils: with_tracers
using Oceananigans.TurbulenceClosures: with_tracers

"""
    SolutionBoundaryConditions(tracers, proposal_bcs)

Construct a `NamedTuple` of `FieldBoundaryConditions` for a model with
fields `u`, `v`, `w`, and `tracers` from the proposal boundary conditions
`proposal_bcs`, which must contain the boundary conditions on `u`, `v`, and `w`
and may contain some or all of the boundary conditions on `tracers`.
"""
SolutionBoundaryConditions(tracers, proposal_bcs) =
    with_tracers(tracers, proposal_bcs, default_tracer_bcs, with_velocities=true)

"""
    HorizontallyPeriodicSolutionBCs(u=HorizontallyPeriodicBCs(), ...)

Construct `SolutionBoundaryConditions` for a horizontally-periodic model
configuration with solution fields `u`, `v`, `w`, `T`, and `S` specified by keyword arguments.

By default `HorizontallyPeriodicBCs` are applied to `u`, `v`, `T`, and `S`
and `HorizontallyPeriodicBCs(top=NoPenetrationBC(), bottom=NoPenetrationBC())` is applied to `w`.

Use `HorizontallyPeriodicBCs` when constructing non-default boundary conditions for `u`, `v`, `w`, `T`, `S`.
"""
function HorizontallyPeriodicSolutionBCs(;
    u = HorizontallyPeriodicBCs(),
    v = HorizontallyPeriodicBCs(),
    w = HorizontallyPeriodicBCs(top=NoPenetrationBC(), bottom=NoPenetrationBC()),
    tracerbcs...)

    return merge((u=u, v=v, w=w), tracerbcs)
end

"""
    ChannelSolutionBCs(u=ChannelBCs(), ...)

Construct `SolutionBoundaryConditions` for a reentrant channel model
configuration with solution fields `u`, `v`, `w`, `T`, and `S` specified by keyword arguments.

By default `ChannelBCs` are applied to `u`, `v`, `T`, and `S`
and `ChannelBCs(top=NoPenetrationBC(), bottom=NoPenetrationBC())` is applied to `w`.

Use `ChannelBCs` when constructing non-default boundary conditions for `u`, `v`, `w`, `T`, `S`.
"""
function ChannelSolutionBCs(;
    u = ChannelBCs(),
    v = ChannelBCs(north=NoPenetrationBC(), south=NoPenetrationBC()),
    w = ChannelBCs(top=NoPenetrationBC(), bottom=NoPenetrationBC()),
    tracerbcs...)

    return merge((u=u, v=v, w=w), tracerbcs)
end


text:
module adapters_odbc_mysql_repo_test

using Test
using Octo.Adapters.ODBC

Repo.debug_sql()

#=
Repo.connect(
    adapter  = Octo.Adapters.ODBC,
    database = Octo.DBMS.MySQL,
    dsn      = "MySQL-test",
    username = "root",
    password = "",
)

Repo.disconnect()
=#

end # module adapters_odbc_mysql_repo_test


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: cost_and_usage_report_service
using AWS.Compat
using AWS.UUIDs

"""
    delete_report_definition()
    delete_report_definition(params::Dict{String,<:Any})

Deletes the specified report.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ReportName"`: The name of the report that you want to delete. The name must be unique,
  is case sensitive, and can't include spaces.
"""
delete_report_definition(; aws_config::AbstractAWSConfig=global_aws_config()) = cost_and_usage_report_service("DeleteReportDefinition"; aws_config=aws_config)
delete_report_definition(params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = cost_and_usage_report_service("DeleteReportDefinition", params; aws_config=aws_config)

"""
    describe_report_definitions()
    describe_report_definitions(params::Dict{String,<:Any})

Lists the AWS Cost and Usage reports available to this account.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`:
- `"NextToken"`:
"""
describe_report_definitions(; aws_config::AbstractAWSConfig=global_aws_config()) = cost_and_usage_report_service("DescribeReportDefinitions"; aws_config=aws_config)
describe_report_definitions(params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = cost_and_usage_report_service("DescribeReportDefinitions", params; aws_config=aws_config)

"""
    modify_report_definition(report_definition, report_name)
    modify_report_definition(report_definition, report_name, params::Dict{String,<:Any})

Allows you to programatically update your report preferences.

# Arguments
- `report_definition`:
- `report_name`:

"""
modify_report_definition(ReportDefinition, ReportName; aws_config::AbstractAWSConfig=global_aws_config()) = cost_and_usage_report_service("ModifyReportDefinition", Dict{String, Any}("ReportDefinition"=>ReportDefinition, "ReportName"=>ReportName); aws_config=aws_config)
modify_report_definition(ReportDefinition, ReportName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = cost_and_usage_report_service("ModifyReportDefinition", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ReportDefinition"=>ReportDefinition, "ReportName"=>ReportName), params)); aws_config=aws_config)

"""
    put_report_definition(report_definition)
    put_report_definition(report_definition, params::Dict{String,<:Any})

Creates a new report using the description that you provide.

# Arguments
- `report_definition`: Represents the output of the PutReportDefinition operation. The
  content consists of the detailed metadata and data file information.

"""
put_report_definition(ReportDefinition; aws_config::AbstractAWSConfig=global_aws_config()) = cost_and_usage_report_service("PutReportDefinition", Dict{String, Any}("ReportDefinition"=>ReportDefinition); aws_config=aws_config)
put_report_definition(ReportDefinition, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = cost_and_usage_report_service("PutReportDefinition", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ReportDefinition"=>ReportDefinition), params)); aws_config=aws_config)


text:
@symbol_func function R_CS(cur_reactor::AbstractReactor)
  cur_R = safe_symbol(cur_reactor, :R_0)

  cur_R -= a(cur_reactor)

  cur_R -= b(cur_reactor)

  cur_R -= c(cur_reactor)

  cur_R -= K_DO(cur_reactor)

  cur_R /= ( 1 + K_DR(cur_reactor) )

  cur_R
end


text:
using ELabFTW
using Test

# @testset "GET request" begin
#     api = ELabFTW.BASE_API

#     r = ELabFTW.send_request(api, "experiments/";
#         params=Dict("limit" => 15, "offset" => 0)
#     )

#     @show r
# end

@testset "Methods" begin
    experiments = get_all_experiments()

    id = first(experiments).id

    experiment = get_experiment(id)

    @show experiment

    # @show first(experiments)

    # experiments_lim = get_all_experiments(;limit=1)
    # @test length(experiments_lim) == 1

    # experiments_off = get_all_experiments(;limit=1, offset=1)
    # @test length(experiments_off) == 1
    # @test first(experiments_lim).id > first(experiments_off).id
end


text:
module HistogramPlot

using StatsBase, PyPlot, PyCall, HEP

#draws a familiar 1-D histogram with errors
#ax - PyPlot axes
#hh - Histogram
#kwords:
#do_error - default true, false to skip drawing errorbars
#scaling - multiply bar heights and errors ny this number

function barplot(ax::PyObject, hh, color::ASCIIString;kwargs...)
    kwargsd = [k=>v for (k, v) in kwargs]
    do_error = pop!(kwargsd, :do_error, true)
    scaling = pop!(kwargsd, :scaling, 1.0)
    empty_overflow = pop!(kwargsd, :empty_overflow, false)
    xs = zeros(2 * nbins(hh)[1]-2);
    ys = zeros(2 * nbins(hh)[1]-2);

#     edgs = deepcopy(edges(hh)[1])
#     edgs[1] = edgs[2] - diff(edgs)[2]
#     edgs[end] = edgs[end-1] + diff(edgs)[end-1]
    conts = deepcopy(contents(hh))

    if empty_overflow
        overflow = conts[end]
        conts[end] = 0
        underflow = conts[1]
        conts[1] = 0
        overflow>0 && println("overflow: $overflow")
        underflow>0 && println("underflow: $underflow")
    end
    if scaling == :normed
        scaling = 1.0 / sum(conts)
    end

    for i=1:2*nbins(hh)[1]-2
        xs[i] = edges(hh)[1][1 + floor(i/2)]
        ys[i] = scaling * conts[1 + floor((i-1)/2)]
        #println(xs[i], " ", ys[i])
    end
    unshift!(xs, xs[1])
    unshift!(ys, 0)

    push!(xs, xs[end])
    push!(ys, 0)

    p = ax[:plot](xs, ys; color=color, kwargsd...)
    do_error && ax[:errorbar](
    midpoints(edges(hh)[1]), scaling * conts, scaling * errors(hh), fmt=nothing, ecolor=color
    )
    return p
end

function barplot{T<:Float64, N, E<:(Array{Float64,1},Array{Float64,1},)}(ax::PyObject, hh::ErrorHistogram{T,N,E}, color::ASCIIString;kwargs...)
    kwargsd = [k=>v for (k, v) in kwargs]
    do_error = pop!(kwargsd, :do_error, true)
    scaling = pop!(kwargsd, :scaling, 1.0)

    ws = deepcopy(hh.weights)
    p = ax[:matshow](
        hh.weights[2:end-1,2:end-1],
        interpolation="none",
        cmap=color,
        extent=[hh.edges[2][2],
        hh.edges[2][end-1],
        hh.edges[1][end-1],
        hh.edges[1][2]],
        aspect="auto"
    )
    #ax[:set_aspect](10)
    return p
end

export barplot
end # module


text:
include("./PersistentVectorBench.jl")
include("./PersistentHashMapBench.jl")


text:
module Plotting

using Reexport
@reexport using Plots
@reexport using StructArrays
@reexport using RecursiveArrayTools
@reexport using LaTeXStrings

using Flight.Modeling

export TimeHistory
export plots, save_plots, thplot, thplot!

########################## plots method ###############################

#our entry plotting entry point cannot be a recipe. a recipe is called within
#the plot() pipeline, which creates a single figure. however, a Vector of System
#outputs will typically need to generate multiple plots from its values. for
#this we define a new method plots(). maybe define different detail levels /
#modes for

plots(args...; kwargs...) = println("Not implemented")

function plots(mdl::Model; mode::Symbol = :basic,
    save_path::Union{String,Nothing} = nothing, kwargs...)
    #generate default path tmp/plots/current_date
    save_path = (save_path === nothing ?
        joinpath("tmp", Dates.format(now(), "yyyy_mm_dd_HHMMSS")) : save_path)
    mkpath(save_path)
    plots(mdl.log.t, mdl.log.saveval; mode, save_path, kwargs...)
end

#for Systems without outputs
plots(::AbstractVector{<:Real}, ::AbstractVector{Nothing}; kwargs...) = nothing

#takes care of SystemGroups
function plots(t, data::AbstractVector{<:NamedTuple}; mode, save_path, kwargs...)

    c = data |> StructArray |> StructArrays.components
    for (c_label, c_data) in zip(keys(c), values(c))
        println("Generating plots for $c_label")
        save_path_c = mkpath(joinpath(save_path, String(c_label)))
        plots(t, c_data; mode, save_path = save_path_c, kwargs...)
    end

end

function save_plots(d::Dict{String,Plots.Plot}; save_path, format = :png)
    for (id, p) in zip(keys(d), values(d))
        savefig(p, joinpath(save_path, id*"."*String(format)))
    end
end

############################ TimeHistory #################################

mutable struct TimeHistory{D}
    t::AbstractVector{<:Real}
    data::D
end

#for all the following TimeHistory subtypes, a single figure is enough, so they
#can be handled by recipes
@recipe function f(th::TimeHistory{<:AbstractVector{<:Real}})

    xguide --> L"$t \: (s)$"
    return th.t, th.data

end

@recipe function f(th::TimeHistory{<:AbstractMatrix{<:Real}}; th_split = :none)

    xguide --> L"$t \ (s)$"

    vlength = size(th.data)[2]
    label --> (vlength <= 3 ?  ["x" "y" "z"][:, 1:vlength] : (1:vlength)')
    if th_split === :h
        layout --> (1, vlength)
        link --> :y #alternative: :none
    elseif th_split === :v
        layout --> (vlength, 1)
    else
        layout --> 1
    end

    return th.t, th.data

end

#convert to TimeHistory{Matrix} and return it to the pipeline for dispatching
@recipe function f(th::TimeHistory{<:AbstractVector{<:AbstractVector{<:Real}}})
    return TimeHistory(th.t, Array(VectorOfArray(th.data))')
end

thplot(t, data; kwargs...) = plot(TimeHistory(t, data); kwargs...)
thplot!(t, data; kwargs...) = plot!(TimeHistory(t, data); kwargs...)

end


text:
using Permutations, StaticArrays, FileIO

println("Building of HyperoctahedralGroup.jl started.")

const xyp = SA[0 -1 0 0; 1 0 0 0; 0 0 1 0; 0 0 0 1]
const xym = SA[0 1 0 0; -1 0 0 0; 0 0 1 0; 0 0 0 1]

const xzp = SA[0 0 -1 0; 0 1 0 0; 1 0 0 0; 0 0 0 1]
const xzm = SA[0 0 1 0; 0 1 0 0; -1 0 0 0; 0 0 0 1]

const xwp = SA[0 0 0 -1; 0 1 0 0; 0 0 1 0; 1 0 0 0]
const xwm = SA[0 0 0 1; 0 1 0 0; 0 0 1 0; -1 0 0 0]

const yzp = SA[1 0 0 0; 0 0 -1 0; 0 1 0 0; 0 0 0 1]
const yzm = SA[1 0 0 0; 0 0 1 0; 0 -1 0 0; 0 0 0 1]

const ywp = SA[1 0 0 0; 0 0 0 -1; 0 0 1 0; 0 1 0 0]
const ywm = SA[1 0 0 0; 0 0 0 1; 0 0 1 0; 0 -1 0 0]

const zwp = SA[1 0 0 0; 0 1 0 0; 0 0 0 -1; 0 0 1 0]
const zwm = SA[1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 -1 0]

function create_elements()
    elements = [SA[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]]
    counter = 0
    while length(elements) < 192 || counter == 16
        if true
            new_elements = [xyp*element for element in elements]
            for (i,e) in enumerate(new_elements)
                if e ∉ elements
                    push!(elements,e)
                end
            end
        end

        if true
            new_elements = [xzp*element for element in elements]
            for (i,e) in enumerate(new_elements)
                if e ∉ elements
                    push!(elements,e)
                end
            end
        end

        if true
            new_elements = [xwp*element for element in elements]
            for (i,e) in enumerate(new_elements)
                if e ∉ elements
                    push!(elements,e)
                end
            end
        end

        if false
            new_elements = [yzp*element for element in elements]
            for (i,e) in enumerate(new_elements)
                if e ∉ elements
                    push!(elements,e)
                end
            end
        end

        if false
            new_elements = [ywp*element for element in elements]
            for (i,e) in enumerate(new_elements)
                if e ∉ elements
                    push!(elements,e)
                end
            end
        end

        if false
            new_elements = [zwp*element for element in elements]
            for (i,e) in enumerate(new_elements)
                if e ∉ elements
                    push!(elements,e)
                end
            end
        end
        counter += 1
    end
    return elements
end

function verifyelements(elements)
    if length(unique(elements)) != 192
        return false
    end
    ones_1x4 = SA[1 1 1 1]
    ones_4x1 = reshape(ones_1x4,(4,1))
    for el in elements
        if size(el) != (4,4)
            return false
        end
        if !issubset(unique(el),[0,1,-1])
            return false
        end
        if map(abs,sum(el,dims=1)) != ones_1x4
            return false
        end
        if map(abs,sum(el,dims=2)) != ones_4x1
            return false
        end
    end
    return true
end

let elements_file = joinpath(@__DIR__,"elements.jld2")
    timestarted = time()
    verified = false
    howstring = "loaded"
    if isfile(elements_file)
        println("Loading group elements.")
        elements = load(elements_file,"elements")
        verified = verifyelements(elements)
    end
    if !verified
        howstring = "created"
        println("Elements could not be loaded. Creating group elements.")
        global elements = create_elements()
        save(elements_file,"elements",elements)
    end
    timeelapsed = time()-timestarted
    println("All elements $howstring after $timeelapsed s.")
end

function write_elements()
    f = open(joinpath(@__DIR__,"elements.jl"),"w")
    write(f,"elements = [")
    for el in elements[1:end-1]
        write(f,"SA")
        print(f,el)
        write(f,",\n")
    end
    write(f,"SA")
    print(f,last(elements))
    write(f,"]\n")
    close(f)
end

write_elements()

function order(e)
    r = e*e
    n = 1
    while r != e
        r = r*e
        n+=1
    end
    return n
end

function group_act_slow(a::Integer,b::Integer)
    if 192 >= a > 0 && 192 >= b > 0
        for i = 1:192
            if elements[i] == elements[a]*elements[b]
                return typeof(a)(i)
            end
        end

        println(elements[a])
        println(elements[b])
        println(elements[a]*elements[b])
        error("Not found? How is that even possible?")
    else
        error("Not defined.")
    end
    return zero(typeof(a))
end

const g = [group_act_slow(a,b) for a in 1:192, b in 1:192]

function ⊗(a::Integer,b::Integer)
    if 192 >= a > 0 && 192 >= b > 0
        return g[a,b]
    else
        error("$a o $b Not defined.")
    end
end

function insertsorted!(set,element;startindex = 1)
    if !issorted(set)
        error("Array is not sorted.")
    end
    i = startindex
    while  i <= length(set) && element > set[i]
        i += 1
    end
    i > length(set) ? push!(set, element) : insert!(set, i, element)
    return i
end

function complete!(set)
    added = -1
    queue = vec([(a,b) for a in set, b in set])
    while added != 0
        added = 0
        newElements = empty(set)
        for i in length(queue):-1:1
            a,b = queue[i]
            if ⊗(a,b) ∉ set
                push!(newElements,⊗(a,b))
                added = added+1
            end
            pop!(queue)
        end
        newElements = sort(unique(newElements))
        for a in newElements
            for b in set
                push!(queue,(a,b))
                push!(queue,(b,a))
            end
            for b in newElements
                push!(queue,(a,b))
            end
        end

        i = 1
        for a in newElements
            i = insertsorted!(set,a, startindex = i)
        end
    end
    return set
end

function completewith!(set,toadd)
    toadd = setdiff(toadd,set)
    isempty(toadd) && return set
    added = -1
    queue = vcat(vec([(a,b) for a in set, b in toadd]),vec([(b,a) for a in set, b in toadd]),vec([(a,b) for a in toadd, b in toadd]))
    for a in toadd
        insertsorted!(set,a)
    end
    while added != 0
        added = 0
        newElements = empty(set)
        for i in length(queue):-1:1
            a,b = queue[i]
            if ⊗(a,b) ∉ set
                push!(newElements,⊗(a,b))
                added = added+1
            end
            pop!(queue)
        end
        newElements = sort(unique(newElements))
        for a in newElements
            for b in set
                push!(queue,(a,b))
                push!(queue,(b,a))
            end
            for b in newElements
                push!(queue,(a,b))
            end
        end

        i = 1
        for a in newElements
            i = insertsorted!(set,a, startindex = i)
        end
    end
    return set
end

function getsubgroups(;maxgenerators = 4)
    subgroupdict = Dict()
    subgroupdict[(generators = 0,)] = [[1]]
    for n = 1:maxgenerators
        timestarted = time()
        println("Generating with $n generators.")
        subgroupbase = subgroupdict[(generators = n-1,)]
        subgroupstoadd_duplicates = vcat([[completewith!(copy(s),[a]) for s in subgroupbase] for a in 1:192]...)
        subgroupstoadd = setdiff(unique(subgroupstoadd_duplicates),vcat([subgroupdict[(generators = i,)] for i = 0:n-1]...))
        subgroupdict[(generators = n,)] = subgroupstoadd
        timeelapsed = time()-timestarted
        println("  Generated after $timeelapsed s.")
    end
    return subgroupdict
end

function addgenerator(subgroups,generator)
    newsubgroups = empty(subgroups)
    for subgroup in subgroups
        newsubgroup = copy(subgroup)
        generator ∈ subgroup && continue
        completewith!(newsubgroup,[generator])
        push!(newsubgroups,newsubgroup)
    end
    return newsubgroups
end

function isgroup(set)
    for a in set
        for b in set
            if ⊗(a,b) ∉ set
                return false
            end
        end
    end
    return true
end

function verifysubgroupdict(subgroupdict)
    #TODO: implement subgroupdict verification
    return true
end

let subgroupdict_file = joinpath(@__DIR__,"subgroupdict.jld2")
    timestarted = time()
    verified = false
    howstring = "loaded"
    if isfile(subgroupdict_file)
        println("Loading all subgroups.")
        subgroupdict = load(subgroupdict_file,"subgroupdict")
        verified = verifysubgroupdict(subgroupdict)
    end
    if !verified
        howstring = "generated"
        println("Generating all subgroups.")
        global subgroupdict = getsubgroups()
        for igens = 1:4
            permute!(subgroupdict[(generators = igens,)],sortperm(length.(subgroupdict[(generators = igens,)])))
        end
        subgroupdict[:allsubgroups] = vcat([subgroupdict[(generators = i,)] for i = 0:4]...)
        save(subgroupdict_file,"subgroupdict",subgroupdict)
    end
    timeelapsed = time()-timestarted
    println("All subgroups $howstring after $timeelapsed s.")
end

function write_subgroups()
    f = open(joinpath(@__DIR__,"subgroups.jl"),"w")
    write(f,"subgroupdict = Dict(")
    for key in keys(subgroupdict)
        print(f,key)
        write(f," => ")
        print(f,subgroupdict[key])
        write(f,",\n")
    end
    skip(f,-2)
    write(f,")\n")
    close(f)
end

write_subgroups()

const vertices = vec([SA[a,b,c,d] for a in [-1,1], b in [-1,1], c in [-1,1], d in [-1,1]])

function element_as_vertex_perm(element)
    permuted_vertices = [element*vertex for vertex in vertices]
    vertex_perm = [findfirst([v_ori==v_new for v_ori in vertices]) for v_new in permuted_vertices]
    return Permutation(vertex_perm)
end

elements_as_vertex_perms = [element_as_vertex_perm(e) for e in elements]

function allequal(array)
    previous = first(array)
    result = true
    for val in array
        result = result && isequal(previous,val)
        previous = val
        !result && continue
    end
    return result
end

function isequalcoloring(coloring1,coloring2)
    if length(unique(coloring1)) != length(unique(coloring2))
        return false
    end
    for c in unique(coloring1)
        v1 = coloring1 .== c
        if !allequal(coloring2[v1])
            return false
        end
    end
    return true
end

function get_subgroup_properties()
    subgroup_properties = [Dict() for s in subgroupdict[:allsubgroups]]
    for (i,sdict) in enumerate(subgroup_properties)
        sdict[:subgroup] = subgroupdict[:allsubgroups][i]
        sdict[:order] = length(sdict[:subgroup])
        j = i
        ngenerators = -1
        while j > 0
            ngenerators += 1
            with_ngenerators = length(subgroupdict[(generators = ngenerators,)])
            with_ngenerators == 0 && break
            j -= with_ngenerators
        end
        sdict[:ngenerators] = ngenerators
    end

    for (is,sdict) in enumerate(subgroup_properties)
        sdict[:coloring] = collect(1:16)
        for (iv,v) in enumerate(vertices)
            group_action_on_v = [elements[e]*v for e in sdict[:subgroup]]
            for connected_vertex in group_action_on_v
                i_con_vert = 1
                while vertices[i_con_vert] != connected_vertex && i_con_vert < 17
                    i_con_vert = i_con_vert+1
                end
                #if connected_vertex has lower color than v, leave it be (should not happen)
                sdict[:coloring][i_con_vert] = min(sdict[:coloring][i_con_vert],sdict[:coloring][iv])
            end
        end
    end

    for (is,sdict) in enumerate(subgroup_properties)
        subgroupcoloring = sdict[:coloring]
        newcolorvalue = 1
        newsubgroupcoloring = copy(subgroupcoloring)
        for colorvalue in sort(unique(subgroupcoloring))
            newsubgroupcoloring[subgroupcoloring.==colorvalue] .= newcolorvalue
            newcolorvalue = newcolorvalue+1
        end
        sdict[:coloring_sorted] = newsubgroupcoloring
    end

    for (is,sdict) in enumerate(subgroup_properties)
        sdict[:coloring_equiv] = zeros(Int,16)
    end

    for (is,sdict) in enumerate(subgroup_properties)
        sgcolor = sdict[:coloring_sorted]
        found = false
        for p in elements_as_vertex_perms
            sgcolorperm = sgcolor[[getindex(p,k) for k in axes(sgcolor,1)]]
            for other_sdict in subgroup_properties
                if isequalcoloring(sgcolorperm,other_sdict[:coloring_equiv])
                    found = true
                    sdict[:coloring_equiv] = other_sdict[:coloring_equiv]
                    continue
                end
            end
            found && continue
        end
        found && continue
        sdict[:coloring_equiv] = sgcolor
    end

    for (is,sdict) in enumerate(subgroup_properties)
        sdict[:coloring_equiv_index] = 0
    end

    coloring_equiv_index = 1
    for (is,sdict) in enumerate(subgroup_properties)
        found = false
        for other_is in 1:is-1
            if sdict[:coloring_equiv] == subgroup_properties[other_is][:coloring_equiv]
                found = true
                sdict[:coloring_equiv_index] = subgroup_properties[other_is][:coloring_equiv_index]
                continue
            end
        end
        found && continue
        sdict[:coloring_equiv_index] = coloring_equiv_index
        coloring_equiv_index += 1
    end

    for col_eq_i = 1:coloring_equiv_index-1
        subgs = [sdict[:coloring_equiv_index] == col_eq_i for sdict in subgroup_properties]
        if length(unique([sdict[:ngenerators] for sdict in subgroup_properties[subgs]])) > 1
            n = length(unique(first(subgroup_properties[subgs])[:coloring_equiv]))
            println("Varying number of generators for coloring_equiv_index = $col_eq_i with $n colors")
        end
    end

    mirrorx = [-1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
    elements_as_vertex_perms_mirror = vcat([element_as_vertex_perm(e) for e in elements],[element_as_vertex_perm(e*mirrorx) for e in elements])

    for (is,sdict) in enumerate(subgroup_properties)
        sdict[:coloring_equiv_mirror] = zeros(Int,16)
    end

    for (is,sdict) in enumerate(subgroup_properties)
        println("Processing $is-th subgroup.")
        sgcolor = sdict[:coloring_sorted]
        found = false
        for p in elements_as_vertex_perms_mirror
            sgcolorperm = sgcolor[[getindex(p,k) for k in axes(sgcolor,1)]]
            for other_sdict in subgroup_properties
                if isequalcoloring(sgcolorperm,other_sdict[:coloring_equiv_mirror])
                    found = true
                    sdict[:coloring_equiv_mirror] = other_sdict[:coloring_equiv_mirror]
                    continue
                end
            end
            found && continue
        end
        found && continue
        sdict[:coloring_equiv_mirror] = sgcolor
    end

    for (is,sdict) in enumerate(subgroup_properties)
        sdict[:coloring_equiv_mirror_index] = 0
    end

    coloring_equiv_mirror_index = 1
    for (is,sdict) in enumerate(subgroup_properties)
        found = false
        for other_is in 1:is-1
            if sdict[:coloring_equiv_mirror] == subgroup_properties[other_is][:coloring_equiv_mirror]
                found = true
                sdict[:coloring_equiv_mirror_index] = subgroup_properties[other_is][:coloring_equiv_mirror_index]
                continue
            end
        end
        found && continue
        sdict[:coloring_equiv_mirror_index] = coloring_equiv_mirror_index
        coloring_equiv_mirror_index += 1
    end

    for col_eq_i = 1:coloring_equiv_mirror_index-1
        subgs = [sdict[:coloring_equiv_mirror_index] == col_eq_i for sdict in subgroup_properties]
        if length(unique([sdict[:ngenerators] for sdict in subgroup_properties[subgs]])) > 1
            n = length(unique(first(subgroup_properties[subgs])[:coloring_equiv_mirror]))
            println("Varying number of generators for coloring_equiv_index = $col_eq_i with $n colors")
        end
    end
    return subgroup_properties
end

let subgroup_properties_file = joinpath(@__DIR__,"subgroup_properties.jld2")
    println("Calculating subgroup properties.")
    global subgroup_properties = get_subgroup_properties()
    save(subgroup_properties_file,"subgroup_properties", subgroup_properties)
end

function write_groupmatrix()
    f = open(joinpath(@__DIR__,"groupmatrix.jl"),"w")
    groupmatrixname = :g
    write(f,"const $(string(groupmatrixname)) = [")
    for row in axes(eval(groupmatrixname),1)
        for col in axes(eval(groupmatrixname),2)
            write(f,"$(g[row,col]) ")
        end
        write(f,";\n")
    end
    skip(f,-2)
    write(f,"]\n")
    close(f)
end

write_groupmatrix()


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

using Markdown
using InteractiveUtils

# ╔═╡ 7608d6b8-0dee-4890-9d82-82ceb43e2399
using LinearAlgebra

# ╔═╡ f647cff6-88e1-11eb-3085-d31e4f22c78c
using Yao

# ╔═╡ 45d06384-88d7-11eb-1c35-03e184c98c8b
include("./essentials.jl")

# ╔═╡ 33ab855e-0e9d-49ba-bc05-cf758fbf8674
md"""
A measurement is a central concept in quantum mechanics. An easy way to think about it as a sample from a probability distribution: it is a random variable with a number of outcomes, each outcome is produced with a certain probability.

Measurement connect the quantum world to our classical one: we cannot directly observe the quantum state in nature, we can only gather statistics about it with measurements. It sounds like a harsh boundary between a quantum and a classical system that can only be bridged by measurement. The reality is more subtle: unless a quantum system is perfectly isolated, it interacts with its surrounding environment. This leads to introduction of mixed states, which in one limit recover classical probabilities.

# More on the bra-ket notation

Before we take a deep dive into what measurements are, we need to introduce one more notation to complement the ket: it called a bra and it is denoted by $\langle\psi|$ for some quantum state $|\psi\rangle$. Together they form the bra-ket or Dirac notation. A bra is the conjugate transpose of a ket, and the other way around. This also means that a bra is a row vector. For instance, this is the bra for $|0\rangle$
"""

# ╔═╡ 59f6626a-88d3-11eb-2a06-2982fa6594fc
begin
	zero_ket = hcat([1, 0])
	Text("|0〉 ket: "), zero_ket
end

# ╔═╡ 57c49654-88d0-11eb-3856-1d2468c2a406
Text("〈0| bra: "), (conj(zero_ket'ᵀ))

# ╔═╡ ab6ec464-989d-45b0-b70e-a55e89a705ae
md"""
This makes it very easy to write dot products: if we write a bra followed by a ket, that is exactly what the dot product is. This is so common that we often drop one of the vertical bars, and just write $\langle 0|0\rangle$, for instance. Since quantum states are normalized, the inner product of any quantum state with itself is always one:
"""

# ╔═╡ 695cbee6-6980-462c-a537-9ee803109066
conj(zero_ket'ᵀ) * zero_ket

# ╔═╡ abe6fb03-d5fd-41a3-9349-8383b4f4cdd3
md"""
Similarly, orthogonal vectors always give 0. E.g. $\langle 0|1\rangle$:
"""

# ╔═╡ 33404c49-bdba-4fb7-a4ec-1109d139c656
begin
	one_ket = hcat([0, 1])
	conj(zero_ket'ᵀ) * one_ket
end

# ╔═╡ ff813c16-f62f-4083-849f-033b8b20baf9
md"""
What about a ket and a bra? That is going to be a matrix: essentially the outer product of the two vectors. Here's $|0\rangle\langle 0|$:
"""

# ╔═╡ 6f8bdaed-0dc3-4c50-ac9e-a357e5633350
zero_ket * conj(zero_ket'ᵀ)

# ╔═╡ b916a917-a1ed-472a-b8de-15df3e73a964
md"""
This should look familiar: it is a projection to the first element of the canonical basis. It is true in general that $|\psi\rangle\langle\psi|$ is going to be a projector to $|\psi\rangle$. It is very intuitive: take some other quantum state $|\phi\rangle$ and apply the matrix $|\psi\rangle\langle\psi|$ on it: $|\psi\rangle\langle\psi|\phi\rangle$. Now the right-most two terms are a bra and a ket, so it is a dot product: the overlap between $|\phi\rangle$ and $|\psi\rangle$. Since this is a scalar, it just scales the left-most term, which is the ket $|\psi\rangle$, so in effect, we projected $|\phi \rangle$ on this vector.
"""

# ╔═╡ e7916dcd-991c-4b6f-8f3c-43f52bf51ad5
md"""
# Measurements

A measurement in quantum mechanics is an operator-valued random variable. The theory of measurements is rich and countless questions about them are still waiting to be answered. Most quantum computers that we have today, however, only implement one very specific measurement, which makes our discussion a lot simpler. This measurement is in the canonical basis. In other words, the measurement contains two projections, $|0\rangle\langle 0|$ and $|1\rangle\langle 1|$, and this measurement can be applied to any of the qubits of the quantum computer.

We already saw how applying a projection on a vector works. If we want to make a scalar value of that, we need to add a bra to the left. For instance, for some state $|\psi\rangle$, we get a scalar for $\langle\psi|0\rangle\langle 0|\psi\rangle$. This is called the expectation value of the operator $|0\rangle\langle 0|$. To put this in context, let us apply the projection $|0\rangle\langle 0|$ on the superposition $\frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)$, which is the column vector $\frac{1}{\sqrt{2}}\begin{bmatrix} 1\\ 0\end{bmatrix}$.
"""

# ╔═╡ 4ff913cd-0717-40fd-a1f4-d099678ff933
let
	ψ = hcat([1,0])/√2
	Π_0 = zero_ket * conj(zero_ket'ᵀ)
	conj(ψ'ᵀ) * Π_0 * ψ
end

# ╔═╡ f8bb6198-0c97-4b62-aa77-af71dcfbe3e6
md"""
That is exactly one half, the square of the absolute value of the probability amplitude corresponding to $|0\rangle$ in the superposition! This is the mathematical formalism of what we had said earlier: given a state $|\psi\rangle = a_0|0\rangle + a_1|1\rangle$, we get an output $i$ with probability $|a_i|^2$. This is known as the *Born rule*. Now we have a recipe to extract probabilities with projections. This is exactly what is implemented in the quantum simulator. The measurement in the simulator is what we described here. Let's create an equal superposition with the Hadamard gate (see a later notebook for quantum circuits), apply the measurement, and observe the statistics:
"""

# ╔═╡ c5bbbc1d-f298-4f68-82d4-c8ddec2c9b7b
let
	q = zero_state(1)
	essentials.plot_hist(measure(q |> put(1,1=>H), nshots=100))
end

# ╔═╡ f67899cf-3bf5-4b9b-bdc0-e43576eb28bc
md"""
You see that the outcome is random, with roughly half of the outcomes being 0.

There is something additional happening. The measurement has a random outcome, but once it is performed, the quantum state is in the corresponding basis vector. That is, the superposition is destroyed. This is referred to as the collapse of the wavefunction. It is the subject of many ongoing debates and research results how and why it happens, but what matters to us is that we can easily calculate the quantum state after the measurement. Just projecting it to the basis vector is insufficient, since that would not be normalized, so we have to renormalize it. Mathematically it is expressed by the somewhat convoluted expression $\frac{|i\rangle\langle i|\psi\rangle}{\sqrt{\langle\psi|i\rangle\langle i|\psi\rangle}}$ if we observe the output $i$. For instance, if we observe zero after measuring the superposition $\frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)$, the state after the measurement will be
"""

# ╔═╡ 5e599569-48e9-477d-9020-ba546ae158ce
let
	ψ = hcat([√2/2, √2/2])
	Π_0 = zero_ket * conj(zero_ket'ᵀ)
	probability_0 = conj(ψ'ᵀ) * Π_0 * ψ
	Π_0 * ψ/√probability_0
end

# ╔═╡ 9e433703-d4b9-439d-af4c-50f06fa4e241
md"""
which is just a very long way of saying we get $|0\rangle$.

You can easily see this by putting two measurements in a sequence on the same qubit. The second one will always give the same outcome as the first. The first one is random, but the second one will be determined, since there will be no superposition in the computational basis after the first measurement. Let's simulate this by writing out the results of the two measurements into two different classical registers:
"""

# ╔═╡ b8480514-54a0-44e1-906f-ec891395dc36
let
	q = zero_state(1)
	q = q |> put(1, 1=>H)
	a = measure!(q)
	b = measure(q, nshots=100)
	a, sum(b .== a)
end

# ╔═╡ e63a1c5f-e06b-4880-a735-a09b741c5688
md"""
There is no output like 01 or 10.
"""

# ╔═╡ 31832051-7c32-4f49-8751-2d79fc99076a
md"""
# Measuring multiqubit systems

Most quantum computers implement local measurements, which means that each qubit is measured separately. So if we have a two qubit system where the first qubit is in the equal superposition and the second one is in $|0\rangle$, that is, we have the state $\frac{1}{\sqrt{2}}(|00\rangle + |01\rangle)$, we will observe 0 and 0 as outcomes of the measurements on the two qubits, or 0 and 1.
"""

# ╔═╡ 917037ea-4132-4816-8d0c-0004781d3777
let
	q = zero_state(2)
	circuit = put(2, 1=>H)
	essentials.plot_hist(measure(q |> circuit, nshots=1024))
end

# ╔═╡ daf68859-d5fc-4df0-9587-950b4c543d54
md"""
What happens if we make measurements on an entangled state? Let's look at the statistics again on the $|\phi^+\rangle$ state:
"""

# ╔═╡ 7129a922-2061-4ecc-a247-3db7903c5f84
let
	q = zero_state(2)
	circuit = chain(2, put(1=>H), control(1, 2=>X))
	essentials.plot_hist(measure(q |> circuit, nshots=1024))
end

# ╔═╡ d7de4321-daee-404b-b408-975928463eea
md"""
We only observe 00 and 11. Since the state is $\frac{1}{\sqrt{2}}(|00\rangle+|11\rangle)$, this should not come as a shock. Yet, there is something remarkable going on here. At the end of the last section, we saw the same statistics, but from measurements on the same qubit. Now we have two, spatially separate qubits exhibiting the same behaviour: this is a very strong form of correlations. This means that if we measure just one qubit, and get, say, 0 as the outcome, we *know* with certainty that if we measured the other qubit, we would also get 0, even though the second measurement is also a random variable.

To appreciate this better, imagine that your are tossing two unbiased coins. If you observe heads on one, there is absolutely nothing that you can say about what the other one might be other than a wild guess that holds with probability 0.5. If you play foul and you biased the coins, you might improve your guessing accuracy. Yet you can never say with certainty what the other coin will be based on the outcome you observed on one coin, except for the trivial case when the other coin deterministically gives the same face always.

Remarkable as it is, there is no activation or instantaneous (faster than the speed of light) signalling happening between the qubits, though. Your measurement was local to the qubit and so is your information. If there is somebody else doing the measurement on the other qubit, you would have to inform the person through classical communication channels that you happen to know what the outcome will be. So while we certainly cannot violate the theory of relativity with entanglement, this strong form of correlation is still central to many quantum algorithms.
"""

# ╔═╡ a9c7d4bc-a74d-4a64-8c75-19f1f6fec9f1
md"""
# Mixed states

If all was clear until now, this is where it gets messy. A ket and a bra is a projection, as we explained above. More than that, it is also a density matrix. A density matrix is another way of writing a quantum state, instead of kets. So, for instance we could write $\rho = |\psi\rangle\langle\psi|$, where $\rho$ is the density matrix for $|\psi\rangle$. The Born rule still applies, but now we have to take the trace of the result: for instance, $\mathrm{Tr}[|0\rangle\langle 0|\rho]$ would be the probability of seeing 0. See it in action:
"""

# ╔═╡ ddfb368c-d940-4146-9020-7e22d9f66ef1
begin
	ψ = hcat([1, 1])/√2
	ρ = ψ * conj(ψ'ᵀ)
	Π_0 = zero_ket * conj(zero_ket'ᵀ)
	tr(Π_0 * ρ)
end

# ╔═╡ c50ecff4-5234-49e7-8fa0-b29049443545
md"""
We get one half again. The renormalization after a measurement happens in a similar way: $\frac{|0\rangle\langle 0|\rho|0\rangle\langle 0|}{\mathrm{Tr}[|0\rangle\langle 0|\rho]}$.
"""

# ╔═╡ 26300e25-f24f-4392-b963-97ee1b04666c
begin
	probability_0 = tr(Π_0 * ρ)
	Π_0 * ρ * Π_0/probability_0
end

# ╔═╡ 2b5d8fe2-9222-4193-98c8-390364473a5a
md"""
So why do we need this at all? Every state we have mentioned so far is called a *pure state*: these are kets or a density matrix created as a ket and a bra. There are other states called *mixed states*: these are *classical* probability distributions over pure states. Formally, a mixed state is written as $\sum_i p_i |\psi_i\rangle\langle\psi_i|$, where $\sum_i p_i=1$, $p_i\geq 0$. This reflects our classical ignorance over the underlying quantum states. Compare the density matrix of the equal superposition $\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$ and the mixed state $0.5(|0\rangle\langle 0|+|1\rangle\langle 1|)$:
"""

# ╔═╡ dd88aef3-5954-47e7-8883-7e4c03b9df10
let
	zero_ket = hcat([1, 0])
	one_ket = hcat([0, 1])
	ψ = (zero_ket + one_ket)/√2
	global t2 = "Density matrix of the equally mixed state of |0〉〈0| and |1〉〈1| :", (zero_ket * conj(zero_ket'ᵀ) + one_ket * conj(one_ket'ᵀ))/2
	"Density matrix of the equal superposition: ", (ψ * conj(ψ'ᵀ))
end

# ╔═╡ 68da7bd4-88dd-11eb-0914-4b8279395d1b
t2

# ╔═╡ b365c769-07aa-4df9-b933-3ce55cbb6ab4
md"""
The off-diagonal elements are gone in the second case. The off-diagonal elements are also called coherences: their presence indicates that the state is quantum. The smaller these values are, the closer the quantum state is to a classical probability distribution.

The second density matrix above has only diagonal elements and they are equal: this is the equivalent way of writing a uniform distribution. We know that the uniform distribution has maximum entropy, and for this reason, a density matrix with this structure is called a maximally mixed state. In other words, we are perfectly ignorant of which elements of the canonical basis constitute the state.

We would like a quantum state to be perfectly isolated from the environment, but in reality, the quantum computers we have today and for the next couple of years cannot achieve a high degree of isolation. So coherences are slowly lost to the environment -- this is a process called decoherence. The speed at which this happens determines the length of the quantum algorithms we can run on the quantum computer: if it happens fast, we have time to apply a handful gates or do any other form calculation, and then we quickly have to pull out (measure) the results.
"""

# ╔═╡ Cell order:
# ╟─33ab855e-0e9d-49ba-bc05-cf758fbf8674
# ╠═7608d6b8-0dee-4890-9d82-82ceb43e2399
# ╠═59f6626a-88d3-11eb-2a06-2982fa6594fc
# ╠═57c49654-88d0-11eb-3856-1d2468c2a406
# ╟─ab6ec464-989d-45b0-b70e-a55e89a705ae
# ╠═695cbee6-6980-462c-a537-9ee803109066
# ╟─abe6fb03-d5fd-41a3-9349-8383b4f4cdd3
# ╠═33404c49-bdba-4fb7-a4ec-1109d139c656
# ╟─ff813c16-f62f-4083-849f-033b8b20baf9
# ╠═6f8bdaed-0dc3-4c50-ac9e-a357e5633350
# ╟─b916a917-a1ed-472a-b8de-15df3e73a964
# ╟─e7916dcd-991c-4b6f-8f3c-43f52bf51ad5
# ╠═4ff913cd-0717-40fd-a1f4-d099678ff933
# ╟─f8bb6198-0c97-4b62-aa77-af71dcfbe3e6
# ╠═45d06384-88d7-11eb-1c35-03e184c98c8b
# ╠═f647cff6-88e1-11eb-3085-d31e4f22c78c
# ╠═c5bbbc1d-f298-4f68-82d4-c8ddec2c9b7b
# ╟─f67899cf-3bf5-4b9b-bdc0-e43576eb28bc
# ╠═5e599569-48e9-477d-9020-ba546ae158ce
# ╟─9e433703-d4b9-439d-af4c-50f06fa4e241
# ╠═b8480514-54a0-44e1-906f-ec891395dc36
# ╟─e63a1c5f-e06b-4880-a735-a09b741c5688
# ╟─31832051-7c32-4f49-8751-2d79fc99076a
# ╠═917037ea-4132-4816-8d0c-0004781d3777
# ╟─daf68859-d5fc-4df0-9587-950b4c543d54
# ╠═7129a922-2061-4ecc-a247-3db7903c5f84
# ╟─d7de4321-daee-404b-b408-975928463eea
# ╟─a9c7d4bc-a74d-4a64-8c75-19f1f6fec9f1
# ╠═ddfb368c-d940-4146-9020-7e22d9f66ef1
# ╟─c50ecff4-5234-49e7-8fa0-b29049443545
# ╠═26300e25-f24f-4392-b963-97ee1b04666c
# ╟─2b5d8fe2-9222-4193-98c8-390364473a5a
# ╠═dd88aef3-5954-47e7-8883-7e4c03b9df10
# ╠═68da7bd4-88dd-11eb-0914-4b8279395d1b
# ╟─b365c769-07aa-4df9-b933-3ce55cbb6ab4

max_stars_repo_name: VarLad/qml-mooc

max_issues_repo_name: VarLad/qml-mooc

max_forks_repo_name: VarLad/qml-mooc


text:
module CachedCalls

import Base.show

export AbstractCachedCall,
        CachedCall,
        description,
        func,
        func_args,
        result,
        call!,
        uncall!,
        is_executed

abstract type AbstractCachedCall{F<:Function, A<:Tuple, R} end

description(cc::T) where {T<:AbstractCachedCall} = cc.description
func(cc::T) where {T<:AbstractCachedCall}  = cc.f
func_args(cc::T) where {T<:AbstractCachedCall} = cc.args
result(cc::T) where {T<:AbstractCachedCall} =
    is_executed(cc) ? cc.result : error("Call has not been executed.")

is_executed(cc::T) where {T<:AbstractCachedCall} = !(cc.result === nothing)

#uncall! and call! are also applied recursively to any args that are AbstractCachedCalls
function call!(cc::AbstractCachedCall)
    if !is_executed(cc)
        cc.result = func(cc)(call!(func_args(cc)...)...)
    end
    return result(cc)
end

call!() = ()
call!(a::CC, args...) where CC <: AbstractCachedCall = (call!(a), call!(args...)...)
call!(a, args...) = (a, call!(args...)...)

function uncall!(cc::AbstractCachedCall{F,A,R}) where {F, A, R}
    if is_executed(cc)
        cc.result = nothing
        uncall!(func_args(cc)...)
    end
    return cc
end

uncall!() = ()
uncall!(a::AbstractCachedCall{F,A,R}, args...) where {F, A, R} =
    begin uncall!(a); uncall!(args...) end
uncall!(a, args...) = uncall!(args...)

function Base.show(s::IO, cc::AbstractCachedCall{F, A, R}) where {F, A, R}
    stat = "Unexecuted"
    if is_executed(cc)
        stat = "Executed"
    end
    write(s, "$stat cached call of function $(func(cc)) with arguments:")
    for a in func_args(cc)
        write(s, "\n$a")
    end
    if is_executed(cc)
        write(s, "\nResult: $(result(cc))")
    end
end

mutable struct CachedCall{F, A, R} <: AbstractCachedCall{F, A, R}
    description::String
    f::F
    args::A
    result::Union{R, Nothing}
end

CachedCall(description::String, f::F, args::A, result_type = Any) where {F<:Function, A<:Tuple} =
    CachedCall{F, A, result_type}(description, f, args, nothing)
CachedCall(f::F, args::A, result_type = Any) where {F<:Function, A<:Tuple} = CachedCall("", f, args, result_type)

end # module


text:
module GD
export GDContext, GradientDescent, fit!

using ..Classification: OptimizerContext, Optimizer, σ, logcrossentropy, stopearly!, sigmoid!, logcrossentropy!, postprocess!
import ..Classification: fit!
using LinearAlgebra: mul!
using Statistics: mean!

struct GDContext <: OptimizerContext
    learningrate::Float64
    tolerance::Float64
    maxiterations::Int
end

GDContext(;learningrate=0.1, tolerance=eps(Float64),
          maxiterations=1000) = GDContext(learningrate, tolerance, maxiterations)

mutable struct GradientDescent <: Optimizer
    context::GDContext
    converged::Bool
    earlystopping::Bool
    hasvalidationset::Bool
    iterations::Int
    loss::Vector{Float32}
    # Lower memory usage is better than high precision
    validationset::Tuple{Matrix{Float32}, Vector{Int8}}
    validationloss::Vector{Float32}
end

function GradientDescent(context=GDContext())
    GradientDescent(context, false, false, false, 0, Float32[],
                   (zeros(Float32, 1, 1), Int8[]), Float32[])
end

function fit!(optim::GradientDescent, β::Vector{T},
              X::Matrix{<:Real}, y::Vector{<:Real}) where {T<:Real}
    context = optim.context
    previousloss = -Inf
    optim.loss = zeros(Float32, context.maxiterations)
    # Preallocate everything
    ŷ         = zeros(Float32, length(y))
    residuals = zeros(Float32, length(y))
    ∇array    = zeros(T, (1, length(β)))
    ∇         = @view ∇array[1, :]
    tmp1      = zeros(T, length(y))
    tmp       = zeros(T, size(X))
    loss      = [0.0]
    tmp2      = zeros(Float32, length(y))

    for i in 1:context.maxiterations
        optim.iterations += 1
        # ŷ = σ(X*β)
        mul!(tmp1, X, β)
        sigmoid!(ŷ, tmp1)
        logcrossentropy!(loss, tmp2, y, ŷ)
        optim.loss[i] = loss[1]

        if abs(previousloss - loss[1]) < context.tolerance
            optim.converged = true
            break
        else
            previousloss = loss[1]
        end

        # Gradient descent
        residuals .= ŷ .- y
        tmp .= X.*residuals
        mean!(∇array, tmp)

        β .-= context.learningrate .* ∇

        if stopearly!(optim, β, i)
            break
        end
    end
    postprocess!(optim)
end

end


text:
module LsqFit

    export curve_fit,
           estimate_errors,
           estimate_covar

    using Optim
    using Calculus
    using Distributions

    include("curve_fit.jl")

end


text:
using TextAnalysis

model = TextAnalysis.SentimentAnalyzer()

model(StringDocument("hello world"))

model(StringDocument("an incredibly boring film"))

model(StringDocument("a highly enjoyable ride"))

model(StringDocument("a great film"))

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
for func in [:clblasSrotm, :clblasDrotm]
    @eval @blasfun $func(N::Csize_t,
                              X::cl.CL_mem, offx::Csize_t, incx::Cint,
                              Y::cl.CL_mem, offy::Csize_t, incy::Cint,
                              SPARAM::cl.CL_mem, offSparam::Csize_t,
                              n_queues::cl.CL_uint,
                              queues::Ptr{cl.CL_command_queue},
                              n_events_in_wait_list::cl.CL_uint,
                              event_wait_list::Ptr{cl.CL_event},
                              events::Ptr{cl.CL_event})

    @eval @blasfun2 $func(N::Csize_t,
                               X::cl.CL_mem, offx::Csize_t, incx::Cint,
                               Y::cl.CL_mem, offy::Csize_t, incy::Cint,
                               SPARAM::cl.CL_mem, offSparam::Csize_t)
end


text:
###
# Const-constrained QR optimizer for sensor placement.
###

mutable struct CostQRPivot
    Ψ::AbstractArray
    pivots::AbstractArray
    sensor_costs::AbstractArray #Costs/weights associated with each sensor. Postsive is bad, negative good.
end

function CostQRPivot(qrpivot::QRPivot,sensor_costs::Vector{Float64})
    return CostQRPivot(qrpivot.Ψ,qrpivot.pivots,sensor_costs)
end

function CostQRPivot(Ψ::AbstractArray,pivots::Vector,sensor_costs::Vector{Float64})
    return CostQRPivot(Ψ,pivots,sensor_costs)
end

function fit(cost_qr_pivot::CostQRPivot)
    n,m = size(cost_qr_pivot.Ψ)

    if cost_qr_pivot.sensor_costs !== n
        throw(DomainError(cost_qr_pivot.sensor_costs,"The specified sensor costs are inconsistent with the basis dimensions"))
    end
    #Initialize helper variables
    R,p,k = make_helper_variables(cost_qr_pivot.Ψ)

    for j in 1:k
        u,i_piv = qr_reflector(R[j:n,j:m],cost_qr_pivot.sensor_costs)
        #Track column pivots
        i_piv += j
        p[[j,i_piv]] = p[[i_piv,j]]
        #Switch column
        R[:,[j,i_piv]] = R[:,[i_piv,j]]
        #Apply reflector
        R[j:n,j:m] -= outer(u,dot(u,R[j:n,j:m]))
        R[j+1:n,j] = 0
    end
    cost_qr_pivot.pivots = p
end

function make_helper_variables(Ψ)
    """Constructs the reflector helper variables"""
    R = conj(Ψ)
    p = [i for i in 1:n]
    k = min(m,n)
    return R,p,k
end

function qr_reflector(r,costs)
    dlens = sqrt(sum(abs(eachrow(r))^2)) #Norm of each column
    i_piv = argmax(dlens-costs) # choose pivot
    dlen = dlens[i_piv]

    if dlen > 0
        u = r[:,i_piv]/dlen
        u[1] += sign(u[1])+(u[1]==0)
        u/=sqrt(abs(u[1]))
    else
        u = r[:,i_piv]
        u[1] = sqrt(2)
    end
    return u,i_piv
end