text:
"""
    Video

Defines the video canvas for an animation.

# Fields
- `width::Int` the width in pixel
- `height::Int` the height in pixel
- `objects::Vector{AbstractObject}` the objects defined in this video
- `background_frames::Vector{Int}` saves for which frames a background is defined
- `defs::Dict{Symbol, Any}` Some definitions which should be accessible throughout the video.
"""
mutable struct Video
    width::Int
    height::Int
    objects::Vector{AbstractObject}
    background_frames::Vector{Int}
    defs::Dict{Symbol,Any}
end

"""
    CURRENT_VIDEO

holds the current video in an array to be declared as a constant
The current video can be accessed using CURRENT_VIDEO[1]
"""
const CURRENT_VIDEO = Array{Video,1}()

"""
    Video(width, height)

Create a video with a certain `width` and `height` in pixel.
This also sets `CURRENT_VIDEO`.
"""
function Video(width, height)
    # some luxor functions need a drawing ;)
    Drawing()
    video = Video(width, height, AbstractObject[], Int[], Dict{Symbol,Any}())
    if isempty(CURRENT_VIDEO)
        push!(CURRENT_VIDEO, video)
    else
        CURRENT_VIDEO[1] = video
    end
    empty!(CURRENT_OBJECT)
    return video
end


text:
include("misc.jl")
include("findMin.jl")

