text:
include("wrapper/preprocWrapper.jl")
include("../src/preproc.jl")

# using Plots
# using DelimitedFiles
using DataFrames
using CSV

function test1()
    A1 = [1.0 2 1 0;2 3 1 1;3 4 1 2;4 5 1 3;5 6 1 4]
    B1 = [6.0 7 1 5;7 1 -6 13;-4 8 9 -2]

    A2 = [1.0 2 1 0;2 3 1 1;3 4 1 2;4 5 1 3;5 6 1 4]
    B2 = [6.0 7 1 5;7 1 -6 13;-4 8 9 -2]

    A = [1.0 2 1 0;2 3 1 1;3 4 1 2;4 5 1 3;5 6 1 4]
    B = [6.0 7 1 5;7 1 -6 13;-4 8 9 -2]
    m, n = size(A1)
    p = size(B1)[1]

    U1, V1, Q1, k1, l1, A1, B1 = dggsvp3!(A1, B1)
    println("k: ", k1)
    println("l: ", l1)
    U2, V2, Q2, k2, l2, A2, B2 = preproc!(A2, B2)

    e = eps(Float64)
    println("Stability test on dggsvp3")
    println("res_a = ", opnorm(U1*A1*Q1' - A, 1)/(max(m,n)*opnorm(A, 1)*e))
    println("res_b = ", opnorm(V1*B1*Q1' - B, 1)/(max(m,n)*opnorm(B, 1)*e))
    println("orthog_u = ", opnorm(U1'*U1-I, 1)/(m*e))
    println("orthog_v = ", opnorm(V1'*V1-I, 1)/(p*e))
    println("orthog_q = ", opnorm(Q1'*Q1-I, 1)/(n*e))
    println("--------------------------------------")
    println("Stability test on preproc")
    println("res_a = ", opnorm(U2*A2*Q2' - A, 1)/(max(m,n)*opnorm(A, 1)*e))
    println("res_b = ", opnorm(V2*B2*Q2' - B, 1)/(max(m,n)*opnorm(B, 1)*e))
    println("orthog_u = ", opnorm(U2'*U2-I, 1)/(m*e))
    println("orthog_v = ", opnorm(V2'*V2-I, 1)/(p*e))
    println("orthog_q = ", opnorm(Q2'*Q2-I, 1)/(n*e))

end

# Example 3 in the special case note
function test2()
    A = [0.  1  2  2; 1  2  1  2; 1  1  1  2; 1  2  2  3];
    B = [0.  0  1  1; 1  2  2  3; 1  1  1  2; 1  1  1  2];
    n = size(A)[1]
    U, V, Q, k, l, A1, B1 = preproc(A, B)
    println("Stability test on preproc")
    e = eps(Float64)
    println("res_a = ", opnorm(U*A1*Q' - A, 1)/(n*opnorm(A, 1)*e))
    println("res_b = ", opnorm(V*B1*Q' - B, 1)/(n*opnorm(B, 1)*e))
    println("orthog_u = ", opnorm(U'*U-I, 1)/(n*e))
    println("orthog_v = ", opnorm(V'*V-I, 1)/(n*e))
    println("orthog_q = ", opnorm(Q'*Q-I, 1)/(n*e))
    return U, V, Q, k, l, A1, B1
end

# Example 4 in the special case note
function test3()
    A =  [4.  2  2  0  0; 0  0  2  0  0; 4  4  4  0  2; 2  2  2  0  0; 2  4  4  0  2];
    B =  [4.  4  2  0  2; 2  2  0  0  2; 2  2  0  0  2; 4  2  2  0  0; 2  0  0  0  0];
    n = size(A)[1]
    U, V, Q, k, l, A1, B1 = preproc(A, B)
    println("Stability test on preproc")
    e = eps(Float64)
    println("res_a = ", opnorm(U*A1*Q' - A, 1)/(n*opnorm(A, 1)*e))
    println("res_b = ", opnorm(V*B1*Q' - B, 1)/(n*opnorm(B, 1)*e))
    println("orthog_u = ", opnorm(U'*U-I, 1)/(n*e))
    println("orthog_v = ", opnorm(V'*V-I, 1)/(n*e))
    println("orthog_q = ", opnorm(Q'*Q-I, 1)/(n*e))
    return U, V, Q, k, l, A1, B1
end

function test(m, p, n)
    A = randn(m, n)
    B = randn(p, n)
    U, V, Q, k, l, A1, B1 = preproc(A, B)
    e = eps(Float64)
    res_a = opnorm(U*A1*Q' - A, 1)/(n*opnorm(A, 1)*e)
    res_b = opnorm(V*B1*Q' - B, 1)/(n*opnorm(B, 1)*e)
    orthog_u = opnorm(U'*U-I, 1)/(n*e)
    orthog_v = opnorm(V'*V-I, 1)/(n*e)
    orthog_q = opnorm(Q'*Q-I, 1)/(n*e)
    return res_a, res_b, orthog_u, orthog_v, orthog_q
end

function test()
    df = DataFrame(m = Int64[], p = Int64[], n = Int64[], res_a = Float64[],
    res_b = Float64[], orthog_u = Float64[],
    orthog_v = Float64[], orthog_q = Float64[])

    # case 1: m >= n & p >= n
    println("=======================Case 1============================")
    for i = 1:20
        a, b, c, d, e = test(60, 50, 40)
        push!(df, [60 50 40 a b c d e])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e = test(300, 250, 200)
        push!(df, [300 250 200 a b c d e])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e = test(900, 750, 600)
        push!(df, [900 750 600 a b c d e])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e = test(1500, 1250, 1000)
        push!(df, [1500 1250 1000 a b c d e])
        println("---------------------------------------------------------------")
    end

    # case 2: m >= n > p
    println("=======================Case 2================================")
    for i = 1:20
        a, b, c, d, e = test(60, 40, 50)
        push!(df, [60 40 50 a b c d e])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e = test(300, 200, 250)
        push!(df, [300 200 250 a b c d e])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e = test(900, 600, 750)
        push!(df, [900 600 750 a b c d e])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e = test(1500, 1000, 1250)
        push!(df, [1500 1000 1250 a b c d e])
        println("---------------------------------------------------------------")
    end

    # case 3: p >= n > m
    println("=======================Case 3==============================")
    for i = 1:20
        a, b, c, d, e = test(40, 60, 50)
        push!(df, [40 60 50 a b c d e])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e = test(200, 300, 250)
        push!(df, [200 300 250 a b c d e])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e = test(600, 900, 750)
        push!(df, [600 900 750 a b c d e])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e = test(1000, 1500, 1250)
        push!(df, [1000 1500 1250 a b c d e])
        println("---------------------------------------------------------------")
    end

    # case 4: n > m & n > p
    println("=====================Case 4============================")
    for i = 1:20
        a, b, c, d, e = test(20, 30, 60)
        push!(df, [20 30 60 a b c d e])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e = test(200, 300, 600)
        push!(df, [200 300 600 a b c d e])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e = test(400, 600, 1200)
        push!(df, [400 600 1200 a b c d e])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e = test(1000, 1500, 3000)
        push!(df, [1000 1500 3000 a b c d e])
        println("---------------------------------------------------------------")
    end

    CSV.write("preproc_rand_matrix_test.csv", df, append=true)
    return df
end


text:
__precompile__()

module MendelIHT

	import Distances: euclidean, chebyshev, sqeuclidean
	import Base.show
	import GLM: glmvar, mueta, fit, linkinv, Link, GeneralizedLinearModel, devresid, checky, canonicallink
	import SpecialFunctions: digamma, trigamma

	using GLM
	using SnpArrays
	using DataFrames
	using Random
	using LinearAlgebra
	using Distributions
	using SparseArrays
	using Distributed

	export loglikelihood, deviance, score!, L0_reg, iht_run_many_models, update_r!
	export iht_path, simulate_random_snparray, make_bim_fam_files, project_k!
	export IHTVariables, use_A2_as_minor_allele, make_snparray, standardize!
	export std_reciprocal, project_group_sparse!, save_prev!, maf_weights
	export simulate_random_response, adhoc_add_correlation, cv_iht
	export simulate_correlated_snparray, naive_impute, cv_iht_distribute_fold, mle_for_θ

	# IHT will only work on single/double precision floats!
	const Float = Union{Float64,Float32}

	include("data_structures.jl")
	include("utilities.jl")
	include("simulate_utilities.jl")
	include("iht.jl")
	include("cross_validation.jl")

end # end module

ext: jl


text:
export double

"""
    double(x)
    double(x::String)

Doubles things

# Examples
```jldoctest
using JuliaProject

double(2)

# output

4
```
"""
double(x) = 2x
double(x::String) = x * x


text:
module AffineArithmetic

using IntervalArithmetic, StaticArrays, LinearAlgebra, DiffRules
using DiffRules: diffrule

import IntervalArithmetic: interval
import Base: +, -, *, /, ^, ==, exp, exp2, exp10, expm1, log, log2,
             log10, log1p, acosh, zero, one, range, show, promote_rule, convert

abstract type AbstractLinearization end
struct MinRange <: AbstractLinearization end
struct Chebyshev <: AbstractLinearization end
const DEFAULT_LIN = MinRange

include("basic_affine1.jl")
#include("basic_affine2.jl")
#include("mixed_affine.jl")

const Affine = BasicAffine1

export Affine, BasicAffine1, MinRange, Chebyshev
export +, -, *, /, ^, ==, exp, exp2, exp10, expm1, log, log2,
       log10, log1p, acosh, zero, one, range, show, promote_rule, convert

end


text:
# -------------------
# Preamble
# -------------------

using LazySets
import Plots
using Plots: plot, plot!, text, lens!, bbox, font, savefig
using LaTeXStrings
using Plots.PlotMeasures

# custom plot shorthands
p(args...; kwargs...) = plot(xlims=(0, 4), ylims=(0, 4), ratio=1, lab="", xlab="x₁", ylab="x₂", args...; kwargs...)
p!(args...; kwargs...) = plot!(lab="", args...; kwargs...);

# utils for plots
__toL(x, digits) = L"%$(round(x, digits=digits))"
xticklatex(vec, digits) = (vec, __toL.(vec, Ref(digits)))
Plots.GR.setarrowsize(3)  # larger arrow heads


text:
import POMDPs: transition, reward, initialstate_distribution
import POMDPs: generate_sr, generate_o

struct W <: POMDP{Int, Bool, Int} end
println("Warning expected:")
@test_throws MethodError initialstate(W(), Random.GLOBAL_RNG)
println("Warning expected:")
@test_throws MethodError generate_s(W(), 1, true, Random.GLOBAL_RNG)
println("Warning expected:")
@test_throws MethodError generate_sr(W(), 1, true, Random.GLOBAL_RNG)
println("Warning expected:")
@test_throws MethodError generate_o(W(), 1, true, 2, Random.GLOBAL_RNG)
println("Warning expected:")
@test_throws MethodError generate_so(W(), 1, true, Random.GLOBAL_RNG)
println("Warning expected:")
@test_throws MethodError generate_sor(W(), 1, true, Random.GLOBAL_RNG)
println("Warning expected:")
@test_throws MethodError generate_or(W(), 1, true, 2, Random.GLOBAL_RNG)

mutable struct B <: POMDP{Int, Bool, Bool} end

transition(b::B, s::Int, a::Bool) = Int[s+a]
@test implemented(generate_s, Tuple{B, Int, Bool, MersenneTwister})
@test generate_s(B(), 1, false, Random.GLOBAL_RNG) == 1

@test !@implemented generate_sor(::B, ::Int, ::Bool, ::MersenneTwister)
# @test_throws MethodError generate_sor(B(), 1, false, Random.GLOBAL_RNG)

reward(b::B, s::Int, a::Bool, sp::Int) = -1.0
@test generate_sr(B(), 1, false, Random.GLOBAL_RNG) == (1, -1.0)

generate_o(b::B, s::Int, a::Bool, sp::Int, rng::AbstractRNG) = sp
@test @implemented generate_o(b::B, s::Int, a::Bool, sp::Int, rng::AbstractRNG)
@test @implemented generate_so(b::B, s::Int, a::Bool, rng::MersenneTwister)
@test @implemented generate_sor(b::B, s::Int, a::Bool, rng::MersenneTwister)
@test generate_sor(B(), 1, true, Random.GLOBAL_RNG) == (2, 2, -1.0) # should throw sor error

initialstate_distribution(b::B) = Int[1,2,3]
@test initialstate(B(), Random.GLOBAL_RNG) in initialstate_distribution(B())

mutable struct C <: POMDP{Nothing, Nothing, Nothing} end
generate_sr(c::C, s::Nothing, a::Nothing, rng::AbstractRNG) = nothing, 0.0
generate_o(c::C, s::Nothing, a::Nothing, sp::Nothing, rng::AbstractRNG) = nothing
@test @implemented generate_sor(::C, ::Nothing, ::Nothing, ::MersenneTwister)
@test generate_sor(C(), nothing, nothing, Random.GLOBAL_RNG) == (nothing, nothing, 0.0)


text:
# Chipmunk's rigid body type. Rigid bodies hold the physical properties of an object like
# it's mass, and position and velocity of it's center of gravity. They don't have an shape on their own.
# They are given a shape by creating collision shapes (cpShape) that point to the body.

@enum(BodyType,
    # A dynamic body is one that is affected by gravity, forces, and collisions.
    # This is the default body type.
dynamic,
    # A kinematic body is an infinite mass, user controlled body that is not affected by gravity, forces or collisions.
    # Instead the body only moves based on it's velocity.
    # Dynamic bodies collide normally with kinematic bodies, though the kinematic body will be unaffected.
    # Collisions between two kinematic bodies, or a kinematic body and a static body produce collision callbacks, but no collision response.
kinematic,
    # A static body is a body that never (or rarely) moves. If you move a static body, you must call one of the cpSpaceReindex*() functions.
    # Chipmunk uses this information to optimize the collision detection.
    # Static bodies do not produce collision callbacks when colliding with other static bodies.
static)

# Allocate a cpBody.
function Body()
    Body(ccall(dlsym(libchipmunk, :cpBodyAlloc), Ptr{Void}, ()))
end

# Allocate and initialize a cpBody.
function Body(mass::Real, moment::Real)
    Body(ccall(dlsym(libchipmunk, :cpBodyNew), Ptr{Void}, (Cdouble, Cdouble,), mass, moment))
end

# Allocate and initialize a cpBody, and set it as a kinematic body.
function KinematicBody()
    Body(ccall(dlsym(libchipmunk, :cpBodyNewKinematic), Ptr{Void}, ()))
end

# Allocate and initialize a cpBody, and set it as a static body.
function StaticBody()
    Body(ccall(dlsym(libchipmunk, :cpBodyNewStatic), Ptr{Void}, ()))
end

# Initialize a cpBody.
function init(body::Body, mass::Real, moment::Real)
    Body(ccall(dlsym(libchipmunk, :cpBodyInit), Ptr{Void}, (Ptr{Void}, Cdouble, Cdouble,), body.ptr, mass, moment))
end

# Destroy and free a cpBody.
function free(body::Body)
    ccall(dlsym(libchipmunk, :cpBodyFree), Void, (Ptr{Void},), body.ptr)
end

# Wake up a sleeping or idle body.
function activate(body::Body)
    ccall(dlsym(libchipmunk, :cpBodyActivate), Void, (Ptr{Void},), body.ptr)
end

# Wake up any sleeping or idle bodies touching a static body.
function activate_static(body::Body, filter::Shape)
    ccall(dlsym(libchipmunk, :cpBodyActivateStatic), Void, (Ptr{Void}, Ptr{Void},), body.ptr, filter.ptr)
end

# Force a body to fall asleep immediately.
function sleep(body::Body)
    ccall(dlsym(libchipmunk, cpBodySleep), Void, (Ptr{Void},), body.ptr)
end

# Force a body to fall asleep immediately along with other bodies in a group.
function sleep_with_group(body::Body, group::Body)
    ccall(dlsym(libchipmunk, :cpBodySleepWithGroup), Void, (Ptr{Void}, Ptr{Void},), body.ptr, group.ptr)
end

# Returns true if the body is sleeping.
function is_sleeping(body::Body)
    Bool(ccall(dlsym(libchipmunk, :cpBodyIsSleeping), Uint8, (Ptr{Void},), body.ptr))
end

# Get the type of the body.
function get_type(body::Body)
    BodyType(ccall(dlsym(libchipmunk, :cpBodyGetType), Int32, (Ptr{Void},), body.ptr))
end

# Set the type of the body.
function set_type(body::Body, body_type::BodyType)
    ccall(dlsym(libchipmunk, :cpBodySetType), Void, (Ptr{Void}, Int32,), body.ptr, body_type)
end

# function accumulate_mass_from_shapes(body::Body)
#   ccall(dlsym(libchipmunk, :cpBodyAccumulateMassFromShapes), Void, (Ptr{Void},), body.ptr)
# end

# Get the space this body is added to.
function get_space(body::Body)
    Space(ccall(dlsym(libchipmunk, :cpBodyGetSpace), Ptr{Void}, (Ptr{Void},), body.ptr))
end

# Get the mass of the body.
function get_mass(body::Body)
    Real(ccall(dlsym(libchipmunk, :cpBodyGetMass), Cdouble, (Ptr{Void},), body.ptr))
end

# Set the mass of the body.
function set_mass(body::Body, mass::Real)
    ccall(dlsym(libchipmunk, :cpBodySetMass), Void, (Ptr{Void}, Cdouble,), body.ptr, mass)
end

# Get the moment of inertia of the body.
function get_moment(body::Body)
    Real(ccall(dlsym(libchipmunk, :cpBodyGetMoment), Cdouble, (Ptr{Void},), body.ptr))
end

# Set the moment of inertia of the body.
function set_moment(body::Body, moment::Real)
    ccall(dlsym(libchipmunk, :cpBodySetMoment), Void, (Ptr{Void}, Cdouble,), body.ptr, moment)
end

# Get the position of a body.
function get_position(body::Body)
    ccall(dlsym(libchipmunk, :cpBodyGetPosition), Vect, (Ptr{Void},), body.ptr)
end

# Set the position of a body.
function get_center_of_gravity(body::Body)
    ccall(dlsym(libchipmunk, :cpBodyGetCenterOfGravity), Vect, (Ptr{Void},), body.ptr)
end

# Get the offset of the center of gravity in body local coordinates.
function set_center_of_gravity(body::Body, cog::Vect)
    ccall(dlsym(libchipmunk, :cpBodySetCenterOfGravity), Void, (Ptr{Void}, Vect,), body.ptr, cog)
end

# Set the offset of the center of gravity in body local coordinates.<F37>
function set_position(body::Body, pos::Vect)
    ccall(dlsym(libchipmunk, :cpBodySetPosition), Void, (Ptr{Void}, Vect,), body.ptr, pos)
end

# Get the velocity of the body.
function get_velocity(body::Body)
    ccall(dlsym(libchipmunk, :cpBodyGetVelocity), Vect, (Ptr{Void},), body.ptr)
end

# Set the velocity of the body.
function set_velocity(body::Body, velocity::Vect)
    ccall(dlsym(libchipmunk, :cpBodySetVelocity), Void, (Ptr{Void}, Vect), body.ptr, velocity)
end

# Get the force applied to the body for the next time step.
function get_force(body::Body)
    ccall(dlsym(libchipmunk, :cpBodyGetForce), Vect, (Ptr{Void},), body.ptr)
end

# Set the force applied to the body for the next time step.
function set_force(body::Body, force::Vect)
    ccall(dlsym(libchipmunk, :cpBodySetForce), Void, (Ptr{Void}, Vect,), body.ptr, force)
end

# Get the angle of the body.
function get_angle(body::Body)
    ccall(dlsym(libchipmunk, :cpBodyGetAngle), Cdouble, (Ptr{Void},), body.ptr)
end

# Set the angle of a body.
function set_angle(body::Body, angle::Real)
    ccall(dlsym(libchipmunk, :cpBodySetAngle), Void, (Ptr{Void}, Cdouble), body.ptr, angle)
end

# Get the angular velocity of the body.
function get_angular_velocity(body::Body)
    ccall(dlsym(libchipmunk, :cpBodyGetAngularVelocity), Cdouble, (Ptr{Void},), body.ptr)
end

# Set the angular velocity of the body.
function set_angular_velocity(body::Body, angular_velocity::Real)
    ccall(dlsym(libchipmunk, :cpBodySetAngularVelocity), Void, (Ptr{Void}, Cdouble,), body.ptr, angular_velocity)
end

# Get the torque applied to the body for the next time step.
function get_torque(body::Body)
    ccall(dlsym(libchipmunk, :cpBodyGetTorque), Cdouble, (Ptr{Void},), body.ptr)
end

# Set the torque applied to the body for the next time step.
function set_torque(body::Body, torque::Real)
    ccall(dlsym(libchipmunk, :cpBodySetTorque), Void, (Ptr{Void}, Cdouble,), body.ptr)
end

# Get the rotation vector of the body. (The x basis vector of it's transform.)
function get_rotation(body::Body)
    ccall(dlsym(libchipmunk, :cpBodyGetRotation), Vect, (Ptr{Void},), body.ptr)
end

# Get the user data pointer assigned to the body.
function get_userdata(body::Body)
    ccall(dlsym(libchipmunk, :cpBodyGetUserData), Ptr{Void}, (Ptr{Void},), body.ptr)
end

# Set the user data pointer assigned to the body.
function set_userdata(body::Body, userdata::Ptr{Void})
    ccall(dlsym(libchipmunk, :cpBodySetUserData), Void, (Ptr{Void}, Ptr{Void},), body.ptr, userdata)
end

# Default velocity integration function.
function update_velocity(body::Body, gravity::Vect, damping::Real, dt::Real)
    ccall(dlsym(libchipmunk, :cpBodyUpdateVelocity), Void, (Ptr{Void}, Vect, Cdouble, Cdouble,), body.ptr, gravity, damping, dt)
end

# Default position integration function.
function update_position(body::Body, dt::Real)
    ccall(dlsym(libchipmunk, :cpBodyUpdatePosition), Void, (Ptr{Void}, Cdouble,), body.ptr, dt)
end

# Convert body relative/local coordinates to absolute/world coordinates.
function local_to_world(body::Body, point::Vect)
    ccall(dlsym(libchipmunk, :cpBodyLocalToWorld), Vect, (Ptr{Void}, Vect,), body.ptr, point)
end

# Convert body absolute/world coordinates to  relative/local coordinates.
function world_to_local(body::Body, point::Vect)
    ccall(dlsym(libchipmunk, :cpBodyWorldToLocal), Vect, (Ptr{Void}, Vect,), body.ptr, point)
end

# Apply a force to a body. Both the force and point are expressed in world coordinates.
function apply_force_at_world_point(body::Body, force::Vect, point::Vect)
    ccall(dlsym(libchipmunk, :cpBodyApplyForceAtWorldPoint), Void, (Ptr{Void}, Vect, Vect,), body.ptr, force, point)
end

# Apply a force to a body. Both the force and point are expressed in body local coordinates.
function apply_force_at_local_point(body::Body, force::Vect, point::Vect)
    ccall(dlsym(libchipmunk, :cpBodyApplyForceAtLocalPoint), Void, (Ptr{Void}, Vect, Vect,), body.ptr, force, point)
end

# Apply an impulse to a body. Both the impulse and point are expressed in world coordinates.
function apply_impulse_at_world_point(body::Body, impulse::Vect, point::Vect)
    ccall(dlsym(libchipmunk, :cpBodyApplyImpulseAtWorldPoint), Void, (Ptr{Void}, Vect, Vect,), body.ptr, impulse, point)
end

# Apply an impulse to a body. Both the impulse and point are expressed in body local coordinates.
function apply_impulse_at_local_point(body::Body, impulse::Vect, point::Vect)
    ccall(dlsym(libchipmunk, :cpBodyApplyImpulseAtLocalPoint), Void, (Ptr{Void}, Vect, Vect,), body.ptr, impulse, point)
end

# Get the velocity on a body (in world units) at a point on the body in world coordinates.
function get_velocity_at_world_point(body::Body, point::Vect)
    ccall(dlsym(libchipmunk, :cpBodyGetVelocityAtWorldPoint), Vect, (Ptr{Void}, Vect,), body.ptr, point)
end

# Get the velocity on a body (in local units) at a point on the body in local coordinates.
function get_velocity_at_local_point(body::Body, point::Vect)
    ccall(dlsym(libchipmunk, :cpBodyGetVelocityAtLocalPoint), Vect, (Ptr{Void}, Vect,), body.ptr, point)
end

# Get the amount of kinetic energy contained by the body.
function kinetic_energy(body::Body)
    ccall(dlsym(libchipmunk, :cpBodyKineticEnergy), Cdouble, (Ptr{Void},), body.ptr)
end

export Body, KinematicBody, StaticBody, BodyType, init, free, is_sleeping, get_type, set_type,
get_space, get_mass, set_mass, get_moment, set_moment, set_position, get_position, get_velocity,
set_velocity, activate, activate_static, sleep, get_angle, set_angle, get_center_of_gravity,
set_center_of_gravity, get_force, set_force, get_angular_velocity, set_angular_velocity,
get_torque, set_torque, get_rotation, get_userdata, set_userdata, update_velocity, update_position,
local_to_world, world_to_local, apply_impulse_at_world_point, apply_force_at_local_point,
apply_impulse_at_world_point, apply_impulse_at_local_point, get_velocity_at_world_point,
get_velocity_at_local_point, kinetic_energy


text:
#*************************************************************************#
#
# Bias Type
#
#*************************************************************************#
"""
    Bias(n)

Initialize `ParametricModel` representing bias operation `z = x + Θ`.
`n` is the input vector dimension.
"""
struct Bias <: ParametricModel
    ninputs::Int
    noutputs::Int
    nparams::Int
end

function Bias(n::Int)
    Bias(n, n, n)
end

#*************************************************************************#
#
# Generate Buffers
#
#*************************************************************************#
function dx_buffer(bias::Bias)
    return Diagonal{Float64}(zeros(bias.noutputs))
end

function dΘ_buffer(bias::Bias)
    return Diagonal{Float64}(zeros(bias.nparams))
end

#*************************************************************************#
#
# Compute Output and its Derivatives
#
#*************************************************************************#
function _evaluate!(bias::Bias,
                    x::AbstractVector{<:Real},
                    Θ::AbstractVector{<:Real},
                    z::AbstractVector{Float64},
                    dx::AbstractMatrix{Float64},
                    dΘ::AbstractMatrix{Float64})

    # Compute z
    z .= x + Θ

    # Compute dx
    if !isempty(dx)
        dx.diag .= 1
    end

    # Compute dΘ
    if !isempty(dΘ)
        dΘ.diag .= 1
    end
end


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

function message(c::GitCommit, raw::Bool=false)
    local msg_ptr::Cstring
    msg_ptr = raw? ccall((:git_commit_message_raw, :libgit2), Cstring, (Ptr{Void},), c.ptr) :
                   ccall((:git_commit_message, :libgit2), Cstring, (Ptr{Void},), c.ptr)
    if msg_ptr == Cstring_NULL
        return nothing
    end
    return unsafe_string(msg_ptr)
end

function author(c::GitCommit)
    ptr = ccall((:git_commit_author, :libgit2), Ptr{SignatureStruct}, (Ptr{Void},), c.ptr)
    @assert ptr != C_NULL
    return Signature(ptr)
end

function committer(c::GitCommit)
    ptr = ccall((:git_commit_committer, :libgit2), Ptr{SignatureStruct}, (Ptr{Void},), c.ptr)
    @assert ptr != C_NULL
    return Signature(ptr)
end

""" Wrapper around `git_commit_create` """
function commit(repo::GitRepo,
                refname::AbstractString,
                msg::AbstractString,
                author::GitSignature,
                committer::GitSignature,
                tree::GitTree,
                parents::GitCommit...)
    commit_id_ptr = Ref(Oid())
    nparents = length(parents)
    parentptrs = Ptr{Void}[c.ptr for c in parents]
    @check ccall((:git_commit_create, :libgit2), Cint,
                 (Ptr{Oid}, Ptr{Void}, Ptr{UInt8},
                  Ptr{SignatureStruct}, Ptr{SignatureStruct},
                  Ptr{UInt8}, Ptr{UInt8}, Ptr{Void},
                  Csize_t, Ptr{Ptr{Void}}),
                 commit_id_ptr, repo.ptr, isempty(refname) ? C_NULL : refname,
                 author.ptr, committer.ptr,
                 C_NULL, msg, tree.ptr,
                 nparents, nparents > 0 ? parentptrs : C_NULL)
    return commit_id_ptr[]
end

"""Commit changes to repository"""
function commit(repo::GitRepo, msg::AbstractString;
                refname::AbstractString=Consts.HEAD_FILE,
                author::Signature = Signature(repo),
                committer::Signature = Signature(repo),
                tree_id::Oid = Oid(),
                parent_ids::Vector{Oid}=Oid[])
    # Retrieve tree identifier
    if iszero(tree_id)
        tree_id = with(GitIndex, repo) do idx; write_tree!(idx) end
    end

    # Retrieve parents from HEAD
    if isempty(parent_ids)
        try # if throws then HEAD not found -> empty repo
            push!(parent_ids, Oid(repo, refname))
        end
    end

    # return commit id
    commit_id  = Oid()

    # get necessary objects
    tree = get(GitTree, repo, tree_id)
    auth_sig = convert(GitSignature, author)
    comm_sig = convert(GitSignature, committer)
    parents = GitCommit[]
    try
        for parent in parent_ids
            push!(parents, get(GitCommit, repo, parent))
        end
        commit_id = commit(repo, refname, msg, auth_sig, comm_sig, tree, parents...)
    finally
        for parent in parents
            finalize(parent)
        end
        finalize(tree)
        finalize(auth_sig)
        finalize(comm_sig)
    end
    return commit_id
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
abstract type AbstractBinningAlgorithm end

mutable struct Bin{T}
    """Number of samples."""
    nsamples::Int
    """Mean of predictions."""
    mean_predictions::T
    """Proportions of targets."""
    proportions_targets::T

    function Bin{T}(
        nsamples::Int, mean_predictions::T, proportions_targets::T
    ) where {T<:Real}
        nsamples ≥ 0 || throw(ArgumentError("the number of samples must be non-negative"))
        return new{T}(nsamples, mean_predictions, proportions_targets)
    end

    function Bin{T}(
        nsamples::Int, mean_predictions::T, proportions_targets::T
    ) where {T<:AbstractVector{<:Real}}
        nsamples ≥ 0 || throw(ArgumentError("the number of samples must be non-negative"))
        nclasses = length(mean_predictions)
        nclasses > 1 || throw(ArgumentError("the number of classes must be greater than 1"))
        nclasses == length(proportions_targets) || throw(
            DimensionMismatch(
                "the number of predicted classes has to be equal to the number of classes",
            ),
        )

        return new{T}(nsamples, mean_predictions, proportions_targets)
    end
end

function Bin(nsamples::Int, mean_predictions::T, proportions_targets::T) where {T}
    return Bin{T}(nsamples, mean_predictions, proportions_targets)
end

"""
    Bin(predictions, targets)

Create bin of `predictions` and corresponding `targets`.
"""
function Bin(predictions::AbstractVector{<:Real}, targets::AbstractVector{Bool})
    # compute mean of predictions
    mean_predictions = mean(predictions)

    # compute proportion of targets
    proportions_targets = mean(targets)

    return Bin(length(predictions), mean_predictions, proportions_targets)
end
function Bin(
    predictions::AbstractVector{<:AbstractVector{<:Real}},
    targets::AbstractVector{<:Integer},
)
    # compute mean of predictions
    mean_predictions = mean(predictions)

    # compute proportion of targets
    nclasses = length(predictions[1])
    proportions_targets = StatsBase.proportions(targets, nclasses)

    return Bin(length(predictions), mean_predictions, proportions_targets)
end

"""
    Bin(prediction, target)

Create bin of a single `prediction` and corresponding `target`.
"""
function Bin(prediction::Real, target::Bool)
    # compute mean of predictions
    mean_predictions = prediction / 1

    # compute proportion of targets
    proportions_targets = target / 1

    return Bin(1, mean_predictions, proportions_targets)
end

function Bin(prediction::AbstractVector{<:Real}, target::Integer)
    # compute mean of predictions
    mean_predictions = prediction ./ 1

    # compute proportion of targets
    proportions_targets = similar(mean_predictions)
    for i in 1:length(proportions_targets)
        proportions_targets[i] = i == target ? 1 : 0
    end

    return Bin(1, mean_predictions, proportions_targets)
end

"""
    adddata!(bin::Bin, prediction, target)

Update running statistics of the `bin` by integrating one additional pair of `prediction`s
and `target`.
"""
function adddata!(bin::Bin, prediction::Real, target::Bool)
    @unpack mean_predictions, proportions_targets = bin

    # update number of samples
    nsamples = (bin.nsamples += 1)

    # update mean of predictions
    mean_predictions += (prediction - mean_predictions) / nsamples
    bin.mean_predictions = mean_predictions

    # update proportions of targets
    proportions_targets += (target - proportions_targets) / nsamples
    bin.proportions_targets = proportions_targets

    return nothing
end
function adddata!(bin::Bin, prediction::AbstractVector{<:Real}, target::Integer)
    @unpack mean_predictions, proportions_targets = bin

    # update number of samples
    nsamples = (bin.nsamples += 1)

    # update mean of predictions
    @. mean_predictions += (prediction - mean_predictions) / nsamples

    # update proportions of targets
    nclasses = length(proportions_targets)
    @. proportions_targets += ((1:nclasses == target) - proportions_targets) / nsamples

    return nothing
end


text:
struct Error <: Exception
    msg::String
end

struct ES
    url::String
    ES(;
        host::String = "localhost",
        port::Integer = 9200,
        user::String = "",
        password::String = "",
        protocol::String = "http",
    ) = new("$protocol://$user:$password@$host:$port")
    ES(url) = new(url)
end


text:
module QuantumOpt

# """
# QuantumOpt:
#
# Module for optimizing quantum measures.
# """

using Convex, SCS, LinearAlgebra

using ..QuantumBehavior

# """
# cvx_prepare_and_measure_max_facet_violation_states(scenario, facet, states):
#
#     Performs semidefinite programming to find the set of quantum measurements
#     which optimally violate the provided facet inequality for the specified
#     quantum states.
#
#   ` This uses cvx optimizer which handles complex-valued hermitian matrices.
#
# Input:
#     num_in/out: Integer, number of inputs/outputs in the scenario
#     facet: Row vector, represents a facet inequality
#     states: Array, contains valid density matrices representing encoded states
#
# Output: Dictionary with keys,
#     violation: Float, the maximal value achieved by quantum strategies. If
#                violation > 0, quantum strategy is nonlocal, violation == 0,
#                strategy lies on facet, violation < 0, strategy lies in polytope
#                volume.
#     povm: array of matrices, the optimal povm.
# """
function cvx_prepare_and_measure_max_facet_violation_states(num_in, num_out, facet, states)
    bound = facet[1]
    dits = size(states[1])[1]

    state_groups = map(i -> sum(facet[(1+1+num_in*(i-1)):(1+num_in+num_in*(i-1))].*states), 1:(num_out - 1))

    # add povm variables and constraints
    povm_variables = map(i -> HermitianSemidefinite(dits), 1:num_out)
    constraints = (sum(map(pv -> real(pv), povm_variables)) == [1 0;0 1])
    constraints += (sum(map(pv -> imag(pv), povm_variables)) == [0 0;0 0])

    # add the objective
    objective = maximize(real(tr(sum(state_groups.*povm_variables[1:(end-1)]))), constraints)

    # optimize model
    solve!(objective, SCS.Optimizer(verbose=0))

    # parse/return results
    violation = bound + objective.optval

    Dict("violation" => violation, "povm" => map(pv -> pv.value, povm_variables))
end

# """
# cvx_prepare_and_measure_max_facet_violation_states2(scenario, facet, states):
#
#     Performs semidefinite programming to find the set of quantum measurements
#     which optimally violate the provided facet inequality for the specified
#     quantum states. This method solves for two povms used.
#
#   ` This uses cvx optimizer which handles complex-valued hermitian matrices.
#
# Input:
#     num_in/out: Integer, number of inputs/outputs in the scenario
#     facet: Row vector, represents a facet inequality
#     states: Array, contains valid density matrices representing encoded states
#
# Output: Dictionary with keys,
#     violation: Float, the maximal value achieved by quantum strategies. If
#                violation > 0, quantum strategy is nonlocal, violation == 0,
#                strategy lies on facet, violation < 0, strategy lies in polytope
#                volume.
#     povms: array of povms, the optimal povm pair.
# """
function cvx_prepare_and_measure_max_facet_violation_states2(α, β, facet, states)
    bound = facet[1]
    dits = size(states[1])[1]

    (α_in, α_out) = α
    (β_in, β_out) = β

    state_groups1 = map(i -> sum(facet[(1+1+α_in*β_in*(i-1)):β_in:(1+α_in*β_in+α_in*β_in*(i-1))].*states), 1:(β_out - 1))
    state_groups2 = map(i -> sum(facet[(1+2+α_in*β_in*(i-1)):β_in:(2+α_in*β_in+α_in*β_in*(i-1))].*states), 1:(β_out - 1))

    povm_variables1 = map(i -> HermitianSemidefinite(dits), 1:β_out)
    povm_variables2 = map(i -> HermitianSemidefinite(dits), 1:β_out)

    constraints = (sum(map(pv -> real(pv), povm_variables1)) == [1 0;0 1])
    constraints += (sum(map(pv -> imag(pv), povm_variables1)) == [0 0;0 0])
    constraints += (sum(map(pv -> real(pv), povm_variables2)) == [1 0;0 1])
    constraints += (sum(map(pv -> imag(pv), povm_variables2)) == [0 0;0 0])

    # add the objective
    objective = maximize(real(tr(sum(state_groups1.*povm_variables1[1:(end-1)] + state_groups2.*povm_variables2[1:(end-1)]))), constraints)

    # optimize model
    solve!(objective, SCS.Optimizer(verbose=0))

    # parse/return results
    violation = bound + objective.optval

    Dict("violation" => violation, "povms" => [map(pv -> pv.value, povm_variables1), map(pv -> pv.value, povm_variables2)])
end

# """
# cvx_prepare_and_measure_max_facet_violation_povm(scenario, facet, povm):
#
#     Performs semidefinite programming to find the set of quantum states
#     which optimally violate the provided facet inequality for the specified
#     quantum measurement.
#
#   ` This uses cvx optimizer which handles complex-valued hermitian matrices.
#
# Input:
#     num_in/out: Integer, number of inputs/outputs in the scenario
#     facet: Row vector, represents a facet inequality
#     povm: Array, contains valid povm to optimize against
#
# Output: Dictionary with keys
#     violation: Float, the maximal value achieved by quantum strategies. If
#                violation > 0, quantum strategy is nonlocal, violation == 0,
#                strategy lies on facet, violation < 0, strategy lies in polytope
#                volume.
#     states: Array of density matrices, the optimal density matrices.
# """
function cvx_prepare_and_measure_max_facet_violation_povm(num_in, num_out, facet, povm)
    bound = facet[1]
    dits = size(povm[1])[1]

    povm_groups = map(i -> sum(facet[(1+i):num_in:end].*povm[1:(end-1)]), 1:num_in)

    # add povm variables and constraints
    ρ_variables = map(i -> HermitianSemidefinite(dits), 1:num_in)
    constraints = map(ρ -> real(tr(ρ)) == 1, ρ_variables)

    # add the objective
    objective = maximize(real(tr(sum(povm_groups.*ρ_variables))), constraints)

    # optimize model
    solve!(objective, SCS.Optimizer(verbose=0))

    # parse/return results
    violation = bound + objective.optval

    Dict("violation" => violation, "states" => map(ρv -> ρv.value, ρ_variables))
end

end


text:
module UnicodeNames

include("./unicodedata.jl")

function lookupchar(str::ASCIIString)
    res = get(unicode16bit, str, 0xd800)
    res != 0xd800 && return UInt32(res)
    get(unicode32bit, str, typemax(UInt32))
end

function loadnames(filename)
    v = Vector{ASCIIString}()
    open(filename, "r") do f
        while true
            l = readline(f)
            l == "" && break
            push!(v, l[1:end-1])
        end
    end
    v
end

function createdicts(filename, v::Vector{ASCIIString})
    open(filename, "w+") do f
        println(f, "const unicode16bit = Dict{ASCIIString, UInt16}(")
        flg16 = true
        for l in v
            vec = split(l, ";")
            num = vec[1]
            str = vec[2]
            alias = vec[11]
            if flg16 && parse(UInt32, num, 16) > 0x0ffff
                println(f, ")\nconst unicode32bit = Dict{ASCIIString, UInt32}(")
                flg16 = false
            end
            if str[1] != '<'
                # Don't save names that simply contain hex representation
                len = length(num)
                if !(str[end-len-1] == '-' && str[end-len:end] == num)
                    println(f, "    \"",str,"\" => 0x",num,",")
                end
            end
            if alias != ""
                println(f, "    \"",alias,"\" => 0x",num,",")
            end
        end
        println(f, ")")
    end
end

end


text:
module CovidDash

using Matte

const title = "Covid evolution :: Live dashboard"
south_america = "['Brazil','Uruguay', 'Argentina','Chile','Paraguay','Peru','Bolivia','Ecuador','Colombia','Venezuela','Suriname']"

function ui()
    sidebar_layout(
        side_panel(
            h1("Choose settings"),
	    p("Source"),
            selector("country", "Country", south_america),
	    selector("all_opt","All","['Yes','No']"),
            br()),
        main_panel(
            plots_output("output_plot")
        )
	)
end

module Server

using StatsPlots , Plots
using DataFrames, CSV , Query

covid_data = CSV.read("data/covid-long16Mar.csv")

function output_plot(country,all_opt)
    if all_opt == "Yes"
	south_america_l = ["Brazil","Uruguay", "Argentina","Chile","Paraguay","Peru","Bolivia","Ecuador","Colombia","Venezuela","Suriname"]
	covid_plot = covid_data |> @filter(_.Country ∈ south_america_l) |> DataFrame
        @df covid_plot scatter(:date, :val, group = :Country)
        @df covid_plot plot!(:date, :val, group=:Country,seriestype = :steppre,legend=:bottomleft)
    else
	if typeof(country) <: AbstractString
		covid_plot = covid_data |> @filter(_.Country ∈ [country]) |> DataFrame
        	@df covid_plot scatter(:date, :val, group = :Country)
 		@df covid_plot plot!(:date, :val, group=:Country,seriestype = :steppre,legend=:bottomleft)
    end
end

end

end
end

max_stars_repo_name: rafaelfrc/covid-dash

max_issues_repo_name: rafaelfrc/covid-dash

max_forks_repo_name: rafaelfrc/covid-dash


text:
@enum ShaderLanguage begin
    SPIR_V = 1
    GLSL   = 2
    HLSL   = 3
end

function file_ext(language::ShaderLanguage, stage::Vk.ShaderStageFlag)
    language in [GLSL, HLSL] || error("Cannot get file extension for language $language")
    @match stage begin
        &Vk.SHADER_STAGE_VERTEX_BIT => "vert"
        &Vk.SHADER_STAGE_FRAGMENT_BIT => "frag"
        &Vk.SHADER_STAGE_TESSELLATION_CONTROL_BIT => "tesc"
        &Vk.SHADER_STAGE_TESSELLATION_EVALUATION_BIT => "tese"
        &Vk.SHADER_STAGE_GEOMETRY_BIT => "geom"
        &Vk.SHADER_STAGE_COMPUTE_BIT => "comp"
        &Vk.SHADER_STAGE_RAYGEN_BIT_KHR => "rgen"
        &Vk.SHADER_STAGE_INTERSECTION_BIT_KHR => "rint"
        &Vk.SHADER_STAGE_ANY_HIT_BIT_KHR => "rahit"
        &Vk.SHADER_STAGE_CLOSEST_HIT_BIT_KHR => "rchit"
        &Vk.SHADER_STAGE_MISS_BIT_KHR => "rmiss"
        &Vk.SHADER_STAGE_CALLABLE_BIT_KHR => "rcall"
        &Vk.SHADER_STAGE_MESH_BIT_NV => "mesh"
        &Vk.SHADER_STAGE_TASK_BIT_NV => "task"
        _ => error("Unknown stage $stage")
    end
end

"""
    shader_stage(file_ext)

Automatically detect a shader stage from a file extension.
Can only be used with [`GLSL`](@ref) and [`HLSL`](@ref).

## Examples

```jldoctest
julia> shader_stage("my_shader.frag", GLSL) == Vk.SHADER_STAGE_FRAGMENT_BIT
true

julia> shader_stage("my_shader.geom", HLSL) == Vk.SHADER_STAGE_GEOMETRY_BIT
true
```
"""
function shader_stage(file::AbstractString, language::ShaderLanguage)
    _, file_ext = splitext(file)
    language in [GLSL, HLSL] || error("Cannot retrieve shader stage from shader language $language")
    @match file_ext begin
        ".vert" => Vk.SHADER_STAGE_VERTEX_BIT
        ".frag" => Vk.SHADER_STAGE_FRAGMENT_BIT
        ".tesc" => Vk.SHADER_STAGE_TESSELLATION_CONTROL_BIT
        ".tese" => Vk.SHADER_STAGE_TESSELLATION_EVALUATION_BIT
        ".geom" => Vk.SHADER_STAGE_GEOMETRY_BIT
        ".comp" => Vk.SHADER_STAGE_COMPUTE_BIT
        ".rgen" => Vk.SHADER_STAGE_RAYGEN_BIT_KHR
        ".rint" => Vk.SHADER_STAGE_INTERSECTION_BIT_KHR
        ".rahit" => Vk.SHADER_STAGE_ANY_HIT_BIT_KHR
        ".rchit" => Vk.SHADER_STAGE_CLOSEST_HIT_BIT_KHR
        ".rmiss" => Vk.SHADER_STAGE_MISS_BIT_KHR
        ".rcall" => Vk.SHADER_STAGE_CALLABLE_BIT_KHR
        ".mesh" => Vk.SHADER_STAGE_MESH_BIT_NV
        ".task" => Vk.SHADER_STAGE_TASK_BIT_NV
        _ => error("Unknown file extension $file_ext")
    end
end

"""
    shader_language(file_ext)

Automatically detect a [`ShaderLanguage`](@ref) from the file extension.
Currently, only .spv, .glsl and .hlsl are recognized.

## Examples

```julia
julia> shader_language("my_shader.glsl")
GLSL
julia> shader_language("my_shader.spv")
SPIR_V
```
"""
function shader_language(file::AbstractString)
    _, file_ext = splitext(file)
    @match file_ext begin
        ".spv" => SPIR_V
        ".glsl" => GLSL
        ".hlsl" => HLSL
        _ => error("Unknown file extension $file_ext")
    end
end


text:
function rocket100()
    model = JuMP.Model()
    # Parsing GAMS Files
    lexed = GAMSFiles.lex(OCTHaGOn.GAMS_DIR * "rocket100.gms")
    gams = GAMSFiles.parsegams(OCTHaGOn.GAMS_DIR * "rocket100.gms")
    GAMSFiles.parseconsts!(gams)

    vars = GAMSFiles.getvars(gams["variables"])
    sets = Dict{String, Any}()
    preexprs, bodyexprs = Expr[], Expr[]

    # Getting variables
    vardict, constdict = generate_variables!(model, gams) # Actual JuMP variables
    objvar = model[:objvar]
    x = [objvar] # :x1 is the objective variable
    symbs = [Symbol("x$(i)") for i=2:608]
    append!(x, [model[symb] for symb in symbs])
    model[:x] = x
    @objective(model, Min, objvar)

    # Initializing GM
    gm = GlobalModel(model = model, name = "rocket100.gms")

    # (Re)writing constraints
    add_nonlinear_constraint(gm, :(x -> -310*sqrt(x[3])*exp(500 - 500*x[104]) + x[508]),
                             name = "e1", vars = [x[3], x[104], x[508]], equality = true)
    for i=1:100
        add_linked_constraint(gm, gm.bbls[end], [x[3+i], x[104+i], x[508+i]])
    end

    add_nonlinear_constraint(gm, :(x -> -sqrt(1/x[104]) + x[205]),
                             name = "e102", vars = [x[104], x[205]], equality = true)
    for i=1:100
        add_linked_constraint(gm, gm.bbls[end], [x[104+i], x[205+i]])
    end

    @constraint(gm.model, -objvar >= x[204])
    @objective(gm.model, Min, objvar)

    add_nonlinear_constraint(gm, :(x -> -0.5*x[2]*(x[3] + x[4]) - x[104] + x[105]),
                             name = "e204", vars = [x[2], x[3], x[4], x[104], x[105]], equality = true)
    for i=1:99
        add_linked_constraint(gm, gm.bbls[end], [x[2], x[3+i], x[4+i], x[104+i], x[105+i]])
    end

    add_nonlinear_constraint(gm, :(x -> -0.5*((x[408] - x[307]*x[206] - x[509])/x[307] + (x[407] - x[306]*x[205] - x[508])/x[306]) *
                                   x[2] - x[3] + x[4]),
                                   name = "e304", vars = [x[2], x[3], x[4], x[205], x[206], x[306],
                                                          x[307], x[407], x[408], x[508], x[509]], equality = true)

    for i=1:99
        add_linked_constraint(gm, gm.bbls[end], [x[2], x[3+i], x[4+i], x[205+i], x[206+i], x[306+i],
                                                 x[307+i], x[407+i], x[408+i], x[508+i], x[509+i]])
    end

    add_nonlinear_constraint(gm, :(x -> x[2]*(x[407] + x[408]) - x[306] + x[307]),
                                   name = "e404", vars = [x[2], x[306], x[307], x[407], x[408]], equality = true)

    for i=1:99
        add_linked_constraint(gm, gm.bbls[end], [x[2], x[306+i], x[307+i], x[407+i], x[408+i]])
    end

    return gm
end


text:
using SeisPlot, PyPlot, SeisAcoustic

# create a simple two layers velocity model with constant density,
nz = 100; nx = 301;
vel = 2800 * ones(nz, nx);  # m/s
vel[31:end,:] .= 3500;      # m/s
rho = 2.0 * ones(nz, nx);   # kg/m^3
rho[61:end,:] .= 2.5;       # m/s

# plotting
SeisPlotTX(vel, wbox=9, hbox=5, title="velocity", cmap="rainbow", vmin=minimum(vel), vmax=maximum(vel));
SeisPlotTX(rho, wbox=9, hbox=5, title="density", cmap="rainbow", vmin=minimum(rho), vmax=maximum(rho));

# number of PML layers
npml = 20;

# top boundary condition
free_surface = true;        #(pml or free_surface)

# vertical and horizontal grid size
dz = 10; dx = 10;

# time step size and maximum modelling length
dt = 0.001; tmax = 1.25;     # use second as unit

# organize these parameters into a structure
# shows the default value for the keyword parameters
# data_format = (Float32 or Float64)
# fd_flag     = ("taylor" or "ls")
# order       = 2 - 10 if we use "ls" to compute the FD coefficients
#             = 2 - n  if we use "taylor" expansion to compute FD coefficients
params = TdParams(rho, vel, free_surface, dz, dx, dt, tmax;
         data_format=Float64, fd_flag="taylor", order=3, npml=20, apml=900.);
# type "?TdParams" to get more information about this structure

# ==============================================================================
#                          single source simulation
# ==============================================================================
isz = 2; isx = floor(Int64, params.nx/2);
# type_flag = "ricker" or "miniphase" or a vector
# location_flag = "index" or "distance"
src = Source(isz, isx, params; ot=0.0, fdom=15.0,
      type_flag="miniphase", amp=100000, location_flag="index");

# receiver location
irx = collect(1:2:params.nx);
irz = 2 * ones(Int64, length(irx));

# get one recordings
rec = multi_step_forward(src, params; rz=irz, rx=irx);
SeisPlotTX(rec.p, dy=0.001, title="one shot")
# type "multi_step_forward" to get more information about forward modelling

# save wavefield every 10 steps
# boundary of wavefield which are used for imaging
# snapshot(vz, vx, pz, px), strength of source side wavefield and so on
path_pre = joinpath(homedir(), "Desktop/pressure.rsf");
path_bnd = joinpath(homedir(), "Desktop/boundary.rsf");
multi_step_forward(src, params; path_pre=path_pre, path_bnd=path_bnd, interval=10);

# read the wavefield
(hdr, pre) = read_RSdata(path_pre);
path_anim  = joinpath(homedir(), "Desktop/wave_propagation");
vmin = minimum(pre[:,:,20]); vmax = maximum(pre[:,:,20]);
SeisAnimation(path_anim, pre; wbox=9, hbox=5, vmin=vmin, vmax=vmax, dx=10, dy=10)

# save every thing generated during forward modelling
# path_shot= joinpath(homedir(), "Desktop/recordings.bin" );
# path_spt = joinpath(homedir(), "Desktop/snapshot.rsf" );
# path_wfd = joinpath(homedir(), "Desktop/wavefield.rsf");
# path_pre = joinpath(homedir(), "Desktop/pressure.rsf" );
# path_bnd = joinpath(homedir(), "Desktop/boundary.rsf" );
# path_lwfd= joinpath(homedir(), "Desktop/last_wfd.rsf" );
# path_sws = joinpath(homedir(), "Desktop/strength.rsf" );
# multi_step_forward(src, params; rz=irz, rx=irx, path_shot=path_shot, path_spt=path_spt, path_wfd=path_wfd, path_pre=path_pre, interval=1,
#                    path_bnd=path_bnd, path_lwfd=path_lwfd, path_sws=path_sws);

# ==============================================================================
#                          simultaneouse source
# ==============================================================================
isx = collect(5:40:params.nx); ns=length(isx); isz = 2*ones(ns);
ot  = 0.3*rand(ns);
# ricker wavelet as the default one
src = get_multi_sources(isz, isx, params; amp=100000, ot=ot, fdom=15);

# forward modeling of simultaneous sources
path_pre = joinpath(homedir(), "Desktop/pressure_semi.rsf");
rec1 = multi_step_forward(src, params;rz=irz, rx=irx, path_pre=path_pre, interval=10);
SeisPlotTX(rec1.p, dy=0.001, title="semitaneous source")

(hdr, pre) = read_RSdata(path_pre);
path_anim  = joinpath(homedir(), "Desktop/wave_propagation_semi");
vmin = minimum(pre[:,:,50]); vmax = maximum(pre[:,:,50]);
SeisAnimation(path_anim, pre; wbox=9, hbox=5, vmin=vmin, vmax=vmax, dx=10, dy=10)

# ==============================================================================
#                multiple shot simulation
# ==============================================================================
isx = collect(120:30:params.nx); ns=length(isx); isz = 2*ones(ns);
src = get_multi_sources(isz, isx, params; amp=100000, fdom=15);

# simulate OBN acquisition with fixed receiver location
irx = collect(1:2:params.nx);
irz = 2 * ones(Int64, length(irx));

# do parallel simulation
dir_obs = joinpath(homedir(), "Desktop/shotgather")
get_shotgather(dir_obs, irz, irx, src, params);

path_shot = joinpath(dir_obs, "recordings_3.bin");
rec       = read_recordings(path_shot)
SeisPlotTX(rec.p, dy=0.001, title="OBN shot");

# simulate towed streamer acquisition with moving receiver
ns  = length(src)
irx = Vector{Vector{Int64}}(undef, ns)
irz = Vector{Vector{Int64}}(undef, ns)
for i = 1 : ns
    irx[i]= collect(isx[i]-102: isx[i]-2);
    irz[i]= 2 * ones(Int64, length(irx[i]));
end

# do parallel simulation
dir_obs = joinpath(homedir(), "Desktop/shotgather_tow")
get_shotgather(dir_obs, irz, irx, src, params);

path_shot = joinpath(dir_obs, "recordings_7.bin");
rec       = read_recordings(path_shot)
SeisPlotTX(rec.p, dy=0.001, title="towed streamer");

# end of documentation


text:
@testset "Filtering" begin
  @testset "UniqueCoordsFilter" begin
    X = rand(3,100)
    z = rand(100)
    n = [string(i) for i in 1:100]
    Xd = hcat(X, X[:,1:10])
    zd = vcat(z, z[1:10])
    nd = vcat(n, n[1:10])
    sdata = PointSetData(OrderedDict([:z=>zd, :n=>nd]), Xd)
    ndata = filter(sdata, UniqueCoordsFilter())
    @test npoints(ndata) == 100
  end
end


text:
using Dates
using PyPlot
using Colors

# Little tools

function makehomesavepath(foldername)
    rightnow = Dates.format(Dates.now(), "yyyy-mm-dd_HH-MM-SS")
    savepath = joinpath(homedir(), foldername * "_" * rightnow)

    while isdir(savepath)
        sleep(1.)
        rightnow = Dates.format(Dates.now(), "yyyy-mm-dd_HH-MM-SS")
        savepath = joinpath(homedir(), foldername * "_" * rightnow)
    end
    savepath
end

# Return a tuple with the values of x at the positions in i and the values of param elsewhere
# i: Array{Int64,1} index of parameter to insert the value x
# x: Float64
# param: tuple with learner's parameters (learnerspecs.param)
function parameterset(param, x, i)
    # tuple([j == i ? x : param[j] for j in 1:length(param)]...)
    tuple([any(in.(j, i)) ? x[findall(in.(j,i))...] : param[j] for j in 1:length(param)]...)
end

searchdirendings(path, key) = filter(x->endswith(x, key), readdir(path))

searchoccurence(a, key) = filter(x->occursin(key, x), a)

cmap(i) = get_cmap("jet")(i)

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# --------- Not used eventually for the paper, only for inspection -------------
function clisthardcodedG()
    clist = [(0., 0., 0.), # bayes
            (1., 0.4, 0.7), # gnassarJN
            (0.54, 0.17, 0.88), # gnassarNatN # (1., 0.08, 0.57),
            (0., 0., 1.), # gsor1
            (0.37, 0.59, 0.63), # gsor10
            (0., 0.75, 1.), # gsor20
            (0.4, 0.4, 0.4), # leaky
            (0.68, 1.0, 0.18), # pf1
            (0.19, 0.8, 0.19), # pf10
            (0., 1., 0.), # pf20
            (1., 0.5, 0.), # Smile - orange
            (0.6, 0., 0.) # Var Smile
    ]
end
function clisthardcodedM()
    clist = [(0., 0., 0.), # bayes
            (0.4, 0.4, 0.4), # leaky
            (0.68, 1.0, 0.18), # pf1
            (0.19, 0.8, 0.19), # pf10
            (0., 1., 0.), # pf20
            (1., 0.5, 0.), # Smile - orange
            (0.6, 0., 0.), # Var Smile
            (0., 0., 1.), # tsor1
            (0.37, 0.59, 0.63), # tsor10
            (0., 0.75, 1.) # tsor20
    ]
end
function clisthardcoded_robM()
    clist = [(0.4, 0.4, 0.4), # leaky
            (0.19, 0.8, 0.19), # pf10
            (0.68, 1.0, 0.18), # pf1
            (0., 1., 0.), # pf20
            (1., 0.5, 0.), # Smile - orange
            (0.6, 0., 0.), # Var Smile
            (0., 0., 0.), # bayes
            (0.37, 0.59, 0.63), # tsor10
            (0., 0., 1.), # tsor1
            (0., 0.75, 1.), # tsor20
        ]
end
function clisthardcoded_robG()
    clist = [(0., 0., 0.), # bayes
            (1., 0.4, 0.7), # gnassarJN
            (0.54, 0.17, 0.88), # gnassarNatN # (1., 0.08, 0.57),
            (0.37, 0.59, 0.63), # gsor10
            (0., 0., 1.), # gsor1
            (0., 0.75, 1.), # gsor20
            (0.4, 0.4, 0.4), # leaky
            (0.19, 0.8, 0.19), # pf10
            (0.68, 1.0, 0.18), # pf1
            (0., 1., 0.), # pf20
            (1., 0.5, 0.), # Smile - orange
            (0.6, 0., 0.) # Var Smile
        ]
end


text:
function swap(array, i, j)
  tmp = array[i]
  array[i] = array[j]
  array[j] = tmp
end

function partition(array, left, right)
  pivotIndex = rand(left:right)
  pivot = array[pivotIndex]
  swap(array, pivotIndex, right)

  pos = left
  for i = left:right-1
    if array[i] < pivot
      swap(array, i, pos)
      pos = pos+1
    end
  end

  swap(array, pos, right)

  return pos
end

function quick_sort(array, left, right)
  if left < right
    pivotIndex = partition(array, left, right)
    quick_sort(array, left, pivotIndex-1)
    quick_sort(array, pivotIndex+1, right)
  end
end

function is_sorted(array)
  max = 0
  i = 1
  while i < length(array)
    if array[i] < max
      return false
    end
    max = array[i]
    i = i+1
  end
  return true
end

function generate_random_array(length)
  rand(0:length, length)
end

for i = 1:100
  array = generate_random_array(1000)
  quick_sort(array, 1, length(array))
  if !is_sorted(array)
    println("Error! Array is not sorted!")
    exit()
  end
end

println("Quick sort works fine!")

max_stars_repo_name: avi-pal/al-go-rithms


text:
module OPdist
#
# Module OPdist contains functions that create and manipulate probability
# distributions.
#

export binit, moments, tdist

using SpecialFunctions
using LinearAlgebra

using ..OPccmt
using ..OPpoly

# ===================================================================
function binit(dat,wgt,Nb,mflag)
#
# Bins it.  Creates a distribution from input data.
#
# Version:        Changes:
# --------        -------------
# 21.05.2020      Original code.
#
# Version:        Verification:
# --------        -------------
# 21.05.2020
#
# Inputs:
# -------
# dat             : Raw data, input as an N-by-M array, M coordinates of
#                   N dimensions.
# wgt             : Weights to attach to each point; in a typical
#                   application this will be 1, and this can be input as
#                   a single scalar.
# Nb              : Number of bins along each dimension.
# mflag           : Method flag.  1=nearest bin, 2=linear weight.
#
# Outputs:
# --------
# dis             : The binned distribution N-by-prod(Nb).
# xb              : N-by-3 array containing xmin, dx, xmax for each dimension.

   M = size(dat,2)
   N = size(dat,1)

   prodNb = prod(Nb)

   # Get the bins.
   xb = zeros(N,3)
   for idim = 1:N
      mx = maximum(dat[idim,:])
      mn = minimum(dat[idim,:])
      mxmn = mx - mn
      rng = 1.1*mxmn
      dx = rng/(Nb[idim]-1)
      xb[idim,3] = mx + 0.05*mxmn
      xb[idim,1] = mn - 0.05*mxmn
      xb[idim,2] = dx
   end

#print(dat)
#print("\n")
#print(xb)
#print("\n")

   # Bin the data.
   ic, wc = nearestCells(dat,xb)

#print(ic)
#print("\n")
#print(wc)

   dis = zeros(prodNb)
   if (size(wgt,1) == size(dat,2))
      if (mflag == 1)  # Nearest neighbor.
         for ip = 1:M
            jnk, jc = findmax(wc[:,ip])
            dis[ic[jc,ip]] = dis[ic[jc,ip]] + wgt[ip]
         end
      elseif (mflag == 2)  # Linear weight.
         for ip = 1:M
            dis[ic[:,ip]] = dis[ic[:,ip]] + wc[:,ip]*wgt[ip]
         end
      end
   elseif (size(wgt,1) == 1)
      if (mflag == 1)  # Nearest neighbor.
         for ip = 1:M
            jnk, jc = findmax(wc[:,ip])
#print([ip jnk jc ic[jc] dis[ic[jc]]])
#print("\n")
            dis[ic[jc,ip]] = dis[ic[jc,ip]] + wgt
         end
      elseif (mflag == 2)  # Linear weight.
         for ip = 1:M
            dis[ic[:,ip]] = dis[ic[:,ip]] + wc[:,ip]*wgt
         end
      end
   end

   return dis, xb

end  # binit

# ===================================================================
function moments(xc,phi,n)
#
# Compute the nth degree moment of the binned input distribution.
#
# Version:        Changes:
# --------        -------------
# 08.01.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 08.01.2021
#
# Inputs:
# -------
# xc              : Ndim-by-Np array of bin coordinates.
# phi             : Probability within each bin.  This does not
#                   need to be normalized, it may be raw counts.
# n               : Order of the moment to be computed.
#
# Outputs:
# --------
# m               : The order-n moments.  Vector of length Ndim^n.

   Ndim = size(xc,1)
   Np = size(xc,2)

   tot = sum(phi)

   # Calculate the mean.
   mu = zeros(Ndim)
   for ip = 1:Np
      mu = mu + xc[:,ip] .* phi[ip]
   end
   mu = mu / tot

#print(mu,"\n")

   if (n == 0)
      m = 1.0
   elseif (n == 1)
      m = zeros(Ndim)
   else
      # Calculate the desired moments.
      xb = zeros(n,3)
      xb[:,1] .= 1
      xb[:,2] .= 1
      xb[:,3] .= Ndim
      Ncell = Ndim^n
      ic = range(1,length=Ncell)
      ks = xFromCell(ic,xb)
      ks = convert.(Int,round.(ks))

#print(ks,"\n")

      m = zeros(Ncell)
      for ip = 1:Np
         for jc = 1:Ncell
            mul = 1.0
#print("-----------------------\n")
            for jord = 1:n
               mul = mul*(xc[ks[jord,jc],ip] - mu[ks[jord,jc]])
#print(ip," ",jc," ",jord," ",xc[ks[jord,jc],ip]," ",mu[ks[jord,jc]]," ",mul,"\n")
            end
            m[jc] = m[jc] + mul*phi[ip]
         end
      end
      m = m / tot
   end

   return m

end  # moments

# ===================================================================
function tdist(x,mu,Sig,nu)
#
# Compute a t distribution at the specified points.
#
# Version:        Changes:
# --------        -------------
# 18.01.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 18.01.2021
#
# Inputs:
# -------
# x               : The points at which to compute the t distribution.
#                   Ndim rows by Np columns.
# mu              : The mean vector.
# Sig             : The matrix Sigma, noting that for nu > 2 the
#                   covariance is (nu/(nu-2)) Sigma.
# nu              : The t distribution DOF parameter.  Approaches
#                   Gaussian as nu -> Inf.  Heavier tails for small nu.
#
# Outputs:
# --------
# t               : The t distribution values.

   d = size(x,1)
   Np = size(x,2)
   detS = det(Sig)
   invS = inv(Sig)
   xc = x .- mu

   gg = loggamma(0.5*(nu + d)) - loggamma(0.5*nu)
   c1 = (((pi*nu)^d)*detS)^(-0.5)
   c2 = -0.5*(nu + d)

   t = zeros(typeof(x[1]),Np)
   for ip in 1:Np
      t[ip] = exp(gg)*c1*
              ((1.0 + (transpose(xc[:,ip])*invS*xc[:,ip])/nu)^c2)
   end

   return t

end  # tdist

end  # module

max_stars_repo_name: SINTEF-Energy-Wind/OutliersPrimer

max_issues_repo_name: SINTEF-Energy-Wind/OutliersPrimer

max_forks_repo_name: SINTEF-Energy-Wind/OutliersPrimer


text:
struct NotOffloadableError
    ir
    sv
    reason
end

function Base.showerror(io::IO, e::NotOffloadableError)
    println(io, "The specified function is not offloadable. The offending IR was:")
    Base.IRShow.show_ir(io, e.ir; verbose_linetable=false)
    if isa(e.reason, String)
        println(io, e.reason)
    elseif isa(e.reason, Function)
        e.reason(io, e.ir, e.sv)
    end
end

function synthesize_bt_for_line(ir, line)
    entry = ir.lines[line]
    frames = Base.StackFrame[]
    while entry !== 0
        lin = ir.linetable[entry]
        push!(frames, Base.StackFrame(lin.method,  lin.file, lin.line, nothing, false, lin.inlined_at != 0, C_NULL))
        entry = lin.inlined_at
    end
    frames
end

function analyze_unreachable(io::IO, ir, sv)
    print(io, "The called function was guaranteed to error. Analyzing...\n")
    rv = ir.stmts[end]
    @assert isa(rv, ReturnNode)
    # Generally the call before the ReturnNode will be the offending one
    before = ir.stmts[end - 1]
    if ir.types[end - 1] !== Bottom || !isexpr(before, :call)
        println(io, "This IR looks odd. Prior statement was not a call that was guaranteed to error")
        return
    end
    # See if this is a method error
    f = argextype(before.args[1], ir, sv.sp)
    isa(f, Const) || return println(io, "Expected function to be a constant")
    f = f.val
    args = Tuple{[widenconst(argextype(before.args[i], ir, sv.sp)) for i = 2:length(before.args)]...}
    args = Base.to_tuple_type(args)
    tt = Base.signature_type(f, args)
    m = ccall(:jl_gf_invoke_lookup, Any, (Any, UInt), tt, sv.params.world)
    if m === nothing
        synthetic_bt = synthesize_bt_for_line(ir, length(ir.stmts)-1)
        showerror(io, MethodError(f, args, sv.params.world), synthetic_bt)
    end
end

function typeinf_code_slottypes(mi::Core.Compiler.MethodInstance, run_optimizer::Bool, params::Core.Compiler.Params)
    ccall(:jl_typeinf_begin, Cvoid, ())
    result = Core.Compiler.InferenceResult(mi)
    frame = Core.Compiler.InferenceState(result, false, params)
    frame === nothing && return (nothing, Any)
    Core.Compiler.typeinf(frame)
    ccall(:jl_typeinf_end, Cvoid, ())
    frame.inferred || return (nothing, Any)
    return (frame.src, result.result, frame.slottypes)
end

function explain_suboptimal_inference(sig, bt=Base.StackFrame[])
    params = Core.Compiler.CustomParams(typemax(UInt); aggressive_constant_propagation=true)
    methds = Base._methods_by_ftype(sig, 1, params.world)
    x = methds[1]
    meth = Base.func_for_method_checked(x[3], Tuple{sig.parameters[2:end]...})
    mi = Compiler.code_for_method(meth, sig, x[2], params.world)
    (ci, ty, slottypes) = typeinf_code_slottypes(mi, false, params)
    ir = Compiler.inflate_ir(ci, Compiler.spvals_from_meth_instance(mi),
                             Compiler.matching_cache_argtypes(mi, nothing)[1])
    empty!(ir.argtypes); append!(ir.argtypes, slottypes)
    sv = Compiler.OptimizationState(mi, params)
    for (idx, (stmt, type)) in enumerate(zip(ir.stmts, ir.types))
        if type === Any || isa(type, UnionAll)
            if isexpr(stmt, :(=))
                stmt = stmt.args[2]
            end
            if isexpr(stmt, :call)
                println("Found uninferred call to $(stmt.args[1]), recursing")
                append!(bt, synthesize_bt_for_line(ir, idx))
                atys = Any[]
                for i = 1:length(stmt.args)
                    ty = widenconst(argextype(stmt.args[i], ir, sv.sp, slottypes))
                    if ty === Any
                        if isa(stmt.args[i], GlobalRef)
                            printstyled("Global variable ($(stmt.args[i])) was referenced.", color=:bold)
                            Base.show_backtrace(stdout, reverse(bt))
                            return
                        end
                        printstyled("Call argument was `Any` (analysis was unable to determine a precise reason. This should be improved).", color=:bold)
                    end
                    push!(atys, ty)
                end
                if atys[1] == typeof(Core.getfield)
                    if atys[2] == Core.Box
                        printstyled("A closure `Box` was introduced. See Julia issue #15276.", color=:bold)
                        Base.show_backtrace(stdout, reverse(bt))
                        return
                    end
                    printstyled("Unknown destructuring.", color=:bold)
                    Base.show_backtrace(stdout, reverse(bt))
                    return
                end
                nsig = Tuple{atys...}
                return explain_suboptimal_inference(nsig, bt)
            end
        end
    end
end


text:
for N in [Float64, Rational{Int}, Float32]
    B1 = BallInf(N[0, 0], N(1))
    B2 = BallInf(N[4, -4], N(1))
    B3 = BallInf(N[1, -1], N(1))
    C = Complement(B1)

    # double-complement is the identity
    @test Complement(C) == B1

    # dimension
    @test dim(C) == 2

    # membership
    @test N[1, 1] ∉ C && N[2, 2] ∈ C

    # emptiness
    @test !isempty(C) && isempty(Complement(Universe(2)))

    # inclusion
    subset, point = ⊆(B2, C, true)
    @test B2 ⊆ C && subset && point == N[]
    for X in [B1, B3]
        subset, point = ⊆(X, C, true)
        @test X ⊈ C && !subset && point ∈ X && point ∉ C
    end

    # isdisjoint
    for X in [B2, B3]
        res, w = isdisjoint(X, C, true)
        @test !isdisjoint(X, C) && !res && w ∈ X && w ∈ C
        res, w = isdisjoint(C, X, true)
        @test !isdisjoint(C, X) && !res && w ∈ X && w ∈ C
    end
    res, w = isdisjoint(B1, C, true)
    @test isdisjoint(B1, C) && res && w == N[]
    res, w = isdisjoint(C, B1, true)
    @test isdisjoint(C, B1) && res && w == N[]

    # test convexity from the type
    @test isconvextype(typeof(Complement(Universe{N}(2))))
    @test isconvextype(typeof(Complement(EmptySet{N}(2))))
    @test isconvextype(typeof(Complement(HalfSpace(N[1], N(0)))))
end


text:
#function testNHHill_L()

using CGDycore

# Physical parameters
Param=CGDycore.PhysParameters();

nx=60;
ny=2;
lx=6000000;
ly=200000;
x0=-3000000;
y0=0;

Boundary = (;WE="Period", BT="Period")
Param.Grid=CGDycore.CartGrid(nx,ny,lx,ly,x0,y0,CGDycore.OrientFaceCart,Boundary,Param);
Param.TopoS="AgnesiCartX";

Param.H=15600;
nz=40;
zP=zeros(nz,1);
z=zeros(nz+1,1);
dz=Param.H/nz;
zP[1]=dz/2;
for i=2:nz
  zP[i]=zP[i-1]+dz;
end
for i=2:nz+1
  z[i]=z[i-1]+dz;
end

Param.Grid.nz=nz;
Param.Grid.z=z;
Param.Grid.zP=zP;
Param.Grid.dz=dz;

# Model
Param.ModelType="Curl";
Param.Coriolis=false;
Param.Thermo="";
Param.Source=false;
Param.StrideDamp=6000;
Param.Relax=1.e-1;
Param.Damping=true;
Param.Flat=false;
Param.Coriolis=false;
Param.Buoyancy=true;
Param.hC=400;
Param.x0C=0;
Param.aC=100000;
Param.Th0=300;
Param.uMax=10;
Param.vMax=0;
Param.NBr=1.e-2;
Param.Equation="Compressible";
Param.RefProfile=false

# Discretization
OrdPoly=4;
OrdPolyZ=1;
(CG,Param)=CGDycore.Discretization(OrdPoly,OrdPolyZ,CGDycore.JacobiDG3,Param);
Param.HyperVisc=true;
Param.HyperDCurl=1.e4;
Param.HyperDGrad=1.e4;
Param.HyperDDiv=1.e4;
Param.Upwind=false

# Initial conditions
Param.NumV=5;
U=zeros(CG.NumG,nz,Param.NumV);
Param.ProfRho="GravityHill";
Param.ProfTheta="GravityHill";
Param.ProfVel="Const";
Param.RhoPos=1;
Param.uPos=2;
Param.vPos=3;
Param.wPos=4;
Param.ThPos=5;
U[:,:,Param.RhoPos]=CGDycore.Project(CGDycore.fRho,CG,Param);
(U[:,:,Param.uPos],U[:,:,Param.vPos])=CGDycore.ProjectVec(CGDycore.fVel,CG,Param);
U[:,:,Param.ThPos]=CGDycore.Project(CGDycore.fTheta,CG,Param).*U[:,:,Param.RhoPos];

# Output
Param.vtkFileName="Hill_L";
Param.vtk=0;
vtkGrid=CGDycore.vtkCGGrid(CG,CGDycore.TransCart,CGDycore.Topo,Param);
Param.cNames = [
  "Rho",
  "u",
  "v",
  "w",
  "Th"
]

IntMethod="Rosenbrock";
#IntMethod="RungeKutta";
if strcmp(IntMethod,"Rosenbrock")
  dtau=30;
else
  dtau=.4;
end
nIter=6000;
Param.RK=CGDycore.RungeKuttaMethod("RK4");
Param.ROS=CGDycore.RosenbrockMethod("SSP-Knoth");
CFL=0.125;
time=0;
Param.EndTime=216000;
nIter=Param.EndTime/dtau;
PrintTime=10000;
PrintInt=floor(PrintTime/dtau);
Param.vtk=CGDycore.vtkOutput(U,vtkGrid,CG,Param);
nIter=20

# Print initial conditions
Param.vtk=CGDycore.vtkOutput(U,vtkGrid,CG,Param);
#

str = IntMethod
if str == "Rosenbrock"
    @time begin
      for i=1:nIter
        @info "Iteration: $i"
        @show dtau
        U .= CGDycore.RosenbrockSchur(U,dtau,CGDycore.FcnNHCurlVec,CGDycore.JacSchur,CG,Param);
        if mod(i,PrintInt)==0
          Param.vtk=CGDycore.vtkOutput(U,vtkGrid,CG,Param);
        end
      end
    end

elseif str == "RungeKutta"
    for i=1:nIter
      @info "Iteration: $i"

      U .= CGDycore.RungeKuttaExplicit(U,dtau,CGDycore.FcnNHCurlVec,CG,Param);

      time[1] += dtau;
      if mod(i,PrintInt)==0
        Param.vtk=CGDycore.vtkOutput(U,vtkGrid,CG,Param);
      end
    end
else
  error("Bad str")
end


text:
#!/usr/bin/env julia

import GR

const theta = 70.0   # initial angle
const gamma = 0.1    # damping coefficient
const L = 0.2        # wire length

w, h, ball = GR.readimage("ball.png")

function rk4(x, h, y, f)
  k1 = h * f(x, y)
  k2 = h * f(x + 0.5 * h, y + 0.5 * k1)
  k3 = h * f(x + 0.5 * h, y + 0.5 * k2)
  k4 = h * f(x + h, y + k3)
  return x + h, y + (k1 + 2 * (k2 + k3) + k4) / 6.0
end

function deriv(t, state)
  theta, omega = state
  return [omega, -gamma * omega - 9.81 / L * sin(theta)]
end

function draw_cradle(theta)
  GR.clearws()
  GR.setviewport(0, 1, 0, 1)
  GR.setcolorrep(1, 0.7, 0.7, 0.7)
  # draw pivot point
  GR.fillarea([0.3, 0.7, 0.7, 0.3], [0.79, 0.79, 0.81, 0.81])
  # draw balls
  for i = -2:2
    x = [0.5, 0.5] .+ i * 0.06
    y = [0.8, 0.4]
    if (theta < 0 && i == -2) || (theta > 0 && i == 2)
      x[2] += sin(theta) * 0.4
      y[2] = 0.8 - cos(theta) * 0.4
    end
    GR.polyline(x, y)               # draw wire
    GR.drawimage(x[2]-0.03, x[2]+0.03, y[2]-0.03, y[2]+0.03, 50, 50, ball)
  end
  GR.updatews()
end

function main()
  t = 0.0
  dt = 0.01
  state = [theta * pi / 180, 0]

  start = refresh = time_ns()

  while t < 30
    t, state = rk4(t, dt, state, deriv)
    theta, omega = state

    if time_ns() - refresh > 20 * 1000000   # 20ms
      draw_cradle(theta)
      refresh = time_ns()
    end

    now = (time_ns() - start) / 1000000000
    if t > now
        sleep(t - now)
    end
  end
end

main()


text:
module PendulumWithDamper

using ModiaLang

# ModiaLang models
include("$(ModiaLang.path)/models/Blocks.jl")
include("$(ModiaLang.path)/models/Electric.jl")
include("$(ModiaLang.path)/models/Rotational.jl")

import Modia3D
using  Modia3D.ModiaInterface

Pendulum = Model(
    m = 1.0,
    g = 9.81,
    vmat1 = VisualMaterial(color="Sienna", transparency=0.5),
    vmat2 = VisualMaterial(color="Red"),
    world = Object3D(feature=Scene()),
    worldFrame = Object3D(parent=:world,
                          feature=Visual(shape=CoordinateSystem(length=0.5))),

    frame0 = Object3D(feature=Solid(shape=Beam(axis=1, length=1.0, width=0.2, thickness=0.2),
                                    massProperties=MassProperties(mass=:m),
                                    visualMaterial=:(vmat1))),
    frame1 = Object3D(parent=:(frame0),
                      translation=[-0.5, 0.0, 0.0]),
    cyl    = Object3D(parent=:(frame1),
                      feature=Visual(shape=Cylinder(axis=3, diameter=0.2/2, length=1.2*0.2),
                                     visualMaterial=:(vmat2))),
    rev    = RevoluteWithFlange(obj1=:world, obj2=:frame1)
)

PendulumWithDamp = Model(
    pendulum = buildModia3D(Pendulum | Map(m=2.0, rev=Map(phi=Var(init=1.0)))),

    damper = Damper | Map(d=0.5),
    fixed = Fixed,
    connect = :[(damper.flange_b, pendulum.rev.flange),
                (damper.flange_a, fixed.flange)]
)

#@showModel PendulumWithDamp

pendulumWithDamper = @instantiateModel(PendulumWithDamp, unitless=true, log=false, logStateSelection=false, logCode=true)

stopTime = 10.0
requiredFinalStates = [-1.578178283450938, 0.061515170100766486]
simulate!(pendulumWithDamper, stopTime=stopTime, log=true, logStates=true, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
plot(pendulumWithDamper, ["pendulum.rev.flange.phi", "pendulum.rev.variables[1]"], figure=1)

end


text:
# Practical Machine learning
# Clustering based learning - L-mena clustering
# Chapter 8

using Images
using ImageView
using Color

# Run non-interactively, accept either an image or a directory
# If directory, iterate through image filetypes

# Enter filepath and number of dominant colors wanted, k
function dominant_colors (filename, k)
	if isfile(filename)
		img = imread(filename,RGB)
		init(img, filename, k)
	elseif isdir(filename)
		files = filter!(abspath(r"\.(?:jpe?g|gif|png|tiff)$")), readdir()) #creates an array of filenames by filtering out only files that in file extensions.
		for i in files
			dominant_colors(i)
		end
	else
		error("No image found.")
	end
end

type Point
	coords::Array # Color associated with pt, a 3D array.
	ct:Int # Count
end

type Cluster
	points::Array # Points associated with cluster, C_k
	centroid::Point # Center of cluster, assumed mean of pt values
	k::Int # Cluster count
end

function init(img, filename, k)
# Convert color space from sRGB (linear) to CIEXYZ to CIELAB
	run(`convert $filename -thumbnail 200x200 $filename`) #convert to thumbnail via ImageMagick CLI
	img = convert(Image{LAB}, img) #use Color to convert to LAB automagically
	points = getpoints(img)
	randclusters(points, k)
	kmeans(points, k)

function getpoints(img)
	points = []
	count = 0
	for count, color in img[1:width(img)]
		for count, color in img[2:height(img)]
			count += 1
			points.append(Point(color, count))
		end
	return points
end

# Sq. euclidean distance
function distance (pt1, pt2)
	return mapreduce((pt1.coords[i]-pt2.coords[i])**2, +, 1:length(pt1.coords))
end

# Randomly assign pixels to represent intial centroid/clusters
function randclusters(pts::Array, k)
	kclusters = []
	for n = 1:k
		kclusters.append(Cluster(pts, pts[rand(1:end)], n)
	return kclusters
end

# Recalculate Centroid
function recenter(points, )

# K-Means Algorithm (alternate between resigning points to a cluster based on similarity and cluster centroid based on the points assigned)
function kmeans()
# Repeat n number of times

# Optionally convert returned clusters back to sRGB


text:
using Colors, ColorSchemes
using DataFrames
using StatsPlots
using Statistics
using GLM
using Measures
using LsqFit

pyplot(grid=false)
fnt = Plots.font("Arial")
default(titlefont=fnt, guidefont=fnt, tickfont=fnt, legendfont=fnt)
#  Plots.resetfontsizes()

noise_level = 0.003544

###############################################################################
### FIG S3A
###############################################################################

function get_pct_off(t_recruit, t_release, k)
    ks = k['s']
    ka = k['a']
    ki = k['i']
    y1 = (ks+ka+ki - sqrt((ks+ka+ki)^2 - 4*ks*ki))/2
    y2 = (ks+ka+ki + sqrt((ks+ka+ki)^2 - 4*ks*ki))/2

    function c_r(t_recruit)
        return ks/(y2-y1) * (exp(-y1*t_recruit) - exp(-y2*t_recruit))
    end

    function c_i(t_recruit)
        if (ki==0)
            return 0
        else
            term1 = (ks*ki)/(y1*(y2-y1)) * exp(-y1*t_recruit)
            term2 = (ks*ki)/(y2*(y2-y1)) * exp(-y2*t_recruit)
            return 1 - term1 + term2
        end
    end

    return c_r(t_recruit) * exp(-ka * t_release) + c_i(t_recruit)
end

# see https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1023694
function get_single_cell_cc(t_recruit, t_release, k)
    p = 1 - get_pct_off(t_recruit, t_release, k) # error probability
    return log2(1 + (1-p)^(p/(1-p)))
end

function get_population_cc(t_recruit, t_release, k, variance=0.05)
    power = get_pct_off(t_recruit, t_release, k)
    return maximum([0, 0.5 * log2(1 + power/variance)])
end

###############################################################################
### FIG 5D
###############################################################################

t_recruits = collect(0:0.001:1)
t_release = 0
k_p_ccs = [get_population_cc(t_recruit, t_release, k_krab, noise_level) for t_recruit in t_recruits]
h_p_ccs = [get_population_cc(t_recruit, t_release, k_hdac, noise_level) for t_recruit in t_recruits]

p1 = plot(t_recruits, k_p_ccs,
         label="KRAB",
         color=color_map["KRAB"],
         linewidth=2)
p1 = plot!(t_recruits, h_p_ccs,
          label="HDAC4",
          color=color_map["HDAC4"],
          linewidth=2,
          xlabel="Days Recruited",
          ylabel="Channel Capacity (bits)", yticks=[0, 1, 2, 3, 4],
          xticks = [0, 0.5, 1],
          legend=:right,
          legendfont = font(10),
          size=(600, 300))

savefig(p1, "fig_s3a.pdf")

###############################################################################
### FIG 5E
###############################################################################

t_recruits = 5
t_release = 1
sigmas = collect(0.00001:0.00001:0.25)
k_p_ccs = [get_population_cc(t_recruit, t_release, k_krab, sigma) for sigma in sigmas]
h_p_ccs = [get_population_cc(t_recruit, t_release, k_hdac, sigma) for sigma in sigmas]

p1 = plot(sigmas, k_p_ccs,
         label="KRAB",
         color=color_map["KRAB"],
         linewidth=2)
p1 = plot!(sigmas, h_p_ccs,
          label="HDAC4",
          color=color_map["HDAC4"],
          linewidth=2,
          xlabel="Output Noise Variance",
          # xscale=:log10,
          # legend=:outerright,
          ylabel="Channel Capacity (bits)",
          ylims=(0,8.5),
          legendfont = font(10),
          size=(600, 300))

println("KRAB: $(get_population_cc(t_recruit, t_release, k_krab, noise_level))")
println("HDAC4: $(get_population_cc(t_recruit, t_release, k_hdac, noise_level))")

p1 = vline!(p1, [noise_level], color=:black, linestyle=:dash, label=false)
savefig(p1, "fig_s3b.pdf")

###############################################################################
### FIG 5F
###############################################################################

function get_population_cc_input_noise(t_recruit, t_release, k, input_variance = 0.05, output_variance=0.05)
    target_poff = get_pct_off(t_recruit, t_release, k)

    # input b/w zero and 1 (power=1), and noise is input variance
    input_variance_scaled = target_poff * input_variance

    overall_cc = 0.5*log2(1 + target_poff/(input_variance_scaled + output_variance))

    return overall_cc
end

k_hdac = Dict('a' => 0.76,
              's' => 8.6,
              'i' => 0)

k_krab = Dict('a' => 0.31,
              's' => 11,
              'i' => 0.13)

t_recruit = 5
t_release = 5

sigmas = collect(0.00001:0.00001:0.25)

input_noise_sigmas = [0, 0.001, 0.01, 0.1, 1]
noise_colors = reverse([colorant"purple4", colorant"purple1", colorant"mediumorchid2", colorant"violet", colorant"plum1", colorant"thistle1"])
noise_colors = reverse(noise_colors)
noise_colors = get(ColorSchemes.viridis, reverse(range(0.0, stop=0.9, length=length(noise_colors))))

krab_cc_list = [[get_population_cc_input_noise(1, t_release, k_krab, input_noise_sigmas[i], sigma) for sigma in sigmas] for i in 1:length(input_noise_sigmas)]
hdac_cc_list = [[get_population_cc_input_noise(1, t_release, k_hdac, input_noise_sigmas[i], sigma) for sigma in sigmas] for i in 1:length(input_noise_sigmas)]

# t1_p_ccs = [get_population_cc_input_noise(1, t_release, k_krab, sigma) for sigma in sigmas]
# t2_p_ccs = [get_population_cc_input_noise(2, t_release, k_krab, sigma) for sigma in sigmas]
# t3_p_ccs = [get_population_cc_input_noise(3, t_release, k_krab, sigma) for sigma in sigmas]
# t4_p_ccs = [get_population_cc_input_noise(4, t_release, k_krab, sigma) for sigma in sigmas]
# t5_p_ccs = [get_population_cc_input_noise(5, t_release, k_krab, sigma) for sigma in sigmas]
#
# cc_list = [t1_p_ccs, t2_p_ccs]
p1 = plot(xlabel="Output Noise Variance",
          ylabel="Channel Capacity (bits)",
          title="KRAB",
          legendtitle="Input Noise Variance",
          legend=false,
          legendfont = font(10),
          xscale=:log10,
          ylims=(-0.25,8.5),
          size=(600, 300))
for i in 1:length(krab_cc_list)
    p1 = plot!(sigmas, krab_cc_list[i], label="$(input_noise_sigmas[i])",
               color=noise_colors[i],
               linewidth=2)
end

p2 = plot(xlabel="Output Noise Variance",
          ylabel="Channel Capacity (bits)",
          title="HDAC4",
          legendtitle="Input Noise Variance",
          legendfont = font(10),
          xscale=:log10,
          legend=:outerright,
          ylims=(-0.25,8.5),
          size=(1200, 300))
for i in 1:length(hdac_cc_list)
    p2 = plot!(sigmas, hdac_cc_list[i], label="$(input_noise_sigmas[i])",
               color=noise_colors[i],
               linewidth=2)
end

p1 = vline!(p1, [noise_level], color=:black, linestyle=:dash, label=false)
p2 = vline!(p2, [noise_level], color=:black, linestyle=:dash, label=false)

lay = @layout[a  b]
p = plot(p1, p2, layout=lay)
savefig(p, "fig_s3c.pdf")


text:
# # [Custom Title 2](@id custom_id_2)

# This is the content


text:
using Documenter, MinimaxCalibratedEBayes
const MCEB = MinimaxCalibratedEBayes

makedocs(;
    modules=[MinimaxCalibratedEBayes],
    format=  Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages=[
        "Home" => "index.md",
        "Manual" => Any[
            "man/targets.md",
            "man/prior_classes.md",
            "man/samples.md",
            "man/pilot.md",
            "man/quasiminimax_linear.md",
            "man/inference.md"],
         "Tutorials" => "tutorials.md",
    ],
    repo="https://github.com/nignatiadis/MinimaxCalibratedEBayes.jl/blob/{commit}{path}#L{line}",
    sitename="MinimaxCalibratedEBayes.jl",
    authors="Nikos Ignatiadis <nikos.ignatiadis01@gmail.com>",
    assets=String[],
)

if get(ENV, "CI", nothing) == "true"
    let start_dir = pwd()
        cd(@__DIR__)
        mkpath("build/tutorials")
        tutorial_names = ["linear_estimation.html"; "data_analysis.html"]
        in_path_names = joinpath.("tutorials", tutorial_names)
        out_path_names = joinpath.("build/tutorials",tutorial_names)
        cp.(in_path_names, out_path_names, force=true)
        cd(start_dir)
    end
end

deploydocs(;
    repo="github.com/nignatiadis/MinimaxCalibratedEBayes.jl",
)


text:
using LabelledArrays

export verdes

"""
    eom_verdes(u, p, t) → SVector{3}

Equations of motion for a 3D system where the response X is a highly
nonlinear combination of Y and Z [1]. The forcings Y and Z involve
sines and cosines, respectively, and have different periods.

The implementation here allows for tuning the period of the
forcing signals.

# References
Verdes, P. F. "Assessing causality from multivariate time series." Physical Review E 72.2 (2005): 026222.
"""
function eom_verdes(u, p, t)
    x, y, z = (u...,)
    ωy, ωz, σx, σy, σz = (p...,)

    ηx = σx == 0 ? 0 : rand(Normal(0, σx))
    ηy = σy == 0 ? 0 : rand(Normal(0, σy))
    ηz = σz == 0 ? 0 : rand(Normal(0, σz))

    dx = y*(18y - 27y^2 + 10)/2 + z*(1-z) + ηx
    dy = (1 - cos((2*pi/ωy) * t))/2 + ηy
    dz = (1 - sin((2*pi/ωz) * t))/2 + ηz
    return SVector{3}(dx, dy, dz)
end

function verdes(u₀, ωy, ωz, σx, σy, σz)
    p = LVector(ωy = ωy, ωz = ωz, σx = σx, σy = σy, σz = σz)

    DiscreteDynamicalSystem(eom_verdes, u₀, p)
end

"""
    verdes(;u₀ = rand(3), ωy = 315, ωz = 80,
        σx = 0.0, σy = 0.0, σz = 0.0) → DiscreteDynamicalSystem

Intitialise a 3D system where the response X is a highly nonlinear combination
of Y and Z. The forcings Y and Z involve sines and cosines, and
have different periods, which controlled by `ωy` and `ωz`.

The equations of motion are

```math
\\begin{aligned}
x(t+1) &= \\dfrac{y(t)(18y(t) - 27y(t)^2 + 10)}{2} + z(t)(1-z(t)) + ηx \\
y(t+1) &= \\dfrac{(1 - \\dfrac{\\cos(2\\pi)}{\\omega y}t)}{2} + ηy \\
z(t+1) &= \\dfrac{(1 - \\dfrac{\\sin(2\\pi)}{\\omega z}t)}{2} + ηz
\\end{aligned}
```
where ηx, ηy, ηz is gaussian noise with mean 0 and standard deviation `σx`, `σy`
and `σz`.

# References
Verdes, P. F. "Assessing causality from multivariate time series." Physical
Review E 72.2 (2005): 026222.
"""
verdes(;u₀ = rand(3),
    ωy = 315, ωz = 80,
    σx = 0.01, σy = 0.01, σz = 0.01) =
    verdes(u₀, ωy, ωz, σx, σy, σz)


text:
# Simple example softmax function with associated utils
# Doesn't do a good job of making util function reusable
using cg
using MNIST

using ImageMagick
using Images

function softmax_regression(input_dim::Int, output_dim::Int, step_size::Float64)
    # Wow that's verbose - start overloading so constants become nodes automatically
    W1 = cg.variable(cg.zeros(cg.constant([output_dim, input_dim])), "W1")
    b1 = cg.variable(cg.zeros(cg.constant(output_dim)), "b1")

    input = cg.placeholder([input_dim], "input")
    label = cg.placeholder([output_dim], "label")
    unnormed = cg.broadcast("+", cg.matmul(W1, input), b1)
    prediction = cg.softmax_stable(unnormed)
    loss = cg.crossentropy(label, prediction)

    train = cg.sgd_optimizer(loss, [W1, b1], cg.constant(step_size))

    (input, label, prediction, loss, train)
end

function onehot(labels::AbstractArray; min::Int = 0, max::Int = 9)
    labels = round(Int, labels)
    result = zeros(Float64, max - min + 1, length(labels))
    for (i, label) in enumerate(labels)
        result[label - min + 1, i] = 1.0
    end
    result
end

function shuffle_data(data, labels)
    perm = shuffle(collect(1:size(labels, 2)))
    data[:, perm], labels[:, perm]
end

# labels = 1hot encoded
function train_steps(sess, batch_size, steps, data, labels)
    @assert size(data, 2) == size(labels, 2)
    @assert batch_size <= size(labels, 2)
    sdata, slabel = shuffle_data(data, labels)
    tsteps = 0
    bs = 1
    while tsteps < steps
        be = bs + batch_size - 1
        # Reshuffle and startover
        if be > size(labels, 2)
            sdata, slabel = shuffle(data, labels)
            bs = 1
            continue
        end
        tsteps += 1
        sess.values[input] = sdata[:, bs:be]
        sess.values[label] = slabel[:, bs:be]
        cg.interpret(sess, sgd)
    end
end

function compute_loss(session, loss, data, labels)
    session.values[input] = data
    session.values[label] = labels
    cg.interpret(session, loss)
end

function visualize_filters(filters, w, h)
    output = zeros(h, w * size(filters,1))
    for i in 1:size(filters,1)
        start = (i-1) * h
        filter = filters[i, :]
        high = maximum(filter)
        low = minimum(filter)
        scaled = (filter - low) / (high - low)
        output[:, start+1:start+w] = scaled
    end
    convert(Image, output)
end

function accuracy(data, labels)
    sess.values[input] = data
    result = cg.interpret(sess, output)
    correct = findmax(labels, 1)[2]
    prediction = findmax(result, 1)[2]
    sum(correct .== prediction) / length(correct)
end

# Scale the data to [0,1] to avoid overflow
test_raw = testdata()
train_raw = traindata()
test = (test_raw[1] / 255.0, test_raw[2])
train = (train_raw[1] / 255.0, train_raw[2])

input, label, output, loss, sgd = softmax_regression(784, 10, 0.01)
sess = cg.Session(sgd)

if !isinteractive()
    @show compute_loss(sess, loss, train[1], onehot(train[2]))
    @show compute_loss(sess, loss, test[1], onehot(test[2]))
    train_steps(100, 1000, train[1], onehot(train[2]))
    @show compute_loss(sess, loss, train[1], onehot(train[2]))
    @show compute_loss(sess, loss, test[1], onehot(test[2]))
end


text: include("GatesTest.jl")


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie
CairoMakie.activate!() # HIDE
let
    x = LinRange(-10,10,200)
    fg(x,μ,σ) = exp.(.-(x.-μ).^2 ./(2σ^2))./(σ*√(2π))
    fig = Figure(resolution = (700, 450))
    ax = Axis(fig, xlabel = "x", ylabel = "y")
    # filled curve 1
    linea1 = lines!(x, fg(x, 0.25, 1.5), color = "#E69F00")
    fillB1 = band!( x, fill(0,length(x)), fg(x, 0.25, 1.5); color =("#E69F00", 0.1))
    # filled curve 2
    linea2 = lines!(x, fg(x, 2, 1), color = "#56B4E9")
    fillB2 = band!( x, fill(0,length(x)), fg(x, 2, 1); color = ("#56B4E9", 0.1))
    # filled curve 3
    linea3 = lines!(x, fg(x, -1, 2), color = "#009E73")
    fillB3 = band!( x, fill(0,length(x)), fg(x, -1, 2); color = ("#009E73", 0.1))
    leg = Legend(fig, [[linea1, fillB1], [linea2, fillB2], [linea3, fillB3]],
    ["μ = 0.25, σ = 1.25", "μ = 2, σ = 1", "μ = -1, σ = 2"], halign = :left,
    valign = :top, tellheight = false, tellwidth = false, margin = (10, 10, 10, 10))
    fig[1, 1] = ax
    fig[1, 1] = leg
    fig
    #save("FilledLines.png"), fig, px_per_unit = 2.0)
    save(joinpath(@__DIR__, "output", "FilledLines.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status("CairoMakie") # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
##########################################################################
# Spherical Harmonics phase retrieval
##########################################################################

"Contraint in real space for spherical harmonics RAAR"
function A(density::SurfaceVolume, Koff::Integer)
	newdensity = deepcopy(density)

	#enforce support and positivity
	for k = 1:density.KMAX
		newdensity.surf[k] = k > Koff ?  zeros(num_val(density.LMAX)) : max.(real(density.surf[k]), 0.0)
	end
	return newdensity
end

"applies projections in Fourier space"
function B(density::SurfaceVolume, amplitudes::SurfaceVolume)
	# Fourier Transformation
	fourier = forward(density)

	# Extract the phases from the Fourier transformation
	for k = 1:fourier.KMAX
		fourier.surf[k] = amplitudes.surf[k] .* exp.(1im*angle.(fourier.surf[k]))
	end

	# Inverse Fourier Transform
	return backward(fourier)
end

"Spherical Harmonics RAAR implementation"
function sRAAR(intensity::SphericalHarmonicsVolume, iterations::Integer, beta0::Float64 = 0.75, beta_max::Float64 = 0.75, tau::Float64 = 100.0, plotting::Bool=false; cutoff_factor::Float64=0.5)

	Koff = ceil(Int64, cutoff_factor*intensity.KMAX)

	amplitudes = getSurfaceVolume(intensity)
	amplitudes.surf = map((y)->sqrt.(y),  map((x)-> max.(x, 0.0),  real(amplitudes.surf)) )

	# #Calculating the phases of a ball
	# ball = deepcopy(intensity)
	# for k = 1:ball.KMAX
	#     ball.coeff[k] = zeros(Complex{Float64}, Base.size(ball.coeff[k]))
	#     if k <= Integer(ball.KMAX/2) setc(ball, k, 0, 0, 1.0) end
	# end
	# ball_phases = getSurfaceVolume(forward(ball))

	# Setting up the starting amplitudes with random initial phases
	fourier = deepcopy(amplitudes)
	for k = 1:fourier.KMAX
		phases = pi * 2 * rand(num_val(fourier.LMAX))
		# phases = angle(ball_phases.surf[k])
		fourier.surf[k] = amplitudes.surf[k] .* exp.(1im * phases)
	end

	density = backward(fourier)
	tmp1 = 2.0*B(density, amplitudes) - density #reflection
	diffa = deepcopy(tmp1)
	tmp2 = deepcopy(tmp1)

	# progress = Progress(iterations, 1, "Phase retrieval...", 50)
	diff_dens_list = Float64[]
	diff_amp_list = Float64[]

	for n = 1:iterations
		#Alternate reflections
		beta = exp((-n/tau)^3.0)*beta0 + (1.0-exp((-n/tau)^3.0))*beta_max
		tmp3 = A(tmp1, Koff)
		tmp_u = 0.5*(beta*(2.0*tmp3 -tmp1) + (1.0-beta)*tmp1 + density)
		tmp2 = B(tmp_u, amplitudes)

		#Difference tracking in real space
		diff_dens = sum(abs, density - tmp_u)/sum(abs, density)
		push!(diff_dens_list, diff_dens)

		diff_amp = 0.0
		#Difference tracking in Fourier space
		# d_fourier = forward(tmp_u)
		# diff_amp = sum(abs, [sum(abs, amplitudes.surf[k]) > 1e-16 ? sum(abs,  abs(d_fourier.surf[k]) - amplitudes.surf[k])/sum(abs, amplitudes.surf[k]) : 0.0 for k = 1:density.KMAX])
		# push!(diff_amp_list, diff_amp)

		#Last step
		density = deepcopy(tmp_u)
		tmp1 = 2.0*tmp2 - tmp_u

		#show progress
		# progress.desc = "$(diff_dens) ($(diff_amp)): "
		# next!(progress)
	end
	if plotting
		plot(collect(1:iterations), diff_dens_list, label="diff dens")
		plot(collect(1:iterations), diff_amp_list, label="diff amp")
	end

	density = A(tmp2, Koff)
	return getSphericalHarmonicsVolume(density)
end

##########################################################################
# Cubic phase retrieval
##########################################################################

"Creates a cube of 1.0s and 0.0s with the non-zero region in the center of the cube"
function createCenterCube(cubesize, ratio)
	core = ceil(cubesize/ratio)
	edge = floor( (cubesize - core)/2 )
	center = cubesize/2
	println("Creating support cube with cubesize:$cubesize,ratio:$ratio,edge:$edge,core:$core")
	c = [x > edge && x <= cubesize-edge &&
	y > edge && y <= cubesize-edge &&
	z > edge && z <= cubesize-edge ? 1.0 : 0.0 for x=1:cubesize, y=1:cubesize, z=1:cubesize]
	# c = [ norm([x,y,z] - [center, center, center]) < cubesize/4 ? 1.0 : 0.0 for x=1:cubesize, y=1:cubesize, z=1:cubesize  ]
	# saveCube(c, cubesize, "supportCube.mrc")
	return c
end
#
#
# "Creates a cube with a gaussian kernel of defined sigma"
# function gaussian_kernel(size, sigma)
#     sigma = max(sigma, 0.0)
#     local center = Float64[ceil(size/2),ceil(size/2),ceil(size/2)]
#     kern = [exp(-norm([x,y,z]-center)^2/sigma) for x=1:size, y=1:size, z=1:size]
# end
#
# "Applies a gaussian kernel to a data cube via FFT"
# function gaussian_filter(data, kernel)
#     fkern = fft(ifftshift(kernel))
#     fdata = fft((data))
#     return real((ifft(fkern.*fdata)))
# #     return ones(10,10,10)
# end
#
# "Symmetrizes the cube and multiplies a factor 1.2 to the center voxel as preparation for phase retrieval (this is emprical)"
# function symmetrize_cube(cube)
# 	size = Base.size(cube)[1]
# 	mid = ceil(Int,cubesize/2)
# 	newcube = [(cube[x,y,z] + cube[(size-x)+1,(size-y)+1,(size-z)+1])*0.5 for x = 1:size, y=1:size, z=1:size]
# 	newcube[mid,mid,mid] = newcube[mid,mid,mid]*1.15
# 	return newcube
# end
#
# "Cubic RAAR phase retrieval"
# function craar(coeff,iter=501, b1 = 0.85, b2 = 0.99, tau = 350, K=KMAX, cubesize=cubesize, outputfile="density.mrc")
# 	supportCube = createCenterCube(cubesize, 2)
#
# 	local newIntensCube
# 	if length(coeff) == KMAX
# 		newIntensCube = getCube(coeff, cubesize, K)
# 	else
# 		newIntensCube = coeff
# 	end
#
# 	# newIntensCube = symmetrize_cube(newIntensCube)
# 	gkernel = gaussian_kernel(cubesize, 1.0)
# 	# newIntensCube = gaussian_filter(newIntensCube, gkernel)
#
# 	# newIntensCube = deepcopy(intensCube)
# 	# newIntensCube[floor(cubesize/2+1),floor(cubesize/2+1),floor(cubesize/2+1)] = 0.0#intensCube[cubesize/2+1,cubesize/2+1,cubesize/2+1]
# 	# newIntensCube[cubesize/2+1,cubesize/2+1,cubesize/2+1] = reduce(max, newIntensCube)
# 	amplitudes = sqrt(max(real(newIntensCube), 0.0))
# 	intensStart = amplitudes.* exp( 2.0*pi*rand(cubesize,cubesize,cubesize)*1im)
#
# 	amplitudes = ifftshift(amplitudes)
# 	intensStart = ifftshift(intensStart)
# 	supportCube = ifftshift(supportCube)
#
# 	res = RAAR(iter, b1, b2, tau, amplitudes, intensStart,supportCube, gkernel)
# 	res = fftshift(res)
# 	saveCube(res, cubesize, outputfile)
# 	# res = [res[cubesize-i+1,j,k] for i=1:cubesize,j=1:cubesize,k=1:cubesize]
# 	# saveCube(res, cubesize,"res_m.mrc")
# 	# say("Done")
#
# end
#
# "Magnitude projection for RAAR"
# function MagProj(u, M)
#
# 	epsilon = 3e-16;
#
# 	modsq_u = abs(u).^2
# 	denom2 = modsq_u + epsilon
# 	denom = sqrt(denom2)
# 	r_eps = modsq_u ./ denom - M
# 	dr_eps = (denom2 + epsilon)./(denom2 .* denom)
# 	unew = (1- dr_eps.*r_eps) .* u
# 	return unew
# end
#
# # function RA(data, gkernel, support)
# # 	filtered = gaussian_filter(data, gkernel)
# # 	threshold = maximum(filtered) * 0.1
# # 	new_support = map((x) -> x > threshold ? 1.0: 0.0, filtered)
# # 	u = data.*support.*new_support
# # 	return max(u, 0.0)
# # end
#
# "Real space constraints"
# function RA(data, gkernel, support)
# 	u = data.*support
# 	return max(real(u), 0.0)
# end
#
# "Actual RAAR implementation"
# function RAAR(maxit, beta0, beta_max, tau, amplitudes, intens, support, gkernel)
#
# 	# g = 0.0
# 	beta = beta0
# 	u = ifft(intens)
# 	m_proj = MagProj(intens,amplitudes)
# 	newu = ifft(m_proj)
#
# 	tmp1 = real(2*newu - u)
#
# 	@showprogress 1 "Phasing..." 50 for iter=1:maxit
# 		beta = exp((-iter/tau)^3.0)*beta0 + (1-exp((-iter/tau)^3.0))*beta_max
#
# 		r_sp = 2*RA(tmp1, gkernel, support)  - tmp1
# 		tmp_u = 0.5*(beta*r_sp + (1-beta)*tmp1 +u)
# 		u = real(tmp_u)
# 		intens = fft(u)
# 		m_proj = MagProj(intens,amplitudes)
# 		newu = ifft(m_proj)
# 		tmp1 = real(2*newu - u)
# 	end
# 	return RA(real(newu), gkernel, support)
# end


text:
function unwrap(err)
    if !(:error in fieldnames(typeof(err)))
        return err
    end
    unwrap(err.error)
end

macro m_throw(mac)
    quote
        try
            eval($(esc(Expr(:inert, mac))))
        catch _e
            throw(unwrap(_e))
        end
    end
end