function softmaxObj(w,X,y, k)
        (n, d) = size(X)
        f = 0
        for i = 1:n
          b = 0
          for c = 1:k
            b += exp(w[:, c]' * X[i, :])
          end
          f += -w[:, y[i]]' * X[i, :] + log(b)
        end
	g = zeros(d,k)
        for i = 1:n
          b = 0
          for c = 1:k
            b += exp(w[:, c]' * X[i, :])
          end
          for c = 1:k
            if y[i] == c
	      g[:,c] += -X[i, :]
            end
            g[:,c] += 1/b * exp(w[:, c]' * X[i, :]) * X[i, :]
          end
        end
        #a = -w[:, y]*X
        #b = w[:, ones(k, size(y,2)) * (1:Int(k))']*X
	#f = sum(a + log(b))
	display(f)
	display(g)
	return (f,g)
end

# Multi-class softmax classifier
function softmaxClassifier(X,y)
	(n,d) = size(X)
	k = maximum(y)

	W = zeros(d,k)

	display(size(W))

	# Each binary objective has the same features but different lables
	funObj(w) = softmaxObj(w,X,y, k)

	W = findMin(funObj, W, verbose=true)#, derivativeCheck=true)

	# Make linear prediction function
	#predict(Xhat) = findmax(Xhat*W)[2]
	predict(Xhat) = mapslices(indmax,Xhat*W,2)

	return LinearModel(predict,W)
end

ext: jl


text:
struct DriftPath{T <: SSDFloat}
    e_path::Vector{<:AbstractCoordinatePoint{T}}
    h_path::Vector{<:AbstractCoordinatePoint{T}}
end

function get_velocity_vector(interpolation_field::Interpolations.Extrapolation{<:StaticVector{3}, 3}, point::CartesianPoint{T})::CartesianVector{T} where {T <: SSDFloat}
    return CartesianVector{T}(interpolation_field(point.x, point.y, point.z))
end

@inline function get_velocity_vector(interpolated_vectorfield, point::CylindricalPoint{T}) where {T <: SSDFloat}
    return CartesianVector{T}(interpolated_vectorfield(point.r, point.φ, point.z))
end

include("ChargeDriftCylindrical.jl")
include("ChargeDriftCartesian.jl")

function _drift_charges( detector::SolidStateDetector{T}, grid::Grid{T, 3}, point_types::PointTypes{T, 3},
                        starting_points::Vector{CartesianPoint{T}},
                        velocity_field_e::Interpolations.Extrapolation{SVector{3, Float64}, 3},
                        velocity_field_h::Interpolations.Extrapolation{SVector{3, Float64}, 3};
                        Δt::T = T(1f-9), n_steps::Int = 2000, verbose::Bool = true)::Vector{DriftPath{T}} where {T <: SSDFloat}

    drift_paths::Vector{DriftPath{T}} = Vector{DriftPath{T}}(undef, length(starting_points))

    for i in eachindex(starting_points)
        drift_path_e::Vector{CartesianPoint{T}} = zeros(CartesianPoint{T}, n_steps )#Vector{CartesianPoint{T}}(undef, n_steps)
        drift_path_h::Vector{CartesianPoint{T}} = zeros(CartesianPoint{T}, n_steps )#Vector{CartesianPoint{T}}(undef, n_steps)
        drift_charge!(drift_path_e, detector, point_types, grid, starting_points[i], Δt, velocity_field_e, verbose = verbose)
        drift_charge!(drift_path_h, detector, point_types, grid, starting_points[i], Δt, velocity_field_h, verbose = verbose)
        drift_paths[i] = DriftPath{T}( drift_path_e, drift_path_h )
    end

    return drift_paths
end

function _drift_charge( detector::SolidStateDetector{T}, grid::Grid{T, 3}, point_types::PointTypes{T, 3},
                       starting_point::CartesianPoint{T},
                       velocity_field_e::Interpolations.Extrapolation{SVector{3, Float64}, 3},
                       velocity_field_h::Interpolations.Extrapolation{SVector{3, Float64}, 3};
                       Δt::T = T(1f-9), n_steps::Int = 2000, verbose::Bool = true)::Vector{DriftPath{T}} where {T <: SSDFloat}
    return _drift_charges(detector, grid, point_types, [starting_point], velocity_field_e, velocity_field_h, Δt = Δt, n_steps = n_steps, verbose = verbose)
end


text:
# https://www.philipzucker.com/egraph-2/
# https://github.com/philzook58/EGraphs.jl/blob/main/src/matcher.jl
# https://www.hpl.hp.com/techreports/2003/HPL-2003-148.pdf
# TODO support destructuring
# DONE support type assertions

# ematching seems to be faster without spawning tasks

# we keep a pair of EClass, Any in substitutions because
# when evaluating dynamic rules we also want to know
# what was the value of a matched literal
const Sub = Base.ImmutableDict{Any, Tuple{EClass, Any}}

# https://www.hpl.hp.com/techreports/2003/HPL-2003-148.pdf
# page 48
"""
From [https://www.hpl.hp.com/techreports/2003/HPL-2003-148.pdf](https://www.hpl.hp.com/techreports/2003/HPL-2003-148.pdf)
The iterator `ematchlist` matches a list of terms `t` to a list of E-nodes by first finding
all substitutions that match the first term to the first E-node, and then extending
each such substitution in all possible ways that match the remaining terms to
the remaining E-nodes. The base case of this recursion is the empty list, which
requires no extension to the substitution; the other case relies on Match to find the
substitutions that match the first term to the first E-node.
"""
function ematchlist(e::EGraph, t::Vector{Any}, v::Vector{Int64}, sub::Sub)::Vector{Sub}
    c = Vector{Sub}()

    if length(t) != length(v) || length(t) == 0 || length(v) == 0
        push!(c, sub)
    else
        for sub1 in ematch(e, t[1], v[1], sub)
            for sub2 in ematchlist(e, t[2:end], v[2:end], sub1)
                push!(c, sub2)
            end
        end
    end
    return c
end

# sub should be a map from pattern variables to Id
function ematch(e::EGraph, t::Symbol, v::Int64, sub::Sub; lit=nothing)::Vector{Sub}
    if haskey(sub, t)
        return find(e, first(sub[t])) == find(e, v) ? [sub] : []
    else
        return [ Base.ImmutableDict(sub, t => (EClass(find(e, v)), lit)) ]
    end
end

function ematch(e::EGraph, t, v::Int64, sub::Sub; lit=nothing)::Vector{Sub}
    c = Vector{Sub}()
    id = find(e,v)
    for n in e.M[id]
        if (t isa QuoteNode ? t.value : t) == n
            if haskey(sub, t)
                union!(c, find(e, first(sub[t])) == id ? [sub] : [])
            else
                union!(c, [ Base.ImmutableDict(sub, t => (EClass(id), n))])
            end
            # union!(c, ematchlist(e, t.args[start:end], n.args[start:end] .|> x -> x.id, sub))
        end
    end
    return c
end

function ematch(e::EGraph, t::Expr, v::Int64, sub::Sub; lit=nothing)::Vector{Sub}
    c = Vector{Sub}()

    for n in e.M[find(e,v)]
        if isexpr(t, :(::))
            # right hand of type assertion
            # tr = t.args[2]
            typ = t.args[2]

            # println(n)
            # println(typ, " ", typeof(typ))
            # if tr isa Type
            #     typ = tr
            # elseif tr isa Symbol
            #     if haskey(sub, tr)
            #         typ = sub[tr][2]
            #     else
            #         # add the type to the egraph
            #         type_id = add!(e, typeof(n))
            #         sub =  Base.ImmutableDict(sub, t.args[2] => (type_id, typeof(n)))
            #         typ = typeof(n)
            #         # union!(c, ematch(e, t, v, sub))
            #         # continue
            #     end
            # # elseif isexpr(tr, :curly)
            # # TODO allow for parametric type variables
            # # see https://dl.acm.org/doi/pdf/10.1145/3276483
            # else
            #     error("Unsupported type assertion '", t, "'")
            # end

            !(typeof(n) <: typ) && continue
            union!(c, ematch(e, t.args[1], v, sub; lit=n))
            continue
        end

        # otherwise ematch on an expr
        (!(n isa Expr) || getfunsym(n) != getfunsym(t)) && continue
        union!(c, ematchlist(e, getfunargs(t), getfunargs(n) .|> x -> x.id, sub))
    end
    return c
end


text:
function dive_in(mod)
    for n in names(mod, true)
        @eval import $(Symbol(mod)): $n
    end
end


text:
# The implementation of caches is influenced by the SparseDiffTools.jl. We implement this
# custom code since the structure of the system models is not compatible with the functionalities
# in SparseDiffTools

abstract type Cache end

struct JacobianCache{F, U <: ForwardDiff.Dual} <: Cache
    f!::F
    bus_count::Int
    ode_output::Vector{Float64}
    ode_output_dual::Vector{U}
    branches_ode::Vector{Float64}
    branches_ode_dual::Vector{U}
    current_balance::Vector{Float64}
    current_balance_dual::Vector{U}
    inner_vars::Vector{Float64}
    inner_vars_dual::Vector{U}
    global_vars::Vector{Float64}
    global_vars_dual::Vector{U}
end

function JacobianCache{U}(F, inputs::SimulationInputs) where {U <: ForwardDiff.Dual}
    n_inj = get_injection_n_states(inputs)
    n_branches = get_branches_n_states(inputs)
    @debug "injection ode size = $n_inj branches ode size = $n_branches"
    bus_count = get_bus_count(inputs)
    inner_vars_count = get_inner_vars_count(inputs)
    n_global_vars = length(keys(get_global_vars_update_pointers(inputs)))
    return JacobianCache{typeof(F), U}(
        F,
        bus_count,
        zeros(Float64, n_inj),
        zeros(U, n_inj),
        zeros(Float64, n_branches),
        zeros(U, n_branches),
        zeros(Float64, 2 * bus_count),
        zeros(U, 2 * bus_count),
        zeros(Float64, inner_vars_count),
        zeros(U, inner_vars_count),
        setindex!(zeros(Float64, n_global_vars), 1.0, GLOBAL_VAR_SYS_FREQ_INDEX),
        setindex!(zeros(U, n_global_vars), convert(U, 1.0), GLOBAL_VAR_SYS_FREQ_INDEX),
    )
end

get_current_injections_r(jc::JacobianCache, ::Type{Float64}) =
    view(jc.current_balance, 1:(jc.bus_count))
get_current_injections_i(jc::JacobianCache, ::Type{Float64}) =
    view(jc.current_balance, (jc.bus_count + 1):(2 * jc.bus_count))
get_ode_output(jc::JacobianCache, ::Type{Float64}) = jc.ode_output
get_branches_ode(jc::JacobianCache, ::Type{Float64}) = jc.branches_ode
get_current_balance(jc::JacobianCache, ::Type{Float64}) = jc.current_balance
get_inner_vars(jc::JacobianCache, ::Type{Float64}) = jc.inner_vars
get_global_vars(jc::JacobianCache, ::Type{Float64}) = jc.global_vars

get_current_injections_r(jc::JacobianCache, ::Type{T}) where {T <: ForwardDiff.Dual} =
    reinterpret(T, jc.current_balance_dual[1:(jc.bus_count)])
get_current_injections_i(jc::JacobianCache, ::Type{T}) where {T <: ForwardDiff.Dual} =
    reinterpret(T, jc.current_balance_dual[(jc.bus_count + 1):end])
get_ode_output(jc::JacobianCache, ::Type{T}) where {T <: ForwardDiff.Dual} =
    reinterpret(T, jc.ode_output_dual)
get_branches_ode(jc::JacobianCache, ::Type{T}) where {T <: ForwardDiff.Dual} =
    reinterpret(T, jc.branches_ode_dual)
get_current_balance(jc::JacobianCache, ::Type{T}) where {T <: ForwardDiff.Dual} =
    reinterpret(T, jc.current_balance_dual)
get_inner_vars(jc::JacobianCache, ::Type{T}) where {T <: ForwardDiff.Dual} =
    reinterpret(T, jc.inner_vars_dual)
get_global_vars(jc::JacobianCache, ::Type{T}) where {T <: ForwardDiff.Dual} =
    reinterpret(T, jc.global_vars_dual)

struct SimCache{F} <: Cache
    f!::F
    bus_count::Int
    ode_output::Vector{Float64}
    branches_ode::Vector{Float64}
    current_balance::Vector{Float64}
    inner_vars::Vector{Float64}
    global_vars::Vector{Float64}
end

function SimCache(f!, inputs::SimulationInputs)
    n_inj = get_injection_n_states(inputs)
    n_branches = get_branches_n_states(inputs)
    @debug "injection ode size = $n_inj branches ode size = $n_branches"
    bus_count = get_bus_count(inputs)
    inner_vars_count = get_inner_vars_count(inputs)
    n_global_vars = length(keys(get_global_vars_update_pointers(inputs)))
    return SimCache{typeof(f!)}(
        f!,
        bus_count,
        zeros(Float64, n_inj),
        zeros(Float64, n_branches),
        zeros(Float64, 2 * bus_count),
        zeros(Float64, inner_vars_count),
        setindex!(zeros(Float64, n_global_vars), 1.0, GLOBAL_VAR_SYS_FREQ_INDEX),
    )
end

function get_current_injections_r(sc::SimCache, ::Type{Float64})
    return view(sc.current_balance, 1:(sc.bus_count))
end

function get_current_injections_i(sc::SimCache, ::Type{Float64})
    return view(sc.current_balance, ((sc.bus_count + 1):(sc.bus_count * 2)))
end

get_ode_output(sc::SimCache, ::Type{Float64}) = sc.ode_output
get_branches_ode(sc::SimCache, ::Type{Float64}) = sc.branches_ode
get_current_balance(sc::SimCache, ::Type{Float64}) = sc.current_balance
get_inner_vars(sc::SimCache, ::Type{Float64}) = sc.inner_vars
get_global_vars(sc::SimCache, ::Type{Float64}) = sc.global_vars

get_ω_sys(cache::Cache, T::Type{<:Union{Float64, ForwardDiff.Dual}}) =
    get_global_vars(cache, T)[GLOBAL_VAR_SYS_FREQ_INDEX]


text:
function lgd10()
#
#  test polynomial suggested by Goedecker
#   Legendre polynomial of degree 10
#
    p = lgd(10);
    z = [0.97390652851717
   0.86506336668898
   0.67940956829903
  -0.97390652851717
  -0.86506336668899
  -0.67940956829903
   0.43339539412925
  -0.43339539412925
   0.14887433898163
  -0.14887433898163]
    z = [z ones(10)]
    p, PolyZeros(z)
end


text:
p = [ 10*rand(2) for i in 1:100 ]

ijcell(p,cutoff) = [ trunc(Int64,p[1]/cutoff) + 1, trunc(Int64,p[2]/cutoff) + 1 ]

function cells(p,cutoff,side)
  ncells_x = trunc(Int64,side/cutoff) + 1
  ncells_y = trunc(Int64,side/cutoff) + 1
  celllist = zeros()
  for i in 1:length(p)

  end
  return celllist
end


text:
module FreeType

include("LibFreeType.jl")
using .LibFreeType

include("compat.jl")

end


text:
#!/home/mkschleg/opt/bin/julia
#SBATCH -o test_err.out # Standard output
#SBATCH -e test_out.err # Standard error
#SBATCH --mem-per-cpu=1000M # Memory request of 1 GB
#SBATCH --time=01:00:00 # Running time of 10 minutes
#SBATCH --ntasks=10
#SBATCH --cpus-per-task=1
#SBATCH --account=def-whitem

using Pkg
Pkg.activate(".")

using Reproduce

const save_loc = "test_exp"
const exp_file = "examples/experiment.jl"
const exp_module_name = :Main
const exp_func_name = :main_experiment

function test_experiment()

    arg_list = [
    ["--opt1", "47"]]

    static_args = ["--steps", "102902"]

    args_iterator = ArgLooper(arg_list, static_args, 5:8, "--opt2")

    experiment = Experiment(save_loc,
                            exp_file,
                            exp_module_name,
                            exp_func_name,
                            args_iterator)

    create_experiment_dir(experiment)
    add_experiment(experiment; settings_dir="settings")
    ret = job(experiment; num_workers=6, extra_args=[save_loc])
    post_experiment(experiment, ret)

end

test_experiment()


text:
"""
    gc_windows(positions::Array, refseq::NamedTuple; window=100, step=10, noslide=false)

Calculate the average GC content for genomic regions using a sliding window. Positions are input as an array of positions (1:200-300). Note that a genomic feature may be skipped by inputing a 2D matrix of postions (e.g. 1:100-200  1:230-300).

The positions array is created using the readbed function and the refseq input is generated using the readfasta function.
"""
function gc_windows(positions::Array, refseq::NamedTuple; window=100, step=10, noslide=false)

    ref = refseq
    m = []
    xvals = []

    for i in 1:size(positions,1)

        if size(positions,2) == 2
            s1 = getref(positions[i,1], ref)
            s2 = getref(positions[i,2], ref)
            s = s1 * s2
        else
            s = getref(positions[i], ref)
        end

        if noslide == true
            window = length(s)
            step = length(s)
        end

        gcdata = sliding_window(s, window, step, gc_simple)
        gcdata = round.(gcdata; digits=4)
        push!(m, gcdata)

    end

    return m

end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct KauraPLL <: FrequencyEstimator
        ω_lp::Float64
        kp_pll::Float64
        ki_pll::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
    end

Parameters of a Phase-Locked Loop (PLL) based on Kaura, Vikram, and Vladimir Blasko. "Operation of a phase locked loop system under distorted utility conditions." IEEE Transactions on Industry applications 33.1 (1997): 58-63.

# Arguments
- `ω_lp::Float64`: PLL low-pass filter frequency (rad/sec), validation range: (0, nothing)
- `kp_pll::Float64`: PLL proportional gain, validation range: (0, nothing)
- `ki_pll::Float64`: PLL integral gain, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`
"""
mutable struct KauraPLL <: FrequencyEstimator
    "PLL low-pass filter frequency (rad/sec)"
    ω_lp::Float64
    "PLL proportional gain"
    kp_pll::Float64
    "PLL integral gain"
    ki_pll::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
end

function KauraPLL(ω_lp, kp_pll, ki_pll, ext=Dict{String, Any}(), )
    KauraPLL(ω_lp, kp_pll, ki_pll, ext, [:vd_pll, :vq_pll, :ε_pll, :θ_pll], 4, )
end

function KauraPLL(; ω_lp, kp_pll, ki_pll, ext=Dict{String, Any}(), )
    KauraPLL(ω_lp, kp_pll, ki_pll, ext, )
end

# Constructor for demo purposes; non-functional.
function KauraPLL(::Nothing)
    KauraPLL(;
        ω_lp=0,
        kp_pll=0,
        ki_pll=0,
        ext=Dict{String, Any}(),
    )
end

"""Get KauraPLL ω_lp."""
get_ω_lp(value::KauraPLL) = value.ω_lp
"""Get KauraPLL kp_pll."""
get_kp_pll(value::KauraPLL) = value.kp_pll
"""Get KauraPLL ki_pll."""
get_ki_pll(value::KauraPLL) = value.ki_pll
"""Get KauraPLL ext."""
get_ext(value::KauraPLL) = value.ext
"""Get KauraPLL states."""
get_states(value::KauraPLL) = value.states
"""Get KauraPLL n_states."""
get_n_states(value::KauraPLL) = value.n_states

"""Set KauraPLL ω_lp."""
set_ω_lp!(value::KauraPLL, val::Float64) = value.ω_lp = val
"""Set KauraPLL kp_pll."""
set_kp_pll!(value::KauraPLL, val::Float64) = value.kp_pll = val
"""Set KauraPLL ki_pll."""
set_ki_pll!(value::KauraPLL, val::Float64) = value.ki_pll = val
"""Set KauraPLL ext."""
set_ext!(value::KauraPLL, val::Dict{String, Any}) = value.ext = val
"""Set KauraPLL states."""
set_states!(value::KauraPLL, val::Vector{Symbol}) = value.states = val
"""Set KauraPLL n_states."""
set_n_states!(value::KauraPLL, val::Int64) = value.n_states = val


text:
using LinearAlgebra, SparseArrays
using JuMP
using CSV, DataFrames
using Base.Filesystem
using Printf
using JLD2

function load_sdpa_file(filepath)
    info = CSV.File(joinpath(dirname(filepath), "sdplib_info.csv"))  |> DataFrame
    objectives = Dict(zip(info[1], info[4]))
    filename =  join(split(basename(filepath), ".")[1:end-1], ".")
    objective = haskey(objectives, filename) ? objectives[filename] : NaN
    @printf("Problem: %s with optimal objective: %.4e\n", filename, objective)
    # Load Problem
    jld_filepath = join([split(filepath, ".")[1:end-1]; "jld2"], ".")
    if isfile(jld_filepath) # Attempt to load saved jld2 file
        @load jld_filepath c F A b
    else # Load directly from dat-s file
        print("Loading SDPA file. This might take some time... ")
        t = @elapsed c, F, A, b = _load_sdpa_file(string(jld_filepath[1:end-5], ".dat-s"))
        println(string("Done in ", t, " seconds!"))
        @save jld_filepath c F A b
    end
    return c, F, A, b, objective
end

function _load_sdpa_file(filename)
    #=
    Load problems of the SDPA format from "dat-s" files, i.e. problems of the form

    minimize c'x
    subject to  F1*x1 + ... + Fm*xm - F0 ⪰ 0

    with the real m-dimensional vector x as the decision variable.
    See https://freesoft.dev/program/129366650#sdpa-sparse-format for a detailed description.

    Note that the coefficient matrices Fi have common block diagonal structure. Furthermore, some of
    these blocks are themselves diagonal. Thus we can transform the original problem into

    minimize c'x
    subject to  F1,1*x1 + ... + Fm,1*xm - F0,1 ⪰ 0
                ...
                F1,l*x1 + ... + Fm,l*xm - F0,l ⪰ 0
                A*x ≥ b

    where Fi,j denotes the j-th diagonal block of the i-th coefficient matrix and the linear inequality constraints
    arised from subblocks of the coefficient matrices that have diagonal structure.

    The function returns
    - c::Vector the linear cost
    - F::Array{SparseMatrixCSC,2} F[i, j] is the j-th block of the (i-1)-th coefficient matrix
    - A::SparseMatrixCSC, b::Vector: for the linear inequality constraints
    - opt_val::Float64 the optimal value of the problem. Only supported from problems in the SDPLIB collection.
        These are extracted from a README file of the SDPLIB suite.
        Contains NaN if README is not present or if the extraction of the optimal value failed.
    =#
    file = open(filename);
    lines = readlines(file)
    close(file)

    # Discard all lines that contain comments
    indices = findall(line->length(line) == 0 || (line[1] != '*' && line[1] != '"'), lines)
    lines = lines[indices]

    m = parse(Int, strip(lines[1])) # length of decision variable
    blocks_sizes_string = split(replace(strip(lines[3], ['{','}','(',')']), "+" => ""))
    block_sizes = [parse(Int, str) for str in blocks_sizes_string]
    @assert length(block_sizes) == parse(Int, strip(lines[2]))

    # try first to split by whitespace as delimiter (also remove certain trailing, ending characters)
    c_strings_vector = split(replace(strip(lines[4], ['{','}','(',')']), "+" => ""))
    # otherwise try comma as delimiter
    if length(c_strings_vector) != m
        c_strings_vector = split(replace(strip(lines[4], ['{','}','(',')']), "+" => ""), ",")
    end
    c = [parse(Float64, str) for str in c_strings_vector]

    @assert length(c) == m
    indices = Vector(1:length(block_sizes));
    indices_map = [indices[block_sizes .> 1]; indices[block_sizes .<= 1]]
    diag_sizes = abs.(block_sizes[block_sizes .<= 1])
    block_sizes = block_sizes[block_sizes .> 1]

    F = [spzeros(block_sizes[j], abs(block_sizes[j])) for i = 1:m + 1, j = 1:length(block_sizes)]
    A = spzeros(sum(diag_sizes), m)
    b = zeros(size(A, 1))
    for line in lines[5:end]
        parsed_line = [parse(Float64, str) for str in split(line)]
        matrix_number = Int(parsed_line[1]) + 1
        block_number = findfirst(indices_map .== Int(parsed_line[2]))

        i = Int(parsed_line[3])
        j = Int(parsed_line[4])
        entry = parsed_line[5]

        if block_number <= length(block_sizes)
            F[matrix_number, block_number][i, j] = entry
            F[matrix_number, block_number][j, i] = entry
        else
            # Diagonal block - convert to linear inequalities
            @assert i == j
            idx = i + sum(diag_sizes[1:block_number - length(block_sizes) - 1])

            if matrix_number == 1
                b[idx] = entry
            else
                A[idx, matrix_number - 1] = entry
            end
        end
    end

    return c, F, A, b
end

function extract_upper_triangle(A::SparseMatrixCSC{Tv,Ti}, scaling_factor::Tv = one(Tv)) where {Tv,Ti}
   	result_nnz = Ti(nnz(A) / 2 + nnz(diag(A)) / 2)
   	nzind = zeros(Ti, result_nnz)
   	nzval = zeros(Tv, result_nnz)
   	n = size(A, 1)
    counter = 0
   	for j in 1:n, idx in A.colptr[j]:A.colptr[j + 1] - 1
        i = A.rowval[idx]
      		if i <= j
            counter += 1
            nzind[counter] = Ti(j * (j - 1) / 2 + i)
            if i == j
                nzval[counter] = A.nzval[idx]
            else
                nzval[counter] = scaling_factor * A.nzval[idx]
         			end
      		end
    end
   	return SparseVector{Tv,Ti}(Ti(n * (n + 1) / 2), nzind, nzval)
end

function solve_sdpa_jump(c, F, A, b, solver; kwargs...)
    m = length(c)
    model = Model(with_optimizer(solver; kwargs...))
    @variable(model, x[1:m])
    @constraint(model, A * x .>= b)
    for block_id in 1:size(F, 2)
        A_ = hcat([extract_upper_triangle(F[i + 1, block_id]) for i = 1:m]...)
        b_ = Vector(extract_upper_triangle(F[1, block_id]))
        block_size = size(F[1, block_id], 1)
        @constraint(model, A_ * x - b_ in MOI.PositiveSemidefiniteConeTriangle(block_size))
    end
    @objective(model, Min, dot(c, x))
    JuMP.optimize!(model)
    raw_solver = model.moi_backend.optimizer.model.optimizer.inner
    # raw_solver = MOI.get(model, MOI.RawSolver())
    return value.(x), JuMP.objective_value(model), JuMP.termination_status(model), raw_solver
end

function solve_sdpa_jump_dual(c, F, A, b, solver; kwargs...)
    model = Model(with_optimizer(solver; kwargs...))

    @variable(model, y[1:size(A, 1)] >= 0)
    num_blocks = size(F, 2)
    block_sizes = [size(F[1, block_id], 1) for block_id in 1:num_blocks]
    Y = [@variable(model, [1:block_size, 1:block_size], PSD) for block_size in block_sizes]

    constraint_sum = A' * y
    objective_sum = dot(y, b)
    for block_id in 1:num_blocks
        # Normally, we should also avoid the matrix dot product below, but in most solvers this doesn't matter
        add_to_expression!.(objective_sum, dot(F[1, block_id], Y[block_id]))
        for i = 1:length(c)
            f = SparseVector(F[i + 1, block_id][:])
            add_to_expression!(constraint_sum[i], dot(f.nzval, Y[block_id][f.nzind]))
        end
    end
    @constraint(model, constraint_sum .== c)
    @objective(model, Max, objective_sum)

    JuMP.optimize!(model)
    raw_solver = model.moi_backend.optimizer.model.optimizer.inner
    # raw_solver = MOI.get(model, MOI.RawSolver())

    return [value.(Y[i]) for i = 1:length(block_sizes)], value.(y),
        JuMP.objective_value(model),
        JuMP.termination_status(model),
        raw_solver
end

#=
function solve_sdpa_jump_b(c, F, A, b; solver=SCS.Optimizer, kwargs...)
    # A version of solve_sdpa_jump_b that does not explicitly construct A_, b_
    # It is significantly slower than the original function, solve_sdpa_jump.
    m = length(c)
    model = Model(with_optimizer(solver; kwargs...))
    @variable(model, x[1:m])
    @constraint(model, A*x .>= b)
    for block_id in 1:size(F, 2)
        block_size = size(F[1, block_id], 1)
        X = @variable(model, [1:block_size, 1:block_size], PSD) .- F[1, block_id]
        for i = 1:m
            add_to_expression!.(X, x[i].*F[i + 1, block_id])
        end
        @constraint(model, X .== 0)
    end
    @objective(model, Min, dot(c, x))
    JuMP.optimize!(model)
    return value.(x), JuMP.objective_value(model), JuMP.termination_status(model)
end
=#


text:
module Eigenbroetler

using FileIO
using FFTW
using FITSIO
using FITSIO.Libcfitsio
using Images

export Eigenbrot
export width, height
export size
export isFFT, save
export fft, fftx, ffty, removeDC
export ImageSetting
export LinearScale, LogScale, RootScale
export RealPart, ImagPart, Magn, Phase
export image, pad
export flipver, fliphor, swapxy
export pixel, coords
export scale_chirpz, scale_x_chirpz, scale_y_chirpz
export hilbert_x, hilbert_y
export fft_xshear, fft_yshear, rotate90, fourierRotation

include("colourmaps.jl")
include("eigenbrot.jl")
include("simple.jl")
include("fft.jl")
include("chirpz.jl")
include("transform.jl")

end


text:
#!/bin/bash
# -*- mode: julia -*-
#=
JULIA="${JULIA:-julia --color=yes --startup-file=no}"
exec ${JULIA} "$@" "${BASH_SOURCE[0]}"
=#

JL_PKG = "InteractiveCodeSearch"
VERSION >= v"0.7.0-DEV.5183" && using Pkg

if isfile("Project.toml") || isfile("JuliaProject.toml")
    if VERSION < v"0.7.0-DEV.5183"
        Pkg.clone(pwd())
        Pkg.build(JL_PKG)
    elseif VERSION >= v"1.1.0-rc1"
        Pkg.build(verbose=true)
    else
        Pkg.build()
    end
else
    Pkg.clone(pwd())
    if VERSION >= v"1.1.0-rc1"
        Pkg.build(JL_PKG; verbose=true)
    else
        Pkg.build(JL_PKG)
    end
end

# https://github.com/travis-ci/travis-build/blob/master/lib/travis/build/script/julia.rb


text:
"""
Enum that defines two possible orientations: ```H``` horizontal, ```V``` vertical.
"""
@enum Orientation H V

"""
Abstract type representing a single rectangular tile
"""
abstract type Tile end

"""
Type representing a rectangular tile with two subtiles.
If orientation is set to ```H```, ```subtiles[1]``` is at the left of ```subtiles[2]```, otherwise ```subtiles[1]``` is below ```subtiles[2]```.

Fields:
    - ```subtiles::NTuple{2,Tile}```: subtiles;
    - ```orientation::Oriantation```: orientation of subtiles.
"""
mutable struct TileLeaf <: Tile
    subtiles::NTuple{2,Tile}
    orientation::Orientation

    """
    Inner constructor for the TileLeaf type. Its subtiles must be compatible, i.e. must share their heaight if they are horizonatally oriented, and similarly if vercitally oriented.
    """
    function TileLeaf(subtiles,orientation)
        if orientation == H
            clause = height(subtiles[1]) == height(subtiles[2])
            errorMsg = "Heights of subtiles not compatible!"
        else
            clause = base(subtiles[1]) == base(subtiles[2])
            errorMsg = "Bases of subtiles not compatible!"
        end

        clause ? new(subtiles,orientation) : error(errorMsg)
    end
end

"""
Type representing a rectangular tile with no subtiles.

Fields:
    - ```b::Float64```: base of tile;
    - ```h::Float64```: height of tile;
"""
mutable struct TileNode <: Tile
    b::Float64 # base
    h::Float64 # height
end

"""
Retrieve recursively the total vertical span of ```tile```.
The vertical span is the height of the subtiles if they are stacked horizontally, and it is the sum of the heights of the subtiles it they are stacked vertically.
"""
height(tile::TileLeaf) = tile.orientation == H ? height(tile.subtiles[1]) : height(tile.subtiles[1])+height(tile.subtiles[2])
height(tile::TileNode) = tile.h

"""
Retrieve recursively the total horiznotal span of ```tile```.
See ```height``` for more informations.
"""
base(tile::TileLeaf) = tile.orientation == V ? base(tile.subtiles[1]) : base(tile.subtiles[1])+base(tile.subtiles[2])
base(tile::TileNode) = tile.b

"""
Retrieve both the horizontal and vertical spans of ```tile```.
"""
size(tile::Tile) = (base(tile),height(tile))

"""
Recursively rescale the spans of a ```tile``` by some factor.

```factor``` can be either:
    - a ```Real```, and both the base and height of ```tile``` will be rescaled by ```factor```
    - a ```(Real,Real)```, and the base and height of ```tile``` will be rescaled respectively by ```factor[1]``` and ```factor[2]```.

"""
function rescale!(factor::T, tile::TileNode) where T<:Real
    tile.b = tile.b * factor
    tile.h = tile.h * factor
    return tile
end
function rescale!(factor::T, tile::TileLeaf) where T<:Real
    rescale!(factor, tile.subtiles[1])
    rescale!(factor, tile.subtiles[2])
    return tile
end
function rescale!(factor::Tuple{T1,T2}, tile::TileNode) where {T1<:Real, T2<:Real}
    tile.b = tile.b * factor[1]
    tile.h = tile.h * factor[2]
    return tile
end
function rescale!(factor::Tuple{T1,T2}, tile::TileLeaf) where {T1<:Real, T2<:Real}
    rescale!(factor, tile.subtiles[1])
    rescale!(factor, tile.subtiles[2])
    return tile
end

"""
Analogous to ```rescale!```, but performs a ```deepcopy``` of ```tile``` before rescaling, so that the original ```tile``` is not modified.
"""
function rescale(factor::T, tile::Tile) where T<:Real
    return rescale!(factor, deepcopy(tile))
end
function rescale(factor::Tuple{T1,T2}, tile::Tile) where {T1<:Real, T2<:Real}
    return rescale!(factor, deepcopy(tile))
end


text:
__precompile__(true)

module NBody
  include("nbody/body.jl")
  include("nbody/nbody.jl")
  include("nbody/integrators.jl")

  export EvolveArgs,
          evolve,
          read_nbody
end


text:
const mode_keywords = [
    "normal",
    "multiply",
    "screen",
    "overlay",
    "darken",
    "lighten",
    "color-dodge",
    "color-burn",
    "hard-light",
    "soft-light",
    "difference",
    "exclusion",
    "hue",
    "saturation",
    "color",
    "luminosity",
]

"""
    parse(t::Type{BlendMode}, keyword)
    parse(t::Type{CompositeOperation}, keyword)

Parse a keyword string as a [`BlendMode`](@ref) or [`CompositeOperation`](@ref).
Keywords are case-insensitive, but hyphens cannot be omitted.

# Examples
```jldoctest
julia> parse(BlendMode, "color-burn")
BlendMode{Symbol("color-burn")}()

julia> parse(BlendMode, "soft-light") === BlendSoftLight
true

julia> parse(CompositeOperation, "source-over")
CompositeOperation{Symbol("source-over")}()

julia> parse(CompositeOperation, "source-over") === CompositeSourceOver
true

julia> parse(CompositeOperation, "SourceOver")
ERROR: ArgumentError: invalid keyword: SourceOver
```
"""
function Base.parse(::Type{BlendMode}, keyword::AbstractString)
    keyword in mode_keywords && return BlendMode{Symbol(keyword)}()
    kl = lowercase(strip(keyword))
    kl in mode_keywords && return BlendMode{Symbol(kl)}()
    throw(ArgumentError("invalid keyword: $keyword"))
end

const op_keywords = [
    "clear",
    "copy",
    "destination",
    "source-over",
    "destination-over",
    "source-in",
    "destination-in",
    "source-out",
    "destination-out",
    "source-atop",
    "destination-atop",
    "xor",
    "lighter",
]

function Base.parse(::Type{CompositeOperation}, keyword::AbstractString)
    keyword in op_keywords && return CompositeOperation{Symbol(keyword)}()
    kl = lowercase(strip(keyword))
    kl in op_keywords && return CompositeOperation{Symbol(kl)}()
    throw(ArgumentError("invalid keyword: $keyword"))
end


text:
function giant_component_size(g, present)
    n = nv(g)
    Q = Queue{Int}()
	label = zeros(Int,n)
	cnt = zeros(Int,n)
	maxcnt = 0
    @inbounds for u in vertices(g)
		if present[u]
	        label[u] != 0 && continue
	        label[u] = u
			cnt[u] += 1
	        enqueue!(Q, u)
	        while !isempty(Q)
	            src = dequeue!(Q)
	            for vertex in all_neighbors(g, src)
					if present[vertex]
		                if label[vertex] == 0
		                    enqueue!(Q, vertex)
		                    label[vertex] = u
							cnt[u] += 1
		                end
					end
	            end
	        end
			if cnt[u] > maxcnt
				maxcnt = cnt[u]
			end
		end
    end
    return maxcnt
end

function recover_add_nodes(g, attack_nodes)
	n = nv(g)
	max_comp_sizes = Int[]
	ds = IntDisjointSets(n)
	present = fill(true, n)
	present[attack_nodes] .= false
	comp_sizes = ones(Int, n)
	max_comp_size = 1

	# add edges between present nodes and record component sizes
	for i in vertices(g)
		if present[i]
			for j in neighbors(g,i)
				if present[j]
					newroot = merge_nodes!(ds, comp_sizes, i, j)
					if comp_sizes[newroot] > max_comp_size
						max_comp_size = comp_sizes[newroot]
					end
				end
			end
		end
	end

	for i in reverse(attack_nodes)
		for j in neighbors(g, i)
			if present[j]
				newroot = merge_nodes!(ds, comp_sizes, i, j)
				if comp_sizes[newroot] > max_comp_size
					max_comp_size = comp_sizes[newroot]
				end
			end
		end
		present[i] = true
		push!(max_comp_sizes, max_comp_size)
	end
	reverse(max_comp_sizes)
end


text:
using Documenter, MultipleScattering

makedocs(
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ), # changes urls if built locally
    sitename = "MultipleScattering.jl",
    authors = "Artur L. Gower and Jonathan Deakin",
    source = "src",
    # modules = [MultipleScattering],
    pages = [
        "Home" =>"index.md",
        "Manual" => [
            "manual/intro.md",
            "manual/source.md",
            "manual/shapes.md",
            "manual/particles.md",
            "manual/time_response.md",
            "manual/plot.md",
            "manual/new_types.md"
        ],
        "Library" => [
            "library/base.md",
            "library/acoustics.md",
            "library/random.md"
        ],
        "Theory" => "maths/README.md",
        "Examples" =>
        [
            "example/README.md",
            "example/box_size/README.md",
            "example/hankel_convergence/README.md",
            "example/moments/README.md",
            "example/near_surface_backscattering/README.md",
            "example/random_particles/README.md",
            "example/time_response_single_particle/README.md",
            "example/two_particles/README.md"
        ]
    ]
)

deploydocs(
    branch = "gh-pages",
    target = "build",
    versions = ["stable" => "v^", "v#.#.#"],
    repo = "github.com/JuliaWaveScattering/MultipleScattering.jl.git"
)


text:
@testset "dataitem.jl" begin
    @testset "Standardization of the Variable type" begin

        v = Variable("x", 100 .+ 5*randn(100))

        vstd = Extremes.standardize(v)

        m = mean(vstd.value)
        s = std(vstd.value)

        @test m ≈ 0.0 atol = sqrt(eps())
        @test s ≈ 1.0 atol = sqrt(eps())
        @test vstd.offset ≈ 100 atol = 5
        @test vstd.scale ≈ 5 atol = 1

        ṽ = Extremes.reconstruct(vstd)

        @test ṽ.name == v.name
        @test ṽ.value ≈ v.value atol = sqrt(eps())

        v = Variable("x", [-1; 0 ; 1])
        vstd = Extremes.standardize(v)

        @test vstd.offset == 0 # integer comparison
        @test vstd.scale == 1 # integer comparison

        v = Variable("x", ones(10))
        vstd = Extremes.standardize(v)

        @test vstd.scale == 1 # integer comparison

        vstd = VariableStd("x", [-1; 0 ; 1])
        @test vstd.offset == 0 # integer comparison
        @test vstd.scale == 1 # integer comparison

        @test_throws AssertionError Extremes.VariableStd("x", [-2; 0 ; 2])
        @test_throws AssertionError Extremes.VariableStd("x", [0; 1 ; 2])

        # Test on the transform function
        y = randn(3)
        x₁ = Extremes.standardize(Variable("x₁", [-1; 0; 1]))
        x₂ = Extremes.standardize(Variable("x₂", [0; 1; 2]))
        x₃ = Extremes.standardize(Variable("x₃", [0; 2; 4]))

        model = BlockMaxima(Variable("y", y), locationcov=[x₁; x₂], logscalecov = [x₃])

        θ̃ = collect(1.0:6.0)

        fm_std = MaximumLikelihoodEVA(model, θ̃)

        fm = Extremes.transform(fm_std)

        θ̂ = fm.θ̂

        @test θ̂[1] ≈ θ̃[1] - θ̃[2]*0/1 - θ̃[3]*1/1
        @test θ̂[2] ≈ θ̃[2]/1
        @test θ̂[3] ≈ θ̃[3]/1
        @test θ̂[4] ≈ θ̃[4] - θ̃[5]*2/2
        @test θ̂[5] ≈ θ̃[5]/2
        @test θ̂[6] ≈ θ̃[6]

    end

end


text:
#==========================================================================================
                                Constructors
==========================================================================================#
function DiaSymSemiseparableCholesky(U::AbstractArray, W::AbstractArray, ds::AbstractArray)
    return DiaSymSemiseparableCholesky(size(U,1),size(U,2),U,W,ds)
end
function DiaSymSemiseparableCholesky(U::AbstractArray, V::AbstractArray, σn, σf)
    n, p = size(U)
    W, dbar = dss_create_wdbar(σf*U, σf*V, ones(n)*σn^2)
    return DiaSymSemiseparableCholesky(n, p, σf*U, W, dbar)
end
function DiaSymSemiseparableCholesky(L::DiaSymSemiseparableMatrix)
    W, dbar = dss_create_wdbar(L.Ut, L.Vt, L.d)
    return DiaSymSemiseparableCholesky(L.n, L.p, L.Ut, W, dbar)
end
#==========================================================================================
                    Defining Matrix Properties / Overloading Base
==========================================================================================#
Base.Matrix(K::DiaSymSemiseparableCholesky) = getproperty(K,:L)
Base.size(K::DiaSymSemiseparableCholesky) = (K.n, K.n)
function Base.size(K::DiaSymSemiseparableCholesky,d::Int)
    return (1 <= d && d <=2) ? size(K)[d] : throw(ArgumentError("Invalid dimension $d"))
end
function Base.getindex(K::DiaSymSemiseparableCholesky, i::Int, j::Int)
	i > j && return dot(K.Ut[:,i], K.Wt[:,j])
	i == j && return K.d[i]
	return 0
end
function Base.getproperty(K::DiaSymSemiseparableCholesky, d::Symbol)
    if d === :U
        return UpperTriangular(triu(K.Wt'*K.Ut,1) + Diagonal(K.d))
    elseif d === :L
        return LowerTriangular(tril(K.Ut'*K.Wt,-1) + Diagonal(K.d))
    else
        return getfield(K, d)
    end
end
function Base.propertynames(F::DiaSymSemiseparableCholesky, private::Bool=false)
    return (:U, :L, (private ? fieldnames(typeof(F)) : ())...)
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")},
		 K::DiaSymSemiseparableCholesky{<:Any,<:AbstractArray,<:AbstractArray,<:AbstractArray})
    summary(io, K); println(io)
    show(io, mime, K.L)
end
#==========================================================================================
                            Overloading LinearAlgebra routines
==========================================================================================#
dss_logdet(d) = sum(log,d)
newlogdet(L::DiaSymSemiseparableCholesky) = dss_logdet(L.d)
newlogdet(L::Adjoint{<:Any,<:DiaSymSemiseparableCholesky}) = dss_logdet(L.parent.d)
LinearAlgebra.logdet(L::DiaSymSemiseparableCholesky) = dss_logdet(L.d)
LinearAlgebra.logdet(L::Adjoint{<:Any,<:DiaSymSemiseparableCholesky}) = dss_logdet(L.parent.d)
#==========================================================================================
                        Defining multiplication and inverse
==========================================================================================#
function LinearAlgebra.mul!(y::AbstractArray, L::DiaSymSemiseparableCholesky, b::AbstractArray)
    dss_tri_mul!(y, L.Ut, L.Wt, L.d, b)
    return y
end
function LinearAlgebra.mul!(y::AbstractArray, L::Adjoint{<:Any,<:DiaSymSemiseparableCholesky}, b::AbstractArray)
    dssa_tri_mul!(y, L.parent.Ut, L.parent.Wt, L.parent.d, b)
    return y
end
function LinearAlgebra.ldiv!(y::AbstractArray, L::DiaSymSemiseparableCholesky, b::AbstractArray)
    dss_forward!(y, L.Ut, L.Wt, L.d, b)
    return y
end
function LinearAlgebra.ldiv!(y::AbstractArray, L::Adjoint{<:Any,<:DiaSymSemiseparableCholesky}, b::AbstractArray)
    dssa_backward!(y, L.parent.Ut, L.parent.Wt, L.parent.d, b)
    return y
end
function (Base.:\)(L::DiaSymSemiseparableCholesky, b::AbstractVecOrMat)
    y = similar(b)
    dss_forward!(y, L.Ut, L.Wt, L.d, b)
    return y
end
function (Base.:\)(L::Adjoint{<:Any,<:DiaSymSemiseparableCholesky}, b::AbstractVecOrMat)
    y = similar(b)
    dssa_backward!(y, L.parent.Ut, L.parent.Wt, L.parent.d, b)
    return y
end

#### Inverse of a EGRQSCholesky using the Cholesky factorization ####
function inv(L::DiaSymSemiseparableCholesky)
	return L'\(L\Diagonal(ones(L.n)))
end
function inv(L::DiaSymSemiseparableCholesky, b::AbstractArray)
	return L'\(L\b)
end

#==========================================================================================
                        Relevant Linear Algebra Routines
==========================================================================================#
function fro_norm_L(L::DiaSymSemiseparableCholesky)
	return sum(squared_norm_cols(L.Ut, L.Wt, L.d))
end

function trinv(L::DiaSymSemiseparableCholesky)
	dbar = L.d
	Y, Z = dss_create_yz(L.Ut, L.Wt, dbar)
	return sum(squared_norm_cols(Y,Z, dbar.^(-1)))
end
function tr(L::DiaSymSemiseparableCholesky)
	return sum(L.d)
end
function tr(Ky::DiaSymSemiseparableCholesky, K::SymSemiseparableMatrix)
	p = Ky.p
	c = Ky.d
	U = K.Ut
	V = K.Vt
	Y, Z = dss_create_yz(Ky.Ut, Ky.Wt, Ky.d)
	b = 0.0
	P = zeros(p,p)
	R = zeros(p,p)
	@inbounds for k = 1:Ky.n
		yk = @view Y[:,k]
		zk = @view Z[:,k]
		uk = @view U[:,k]
		vk = @view V[:,k]
		cki = c[k]^(-1)
		b += yk'*P*yk + 2*yk'*R*uk*cki + uk'*vk*(cki^2)
		P += ((uk'*vk)*zk)*zk' + zk*(R*uk)' + (R*uk)*zk'
		R += zk*vk';
	end
	return b
end
#===========================================================================================
                Choleskyesky factorization of Higher-order quasiseparable matrices
===========================================================================================#
"""
    dss_create_wdbar(U, V, d)

Computes `W` and `dbar` such that, `L = tril(UW',-1) + diag(ds)`.
"""
function dss_create_wdbar(U, V, d)
    m, n = size(U)
    P  = zeros(m, m)
    W  = zeros(m, n)
    dbar = zeros(n)
    @inbounds for (u,v,w,i) in zip(eachcol(U),eachcol(V),eachcol(W),eachindex(d))
        w      .= v - P*u
        dbar[i] = sqrt(dot(u,w) + d[i])
        w      .= w/dbar[i]
        add_outer_product!(P,w) #
    end
    return W, dbar
end
#### Multiplying with Ld ####
function dss_tri_mul!(Y,U,W,ds,X)
    p     = size(U,1)
    n_rhs = size(X,2)
    Wbar  = zeros(p,n_rhs)
    @inbounds for (u,w,y,x,i) in zip(eachcol(U),eachcol(W),eachrow(Y),eachrow(X),eachindex(ds))
        add_Y_tri_diag!(y,u,Wbar,x,ds[i]) #Y[i,:] = tmpU'*Wbar + ds[i]*tmpX
        add_product!(Wbar,w,x) # Wbar  += tmpW*tmpX
    end
end
#### Adjoint of Ld ####
function dssa_tri_mul!(Y,U,W,ds,X)
    p     = size(U,1)
    n_rhs = size(X, 2)
    Ubar  = zeros(p,n_rhs)
    for (u,w,y,x,i) in zip(Iterators.reverse(eachcol(U)),
                           Iterators.reverse(eachcol(W)),
                           Iterators.reverse(eachrow(Y)),
                           Iterators.reverse(eachrow(X)),
                           Iterators.reverse(eachindex(ds)))
        add_Y_tri_diag!(y,w,Ubar,x,ds[i])    # Y[i,:] = tmpW'*Ubar + ds[i]*tmpX
        add_product!(Ubar,u,x)                    # Ubar = Ubar + tmpU*tmpX
    end
end
#### Forward substitution ####
function dss_forward!(X,U,W,ds,B)
    p     = size(U,1)
    n_rhs = size(B,2)
    Wbar  = zeros(p,n_rhs)
    for (u,w,x,b,i) in zip(eachcol(U),eachcol(W),eachrow(X),eachrow(B),eachindex(ds))
        x .= (b - Wbar'*u)/ds[i]
        add_inner_plus!(Wbar,w,x)
    end
end
#### Backward substitution ####
function dssa_backward!(X,U,W,ds,B)
    m, n = size(U)
    mx = size(B,2)
    Ubar = zeros(m,mx)
    for (u,w,b,x,i) in zip(Iterators.reverse(eachcol(U)),
                           Iterators.reverse(eachcol(W)),
                           Iterators.reverse(eachrow(B)),
                           Iterators.reverse(eachrow(X)),
                           Iterators.reverse(eachindex(ds)))
        x .= (b - Ubar'*w)/ds[i]
        add_inner_plus!(Ubar,u,x) # Ubar += tmpU .* X[i:i,:]
    end
end

#### Squared norm of columns of L = tril(UW',-1) + diag(dbar) ####
function squared_norm_cols(U,W,dbar)
    m, n = size(U)
    P = zeros(m, m)
    c = zeros(n)
    for (w,u,i) in zip(Iterators.reverse(eachcol(W)),
                       Iterators.reverse(eachcol(U)),
                       Iterators.reverse(eachindex(dbar)))
        c[i]  = dbar[i]^2 + w'*P*w
        add_outer_product!(P,u) # P += tmpU*tmpU'
    end
    return c
end
#### Implicit inverse of  L = tril(UW',-1) + diag(dbar) ####
function dss_create_yz(U, W,dbar)
    m, n = size(U)
    Y = zeros(n,m)
    Z = zeros(n,m)
    dss_forward!(Y, U, W, dbar, U')
    dssa_backward!(Z, U, W, dbar, W')
    # Probably best not to use inv
    return copy(Y'), copy((Z*inv(U*Z - Diagonal(ones(m))))')
end


text:
# TODO: automatically generate by brainflow CICD
function brainflow_url()
    url = "https://github.com/brainflow-dev/brainflow/releases/download/3.7.2/compiled_libs.tar"
    return url
end

max_stars_repo_name: p-riebs/brainflow

max_issues_repo_name: p-riebs/brainflow

max_forks_repo_name: p-riebs/brainflow


text:
include("test_utilities.jl")
include("../src/arnoldi.jl")

@testset "Basic Tests for Iterative Hessenberg" begin

    function Test1(N=10)
        @info "Testing basic full Hessenberg, direct run no checking: "
        A = rand(N, N)
        ih = Arnoldi(A, rand(N, 1))
        for II = 1:N
            ih()
        end
        display(GetHessenberMatrix(ih))
        display(GetOrthogonalMatrix(ih))
        return true
    end

    function Test2(N=10)
        @info "Checking the errors on the Hessenberg decomposition: "
        A = rand(N, N)
        ih = Arnoldi(A, rand(N, 1))
        for II = 1:N-2
            ih()
        end

        H = GetHessenberMatrix(ih)
        print("Matrix H: ")
        display(H)
        Q = GetOrthogonalMatrix(ih)
        print("Matrix Q:")
        display(Q)

        @assert (A*Q[:, 1:end-1] - Q*H).|>abs|>sum|>mean < 1e-10 "Decomposition Errors too big"
        return true
    end

    function Test3(N=10)
        @info "Checking Tridiagonal Decomposition on Hermitian Matrices: "
        A = rand(N, N)
        A = A*A'
        ih = Arnoldi(A, rand(N, 1))
        for II = 1:N-2
            ih()
        end

        H = GetHessenberMatrix(ih)
        println("H Matrix:")
        display(H)
        Q = GetOrthogonalMatrix(ih)
        println("Q Matrix")
        display(Q)
        @info "Assertring Error bounds of decomposition and Tridiagonal Property of the H Matrix"
        @assert (H - tril(H, 1)).|>abs|>sum|>mean < 1e-10 "Not Tridiagonal Enough"
        @assert (A*Q[:, 1:end-1] - Q*H).|>abs|>sum|>mean < 1e-10 "Decomposition Error too big"
        return true
    end

    function Test4(N=10)
        @info "Testing Limited Memory Arnoldi Decomposition"
        A = rand(N, N)
        ih = Arnoldi(A, rand(N, 1); max_k=2)
        for II = 1:N-2
            ih()
        end
        H = GetHessenberMatrix(ih)
        println("Limited Memory H matrix: ")
        display(H)
        return true

    end

    function Test5(N=10)
        @info "Testing Limited Memory Arnoldi Decomposition with a complex matrix"
        A = rand(N, N) + rand(N, N)im
        ih = Arnoldi(A, rand(N, 1)im; max_k=2)
        for II = 1:N-2
            ih()
        end
        H = GetHessenberMatrix(ih)
        println("Limited Memory H matrix: ")
        display(H)
        Q = GetOrthogonalMatrix(ih)
        println("Limited Memory Q Matrix")
        display(Q)
        return true

    end

    @test Test1()
    @test Test2()
    @test Test3()
    @test Test4()
    @test Test5()
end


text:
"""
    setoptions!

This takes `options` arguments for setting DSP options.

# Arguments
- `options`: possible keys are `:param` (path to parameter file),
    `:is_stochastic` (`true` if stochastic model; otherwise, `false`),
    and `:solve_type` (algorithm type; see `Methods`).
"""
function setoptions!(options)
    check_dsp()
    for (optname, optval) in options
        if optname == :param
            readParamFile(dspenv, optval)
        elseif optname == :is_stochastic
            dspenv.is_stochastic = optval
        elseif optname == :solve_type
            if optval in instances(Methods)
                dspenv.solve_type = optval
            else
                @warn("solve_type $optval is not available.")
            end
        else
            @warn("Options $optname is not available.")
        end
    end
end

"""
    load_problem!

Load problem from StructJuMP

# Arguments
- `m`: StructJuMP model
"""
function load_problem!(m::SJ.StructuredModel)

    # set number of blocks (scenarios)
    dspenv.nblocks = SJ.num_scenarios(m)

    if dspenv.is_stochastic
        loadStochasticProblem!(m)
    else
        loadStructuredProblem!(m)
    end

    setBlocks()
end

"""
    loadStochasticProblem!

Load stochastic programming problem from StructJuMP

# Arguments
- `model`: StructJuMP model
"""
function loadStochasticProblem!(model::SJ.StructuredModel)

    nscen = dspenv.nblocks
    ncols1 = length(model.variables)
    nrows1 = length(dspenv.linConstrs[0])
    ncols2 = 0
    nrows2 = 0
    for (id, subm) in SJ.getchildren(model)
        ncols2 = length(subm.variables)
        nrows2 = length(dspenv.linConstrs[id])
        break
    end

    # set scenario indices for each MPI processor
    if dspenv.comm_size > 1
        ncols2 = MPI.Allreduce([ncols2], MPI.MAX, dspenv.comm)[1]
        nrows2 = MPI.Allreduce([nrows2], MPI.MAX, dspenv.comm)[1]
    end

    # set DSPProblem data
    dspenv.numCols[0] = ncols1
    dspenv.numRows[0] = nrows1
    dspenv.colVal[0] = Vector{Float64}(undef, ncols1)
    for (id, blk) in SJ.getchildren(model)
        dspenv.numCols[id] = ncols2
        dspenv.numRows[id] = nrows2
        dspenv.colVal[id] = Vector{Float64}(undef, ncols2)
    end

    setNumberOfScenarios(dspenv, nscen)
    setDimensions(dspenv, ncols1, nrows1, ncols2, nrows2)

    qc_supported = true

    # set problem data
    start, index, value, rlbd, rubd, obj, clbd, cubd, ctype, cname, nqrows, linnzcnt, quadnzcnt, rhs, sense, linstart, linind, linval, quadstart, quadrow, quadcol, quadval = get_model_data(model)
    if nqrows == 0
        loadFirstStage(dspenv, start, index, value, clbd, cubd, ctype, obj, rlbd, rubd)
    else
        if getVersionMajor(dspenv) >= 2
            loadQCQPFirstStage(dspenv, start, index, value, clbd, cubd, ctype, obj, C_NULL, C_NULL, C_NULL, 0, rlbd, rubd, nqrows, linnzcnt, quadnzcnt, rhs, sense, linstart, linind, linval, quadstart, quadrow, quadcol, quadval)
        else
            loadFirstStage(dspenv, start, index, value, clbd, cubd, ctype, obj, rlbd, rubd)
            qc_supported = false
        end
    end

    for (id, blk) in SJ.getchildren(model)
        probability = SJ.getprobability(model)[id]
        start, index, value, rlbd, rubd, obj, clbd, cubd, ctype, cname, nqrows, linnzcnt, quadnzcnt, rhs, sense, linstart, linind, linval, quadstart, quadrow, quadcol, quadval = get_model_data(blk, id)
        if nqrows == 0
            loadSecondStage(dspenv, id-1, probability, start, index, value, clbd, cubd, ctype, obj, rlbd, rubd)
        else
            if getVersionMajor(dspenv) >= 2
                loadQCQPSecondStage(dspenv, id-1, probability, start, index, value, clbd, cubd, ctype, obj, C_NULL, C_NULL, C_NULL, 0, rlbd, rubd, nqrows, linnzcnt, quadnzcnt, rhs, sense, linstart, linind, linval, quadstart, quadrow, quadcol, quadval)
            else
                loadSecondStage(dspenv, id-1, probability, start, index, value, clbd, cubd, ctype, obj, rlbd, rubd)
                qc_supported = false
            end
        end
    end

    if qc_supported == false
        @warn "QCQP is not supported with DSP version $(getVersion(dspenv)). The unsupported objective/constraints are ignored."
    end

    # Set DRO data
    loadDroData(dspenv, dspenv.dro)
end

"""
    loadStructuredProblem!

Load generic block-structured problem from StructJuMP

# Arguments
- `model`: StructJuMP model
"""
function loadStructuredProblem!(model::SJ.StructuredModel)

    ncols1 = length(model.variables)
    nrows1 = length(model.constraints)

    # set DSPProblem data
    dspenv.numCols[0] = ncols1
    dspenv.numRows[0] = nrows1
    dspenv.colVal[0] = Vector{Float64}(undef, ncols1)
    for (id, blk) in SJ.getchildren(model)
        ncols2 = length(blk.variables)
        nrows2 = length(dspenv.linConstrs[id])
        dspenv.numCols[id] = ncols2
        dspenv.numRows[id] = nrows2
        dspenv.colVal[id] = Vector{Float64}(undef, ncols2)
    end

    # load master
    start, index, value, rlbd, rubd, obj1, clbd1, cubd1, ctype1, cname, nqrows, linnzcnt, quadnzcnt, rhs, sense, linstart, linind, linval, quadstart, quadrow, quadcol, quadval = get_model_data(model)
    loadBlockProblem(dspenv, 0, ncols1, nrows1, start[end],
        start, index, value, clbd1, cubd1, ctype1, obj1, rlbd, rubd)

    # Check if the master block has quadratic constraints.
    has_qc = nqrows == 0 ? false : true

    # going over blocks
    for (id, blk) in SJ.getchildren(model)
        ncols2 = length(blk.variables)
        nrows2 = length(dspenv.linConstrs[id])
        start, index, value, rlbd, rubd, obj, clbd, cubd, ctype, cname, nqrows, linnzcnt, quadnzcnt, rhs, sense, linstart, linind, linval, quadstart, quadrow, quadcol, quadval = get_model_data(blk, id)
        loadBlockProblem(dspenv, id, ncols1 + ncols2, nrows2, start[end],
            start, index, value, [clbd1; clbd], [cubd1; cubd], [ctype1; ctype], [obj1; obj], rlbd, rubd)

        # Check if the sub-block has quadratic constraints.
        if !has_qc && nqrows > 0
            has_qc = true
        end
    end

    if has_qc
        @warn "Quadratic constraints are not supported for generic structured programs and thus will be ignored."
    end

    # Finalize loading blocks
    updateBlocks(dspenv)
end

function solve!()
    if dspenv.comm_size == 1
        if dspenv.solve_type == Dual
            if dspenv.is_stochastic
                solveDd(dspenv);
            else
                @warn("Dual decomposition is available for stochastic programming only.")
                return
            end
        elseif dspenv.solve_type == Benders
            if dspenv.is_stochastic
                solveBd(dspenv);
            else
                @warn("Benders decomposition is available for stochastic programming only.")
                return
            end
        elseif dspenv.solve_type == ExtensiveForm
            solveDe(dspenv);
        elseif dspenv.solve_type == DW
            solveDw(dspenv);
        else
            @error("Unexpected error")
            return
        end
    elseif dspenv.comm_size > 1
        if dspenv.solve_type == Dual
            if dspenv.is_stochastic
                solveDdMpi(dspenv);
            else
                @warn("Dual decomposition is available for stochastic programming only.")
                return
            end
        elseif dspenv.solve_type == Benders
            if dspenv.is_stochastic
                solveBdMpi(dspenv);
            else
                @warn("Benders decomposition is available for stochastic programming only.")
                return
            end
        elseif dspenv.solve_type == DW
            solveDwMpi(dspenv);
        elseif dspenv.solve_type == ExtensiveForm
            solveDe(dspenv);
        else
            @error("Unexpected error")
            return
        end
    end

    # solution status
    dspenv.status = getStatus(dspenv)
end

function post_solve!()
    if dspenv.status == 3998
        return
    end

    # get solution time
    dspenv.solve_time = getWallTime(dspenv)

    # primal and dual objective values
    dspenv.primVal = getPrimalBound(dspenv) * dspenv.objective_sense
    dspenv.dualVal = getDualBound(dspenv) * dspenv.objective_sense

    if dspenv.solve_type == Dual
        dspenv.rowVal = getDualSolution(dspenv)
    end

    primVal = dspenv.primVal
    if mysize() > 1
        primVal = MPI.bcast(dspenv.primVal, 0, dspenv.comm)
    end

    if abs(primVal) < 1.0e+20
        primsol = getSolution(dspenv)

        # parse solution to each block
        n_start = 1
        n_end = dspenv.numCols[0]
        dspenv.colVal[0] = primsol[n_start:n_end]
        n_start += dspenv.numCols[0]

        numBlockCols = getNumBlockCols()
        for i in 1:dspenv.nblocks
            n_end += numBlockCols[i]
            dspenv.colVal[i] = primsol[n_start:n_end]
            n_start += numBlockCols[i]
        end
    end
end


text:
# module timeevolution_quantum

export timeevolution_master, timeevolution_mcwf

using QuantumOptics
# using ..quantum

function timeevolution_master(system::MultimodeSystem, T, ρ₀; kwargs...)
    H = Hamiltonian(system)
    Hsparse = SparseOperator(H)
    J = Jump_operators(system)
    Jsparse = map(SparseOperator, J)
    return timeevolution.master(T, ρ₀, Hsparse, Jsparse; kwargs...)
end

function timeevolution_mcwf(system::MultimodeSystem, T, Ψ₀::Ket; kwargs...)
    H = Hamiltonian(system)
    Hsparse = SparseOperator(H)
    J = Jump_operators(system)
    Jsparse = map(SparseOperator, J)
    return timeevolution.mcwf(T, Ψ₀, Hsparse, Jsparse; kwargs...)
end

# end # module


text:
# Deprecate RodriguesVec => RotationVec
Base.@deprecate_binding RodriguesVec RotationVec true

# Deprecate Quat => UnitQuaternion
Base.@deprecate_binding Quat UnitQuaternion true

# Deprecate SPQuat => MRP
Base.@deprecate_binding SPQuat MRP true


text:
# This file reimplements gmsh/tutorial/t1.geo in Julia.

# For all the elementary explanations about the general philosphy of entities in
# Gmsh, see the comments in the .geo file. Comments here focus on the specifics
# of the Julia API.

# The API is entirely defined in the gmsh module
import gmsh

# Before using any functions in the Julia API, Gmsh must be initialized.
gmsh.initialize()

# By default Gmsh will not print out any messages: in order to output messages
# on the terminal, just set the standard Gmsh option "General.Terminal" (same
# format and meaning as in .geo files):
gmsh.option.setNumber("General.Terminal", 1)

# Next we add a new model named "t1" (if gmsh.model.add() is not called a new
# unnamed model will be created on the fly, if necessary):
gmsh.model.add("t1")

# The Julia API provides direct access to the internal CAD kernels. The
# built-in CAD kernel was used in t1.geo: the corresponding API functions have
# the "gmsh.model.geo" prefix. To create geometrical points with the built-in
# CAD kernel, one thus uses gmsh.model.geo.addPoint():
#
# - the first 3 arguments are the point coordinates (x, y, z)
#
# - the next (optional) argument is the target mesh size close to the point
#
# - the last (optional) argument is the point tag
lc = 1e-2
gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
gmsh.model.geo.addPoint(.1, 0,  0, lc, 2)
gmsh.model.geo.addPoint(.1, .3, 0, lc, 3)
gmsh.model.geo.addPoint(0, .3, 0, lc, 4)

# The API to create lines with the built-in kernel follows the same
# conventions: the first 2 arguments are point tags, the last (optional one)
# is the line tag.
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(3, 2, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)

# The philosophy to construct curve loops and surfaces is similar: the first
# argument is now a vector of integers.
gmsh.model.geo.addCurveLoop([4, 1, -2, 3], 1)
gmsh.model.geo.addPlaneSurface([1], 1)

# Physical groups are defined by providing the dimension of the group (0 for
# physical points, 1 for physical curves, 2 for physical surfaces and 3 for
# phsyical volumes) followed by a vector of entity tags. The last (optional)
# argument is the tag of the new group to create.
gmsh.model.addPhysicalGroup(0, [1, 2], 1)
gmsh.model.addPhysicalGroup(1, [1, 2], 2)
gmsh.model.addPhysicalGroup(2, [1], 6)

# Physical names are also defined by providing the dimension and tag of the
# entity.
gmsh.model.setPhysicalName(2, 6, "My surface")

# Before it can be meshed, the internal CAD representation must be synchronized
# with the Gmsh model, which will create the relevant Gmsh data structures. This
# is achieved by the gmsh.model.geo.synchronize() API call for the built-in CAD
# kernel. Synchronizations can be called at any time, but they involve a non
# trivial amount of processing; so while you could synchronize the internal CAD
# data after every CAD command, it is usually better to minimize the number of
# synchronization points.
gmsh.model.geo.synchronize()

# We can then generate a 2D mesh...
gmsh.model.mesh.generate(2)

# ... and save it to disk
gmsh.write("t1.msh")

# Remember that by default, if physical groups are defined, Gmsh will export in
# the output mesh file only those elements that belong to at least one physical
# group. To force Gmsh to save all elements, you can use
#
# gmsh.option.setNumber("Mesh.SaveAll", 1)

# This should be called at the end:
gmsh.finalize()

max_stars_repo_name: Christophe-Foyer/PyFEA

max_issues_repo_name: Christophe-Foyer/pyFEA

max_forks_repo_name: Christophe-Foyer/pyFEA


text:
struct EvenAsphere{T} <: AbstractRotationalSymmetricSurface{T}
    curvature::T
    conic::T
    c4::T
    c6::T
    c8::T
end

function EvenAsphere(curvature, conic, c4, c6, c8, aperture, n, t, id = nothing)
    OpticalSurface(
        EvenAsphere(
            curvature,
            conic,
            c4,
            c6,
            c8,
        ),
        aperture,
        n,
        t,
        id,
    )
end

function sag(radial_dist, s::EvenAsphere)
    z = s.curvature * radial_dist^2 / abs(1 + sqrt(1 - s.curvature^2 * radial_dist^2))
    z = z + s.c4 * radial_dist^4
    z = z + s.c6 * radial_dist^6
    return (z + s.c8 * radial_dist^8)
end

# modelled after us_itera.c example from Zemax user defined surface DLLs
function transfer_to_intersection_evenasphere!(ray, s::EvenAsphere{T}) where {T}
    t = 100.0
    x = ray.x
    y = ray.y
    z = ray.z
    loop = 0
    while abs(t) > 1e-10
        p2 = sqrt(x^2 + y^2)
        alpha = 1.0 - (1.0 + s.conic) * s.curvature^2 * p2
        if alpha < 0.0
            return ray_miss_error()
        end
        sg = sag(p2, s)
        dz = sg - z
        t = dz / ray.cz
        x += ray.cx * t
        y += ray.cy * t
        z += ray.cz * t
        loop += 1
        if loop > 1000
            return intersection_max_iterations_error(1000)
        end
    end

    ray.x = x
    ray.y = y
    ray.z = z
    return ray
end

function refract!(ray, m0, s::EvenAsphere, m1, wavelength)
    if m0 == m1
        return ray
    end

    x = ray.x
    y = ray.y
    z = ray.z
    X = ray.cx
    Y = ray.cy
    Z = ray.cz
    k = s.conic
    cv = s.curvature
    r2 = x^2 + y^2
    if r2 == 0
        ln = 0
        mn = 0
        nn = -1
    else
        alpha0 = 1.0 - (1.0 + k) * cv * cv * r2
        alpha = 1.0 - (1.0 + k) * cv * cv * r2
        if alpha < 0
            return ray_miss_error()
        end
        alpha = sqrt(alpha)
        mm0 =
            (cv / (1.0 + alpha)) *
            (2.0 + (cv * cv * r2 * (1.0 + k)) / (alpha * (1.0 + alpha)))
        mm =
            (cv / (1.0 + alpha)) *
            (2.0 + (cv * cv * r2 * (1.0 + k)) / (alpha * (1.0 + alpha)))
        mm += 4 * s.c4 * r2^(4 / 2 - 1)
        mm += 6 * s.c6 * r2^(6 / 2 - 1)
        mm += 8 * s.c8 * r2^(8 / 2 - 1)

        mm0 += 4 * s.c4 * r2
        mm0 += 6 * s.c6 * r2 * r2
        mm0 += 8 * s.c8 * r2 * r2 * r2

        mx = x * mm
        my = y * mm
        nn = -sqrt(1 / (1 + mx^2 + my^2))
        ln = -mx * nn
        mn = -my * nn

        nn0 = -sqrt(1 / (1 + (mx * mx) + (my * my)))
        ln0 = -mx * nn
        mn0 = -my * nn0
    end
    nr = refractive_index(m0, wavelength) / refractive_index(m1, wavelength)
    cosi = abs(X * ln + Y * mn + Z * nn)
    cosi2 = cosi * cosi
    if cosi2 > 1
        cosi2 = 1
    end
    rad = 1 - ((1 - cosi2) * nr^2)
    if rad < 0
        return total_internal_reflection_error()
    end
    cosr = sqrt(rad)
    gamma = nr * cosi - cosr
    X = nr * X + gamma * ln
    Y = nr * Y + gamma * mn
    Z = nr * Z + gamma * nn

    ray.cx = X
    ray.cy = Y
    ray.cz = Z
    return ray
end

function transfer_to_intersection!(ray, t::T, s::EvenAsphere{T}) where {T}
    # propagiere den Strahl erstmal zur sphärischen Grundfläche,
    # was wahrscheinlich ein ganz guter Startwert ist
    # TODO: missing error handling
    ts = transfer_to_intersection!(ray, t, Sphere(s.curvature))
    if iserror(ts)
        return ts
    else
        return transfer_to_intersection_evenasphere!(ts, s)
    end
end


text:
# ===================== Dependencies ======================
using ArgParse
using Logging
using Distributed
using SharedArrays
#TODO: move to JLD2
using JLD2, FileIO
using DelimitedFiles
using LinearAlgebra, GenericLinearAlgebra
using Combinatorics
using TOML          # used for input
using Printf
using ForwardDiff
using Query
using IntervalArithmetic, IntervalRootFinding
using Roots
using PaddedViews   # used to pad fft arrays
using FFTW          # used for convolutions
using NLsolve

using SeriesAcceleration
using Dispersions

# ======================= Includes ========================
include("$(@__DIR__)/Config.jl")
include("$(@__DIR__)/DataTypes.jl")
include("$(@__DIR__)/helpers.jl")
include("$(@__DIR__)/IO.jl")
include("$(@__DIR__)/GFTools.jl")
include("$(@__DIR__)/GFFit.jl")
include("$(@__DIR__)/ladderDGATools.jl")
include("$(@__DIR__)/lambdaCorrection.jl")
include("$(@__DIR__)/thermodynamics.jl")

# if (myid()==1) && (nprocs()==1)
#   println("activating procs")
#   addprocs(3)
#end
#
# ==================== Argument Parser ====================
s = ArgParseSettings()
@add_arg_table s begin
    "--config", "-c"
        help = "configuration file, default `config.toml`"
        arg_type = String
        default = "config.toml"
end

args = parse_args([], s)
io = stdout
metafmt(level::Logging.LogLevel, _module, group, id, file, line) = Logging.default_metafmt(level, nothing, group,id, nothing, nothing)
logger = ConsoleLogger(io, Logging.Info, meta_formatter=metafmt, show_limited=true, right_justify=0)
global_logger(logger)


text:
## initialize

using BenchmarkTools

###

img1, imgw1, flow1 = gen_init(flow_args=(20,60));
flow_est1, source_reg1, figs1 = sparse_pflap(img1, imgw1)

showflow(flow_est1, figtitle="estim")
showflow(flow1, figtitle="orig")

### porovnani o kolik je single_lap_at_points rychlejsi nez single_lap podle velikosti posunu

img, imgw, flow = gen_init(flow_args=(4,100));
showflow(flow, figtitle="orig")

# vars
fhs = 5
window_size = [11,11]

# classic
clasic_estim = single_lap(img, imgw, fhs, window_size)
showflow(clasic_estim, figtitle="clasic_estim")

# new
spacing = 35
point_count = 30
mask = parent(padarray(trues(size(target).-(2*fhs, 2*fhs)), Fill(false, (fhs, fhs), (fhs, fhs))))
inds = find_edge_points(target, spacing=spacing, number=point_count, mask=mask)
points = LAP_julia.inds_to_points(inds)
new_estim = single_lap_at_points(img, imgw, fhs, window_size, 3, points)
showflow(new_estim, disp_type=:sparse, figtitle="new_estim")

full_new_estim = interpolate_flow(new_estim, inds)
showflow(full_new_estim, figtitle="full_new_estim")
addpoints(inds)

imgshow(img); addpoints(inds)

# # classic
# function lap(img, imgw, fhs, window_size)
#     classic_estim = single_lap(img, imgw, fhs, window_size)
#     LAP_julia.inpaint_nans!(classic_estim)
#     LAP_julia.smooth_with_gaussian!(classic_estim, window_size)
#     return classic_estim
# end

# new
function new_alg(img, imgw, fhs, window_size)
    spacing = 35
    point_count = 35
    mask = parent(padarray(trues(size(img).-(2*fhs, 2*fhs)), Fill(false, (fhs, fhs), (fhs, fhs))))
    inds = find_edge_points(img, spacing=spacing, number=point_count, mask=mask)
    points = LAP_julia.inds_to_points(inds)
    new_estim = single_lap_at_points(img, imgw, fhs, window_size, 3, points)
    full_new_estim = interpolate_flow(new_estim, inds)
    return full_new_estim
end

function get_erroring_params_from_new(maxd)
    # vars
    runs = 0
    while true
        img, imgw, flow = gen_init(flow_args=(maxd, 100));
        fhs = maxd+2
        window_size = [fhs*2+1,fhs*2+1]

        output = new_alg(img, imgw, fhs, window_size)
        if length(output) < 200
            return output
        end
        runs = 1 + runs
        println(runs)
    end
end

out = get_erroring_params_from_new(5)

function compare_speed_dep_on_max_displacement(maxd)

    # vars
    img, imgw, flow = gen_init(flow_args=(maxd, 100));
    fhs = maxd+2
    window_size = [fhs*2+1,fhs*2+1]

    new_flow = new_alg(img, imgw, fhs, window_size)
    classic_flow = lap(img, imgw, fhs, window_size)

    new_bench = @benchmark $new_alg($img, $imgw, $fhs, $window_size)
    classic_bench = @benchmark $lap($img, $imgw, $fhs, $window_size)

    classic_speed = median(classic_bench.times)
    classic_mse = mse(classic_flow, flow)

    new_speed = median(new_bench.times)
    new_mse = mse(new_flow, flow)

    new_mse_pix = mse(warp_img(imgw, -real(new_flow), -imag(new_flow)), flow)
    classic_mse_pix = mse(warp_img(imgw, -real(classic_flow), -imag(classic_flow)), flow)

    println("CLASSIC: speed: ", classic_speed, " mse: ", classic_mse, "mse pix :", classic_mse_pix)
    println("NEW: speed: ", new_speed, " mse: ", new_mse, "mse pix :", new_mse_pix)

    println("mse pix quility increase: ", classic_mse_pix/new_mse_pix)

    speedup = classic_speed / new_speed
    quality_compare = classic_mse / new_mse

    println("speedup depending on points: ", speedup)
    println("quality increase depending on points: ", quality_compare)
    return speedup, quality_compare
end

function run_multiple(fun, run_count; args=[])
    outs = Any[]
    for k in 1:run_count
        println()
        out = fun(args...)
        push!(outs, out)
    end
    # make sure the out is a tuple or array
    if typeof(outs[1]) <: Union{Tuple, Array}
        output_size = length(outs[1])
    else
        output_size = 1
        outs = map(x -> (x), outs)
    end
    # average them
    ret = (LAP_julia.mean([outs[l][k] for l in 1:output_size]) for k in 1:output_size)
    return outs
end

using LinearAlgebra

function end_point_error(u, u_est)
    as, df = LAP_julia.mse.(u, u_est)
end
u = flow

function e_med(u, u_est)
    retrun median((LAP_julia.vec_len.(u) .- LAP_julia.vec_len.(u_est))[:])
end

(real(flow), imag(flow)) .- (real(u_est), imag(u_est)) .^2

as, df = (real(flow), imag(flow)) .- (real(u_est), imag(u_est))
as .^2

function angular_error(u, u_est)
    return 1/cos(
        (1 + dot(u[:], u_est[:])) /
        (sqrt(1 + dot(u[:], u[:])) * sqrt(1 + dot(u_est[:], u_est[:]))))
end

LAP_julia.vec_len(angular_error(flow, u_est))
end_point_error(flow, u_est)

out1, out2 = compare_speed_dep_on_max_displacement(4)

showflow(out1)

run_multiple(compare_speed_dep_on_max_displacement, 4, args=[15])

ar = [[1212, 11], [-1212, 11]]
outs = ar
output_size = 2
[LAP_julia.mean([outs[l][k] for l in 1:output_size]) for k in 1:output_size]

displacements = [2:3:38]
run

### single lap at point showcase:

# vars
maxd = 5
img, imgw, flow = gen_init(flow_args=(maxd, 100));
fhs = maxd+2
window_size = [fhs*2+1,fhs*2+1]
showflow(flow, figtitle="orig")

single_points_full_flow = new_alg(img, imgw, fhs, window_size)
showflow(single_points_full_flow, figtitle="estim using single points")

single_full_flow = lap(img, imgw, fhs, window_size)
showflow(single_full_flow, figtitle="estim using classic single")

# show dif
showflow(single_full_flow .- orig, figtitle="classic single - orig")
showflow(single_points_full_flow .- orig, figtitle="single points - orig")


text:
function leftover(d::DiscreteUnivariateDistribution, q)
    lower = max(minimum(d), 0)
    upper = min(maximum(d), q)
    return sum((q - y) * pdf(d, y) for y in lower:upper; init = 0.0)
end

function leftover(d::ContinuousUnivariateDistribution, q)
    I, _ = quadgk(y -> (q - y) * pdf(d, y), minimum(d), q)
    return I
end

# For normal distributions, we use σ * (ϕ(z) + z*Φ(z)) where z = (q - μ) / σ.
I(z) = pdf(Normal(), z) + z * cdf(Normal(), z)
function leftover(d::Normal, q)
    z = (q - mean(d)) / std(d)
    return std(d) * I(z)
end


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

# these should both use Geom.point by default, see #1062
plot(x=[1,2,3], y=[4,5,6])
plot(layer(x=[1,2,3], y=[4,5,6]))


text:
@testset "map" begin
    m = sprand(5, 5, 0.25)
    n = GBMatrix(m)
    @test map(UnaryOps.LOG, n)[1,1] == map(log, m)[1,1]
    o = map!(>, GBMatrix{Bool}(5, 5),  0.1, n)
    @test o[1,4] == (0.1 > m[1,4])
    @test map(second, n, 1.5)[1,1] == 1.5
    @test (n .* 10)[1,1] == n[1,1] * 10
    # Julia will map over the entire array, rather than just nnz.
    # so just test [1,1]
    @test map((x) -> 1.5, n)[1,1] == map((x) -> 1.5, m)[1,1]
end


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

if isempty(ARGS) || ARGS[1] !== "0"
Sys.__init_build()
# Prevent this from being put into the Main namespace
@eval Module() begin
if !isdefined(Base, :uv_eventloop)
    Base.reinit_stdio()
end
Base.include(@__MODULE__, joinpath(Sys.BINDIR, "..", "share", "julia", "test", "testhelpers", "FakePTYs.jl"))
import .FakePTYs: open_fake_pty

CTRL_C = '\x03'
UP_ARROW = "\e[A"
DOWN_ARROW = "\e[B"

hardcoded_precompile_statements = """
@assert precompile(Tuple{typeof(Base.stale_cachefile), String, String})
@assert precompile(Tuple{typeof(push!), Set{Module}, Module})
@assert precompile(Tuple{typeof(push!), Set{Method}, Method})
@assert precompile(Tuple{typeof(push!), Set{Base.PkgId}, Base.PkgId})
@assert precompile(Tuple{typeof(getindex),  Dict{Base.PkgId,String}, Base.PkgId})
@assert precompile(Tuple{typeof(setindex!), Dict{String,Base.PkgId}, Base.PkgId, String})
@assert precompile(Tuple{typeof(get!), Type{Vector{Function}}, Dict{Base.PkgId,Vector{Function}}, Base.PkgId})
@assert precompile(Tuple{typeof(isassigned), Core.SimpleVector, Int})
@assert precompile(Tuple{typeof(pushfirst!), Vector{Any}, Function})
@assert precompile(Tuple{typeof(Base.parse_cache_header), String})
"""

precompile_script = """
2+2
print("")
@time 1+1
; pwd
? reinterpret
using Ra\t$CTRL_C
\\alpha\t$CTRL_C
\e[200~paste here ;)\e[201~"$CTRL_C
$UP_ARROW$DOWN_ARROW$CTRL_C
123\b\b\b$CTRL_C
\b\b$CTRL_C
f(x) = x03
f(1,2)
[][1]
cd("complet_path\t\t$CTRL_C
"""

julia_exepath() = joinpath(Sys.BINDIR, Base.julia_exename())

have_repl =  haskey(Base.loaded_modules,
                    Base.PkgId(Base.UUID("3fa0cd96-eef1-5676-8a61-b3b8758bbffb"), "REPL"))
if have_repl
    hardcoded_precompile_statements *= """
    @assert precompile(Tuple{typeof(getproperty), REPL.REPLBackend, Symbol})
    """
end

Distributed = get(Base.loaded_modules,
          Base.PkgId(Base.UUID("8ba89e20-285c-5b6f-9357-94700520ee1b"), "Distributed"),
          nothing)
if Distributed !== nothing
    precompile_script *= """
    using Distributed
    addprocs(2)
    pmap(x->iseven(x) ? 1 : 0, 1:4)
    @distributed (+) for i = 1:100 Int(rand(Bool)) end
    """
end

Artifacts = get(Base.loaded_modules,
          Base.PkgId(Base.UUID("56f22d72-fd6d-98f1-02f0-08ddc0907c33"), "Artifacts"),
          nothing)
if Artifacts !== nothing
    precompile_script *= """
    using Artifacts, Base.BinaryPlatforms
    artifacts_toml = abspath($(repr(joinpath(Sys.STDLIB, "Artifacts", "test", "Artifacts.toml"))))
    cd(() -> @artifact_str("c_simple"), dirname(artifacts_toml))
    artifacts = Artifacts.load_artifacts_toml(artifacts_toml)
    platforms = [Artifacts.unpack_platform(e, "c_simple", artifacts_toml) for e in artifacts["c_simple"]]
    best_platform = select_platform(Dict(p => triplet(p) for p in platforms))
    """
end

Pkg = get(Base.loaded_modules,
          Base.PkgId(Base.UUID("44cfe95a-1eb2-52ea-b672-e2afdf69b78f"), "Pkg"),
          nothing)

if Pkg !== nothing
    precompile_script *= Pkg.precompile_script
end

FileWatching = get(Base.loaded_modules,
          Base.PkgId(Base.UUID("7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"), "FileWatching"),
          nothing)
if FileWatching !== nothing
    hardcoded_precompile_statements *= """
    @assert precompile(Tuple{typeof(FileWatching.watch_file), String, Float64})
    @assert precompile(Tuple{typeof(FileWatching.watch_file), String, Int})
    """
end

function generate_precompile_statements()
    start_time = time_ns()
    debug_output = devnull # or stdout

    # Precompile a package
    mktempdir() do prec_path
        push!(DEPOT_PATH, prec_path)
        push!(LOAD_PATH, prec_path)
        pkgname = "__PackagePrecompilationStatementModule"
        mkpath(joinpath(prec_path, pkgname, "src"))
        write(joinpath(prec_path, pkgname, "src", "$pkgname.jl"),
              """
              module $pkgname
              end
              """)
        @eval using __PackagePrecompilationStatementModule
        empty!(LOAD_PATH)
        empty!(DEPOT_PATH)
    end

    mktemp() do precompile_file, precompile_file_h
        # Run a repl process and replay our script
        pts, ptm = open_fake_pty()
        blackhole = Sys.isunix() ? "/dev/null" : "nul"
        if have_repl
            cmdargs = ```--color=yes
                      -e 'import REPL; REPL.Terminals.is_precompiling[] = true'
                      ```
        else
            cmdargs = `-e nothing`
        end
        p = withenv("JULIA_HISTORY" => blackhole,
                    "JULIA_PROJECT" => nothing, # remove from environment
                    "JULIA_LOAD_PATH" => Sys.iswindows() ? "@;@stdlib" : "@:@stdlib",
                    "TERM" => "") do
            sysimg = Base.unsafe_string(Base.JLOptions().image_file)
            run(```$(julia_exepath()) -O0 --trace-compile=$precompile_file --sysimage $sysimg
                   --cpu-target=native --startup-file=no --color=yes
                   -e 'import REPL; REPL.Terminals.is_precompiling[] = true'
                   -i $cmdargs```,
                   pts, pts, pts; wait=false)
        end
        Base.close_stdio(pts)
        # Prepare a background process to copy output from process until `pts` is closed
        output_copy = Base.BufferStream()
        tee = @async try
            while !eof(ptm)
                l = readavailable(ptm)
                write(debug_output, l)
                Sys.iswindows() && (sleep(0.1); yield(); yield()) # workaround hang - probably a libuv issue?
                write(output_copy, l)
            end
            close(output_copy)
            close(ptm)
        catch ex
            close(output_copy)
            close(ptm)
            if !(ex isa Base.IOError && ex.code == Base.UV_EIO)
                rethrow() # ignore EIO on ptm after pts dies
            end
        end
        # wait for the definitive prompt before start writing to the TTY
        readuntil(output_copy, "julia>")
        sleep(0.1)
        readavailable(output_copy)
        # Input our script
        if have_repl
            precompile_lines = split(precompile_script, '\n'; keepempty=false)
            curr = 0
            for l in precompile_lines
                sleep(0.1)
                curr += 1
                print("\rGenerating precompile statements... $curr/$(length(precompile_lines))")
                # consume any other output
                bytesavailable(output_copy) > 0 && readavailable(output_copy)
                # push our input
                write(debug_output, "\n#### inputting statement: ####\n$(repr(l))\n####\n")
                write(ptm, l, "\n")
                readuntil(output_copy, "\n")
                # wait for the next prompt-like to appear
                # NOTE: this is rather inaccurate because the Pkg REPL mode is a special flower
                readuntil(output_copy, "\n")
                readuntil(output_copy, "> ")
            end
            println()
        end
        write(ptm, "exit()\n")
        wait(tee)
        success(p) || Base.pipeline_error(p)
        close(ptm)
        write(debug_output, "\n#### FINISHED ####\n")

        # Extract the precompile statements from the precompile file
        statements = Set{String}()
        for statement in eachline(precompile_file_h)
            # Main should be completely clean
            occursin("Main.", statement) && continue
            push!(statements, statement)
        end

        for statement in split(hardcoded_precompile_statements, '\n')
            push!(statements, statement)
        end

        # Create a staging area where all the loaded packages are available
        PrecompileStagingArea = Module()
        for (_pkgid, _mod) in Base.loaded_modules
            if !(_pkgid.name in ("Main", "Core", "Base"))
                eval(PrecompileStagingArea, :(const $(Symbol(_mod)) = $_mod))
            end
        end

        # Execute the collected precompile statements
        n_succeeded = 0
        include_time = @elapsed for statement in sort(collect(statements))
            # println(statement)
            try
                Base.include_string(PrecompileStagingArea, statement)
                n_succeeded += 1
                print("\rExecuting precompile statements... $n_succeeded/$(length(statements))")
            catch
                # See #28808
                # @error "Failed to precompile $statement"
            end
        end
        println()
        if have_repl
            # Seems like a reasonable number right now, adjust as needed
            # comment out if debugging script
            @assert n_succeeded > 1200
        end

        tot_time = time_ns() - start_time
        include_time *= 1e9
        gen_time = tot_time - include_time
        println("Precompilation complete. Summary:")
        print("Total ─────── "); Base.time_print(tot_time); println()
        print("Generation ── "); Base.time_print(gen_time);     print(" "); show(IOContext(stdout, :compact=>true), gen_time / tot_time * 100); println("%")
        print("Execution ─── "); Base.time_print(include_time); print(" "); show(IOContext(stdout, :compact=>true), include_time / tot_time * 100); println("%")
    end

    return
end

generate_precompile_statements()

# As a last step in system image generation,
# remove some references to build time environment for a more reproducible build.
@eval Base PROGRAM_FILE = ""
@eval Sys begin
    BINDIR = ""
    STDLIB = ""
end
empty!(Base.ARGS)
empty!(Core.ARGS)

end # @eval
end

max_stars_repo_name: mlxd/julia

max_issues_repo_name: mlxd/julia

max_forks_repo_name: mlxd/julia


text:
export gauss_quad_linear,gauss_points
function gauss_quad_linear(N,a,b)
    N=N-1;
    N1=N+1;
    N2=N+2;
    xu=linspace(-1,1,N1)';
    # Initial guess
    y=cos.((2*(0:N)'+1)*pi/(2*N+2))+(0.27/N1)*sin.(pi*xu*N/N2);
    # Legendre-Gauss Vandermonde Matrix
    L=zeros(N1,N2);
    # Derivative of LGVM
    Lp=zeros(N1,N2);
    # Compute the zeros of the N+1 Legendre Polynomial
    # using the recursion relation and the Newton-Raphson method
    y0=2;
    # Iterate until new points are uniformly within epsilon of old points
    while (maximum(abs.(y-y0))>eps(1.0))
        L[:,1]=1;
        Lp[:,1]=0;
        L[:,2]=y;
        #Lp[:,2]=1;
        for k=2:N1
            L[:,k+1]=( (2*k-1)*y'.*L[:,k]-(k-1)*L[:,k-1] )/k;
        end
        Lp=(N2)*( L[:,N1]-y'.*L[:,N2] )./(1-y.^2)';
        y0=y;
        y=y0'-L[:,N2]./Lp;
        y=y'
    end
    y=y'
    # Linear map from[-1,1] to [a,b]
    x=(a*(1-y)+b*(1+y))/2;
    # Compute the weights
    w=(b-a)./((1-y.^2).*Lp.^2)*(N2/N1)^2;
    return x, w
end

function gauss_points(N)
    x,w = gauss_quad_linear(N,-1,1)
    gpt = zeros(N*N,2)
    gpw = zeros(N*N,1)
    for i = 1:N
        u = x[i]
        wu = w[i]
        for j = 1:N
            v = x[j]
            wv = w[j]
            gpt[(i-1)*N + j,1] = u
            gpt[(i-1)*N + j,2] = v
            gpw[(i-1)*N + j] = wu * wv
        end
    end
    return gpt,gpw
end


text:
function cpuv(f::Func, ts...)
  ar, v = infer(f, ts...)
  v = insert_domains(striptypes(v), ar)
  v = tolambda(v, length(ts))
end

cpu(a...) = cpuv(a...) |> syntax |> eval

max_stars_repo_name: JuliaTagBot/Tokamak

max_issues_repo_name: JuliaTagBot/Tokamak

max_forks_repo_name: FluxML/Tokamak


text:
function ripetitore(num::Int, rep::Int, text::AbstractString) :: Any
    """This function makes small easy drawings"""
    for i in 1::num
        println("  "^(i+1), text^rep,  "  "^(i+1), text^rep)
    end
end

#ripetitore(15, 2, "<3 ")

function reverser(str::AbstractString)
    """This easy function reverses the strings"""
    collection = collect(str)
    for char in reverse(collection)
        if char!== "m"
            println(char)
        else
            println("Mausa!")
        end
    end
end

#reverser("pillmmmmaaa")

function dna_complementary(sequence_dna::AbstractString)
    """Simple function to genrate reverese complementary sequence"""
    bases = Dict("a" :: "t", "c" :: "g", "g"::"c", "t"::"a")
    complementary_seq = []
    sequence_low = lowercase(sequence_dna)
    for el in sequence_low
         if el in keys(bases)
            append!(complementary_seq, bases[el])
         else
            println("The base you put is not right, please re-enter the correct sequence and run it again")
            break
         end
    end
    return uppercase(join(complementary_seq, ""))
end

#dna_complementary("gccttgggaatttAAAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGGGGGG")

function traslatore_dna_in_rna(string_dna :: AbstractString) ::AbstractString
    """Simple function to generate RNA starting from DNA sequences"""
    bases = Dict("a" :: "t", "c" :: "g", "g"::"c", "t"::"T")
    println("You provided a sequence of ", length(string_dna))
    sequence_low_rna = lowercase(string_dna)
    complementary_seq = []
    for el in sequence_low_rna
        append!(complementary_seq, bases[el])
    end
    return uppercase(join(complementary_seq, ""))
end

#traslatore_dna_in_rna("gccttgggaattttgcAAAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTTT")

function translator_from_rna(rna_seq :: AbstractString)
    """This function provided a RNA sequence will generate the protein sequence starting from the first encountered AUG"""
    traslation = []
    table = Dict(
        "AUA"::"I", "AUC"::"I", "AUU"::"I", "AUG"::"M",
        "ACA"::"U", "ACC"::"T", "ACG"::"T", "ACU"::"T",
        "AAC"::"N", "AAU"::"N", "AAA"::"K", "AAG"::"K",
        "AGC"::"S", "AGU"::"S", "AGA"::"R", "AGG"::"R",
        "CUA"::"L", "CUC"::"L", "CUG"::"L", "CUU"::"L",
        "CCA"::"P", "CCC"::"P", "CCG"::"P", "CCU"::"P",
        "CAC"::"H", "CAU"::"H", "CAA"::"Q", "CAG"::"Q",
        "CGA"::"R", "CGC"::"R", "CGG"::"R", "CGU"::"R",
        "GUA"::"V", "GUC"::"V", "GUG"::"V", "GUU"::"V",
        "GCA"::"A", "GCC"::"A", "GCG"::"A", "GCU"::"A",
        "GAC"::"D", "GAU"::"D", "GAA"::"E", "GAG"::"E",
        "GGA"::"G", "GGC"::"G", "GGG"::"G", "GGU"::"G",
        "UCA"::"S", "UCC"::"S", "UCG"::"S", "UCU"::"S",
        "UUC"::"F", "UUU"::"F", "UUA"::"L", "UUG"::"L",
        "UAC"::"Y", "UAU"::"Y", "UAA"::"_", "UAG"::"_",
        "UGC"::"C", "UGU"::"C", "UGA"::"_", "UGG"::"W" )

    start_translation = findfirst("AUG", uppercase(rna_seq))[2]
    triplets = 2
    for (num, sequence) in enumerate(start_translation:length(rna_seq))
        if num % 3 == 0 && num < length(rna_seq)
            triplents = rna_seq[sequence:sequence + triplets]
            println(triplents)
            append!(traslation, table[triplents])
        end
    end
    return join(traslation, "")
end

#translator_from_rna("UUUAUGGCCAUGGCGCCCAGAACUGAGAUCAAUAGUACCCGUAUUAACGGGUGAUUCACAUGUAAA")

function translator_from_rna_v2((ARGS[1]) :: AbstractString)
    """This function provided a RNA sequence will generate the protein sequence starting from the first encountered AUG.
    However, in contrast to the previous function this one takes the function argument from the command line"""

    #This will throw an error in case the lenght of the argument is not correct
    length(ARGS) == 1 || error("I need one command line argument")
    traslation = []
    table = Dict(
        "AUA"::"I", "AUC"::"I", "AUU"::"I", "AUG"::"M",
        "ACA"::"U", "ACC"::"T", "ACG"::"T", "ACU"::"T",
        "AAC"::"N", "AAU"::"N", "AAA"::"K", "AAG"::"K",
        "AGC"::"S", "AGU"::"S", "AGA"::"R", "AGG"::"R",
        "CUA"::"L", "CUC"::"L", "CUG"::"L", "CUU"::"L",
        "CCA"::"P", "CCC"::"P", "CCG"::"P", "CCU"::"P",
        "CAC"::"H", "CAU"::"H", "CAA"::"Q", "CAG"::"Q",
        "CGA"::"R", "CGC"::"R", "CGG"::"R", "CGU"::"R",
        "GUA"::"V", "GUC"::"V", "GUG"::"V", "GUU"::"V",
        "GCA"::"A", "GCC"::"A", "GCG"::"A", "GCU"::"A",
        "GAC"::"D", "GAU"::"D", "GAA"::"E", "GAG"::"E",
        "GGA"::"G", "GGC"::"G", "GGG"::"G", "GGU"::"G",
        "UCA"::"S", "UCC"::"S", "UCG"::"S", "UCU"::"S",
        "UUC"::"F", "UUU"::"F", "UUA"::"L", "UUG"::"L",
        "UAC"::"Y", "UAU"::"Y", "UAA"::"_", "UAG"::"_",
        "UGC"::"C", "UGU"::"C", "UGA"::"_", "UGG"::"W" )

    start_translation = findfirst("AUG", uppercase(rna_seq))[2]
    triplets = 2
    for (num, sequence) in enumerate(start_translation:length(rna_seq))
        if num % 3 == 0 && num < length(rna_seq)
            triplents = rna_seq[sequence:sequence + triplets]
            println(triplents)
            append!(traslation, table[triplents])
        end
    end
    return join(traslation, "")
end


text:
using ProjectedGPE
ecut = 40
basis = "Hermite"

ωx = 1
ωy = π
ωz = 5
Ω = [ωx;ωy;ωz]

#check 1D
cinfo = makecinfo(ecut,Ω,basis)
#@step makecinfo("Hermite",ecut,[ωx,ωz])

#test 1D
Ω=[ωx]
M=[5]
x,wx,Tx=makealltrans(M,4,Ω,basis)

#test 2D
Ω=[ωx, ωy]
M=[5,10]
x,wx,Tx,y,wy,Ty=makealltrans(M,4,Ω,basis)

#test 3D
Ω=[ωx, ωy, ωz]
M=[25,17,33]
x,wx,Tx,y,wy,Ty,z,wz,Tz=makealltrans(M,4,Ω,basis)


text:
# The original license on the t-tests from HypothesisTests.jl
# is as follows. We're including it here because we're mostly copying the
# documentation of the functions.
#
# kolmogorov_smirnov.jl
# Kolmogorov–Smirnov
#
# Copyright (C) 2014   Christoph Sawade
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

import HypothesisTests.ExactOneSampleKSTest
import HypothesisTests.ApproximateTwoSampleKSTest

"""
    ExactOneSampleKSTest(uv::AbstractUncertainValue,
        d::UnivariateDistribution, n::Int = 1000) -> ExactOneSampleKSTest

Perform a one-sample exact Kolmogorov–Smirnov test of the null hypothesis that
a draw of `n` realisations of the uncertain value `uv` comes from the
distribution `d` against the alternative hypothesis that the sample is not
drawn from `d`.
"""
function ExactOneSampleKSTest(uv::AbstractUncertainValue,
        d::UnivariateDistribution, n::Int = 1000)
    x = resample(uv, n)
    ExactOneSampleKSTest(x, d)
end

"""
    ExactOneSampleKSTestPooled(ud::UncertainDataset,
        d::UnivariateDistribution, n::Int = 1000) -> ExactOneSampleKSTest

First, draw `n` realisations of each uncertain value in `ud` and pool them
together. Then perform a one-sample exact Kolmogorov–Smirnov test of the null
hypothesis that the pooled values comes from the distribution `d` against the
alternative hypothesis that the sample is not drawn from `d`.
"""
function ExactOneSampleKSTestPooled(ud::UncertainDataset,
        d::UnivariateDistribution, n::Int = 1000)
    x = vcat(resample(ud, n)...,)
    ExactOneSampleKSTest(x, d)
end

"""
    ExactOneSampleKSTestElementWise(ud::UncertainDataset,
        d::UnivariateDistribution, n::Int = 1000) -> Vector{ExactOneSampleKSTest}

First, draw `n` realisations of each uncertain value in `ud`, keeping one
pool of values for each uncertain value.

Then, perform an element-wise (pool-wise) one-sample exact Kolmogorov–Smirnov
test of the null hypothesis that each value pool comes from the distribution `d`
against the alternative hypothesis that the sample is not drawn from `d`.
"""
function ExactOneSampleKSTestElementWise(ud::UncertainDataset,
        d::UnivariateDistribution, n::Int = 1000)
    [ExactOneSampleKSTest(resample(ud[i], n), d) for i = 1:length(ud)]
end

"""
    ApproximateTwoSampleKSTest(d1::AbstractUncertainValue,
        d2::AbstractUncertainValue, n::Int = 1000) -> ApproximateTwoSampleKSTest

Perform an asymptotic two-sample Kolmogorov–Smirnov-test of the null hypothesis
that the distribution furnishing the uncertain value `d1` represent the same
distribution as the distribution furnishing the uncertain value `d2`
against the alternative hypothesis that the furnishing distributions are
different.
"""
function ApproximateTwoSampleKSTest(d1::AbstractUncertainValue,
                                    d2::AbstractUncertainValue,
                                    n::Int = 1000)
    x = resample(d1, n)
    y = resample(d2, n)
    ApproximateTwoSampleKSTest(x, y)
end

"""
    ApproximateTwoSampleKSTestPooled(d1::UncertainDataset,
        d2::UncertainDataset, n::Int = 1000) -> ApproximateTwoSampleKSTest

First, draw `n` realisations of each uncertain value in `d1`, then separately
draw `n` realisations of each uncertain value in `d2`. Then, pool all
realisations for `d1` together and all realisations of `d2` together.

On the pooled realisations, perform an asymptotic two-sample
Kolmogorov–Smirnov-test of the null hypothesis that the distribution furnishing
the `d1` value pool represents the same distribution as the distribution
furnishing the `d2` value pool, against the alternative hypothesis that the
furnishing distributions are different.
"""
function ApproximateTwoSampleKSTestPooled(d1::UncertainDataset,
                                        d2::UncertainDataset,
                                        n::Int = 1000)
    x = vcat(resample(d1, n)...,)
    y = vcat(resample(d2, n)...,)
    ApproximateTwoSampleKSTest(x, y)
end

"""
    ApproximateTwoSampleKSTestElementWise(d1::UncertainDataset,
        d2::UncertainDataset, n::Int = 1000) -> Vector{ApproximateTwoSampleKSTest}

Assuming `d1` and `d2` contain the same number of uncertain observations,
draw `n` realisations of each uncertain value in `d1`, then separately and
separately draw `n` realisations of each uncertain value in `d2`.

Then, perform an asymptotic two-sample Kolmogorov–Smirnov-test of
the null hypothesis that the uncertain values in `d1` and `d2` come from
the same distribution against the alternative hypothesis that the (element-wise)
values in  `d1` and `d2` come from different distributions.

The test is performed pairwise, i.e. ApproximateTwoSampleKSTest(d1[i], d2[i])
with `n` draws for the ``i``-ith pair of uncertain values.
"""
function ApproximateTwoSampleKSTestElementWise(d1::UncertainDataset,
                                                d2::UncertainDataset,
                                                n::Int = 1000)
    N = length(d1)

    [ApproximateTwoSampleKSTest(resample(d1[i], n), resample(d2[i], n)) for i = 1:N]
end

export
ExactOneSampleKSTest,
ExactOneSampleKSTestPooled,
ExactOneSampleKSTestElementWise,

ApproximateTwoSampleKSTest,
ApproximateTwoSampleKSTestPooled,
ApproximateTwoSampleKSTestElementWise