text:
using Plots

# Poisson equation: u_xx = b
# Use centered diff: u_xx = 1/dx^2 (u_i+1 - 2u_i + u_i-1)
true_u(x) = - (1/4π^2) * sin.(2π * x)
u_xx(x) = sin.(2π * x)

Δx = 0.01
x = Δx:Δx:1-Δx # Solve only for the interior: the endpoints are known to be zero!
N = length(x)
B = u_xx(x)
A = zeros(N,N)

# Generate the centered difference matrix manually
for i in 1:N, j in 1:N

  if abs(i-j)<=1
    (A[i,j]+=1)
  end
  if i==j
    (A[i,j]-=3)
  end

end
A = A/(Δx^2)
U = A\B

display(plot([0;x;1],[0;U;0],label="U"))
display(plot!([0; x; 1], [0; true_u(x); 0], label="Analytical"))


text:
# Quelques fonctions utilitaires pour les manipulations de solvers PL ou MIP

# Classe abstraite non utilisé
# abstract type AbstractLpSolver end

"""
Création d'un modèle de solver externe (:cplex, :glpk, ...)

Args:
    solver : symbole :cplex, :glpk, :clp, :cbc

        - :cbc : permet de faire du MIP
        - :clp : resteint aux problème LP mais plus efficace de cbc pour les relaxations.

    mode : (:lp ou :mip) utile pour choisir entre :clp ou :cbc
    log_level (expérimental) : entier 0 pour silence, sinon 1, 2, ...

Autre solver à supporter plus tard Gurobi (concurrent de cplex)

Attention:
    DEPENDANCE EXTERNE : module Args (pour Args.get("external_mip_solver"))
"""
function new_model(;solver = :auto,
                     mode = :lp,
                     log_level = 0, # EXPLOITATION PARTIELLE
                     )
    if !(mode in [:mip, :lp])
        error("ERROR: unknown mode $(mode). Should be :cplex, :glpk, ...")
    end
    if solver == :auto
        solver = Args.get("external_mip_solver")
    end
    if solver in [:clp, :cbc]
        # solver = (mode == :mip ? :cbc : :clp )
        if mode==:mip
            solver = :cbc
        else
            solver = :clp
        end
    end
    if ! (solver in [:cplex, :clp, :cbc, :glpk])
        error("ERROR: unknown solver $(solver). Should be :cplex, :glpk, ...")
    end

    if solver == :glpk
        # voir https://github.com/JuliaOpt/GLPK.jl
        # model = Model(GLPK.Optimizer)  # PAS ENCORE SUPPORTE (05/12/2019)
        model = Model(with_optimizer(GLPK.Optimizer))
    elseif solver == :cplex
        model = Model(with_optimizer(CPLEX.Optimizer))
        # model = Model(CPLEX.Optimizer) # PAS ENCORE SUPPORTE (05/12/2019)
    elseif solver == :clp
        # voir https://github.com/JuliaOpt/Clp.jl
        # model = Model(Clp.Optimizer) # PAS ENCORE SUPPORTE (05/12/2019)
        model = Model(with_optimizer(Clp.Optimizer))
    elseif solver == :cbc
        # Voir https://github.com/JuliaOpt/Cbc.jl pour les options
        # model = Model(Cbc.Optimizer) # PAS ENCORE SUPPORTE (05/12/2019)
        model = Model(with_optimizer(Cbc.Optimizer))
    else
        error("LpTimingSolver: unknown external_mip_solver inconnu : $solver")
    end
    # MOI.set(model, MOI.Silent(), true) # Variante de bas niveau
    if log_level <= 3
        set_silent(model) # unset_silent(model) pour réautoriser affichage
    end
    return model
end

ext: jl


text:
# Base Type
## API:
### initialparameters(rng, l) --> Must return a NamedTuple/ComponentArray
### initialstates(rng, l)     --> Must return a NamedTuple
### parameterlength(l)        --> Integer
### statelength(l)            --> Integer
### l(x, ps, st)

abstract type AbstractExplicitLayer end

"""
    initialparameters(rng::AbstractRNG, l)

Generate the initial parameters of the layer `l`.
"""
initialparameters(::AbstractRNG, ::Any) = NamedTuple()
function initialparameters(rng::AbstractRNG, l::NamedTuple)
    map(Base.Fix1(initialparameters, rng), l)
end

"""
    initialstates(rng::AbstractRNG, l)

Generate the initial states of the layer `l`.
"""
initialstates(::AbstractRNG, ::Any) = NamedTuple()
initialstates(rng::AbstractRNG, l::NamedTuple) = map(Base.Fix1(initialstates, rng), l)

"""
    parameterlength(l)

Return the total number of parameters of the layer `l`.
"""
function parameterlength(l::AbstractExplicitLayer)
    parameterlength(initialparameters(Random.default_rng(), l))
end
function parameterlength(nt::Union{NamedTuple, Tuple})
    length(nt) == 0 ? 0 : sum(parameterlength, nt)
end
parameterlength(a::AbstractArray) = length(a)
parameterlength(x) = 0

"""
    statelength(l)

Return the total number of states of the layer `l`.
"""
statelength(l::AbstractExplicitLayer) = statelength(initialstates(Random.default_rng(), l))
statelength(nt::Union{NamedTuple, Tuple}) = length(nt) == 0 ? 0 : sum(statelength, nt)
statelength(a::AbstractArray) = length(a)
statelength(x::Union{Number, Symbol}) = 1
statelength(x) = 0

"""
    setup(rng::AbstractRNG, l::AbstractExplicitLayer)

Shorthand for getting the parameters and states of the layer `l`. Is equivalent to `(initialparameters(rng, l), initialstates(rng, l))`.
"""
function setup(rng::AbstractRNG, l::AbstractExplicitLayer)
    (initialparameters(rng, l), initialstates(rng, l))
end

"""
    apply(model::AbstractExplicitLayer, x, ps::Union{ComponentArray,NamedTuple}, st::NamedTuple)

Simply calls `model(x, ps, st)`
"""
function apply(model::AbstractExplicitLayer, x, ps::Union{ComponentArray, NamedTuple},
               st::NamedTuple)
    model(x, ps, st)
end

function Base.show(io::IO, x::AbstractExplicitLayer)
    __t = rsplit(string(get_typename(x)), "."; limit=2)
    T = length(__t) == 2 ? __t[2] : __t[1]
    print(io, "$T()")
end

# Abstract Container Layers
abstract type AbstractExplicitContainerLayer{layers} <: AbstractExplicitLayer end

function initialparameters(rng::AbstractRNG,
                           l::AbstractExplicitContainerLayer{layers}) where {layers}
    length(layers) == 1 && return initialparameters(rng, getfield(l, layers[1]))
    return NamedTuple{layers}(initialparameters.(rng, getfield.((l,), layers)))
end

function initialstates(rng::AbstractRNG,
                       l::AbstractExplicitContainerLayer{layers}) where {layers}
    length(layers) == 1 && return initialstates(rng, getfield(l, layers[1]))
    return NamedTuple{layers}(initialstates.(rng, getfield.((l,), layers)))
end

function parameterlength(l::AbstractExplicitContainerLayer{layers}) where {layers}
    sum(parameterlength, getfield.((l,), layers))
end

function statelength(l::AbstractExplicitContainerLayer{layers}) where {layers}
    sum(statelength, getfield.((l,), layers))
end

# Test Mode
"""
    testmode(st::NamedTuple, mode::Bool=true)

Make all occurances of `training` in state `st` `!mode`
"""
testmode(st::NamedTuple, mode::Bool=true) = update_state(st, :training, Val(!mode))

"""
    trainmode(x::Any, mode::Bool=true)

Make all occurances of `training` in state `st` `mode`
"""
trainmode(x::Any, mode::Bool=true) = testmode(x, !mode)

"""
    update_state(st::NamedTuple, key::Symbol, value; layer_check=_default_layer_check(key))

Recursively update all occurances of the `key` in the state `st` with the `value`.
"""
function update_state(st::NamedTuple, key::Symbol, value;
                      layer_check=_default_layer_check(key))
    function _update_state(st, key::Symbol, value)
        return Setfield.set(st, Setfield.PropertyLens{key}(), value)
    end
    return fmap(_st -> _update_state(_st, key, value), st; exclude=layer_check)
end

function _default_layer_check(key)
    _default_layer_check_closure(x) = hasmethod(keys, (typeof(x),)) ? key ∈ keys(x) : false
    return _default_layer_check_closure
end


text:
# used in shortest path calculations

function eccentricity(g::AbstractGraph,
    vs=vertices(g),
    distmx::AbstractMatrix{T}=weights(g)) where T <: Real
    vlen = length(vs)
    eccs = SharedVector{T}(vlen)
    @sync @distributed for i = 1:vlen
        eccs[i] = maximum(LightGraphs.dijkstra_shortest_paths(g, vs[i], distmx).dists)
    end
    d = sdata(eccs)
    maximum(d) == typemax(T) && @warn("Infinite path length detected")
    return d
end

eccentricity(g::AbstractGraph, distmx::AbstractMatrix) =
    eccentricity(g, vertices(g), distmx)

diameter(g::AbstractGraph, distmx::AbstractMatrix=weights(g)) =
    maximum(eccentricity(g, distmx))

periphery(g::AbstractGraph, distmx::AbstractMatrix=weights(g)) =
    LightGraphs.periphery(eccentricity(g, distmx))

radius(g::AbstractGraph, distmx::AbstractMatrix=weights(g)) =
    minimum(eccentricity(g, distmx))

center(g::AbstractGraph, distmx::AbstractMatrix=weights(g)) =
    LightGraphs.center(eccentricity(g, distmx))


text:
using RigidBodyIntegrators
using Test

@testset "RigidBodyIntegrators.jl" begin
    # Write your tests here.
end


text:
# floating point tests used in Julia Base

@testset "flipsign" begin
    x = convert(Single32,-2.0)
    x = flipsign(x,-1.0)
    @test flipsign(x,big(-1.0)) == convert(Single32,-2.0)
end

@testset "maxintfloat" begin
    @test maxintfloat(Single32, Int64) === maxintfloat(Float32, Int64)
    @test maxintfloat(Single32, Int32) === maxintfloat(Float32, Int32)
    @test maxintfloat(Single32, Int16) === maxintfloat(Float32, Int16)
end

@testset "isinteger" begin
    @test !isinteger(Single32(1.2))
    @test isinteger(Single32(12))
    @test isinteger(zero(Single32))
    @test isinteger(-zero(Single32))
    @test !isinteger(nextfloat(zero(Single32)))
    @test !isinteger(prevfloat(zero(Single32)))
    @test isinteger(maxintfloat(Single32))
    @test isinteger(-maxintfloat(Single32))
    @test !isinteger(Single32(Inf))
    @test !isinteger(-Single32(Inf))
    @test !isinteger(Single32(NaN))
end

@testset "round" begin
    x = Single32(rand(Float32))
    A = fill(x,(10,10))
    @test round.(A,RoundToZero) == fill(trunc(x),(10,10))
    @test round.(A,RoundUp) == fill(ceil(x),(10,10))
    @test round.(A,RoundDown) == fill(floor(x),(10,10))
    A = fill(x,(10,10,10))
    @test round.(A,RoundToZero) == fill(trunc(x),(10,10,10))
    @test round.(A,RoundUp) == fill(ceil(x),(10,10,10))
    @test round.(A,RoundDown) == fill(floor(x),(10,10,10))
end

@testset "round2" begin
    x = Single32(rand(Float32))
    for elty2 in (Int32,Int64)
        A = fill(x,(10,))
        @test round.(elty2,A,RoundToZero) == fill(trunc(elty2,x),(10,))
        @test round.(elty2,A,RoundUp) == fill(ceil(elty2,x),(10,))
        @test round.(elty2,A,RoundDown) == fill(floor(elty2,x),(10,))
        A = fill(x,(10,10))
        @test round.(elty2,A,RoundToZero) == fill(trunc(elty2,x),(10,10))
        @test round.(elty2,A,RoundUp) == fill(ceil(elty2,x),(10,10))
        @test round.(elty2,A,RoundDown) == fill(floor(elty2,x),(10,10))
        A = fill(x,(10,10,10))
        @test round.(elty2,A,RoundToZero) == fill(trunc(elty2,x),(10,10,10))
        @test round.(elty2,A,RoundUp) == fill(ceil(elty2,x),(10,10,10))
        @test round.(elty2,A,RoundDown) == fill(floor(elty2,x),(10,10,10))
        @test round.(elty2,A) == fill(round(elty2,x),(10,10,10))
    end
end


text:
#Modify the previous program such that only the users Alice and Bob are greeted with their names.

println("Please Enter your name:")
name = readline()

function greeter(x)
    println("Hello $x, welcome back.")
end

if name == "Bob" || name=="Alice"
    greeter(name)
else
    println("You are not Bob or Alice")
end

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
export BlankEnv

mutable struct BlankEnv <: Environment
  worldAABBmin::Vector
  worldAABBmax::Vector
  keepin_zones::Vector
  keepout_zones::Vector
  obstacle_set::Vector
end

function BlankEnv()
  worldAABBmin = -1000*ones(3)
  worldAABBmax = 1000*ones(3)
  keepin_zones = Vector{GeometryTypes.GeometryPrimitive}()
  keepout_zones = Vector{GeometryTypes.GeometryPrimitive}()
  obstacle_set = Vector{GeometryTypes.GeometryPrimitive}()
  BlankEnv(worldAABBmin, worldAABBmax, keepin_zones, keepout_zones, obstacle_set)
end


text:
function MillerUpdatingRegression(arg0::jint, arg1::jboolean)
    return MillerUpdatingRegression((jint, jboolean), arg0, arg1)
end

function MillerUpdatingRegression(arg0::jint, arg1::jboolean, arg2::jdouble)
    return MillerUpdatingRegression((jint, jboolean, jdouble), arg0, arg1, arg2)
end

function add_observation(obj::MillerUpdatingRegression, arg0::Vector{jdouble}, arg1::jdouble)
    return jcall(obj, "addObservation", void, (Vector{jdouble}, jdouble), arg0, arg1)
end

function add_observations(obj::MillerUpdatingRegression, arg0::Vector{Vector{jdouble}}, arg1::Vector{jdouble})
    return jcall(obj, "addObservations", void, (Vector{Vector{jdouble}}, Vector{jdouble}), arg0, arg1)
end

function clear(obj::MillerUpdatingRegression)
    return jcall(obj, "clear", void, ())
end

function get_diagonal_of_hat_matrix(obj::MillerUpdatingRegression, arg0::Vector{jdouble})
    return jcall(obj, "getDiagonalOfHatMatrix", jdouble, (Vector{jdouble},), arg0)
end

function get_n(obj::MillerUpdatingRegression)
    return jcall(obj, "getN", jlong, ())
end

function get_order_of_regressors(obj::MillerUpdatingRegression)
    return jcall(obj, "getOrderOfRegressors", Vector{jint}, ())
end

function get_partial_correlations(obj::MillerUpdatingRegression, arg0::jint)
    return jcall(obj, "getPartialCorrelations", Vector{jdouble}, (jint,), arg0)
end

function has_intercept(obj::MillerUpdatingRegression)
    return jcall(obj, "hasIntercept", jboolean, ())
end

function regress(obj::MillerUpdatingRegression)
    return jcall(obj, "regress", RegressionResults, ())
end

function regress(obj::MillerUpdatingRegression, arg0::Vector{jint})
    return jcall(obj, "regress", RegressionResults, (Vector{jint},), arg0)
end

function regress(obj::MillerUpdatingRegression, arg0::jint)
    return jcall(obj, "regress", RegressionResults, (jint,), arg0)
end


text:
import StatsBase
square_error2RMSE(sqrerr, n) = sqrt(sqrerr/n)
abs_error2MAE(abserr, n) = sqrerr/n

"""
$(SIGNATURES)

Returns tuples (control-lines, fitting-lines, all-lines):
- statistics of absolute errors (StatsBase.summarystats)
- statistics of square errors
- cumulative h error
- TODO ratio of track-length with errors smaller than:
  - 20%
  - std of h
  - 2std of h
"""
function compare_to_radar(sol::MCMCSol, radarlines_fitting::AbstractVector{<:Integer}=Int[],
                          radarlines_control=nothing; # nothing means all but the radarlines_fitting
                          thin=1, meters_thin=0.0)
    compare_to_radar(sol.theta0.gb.gl, sol.mh2d, sol.sh2d, radarlines_fitting, radarlines_control;
                     thin=thin, meters_thin=meters_thin)
end
function compare_to_radar(gl, hs2d, hs2ds, radarlines_fitting::AbstractVector{<:Integer}=Int[],
                          radarlines_control=nothing; # nothing means all but the radarlines_fitting
                          thin=1, meters_thin=0.0)
    hs2d = hs2d isa Gridded ? hs2d.v : hs2d
    hs2ds = hs2ds isa Gridded ? hs2ds.v : hs2ds
    if gl.h==SomeData(nothing)
        return NaN, NaN, NaN, (), (), ()
    end
    radarlines_all = 1:length(gl.h.vals.splits)
    if radarlines_control==nothing
        radarlines_control = [i for i=1:length(gl.h.vals.splits) if !(i in radarlines_fitting)]
    end
    rls = (radarlines_control, radarlines_fitting, radarlines_all)
#     RMSE = map(x->square_error2RMSE(x...), map(r->thick_err(gl, hs2d, r, thin, sqrerr), rls))
#     max_err = map(r->thick_err(gl, hs2d, r, thin, abserr, +)[1], rls)

# @show    StatsBase.summarystats(thick_err(gl, hs2d, rls[1], thin, abserr, push!, Float64[])[1])
#     @show map(r->isnan(r) StatsBase.summarystats(thick_err(gl, hs2d, r, thin, abserr, push!, Float64[])[1]), rls)
#     MAE = map(x->x[1]/x[2], map(r->thick_err(gl, hs2d, r, thin, abserr, +), rls))
#     quantiles = ()

    SEs = []
    AEs = []
    errs = []
    for r in rls
        err, n = thick_err(gl, hs2d, r, thin, -, push!, Float64[])

        if isequal(err, NaN) || isempty(err)
            push!(AEs, NaN)
            push!(SEs, NaN)
            push!(errs, NaN)
        else
            push!(AEs, StatsBase.summarystats(abs.(err)))
            push!(SEs, StatsBase.summarystats(err.^2))
            push!(errs, StatsBase.summarystats(err))
        end
    end

    # # errors as Traj:
    # err = error_h(hs2d, gl, radarlines_fitting, in_glacier(gl))

    # smaller than
    smaller_20 = () # map(r->/(thick_err(gl, hs2d, r, thin,
                    #              (model,data) -> abs(model-data)/abs(data)<0.2, +)...), rls)
    smaller_std = ()
    smaller_2std = ()

    # remember the other `return`
    return AEs, SEs, errs, smaller_20, smaller_std, smaller_2std
end

"""
$(SIGNATURES)

"""
function compare_to_iv(gl, iv2d, iv2ds, fwdsol)
    error("TODO: Update to same API as compare_to_radar")
    iiv  = Interp.interpolate((dem.x, dem.y), ivs2d, Interp.Gridded(Interp.Linear()) )
    iivs  = Interp.interpolate((dem.x, dem.y), ivs2ds, Interp.Gridded(Interp.Linear()) )

    thin = 1
    if gl.iv.vals isa Gridded{F}
        ivmask = gl.ivmask .& fwdsol.ivs2d_band_mask
        RMSE = square_error2RMSE(iv2d_err(gl, iv2d, ivmask, thin)...)
        max_err = iv2d_err(gl, iv2d, ivmask, thin, abserr, max)
        # produces garbage, would need some abstol:
        # max_rel_err = iv2d_err(gl, iv2d, ivmask, thin, rel_abserr, max)

    elseif gl.iv.vals isa Traj{F}

    end
    quantiles = ()
    return RMSE, max_err, MAE, quantiles, smaller_20, smaller_std, smaller_2std
end

"""
    summarize_mcmc_sols(sols::Vector{<:MCMCSol})

Returns the quartiles of the quartiles of the absolute error
of the solutions ins `sol`.

Return: q1, median, q3

TODO:
- allow using particular fitting lines?
"""
function summarize_mcmc_sols(sols::Vector{<:MCMCSol})
    # do stats on quartiles
    q1 = Float64[]
    q2 = Float64[]
    q3 = Float64[]
    for sol in sols
        ae = compare_to_radar(sol)[1][1]
        push!(q1, ae.q25)
        push!(q2, ae.median)
        push!(q3, ae.q75)
    end
    # now summarize the summaries

    return StatsBase.summarystats.((q1, q2, q3))
end
function summarize_mcmc_sols(gls::Vector, hs2d::Vector, hs2ds::Vector)
    # do stats on quartiles
    q1 = Float64[]
    q2 = Float64[]
    q3 = Float64[]
    for (gl, h, hs) in zip(gls, hs2d, hs2ds)
        ae = compare_to_radar(gl, h, hs)[1][1]
        push!(q1, ae.q25)
        push!(q2, ae.median)
        push!(q3, ae.q75)
    end
    # now summarize the summaries

    return StatsBase.summarystats.((q1, q2, q3))
end


text:
using Plots, Distributions, Random; pyplot()
Random.seed!(0)

n = 2000
data = rand(Normal(),n)
l, m = minimum(data), maximum(data)

delta = 0.3;
bins = [(x,x+delta) for x in l:delta:m-delta]
if last(bins)[2] < m
    push!(bins,(last(bins)[2],m))
end
L = length(bins)

inBin(x,j) = first(bins[j]) <= x && x < last(bins[j])
sizeBin(j) = last(bins[j]) - first(bins[j])
f(j) = sum([inBin(x,j)  for x in data])/n
h(x) = sum([f(j)/sizeBin(j) * inBin(x,j) for j in 1:L])

xGrid = -4:0.01:4
histogram(data,normed=true, bins=L,
    label="Built-in histogram",
    c=:blue, la=0, alpha=0.6)
plot!(xGrid,h.(xGrid), lw=3, c=:red, label="Manual histogram",
    xlabel="x",ylabel="Frequency")
plot!(xGrid,pdf.(Normal(),xGrid),label="True PDF",
    lw=3, c=:green, xlims=(-4,4), ylims=(0,0.5))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
module LanguageFinder

import HTTP
import StatsBase

export LanguageFinder, train_wikipedia_text

include("Wikipedia/wikipedia_corpus_buid.jl")
include("language_finder.jl")

end #module

max_stars_repo_name: nusretipek/LanguageFinder

max_issues_repo_name: nusretipek/LanguageFinder

max_forks_repo_name: nusretipek/LanguageFinder


text:
using PyPlot
using MAT

function replay(fileName)

fileConfig = matread(fileName)

history_W = fileConfig["history_W"]
history_A = fileConfig["history_A"]
history_M = fileConfig["history_M"]
history_F = fileConfig["history_F"]
history_dir = fileConfig["history_dir"]

#history_DAvec = fileConfig["DAvec"]
#history_DMvec = fileConfig["DMvec"]

hold(false)
for t in 1:size(history_M, 3)
    subplot(151)
    pcolormesh(history_M[:,:,t], vmin=0, vmax=0.6)
    title("Mfield")

    subplot(152)
    pcolormesh(history_A[:,:,t], vmin=0, vmax=0.6)
    title("Afield")

    subplot(153)
    pcolormesh(history_F[:,:,t], vmin=0, vmax=1)
    title("Ffield")

    subplot(154)
    pcolormesh(history_W[:,:,t], vmin=0, vmax=1)
    title("Wfield")

    subplot(155)
    title("DirectionField")
    directionfield = history_dir[:,:,t]

    U = cos(directionfield)
    V = sin(directionfield)
    d1, d2 = size(directionfield)
    plt.quiver([1:d1], [1:d2], U, V, linewidth=1.5, headwidth = 0.5)
    yield()
end
end

if length(ARGS) == 1
    replay(first(ARGS))
end


text:
push!(LOAD_PATH, ".")

module GASelection

using Debug
using Population
using Individual
using Random

abstract type Tournament end

function selectParents(individuals::Vector{<: Population.IndFitType{IndType}},  parent_num::Integer, ::Type{Tournament}, tour_size::Integer) where {IndType}
  parents = Vector{Tuple{IndType, IndType}}(undef, parent_num)
  choices = [(individuals[i], i) for i = eachindex(individuals)]

  Debug.ga_debug && println("----- Tournament Start -----\n")

  for i = 1 : parent_num
    tournament = shuffle(choices)[1:tour_size]
    p1 = maximum(tournament)[2]

    if Debug.ga_debug
      println("Selection ", i)
      println("First tournament: ", [(b, a) for ((x, a), b) in tournament], " -> winner ", p1, "\n")
    end

    tournament = eltype(choices)[]
    for (ind, i) in shuffle(choices)
      i == p1 || push!(tournament, (ind, i))
      length(tournament) < tour_size || break
    end
    p2 = maximum(tournament)[2]
    for ((x, a), b) in tournament
      (a <= individuals[p2][2]) || error("LOL")
    end

    Debug.ga_debug && println("Second tournament: ", [(b, a) for ((x, a), b) in tournament], " -> winner ", p2, "\n")

    parents[i] = (individuals[p1][1], individuals[p2][1])
  end

  Debug.ga_debug && println("\n----- Tournament End -----\n")

  return parents
end

end

max_stars_repo_name: Markhyz/Julia-Genetic-Algorithm

max_issues_repo_name: Markhyz/Julia-Genetic-Algorithm

max_forks_repo_name: Markhyz/Julia-Genetic-Algorithm


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

using Markdown
using InteractiveUtils

# ╔═╡ 250252f0-6e61-11ec-2ffe-ebddaa9f0fb3
begin
	using Pkg
	Pkg.activate("../../")
	using Revise
end

# ╔═╡ 0eeb5bb1-bb41-43cc-8c14-eb9c1ebf6fad
using ImgProcMic, FourierTools, IndexFunArrays, ImageShow, FFTW, TestImages, Colors

# ╔═╡ 1b2ebe1a-f3a7-494c-bcca-18df5403ee65
using PlutoTest, Noise, Plots

# ╔═╡ 4b46e310-cc26-4a96-8c2e-0c89307d4e34
md"# Structured Illumination Microscopy

In this homework we are going to work through the full pipeline from capturing a dataset of a SIM to reconstructing the data.

For the demonstration we focus on the 2D case and consider the resolution enhancement only along one direction. That allows us to observe the effect clearly.
"

# ╔═╡ 7054db2c-606e-48eb-ab93-8c0260cb7a81
md"## 1 Simulate PSF"

# ╔═╡ edaacc75-fd91-4d46-a31f-21e738253708
"""
	simulate_PSF(s, r, [T=Float32])

Simulate the incoherent 2D PSF with radius `r` and total output size `s`.
This functions returns the PSF centered around the first entry (top left corner).
Furthermore, it returns the `freq_limit` of the PSF.
The `freq_limit` is the frequency at which the OTF becomes 0.

The procedure is as following:
* create a disc with values of 1 and radius `r` (use `rr2`)
* take the `ift` of it.
* abs2.
* normalize the sum to 1
* center pixcel should be located in the top left
"""
function simulate_PSF(s, r, T=Float32)
	disc = rr2(T, (s)) .< r^2

	# is equivalent of `ifft(disc)` is equivalent to `ifftshift(ift(disc))`
	# but only because
	# we removed the phase by abs2
	disc_ift = ift(disc)

	psf = abs2.(disc_ift)
	psf ./= sum(psf)

	# shift from center to first entry
	psf = psf
	# don't change freq_limit, is already

	freq_limit = r / (s[1] / 2)
	return ifftshift(psf), freq_limit
end

# ╔═╡ b6061e31-e06a-45ad-86f0-006b8309cca7
rr2((4,4)) .< 1.3^2

# ╔═╡ cd5bef41-6246-4d88-91d9-c83b7a47e110
# leave this value equal to 40!
radius = 40

# ╔═╡ 29e6311e-eb5b-4073-81dd-7372ac06498e
psf, freq_limit = simulate_PSF((512,512), radius);

# ╔═╡ 86e0df4f-4dfe-40c3-a6d3-1a394442a9f7
# MTF = magnitude of the optical transfer function (OTF)
plot(fftshift(fftfreq(512)), abs.(ffts(psf)[:, 257]))

# ╔═╡ 963c4f23-e97b-41fd-9423-1080fae412e9
freq_limit

# ╔═╡ c4c8b57b-e09a-43f8-b371-bb1f0761215d
gray_show(psf)

# ╔═╡ 599aabf7-6333-4272-ae07-fec07bdb02aa
img = Float32.(testimage("resolution_test_512"));

# ╔═╡ 920c26f4-5750-4c83-871d-9622b3dca134
gray_show(img);

# ╔═╡ a550fc7f-ff7c-4fd5-9e40-d8ed00da200f
md"### 1 Test"

# ╔═╡ 7cba5a99-21a7-45d6-a92a-7eb1e75e39d3
PlutoTest.@test simulate_PSF((5, 5), 2)[1] ≈ (Float32[0.36 0.10472136 0.01527864 0.01527864 0.10472136; 0.10472136 0.030462677 0.0044444446 0.0044444446 0.030462677; 0.01527864 0.0044444446 0.00064843567 0.00064843567 0.0044444446; 0.01527864 0.0044444446 0.00064843567 0.00064843567 0.0044444446; 0.10472136 0.030462677 0.0044444446 0.0044444446 0.030462677], 0.8)[1]

# ╔═╡ eafd2e8a-6bbd-4fe4-a5fd-d721b78604e1
PlutoTest.@test simulate_PSF((5, 5), 2)[2] == 0.8

# ╔═╡ 75f2ea94-ba30-48d3-b676-8423f14c453c
delta = zeros((5,5));

# ╔═╡ 9b856e9c-b818-4e49-a27d-0c0db623698e
delta[1] = 1;

# ╔═╡ 11009f9c-de4c-4ade-ad30-375ed47c3901
PlutoTest.@test simulate_PSF((5, 5), 3)[1] ≈ delta

# ╔═╡ 4ff25625-e450-47e4-97ae-bb24b903ed98
md"## 2 Structured Illumination

The key idea behind SIM is to illuminate the sample with an illumination pattern of the form

$I_n = 1 + \cos(k_g \cdot x + \theta_n)$

where $k_g$ is the frequency of the illumination pattern (twice the grating constant). $\theta_n$ is a illumination pattern shift which is different for each image. Such an illumination allows to increase the resolution along the first dimension.

We capture three images ($n=\{1,2,3\}$) of the sample which can be expressed as

$Y_n = (I_n \cdot S) * h$
where $h$ is the PSF, $S$ the ideal sample and $Y$ the measurement.

In our example we choose the following 3 different illuminations

$I_{1} = 1 + \cos\left(k_g \cdot x \cdot 2 \pi + 1\cdot \frac{2 \pi}{3} \right)$

$I_2 = 1 + \cos\left(k_g \cdot x  \cdot 2 \pi+ 2 \cdot \frac{2 \pi}{3}\right)$

$I_3 = 1 + \cos\left(k_g \cdot x  \cdot 2 \pi+ 3 \cdot \frac{2 \pi}{3}\right)$

"

# ╔═╡ e62fc2f1-a1e2-4b8a-be08-94444e12b533
begin
	plot(1 .+ cos.(2π/3 .+ 2π .* (0.0:0.01:3)))
	plot!(1 .+ cos.(4π/3 .+ 2π .* (0.0:0.01:3)))
	plot!(1 .+ cos.(2π .* (0.0:0.01:3)))
	plot!(3 .+  cos.(2π/3 .+ 2π .* (0.0:0.01:3)) .+ cos.(4π/3 .+ 2π .* (0.0:0.01:3)) .+ cos.(6π/3 .+ 2π .* (0.0:0.01:3)))
end

# ╔═╡ ca8618e8-8389-4d30-8112-af0e987de417
"""
	illuminate(arr, kg)

Returns a vector of `arr` illuminated three times according to the SIM principle.
Implement the equations from above. Only illumination, no convolution with PSF.

You should return a vector of the three illumination arrays.
`[I_1, I_2, I_3]`.
"""
function illuminate(arr, kg)
	# x coordinates
	x = 1:size(arr, 1)

	# todo fix illumination
	I₁ = (1 .+ cos.(x .* kg .* 2π .+ 1 * 2π / 3)) .* arr
	I₂ = (1 .+ cos.(x .* kg .* 2π .+ 2 * 2π / 3)) .* arr
	I₃ = (1 .+ cos.(x .* kg .* 2π .+ 3 * 2π / 3)) .* arr

	# todo return correct results
	[I₁, I₂, I₃]
end

# ╔═╡ 20aaf9b2-2660-413f-90a8-1f4fdc91d8c7
# toy example
illuminate(ones((4,4)), 0.5)

# ╔═╡ 4898da63-4da0-4538-b08d-cd4fe419584b
Illus = illuminate(ones((512, 512)), freq_limit * 0.8);

# ╔═╡ 8457cfeb-705f-41e8-a203-298c8ee4a7ab
# display the three illuminations
gray_show(hcat(Illus...))

# ╔═╡ 28784e50-285f-47a8-a311-4122ec8f47ce
# check whether illumination still pass through the system.
# if you see the same stripes but with lower contrast, then everything
# is perfect
gray_show(conv(Illus[1], psf))#simulate_PSF(size(Illus[1]), 31)[1]))

# ╔═╡ 6e2814c9-5644-47b6-8a6e-3308e1623d5e
md"### 2 Test"

# ╔═╡ ca0e55ca-ec4f-4cd8-a1f7-a2a63ba9782d
PlutoTest.@test illuminate(ones((4, 4)), 0.5) ≈ [[1.4999999999999993 1.4999999999999993 1.4999999999999993 1.4999999999999993; 0.5000000000000008 0.5000000000000008 0.5000000000000008 0.5000000000000008; 1.4999999999999991 1.4999999999999991 1.4999999999999991 1.4999999999999991; 0.500000000000001 0.500000000000001 0.500000000000001 0.500000000000001], [1.5000000000000004 1.5000000000000004 1.5000000000000004 1.5000000000000004; 0.49999999999999867 0.49999999999999867 0.49999999999999867 0.49999999999999867; 1.5000000000000013 1.5000000000000013 1.5000000000000013 1.5000000000000013; 0.49999999999999845 0.49999999999999845 0.49999999999999845 0.49999999999999845], [0.0 0.0 0.0 0.0; 2.0 2.0 2.0 2.0; 0.0 0.0 0.0 0.0; 2.0 2.0 2.0 2.0]]

# ╔═╡ fa957648-0b4a-48eb-96d5-bb731966693c
Ittt = illuminate(ones((10, 10)), 1/3);

# ╔═╡ 287d85ae-efba-4095-b7ce-58d90bd22492
PlutoTest.@test Ittt[1] .+ Ittt[2] .+ Ittt[3] ≈ 3 .* ones((10, 10))

# ╔═╡ fce697cc-67d1-4373-b2b0-a6a60b8a5893
PlutoTest.@test all(illuminate([1 1; 0 0], 0) .≈ [[0.49999999999999956 0.49999999999999956; 0.0 0.0], [0.5000000000000002 0.5000000000000002; 0.0 0.0],  [2.0 2.0; 0.0 0.0],])

# ╔═╡ b961b3e0-6935-4f29-9f94-be6e21f302e3
PlutoTest.@test  illuminate([1 1; 0 0], 0.2) ≈ [[0.02185239926619431 0.02185239926619431; 0.0 0.0], [1.6691306063588578 1.6691306063588578; 0.0 0.0], [1.3090169943749477 1.3090169943749477; 0.0 0.0]]

# ╔═╡ 6e904b1d-2071-4d48-ad21-ab180e333474
md"## 3 Forward Imaging Process
Now we complete the forward modelling part from the image to the three measurements.

The grating constant should be small enough such that illumination
passes through the system.
Since the optical system has a bandlimit, not arbitrary small illumination patterns pass through the system.
The larger you choose it, the more of the reconstruction is affected by noise.
By  selecting `0.8 * freq_limit` we are definitely below the frequency threshold (e.g. 1.1 would be above)
"

# ╔═╡ e3994f91-a7bc-4811-b6e6-1f67692ff554
# don't change!
kg = 0.8 * freq_limit

# ╔═╡ 71af2b63-475c-429f-a15e-9a212f724a69
# the amount the three OTFs are displaced with respect to each other in Fourier space
# we need the value later for the reconstruction
# this value should stay at 64!
fourier_space_shift = radius * 2 * 0.8

# ╔═╡ 2a78c81d-f8f9-4d99-9c24-1aefb9941294
"""
	forward(img, psf, kg)

From the input `img` create three output images.
The output images are produced by multiplying the three different illumination
to the input image and convolving the three resulting `img`s with the `psf`.
Apply `poisson` to each of the image with a photon number of 1000 (use `N_phot`).
Use the function `poisson` for that.

## Procedure
* Multiplication with illumination
* convolution (use `FourierTools.conv`) with `psf`.
* Add noise with `poisson(img_without_noise, N_phot)`
* Return vectors of the three measuremens,e g. `[meas_1, meas_2, meas_3]`

"""
function forward(img, psf, kg, N_phot=1000)
	# already the illuminated images
	illus = illuminate(img, kg)

	# convolve and add Poisson noise
	Y₁ = poisson(conv(illus[1], psf), N_phot)
	Y₂ = poisson(conv(illus[2], psf), N_phot)
	Y₃ = poisson(conv(illus[3], psf), N_phot)

	# return correct result
	return [Y₁, Y₂, Y₃]
end

# ╔═╡ d68d9f24-bcea-4bcc-b597-3c051e77d99c
Is = forward(img, psf, kg);

# ╔═╡ 23c98291-4c72-4a60-b722-910eadaf7e77
md"Three images which look all slightly different because"

# ╔═╡ 9a75e192-9605-4787-acee-706ee7c90549
gray_show(hcat(Is...))

# ╔═╡ 524b208d-7fc3-4d96-bdc1-2bed72f2a9ee
md"In comparison without any structured illumination"

# ╔═╡ 2c1a5b9e-8b1c-4619-b303-69d29181a5c4
[gray_show(conv(img, psf)) gray_show(Is[1])]

# ╔═╡ bb1e5dd4-0159-47b5-a803-5327952438ee
gray_show(img)

# ╔═╡ 09593275-85bf-42ee-860e-ebbb74866a40
md"## 3 Test"

# ╔═╡ 1103b35d-cb86-4970-beb7-c8463e374f78
# test random example
PlutoTest.@test all(.≈(forward([1 2; 3 4], [0.55 0.15; 0.15 0.15], 0.1, 100_000_000),  [[0.09645833671795113 0.13103484317669054; 0.08808351198799098 0.09683198111576993], [2.6912285614539946 3.133633214436668; 4.25256575017253 4.919802755223418], [2.9120113512125654 3.6367083598655254; 3.7589874298841632 4.283291164312841]], rtol=0.01))

# ╔═╡ ba74ec4a-b38f-4fcb-9b43-2448976153eb
PlutoTest.@test radius == 40

# ╔═╡ ab56310f-671d-4b54-b9c2-9be102c826e3
PlutoTest.@test kg == 1/8

# ╔═╡ db20d906-b22f-4f16-b355-eba9c82809d7
PlutoTest.@test fourier_space_shift == 64

# ╔═╡ b7f11d85-6859-4b5a-814e-003946fbd894
md"# 4 Unmixing
The three measurements are not the different parts of the Fourier spectrum yet.
But they contain the information to high frequencies!
Therefore, we want to unmix the three parts of the Fourier space.
"

# ╔═╡ 96b8f1e6-9159-4616-be08-c88f17e9829e
"""
	extract_components(Is)

Extract the different Fourier components from the 3 measurements (the vector `Is`).
This is achieved by inverting the mixing matrix (`inv(M)`) and multiplying (use a `*` for matrix multiplication, not a pointwise `.*`) it from the left to
`[fft(Is[1]), fft(Is[2]), fft(Is[3])]`.

See the lecture slides for more information.
"""
function extract_components(Is)
	θ₁ = 2π / 3 * 1
	θ₂ = 2π / 3 * 2
	θ₃ = 2π / 3 * 3

	# TODO create correct matrix
	M̂ = [exp(-1im * θ₁) 1 exp(1im * θ₁);
		 exp(-1im * θ₂) 1 exp(1im * θ₂);
		 exp(-1im * θ₃) 1 exp(1im * θ₃)]

	# is the same Ĩ = fft.(Is)
	Ĩ = [fft(Is[1]), fft(Is[2]), fft(Is[3])]

	# TODO return correct result
	return inv(M̂) * Ĩ
end

# ╔═╡ 53af3f77-cbc2-4809-8dda-743e096df10c
Cₙ = extract_components(Is);

# ╔═╡ 4c136f35-e489-4d6e-83e0-e153bfc7040d
md"Here you can see the three different components still. They are not shifted to the correct position, yet.
"

# ╔═╡ bf97497a-e149-4552-9bdd-9fbeb5c24987
gray_show(log1p.(abs.(hcat(fftshift.(Cₙ)...))))

# ╔═╡ b3528818-74f7-471b-8d33-ad7b9fcfc70b
gray_show(log1p.(abs.(ft(img))));

# ╔═╡ 1ff22cf9-ac3b-4849-bd4c-8b0d3a901447
circshift([1,2,3], 1)

# ╔═╡ c091d147-4874-4a19-9855-c78ded59b67a
gray_show(real.(fft(Cₙ[1])))

# ╔═╡ c922f283-3b31-43bd-9855-933be8a9e095
md"### 4 Test"

# ╔═╡ 460a8ea2-d077-48fc-804e-7367ac16eff8
# check if three output array
PlutoTest.@test size(extract_components([randn((2,2)) for i = 1:3])) == (3,)

# ╔═╡ ba0e90f8-a310-40f9-a306-fb8700ea473f
# check if structure and type is ok
PlutoTest.@test typeof(extract_components([randn((2,2)) for i = 1:3])) == Vector{Matrix{ComplexF64}}

# ╔═╡ 5d263158-03b0-4c87-ad90-a81c24493b66
# check a random example for correctness
PlutoTest.@test extract_components([[1 2; 3 4], [1 3; 4 5], [1 10; 20 30]]) ≈ Matrix{ComplexF64}[[16.50000000000001 - 0.8660254037844412im -5.500000000000003 + 0.2886751345948141im; -11.500000000000005 + 0.28867513459481564im 0.5000000000000002 + 0.28867513459481253im], [28.00000000000001 + 1.5597088666161704e-15im -8.000000000000002 - 3.563506943082403e-16im; -16.000000000000007 - 6.01679086153965e-16im 1.1102230246251565e-16 - 1.1102230246251565e-16im], [16.5 + 0.8660254037844446im -5.5 - 0.2886751345948149im; -11.5 - 0.28867513459481736im 0.5 - 0.2886751345948127im]]

# ╔═╡ 515d2174-56e3-4ebc-a475-b84f911a5a2f
md"# 5 Reconstruction
Having the three parts of the Fourier space separated, we can try to combine them in a reasonable way to one Fourier spectrum.
Simply adding them would not be optimal with respect to Signal-to-Noise-Ratio (SNR).
"

# ╔═╡ a97d7c94-e69b-44cd-ae4c-61d8d86808a5
"""
	reconstruct(psf, Cₙ, fourier_space_shift)

`psf` is the PSF of the single image (without SIM part).
`Cₙ` is an vector containing the three components of the Fourier spectrum.
`fourier_space_shift` is the shift the three pictures are shifted with respect to the center.

## Procedure
* First shift the spectrum (`otf = fft(psf)`) by `fourier_space_shift` along the first dimension (`circshift`)
* apply the same shift for the weighting factors
* Combine the Fourier spectra using _weighted averaging_ (see slides)
* also return the effective OTF
"""
function reconstruct(psf, Cₙ, fourier_space_shift)
	# Int shift
	Δ = round(Int, fourier_space_shift)

	otf = fft(psf)
	# TODO fix the weights
	w₋₁ = 0.5 * circshift(otf, +Δ)
	w₀ = otf
	w₁ =  0.5 * circshift(otf, -Δ)

	# TODO fix the mixing
	C₋₁ = circshift(Cₙ[1], +Δ)
	C₀ = Cₙ[2]
	C₁ = circshift(Cₙ[3], -Δ)

	# the small 1f-8 factor is added to prevent division by zero
	res_fourier_space = ((w₋₁ .* C₋₁ .+ w₀ .* C₀ .+ w₁ .* C₁ .+ 1f-8^2) ./
		  (w₋₁ .+ w₀ .+ w₁ .+ 1f-8))

	# todo calculate similarly to res the eff_otf but adapt it according to the slides
	eff_otf = ((w₋₁.^2 .+ w₀.^2 .+ w₁.^2 .+ 1f-8^2) ./
		      (w₋₁ .+ w₀ .+ w₁ .+ 1f-8))

	# todo adapt res_fourier_space to be a real space image
	res = real(ifft(res_fourier_space))

	# don't change return args
	return res, eff_otf
end

# ╔═╡ 9e64def4-b397-4292-8fe9-8a828e08ee0a
res, eff_otf = reconstruct(psf, Cₙ, fourier_space_shift);

# ╔═╡ 54fa3a01-5a8d-4dd3-a521-a60a2f69b541
md"Unfiltered reconstruction"

# ╔═╡ d9c33b8c-6ccf-462f-b5c2-267fd7b93ba3
gray_show(res)

# ╔═╡ 2f14fb93-91bc-44af-8a6d-cebaa1607a51
md"Unfiltered Fourier spectrum"

# ╔═╡ fd53c66d-ba90-4d96-9a11-eabc1dd666dc
gray_show(log1p.(abs.(ft(res))))

# ╔═╡ d7091577-a5f9-4edf-b96f-88b13d9c6f3c
gray_show(log1p.(abs.(ft(img))))

# ╔═╡ 19546bad-d7be-4c09-a248-f7e98c53d42d
"""
	wiener_filter(img, otf, ϵ)

Copy the wiener filter from previous homeworks and adapt it (OTF instead of the PSF is passed).

"""
function wiener_filter(img, otf, ϵ)
	S = 1
	# wiener kernel in fourier space
	W = conj.(otf) .* S ./ (abs2.(otf) .* S .+ ϵ)

	img = real(ifft(fft(img) .* W))
	return img
end

# ╔═╡ 33c156b6-5c18-4377-9920-32c7865898e9
md"### Final Inspection"

# ╔═╡ 0cdf07a5-a87e-449e-ab09-871924cee6d1
md"
If everything is correct, we should definitely see an resolution improvement in the left image (the SIM reconstruction).
The resolution increase is mainly along the first dimension.
"

# ╔═╡ 0a729667-599a-41f1-98f3-ff94608c15bb
[ gray_show(wiener_filter(res, eff_otf, 1e-3)) gray_show(wiener_filter(poisson(conv(img, psf), 300), fft(psf), 1e-3)) gray_show(img)]

# ╔═╡ 0d78c765-1912-4ded-a5e8-aa98afe73cb4
gray_show(log1p.(abs.(ft(wiener_filter(res, eff_otf, 1e-1)))));

# ╔═╡ 5566b0c6-ac96-4342-a62b-0c292469c2eb
[gray_show(wiener_filter(res, eff_otf, 1e-3)) gray_show(wiener_filter(poisson(conv(img, psf), 300), fft(psf), 1e-3))]

# ╔═╡ 3bed0c72-ab5c-4dbf-85a5-6f6d062dec97
md"The left image shows the effective OTF for the SIM system.
On the right we can see the OTF of the same microscope without structured illumination"

# ╔═╡ d32536af-b369-4b39-87f6-4e1172d971ae
[gray_show(fftshift(log1p.(abs.(eff_otf))))  gray_show(log1p.(abs.(ffts(psf))))]

# ╔═╡ 35cdd4f9-87f7-408e-a3fd-7456b7238dd5
md"## 5 Test"

# ╔═╡ 1320b506-d378-4fa1-bbf1-c2f994abfff1
PlutoTest.@test wiener_filter([1, 2, 3, 4, 5], [1, -2, 3, 4, 0], 0.001) ≈ [2.9553001744303087, 3.5816540338359086, 2.9970029970029977, 2.4123519601700862, 3.0387058195756858]

# ╔═╡ Cell order:
# ╠═250252f0-6e61-11ec-2ffe-ebddaa9f0fb3
# ╠═0eeb5bb1-bb41-43cc-8c14-eb9c1ebf6fad
# ╠═1b2ebe1a-f3a7-494c-bcca-18df5403ee65
# ╟─4b46e310-cc26-4a96-8c2e-0c89307d4e34
# ╟─7054db2c-606e-48eb-ab93-8c0260cb7a81
# ╠═edaacc75-fd91-4d46-a31f-21e738253708
# ╠═86e0df4f-4dfe-40c3-a6d3-1a394442a9f7
# ╠═963c4f23-e97b-41fd-9423-1080fae412e9
# ╠═b6061e31-e06a-45ad-86f0-006b8309cca7
# ╠═cd5bef41-6246-4d88-91d9-c83b7a47e110
# ╠═29e6311e-eb5b-4073-81dd-7372ac06498e
# ╠═c4c8b57b-e09a-43f8-b371-bb1f0761215d
# ╠═599aabf7-6333-4272-ae07-fec07bdb02aa
# ╠═920c26f4-5750-4c83-871d-9622b3dca134
# ╟─a550fc7f-ff7c-4fd5-9e40-d8ed00da200f
# ╠═7cba5a99-21a7-45d6-a92a-7eb1e75e39d3
# ╠═eafd2e8a-6bbd-4fe4-a5fd-d721b78604e1
# ╠═75f2ea94-ba30-48d3-b676-8423f14c453c
# ╠═9b856e9c-b818-4e49-a27d-0c0db623698e
# ╠═11009f9c-de4c-4ade-ad30-375ed47c3901
# ╟─4ff25625-e450-47e4-97ae-bb24b903ed98
# ╠═e62fc2f1-a1e2-4b8a-be08-94444e12b533
# ╠═ca8618e8-8389-4d30-8112-af0e987de417
# ╠═20aaf9b2-2660-413f-90a8-1f4fdc91d8c7
# ╠═4898da63-4da0-4538-b08d-cd4fe419584b
# ╠═8457cfeb-705f-41e8-a203-298c8ee4a7ab
# ╠═28784e50-285f-47a8-a311-4122ec8f47ce
# ╟─6e2814c9-5644-47b6-8a6e-3308e1623d5e
# ╠═ca0e55ca-ec4f-4cd8-a1f7-a2a63ba9782d
# ╠═fa957648-0b4a-48eb-96d5-bb731966693c
# ╠═287d85ae-efba-4095-b7ce-58d90bd22492
# ╠═fce697cc-67d1-4373-b2b0-a6a60b8a5893
# ╠═b961b3e0-6935-4f29-9f94-be6e21f302e3
# ╟─6e904b1d-2071-4d48-ad21-ab180e333474
# ╠═e3994f91-a7bc-4811-b6e6-1f67692ff554
# ╠═71af2b63-475c-429f-a15e-9a212f724a69
# ╠═2a78c81d-f8f9-4d99-9c24-1aefb9941294
# ╠═d68d9f24-bcea-4bcc-b597-3c051e77d99c
# ╟─23c98291-4c72-4a60-b722-910eadaf7e77
# ╠═9a75e192-9605-4787-acee-706ee7c90549
# ╟─524b208d-7fc3-4d96-bdc1-2bed72f2a9ee
# ╠═2c1a5b9e-8b1c-4619-b303-69d29181a5c4
# ╠═bb1e5dd4-0159-47b5-a803-5327952438ee
# ╟─09593275-85bf-42ee-860e-ebbb74866a40
# ╠═1103b35d-cb86-4970-beb7-c8463e374f78
# ╠═ba74ec4a-b38f-4fcb-9b43-2448976153eb
# ╠═ab56310f-671d-4b54-b9c2-9be102c826e3
# ╠═db20d906-b22f-4f16-b355-eba9c82809d7
# ╟─b7f11d85-6859-4b5a-814e-003946fbd894
# ╠═96b8f1e6-9159-4616-be08-c88f17e9829e
# ╠═53af3f77-cbc2-4809-8dda-743e096df10c
# ╟─4c136f35-e489-4d6e-83e0-e153bfc7040d
# ╠═bf97497a-e149-4552-9bdd-9fbeb5c24987
# ╠═b3528818-74f7-471b-8d33-ad7b9fcfc70b
# ╠═1ff22cf9-ac3b-4849-bd4c-8b0d3a901447
# ╠═c091d147-4874-4a19-9855-c78ded59b67a
# ╟─c922f283-3b31-43bd-9855-933be8a9e095
# ╠═460a8ea2-d077-48fc-804e-7367ac16eff8
# ╠═ba0e90f8-a310-40f9-a306-fb8700ea473f
# ╠═5d263158-03b0-4c87-ad90-a81c24493b66
# ╟─515d2174-56e3-4ebc-a475-b84f911a5a2f
# ╠═a97d7c94-e69b-44cd-ae4c-61d8d86808a5
# ╠═9e64def4-b397-4292-8fe9-8a828e08ee0a
# ╟─54fa3a01-5a8d-4dd3-a521-a60a2f69b541
# ╠═d9c33b8c-6ccf-462f-b5c2-267fd7b93ba3
# ╟─2f14fb93-91bc-44af-8a6d-cebaa1607a51
# ╠═fd53c66d-ba90-4d96-9a11-eabc1dd666dc
# ╠═d7091577-a5f9-4edf-b96f-88b13d9c6f3c
# ╠═19546bad-d7be-4c09-a248-f7e98c53d42d
# ╟─33c156b6-5c18-4377-9920-32c7865898e9
# ╟─0cdf07a5-a87e-449e-ab09-871924cee6d1
# ╠═0a729667-599a-41f1-98f3-ff94608c15bb
# ╠═0d78c765-1912-4ded-a5e8-aa98afe73cb4
# ╠═5566b0c6-ac96-4342-a62b-0c292469c2eb
# ╟─3bed0c72-ab5c-4dbf-85a5-6f6d062dec97
# ╠═d32536af-b369-4b39-87f6-4e1172d971ae
# ╟─35cdd4f9-87f7-408e-a3fd-7456b7238dd5
# ╠═1320b506-d378-4fa1-bbf1-c2f994abfff1


text:
module GPPerformance

export run_study
export GPNvsTime, GPNvsTimeResult
export GP2NvsTime, GP2NvsTimeResult

using GaussianProcesses, GaussianProcesses2
using CPUTime, Plots, DataFrames
using Parameters

@with_kw struct GPNvsTime
    ds::Vector{Int} = [2,3]
    ns::Vector{Int} = [10,20,50,100,200,500,700,1000,1500,2000]
    seeds::Vector{Int} = collect(1:5)
end

struct GPNvsTimeResult
    data::DataFrame
end

function run_study(study::GPNvsTime)
    mzero = MeanZero()
    kern = SE(0.0,0.0)
    v = 0.02
    logobsnoise = log(sqrt(v))

    ####
    # make sure things are compiled before timing
    X = rand(2, 5)
    y = rand(5)
    gp = GP(X, y, mzero, kern, logobsnoise)
    predict_f(gp, X)
    ####

    df = DataFrame([Int,Int,Int,Float64],[:dims,:n,:seed,:cputime_us],0)
    for d in study.ds, n in study.ns, s in study.seeds
        rng = MersenneTwister(s)
        X = rand(rng, d, n)
        y = rand(rng, n)
        tstart = CPUtime_us()
        gp = GP(X, y, mzero, kern, logobsnoise)
        predict_f(gp, X)
        t = CPUtime_us() - tstart
        push!(df, [d, n, s, t])
    end
    GPNvsTimeResult(df)
end

@with_kw struct GP2NvsTime
    ds::Vector{Int} = [2,3]
    ns::Vector{Int} = [10,20,50,100,200,500,700,1000,1500,2000]
    seeds::Vector{Int} = collect(1:5)
end

struct GP2NvsTimeResult
    data::DataFrame
end

function run_study(study::GP2NvsTime)
    mzero = ZeroMean()
    kern = SquaredExponential(1.0)
    v = 0.02

    ####
    # make sure things are compiled before timing
    X = [rand(2) for _ = 1:5]
    y = rand(5)
    gp = GaussianProcess(mzero, kern, X, y, v)
    GaussianProcesses2.predict(gp, X)
    ####

    df = DataFrame([Int,Int,Int,Float64],[:dims,:n,:seed,:cputime_us],0)
    for d in study.ds, n in study.ns, s in study.seeds
        rng = MersenneTwister(s)
        X = [rand(rng, d) for _ in 1:n]
        y = rand(rng, n)
        tstart = CPUtime_us()
        gp = GaussianProcess(mzero, kern, X, y, v)
        GaussianProcesses2.predict(gp, X)
        t = CPUtime_us() - tstart
        push!(df, [d, n, s, t])
    end
    GP2NvsTimeResult(df)
end

@recipe function plot(result::Union{GPNvsTimeResult,GP2NvsTimeResult})
    xlabel := "number of points"
    ylabel := "cpu time (us)"
    df = aggregate(result.data, [:dims,:n], [mean,std,length])
    for dd in groupby(df, :dims)
        @series begin
            seriestype := :path
            err := dd[:cputime_us_std] / dd[1,:seed_length]
            label := "$(string(result))-$(dd[1,:dims])d"
            dd[:n], dd[:cputime_us_mean]
        end
    end
end
string(result::GPNvsTimeResult) = "GP1"
string(result::GP2NvsTimeResult) = "GP2"

@recipe function plot(result1::GPNvsTimeResult, result2::GP2NvsTimeResult)
   @series begin
       result1
   end
   @series begin
       result2
   end
end

end


text:
# coding: utf-8

# In[1]:

#another fractal script via recursion
#this script heavily involves analytic geometry
#a good material to help you rewind high school geometry
# https://youthconway.com/handouts/analyticgeo.pdf
using Plots

# In[2]:

#simple solution to get coefficients of the equation
function get_line_params(x1,y1,x2,y2)

    slope=(y1-y2)/(x1-x2)
    intercept=y1-slope*x1

    return slope,intercept

end

#compute euclidean distance
function euclidean_distance(point1,point2)
    return √((point1[1]-point2[1])^2+(point1[2]-point2[2])^2)
end

#standard solution to quadratic equation
function solve_quadratic_equation(A,B,C)
    x1=(-B+√(B^2-4*A*C))/(2*A)
    x2=(-B-√(B^2-4*A*C))/(2*A)
    return [x1,x2]
end

#analytic geometry to compute target datapoints
function get_datapoint(pivot,measure,length,direction="inner")

    #for undefined slope
    if pivot[1]==measure[1]
        y1=pivot[2]+length
        y2=pivot[2]-length
        x1=pivot[1]
        x2=pivot[1]

    #for general cases
    else

        #get line equation
        slope,intercept=get_line_params(pivot[1],pivot[2],
                                   measure[1],measure[2],)

        #solve quadratic equation
        A=1
        B=-2*pivot[1]
        C=pivot[1]^2-length^2/(slope^2+1)
        x1,x2=solve_quadratic_equation(A,B,C)

        #get y from line equation
        y1=slope*x1+intercept
        y2=slope*x2+intercept

    end

    if direction=="inner"

        #take the one between pivot and measure points
        if euclidean_distance((x1,y1),measure)<euclidean_distance((x2,y2),measure)
            datapoint=(x1,y1)
        else
            datapoint=(x2,y2)
        end

    else

        #take the one farther away from measure points
        if euclidean_distance((x1,y1),measure)>euclidean_distance((x2,y2),measure)
            datapoint=(x1,y1)
        else
            datapoint=(x2,y2)
        end

    end

    return datapoint

end

#recursively compute the coordinates of koch curve data points
#to effectively connect the data points,the best choice is to use turtle
#it would be too difficult to connect the dots via analytic geometry
function koch_snowflake(base1,base2,base3,n,arr)

    #base case
    if n==0
        return

    else

        #find mid point
        #midpoint between base1 and base2 has to satisfy two conditions
        #it has to be on the same line as base1 and base2
        #assume this line follows y=kx+b
        #the midpoint is (x,kx+b)
        #base1 is (α,kα+b),base2 is (δ,kδ+b)
        #the euclidean distance between midpoint and base1 should be
        #half of the euclidean distance between base1 and base2
        #(x-α)^2+(kx+b-kα-b)^2=((α-δ)^2+(kα+b-kδ-b)^2)/4
        #apart from x,everything else in the equation is constant
        #this forms a simple quadratic solution to get two roots
        #one root would be between base1 and base2 which yields midpoint
        #and the other would be farther away from base2
        #this function solves the equation via (-B+(B^2-4*A*C)^0.5)/(2*A)
        #alternatively,you can use scipy.optimize.root
        #the caveat is it does not offer both roots
        #a wrong initial guess could take you to the wrong root
        midpoint=get_datapoint(base1,base2,euclidean_distance(base1,base2)/2)

        #compute the top point of a triangle
        #the computation is similar to midpoint
        #the euclidean distance between triangle_top and midpoint should be
        #one third of the distance between base3 and midpoint
        triangle_top=get_datapoint(midpoint,base3,
                                   euclidean_distance(midpoint,base3)/3,
                                   "outer")

        #two segment points divide a line into three equal parts
        #the computation is almost the same as midpoint
        #the euclidean distance between segment1 and base1
        #should be one third of the euclidean distance between base2 and base1
        segment1=get_datapoint(base1,base2,euclidean_distance(base1,base2)/3)
        segment2=get_datapoint(base2,base1,euclidean_distance(base1,base2)/3)

        #compute the nearest segment point of the neighboring line
        segment_side_1=get_datapoint(base1,base3,euclidean_distance(base1,base3)/3)
        segment_side_2=get_datapoint(base2,base3,euclidean_distance(base2,base3)/3)

        append!(arr,[segment1,segment2,triangle_top])

        #recursion
        koch_snowflake(base1,segment1,segment_side_1,n-1,arr)
        koch_snowflake(segment1,triangle_top,segment2,n-1,arr)
        koch_snowflake(triangle_top,segment2,segment1,n-1,arr)
        koch_snowflake(segment2,base2,segment_side_2,n-1,arr)
    end
    arr
end

# In[3]:

#set data points
point1=(0.0,0.0)
point2=(3.0,0.0)
point3=(1.5,1.5*√(3))

n=4

# In[4]:

#collect coordinates
coordinates=[point1,point2,point3]
append!(coordinates,koch_snowflake(point1,point2,point3,n,[]))
append!(coordinates,koch_snowflake(point3,point1,point2,n,[]))
append!(coordinates,koch_snowflake(point2,point3,point1,n,[]));

# In[5]:

#viz
#visually u can tell some of the data points are miscalculated
#purely caused by floating point arithmetic
gr(size=(250,250))
scatter([i[1] for i in coordinates],
        [i[2] for i in coordinates],markersize=2,
        legend=false,showaxis=false,ticks=false,grid=false)

max_stars_repo_path: koch snowflake.jl

max_issues_repo_path: koch snowflake.jl

max_forks_repo_path: koch snowflake.jl


text:
include("SWAP.jl")
include("SCAL.jl")
include("SSCAL.jl")
include("COPY.jl")
include("AXPY.jl")
include("DOT.jl")
include("ROTG.jl")
include("ROTMG.jl")
include("ROT.jl")
include("ROTM.jl")
include("NRM2.jl")
include("iAMAX.jl")
include("ASUM.jl")


text:
using OPFInstances

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

@testset "getinstancepath tests" begin
    @test ispath(getinstancepath("Matpower", "m", "case9"))
    @test ispath(getinstancepath("Matpower", "QCQP", "case9"))
    @test ispath(getinstancepath("GOC", "Phase_0_IEEE14", "scenario_4"))
end


text:
module RepeatArrays

export myrepeat

import Base.size

struct RepeatArray{T}
    A::Array{T}

    # Array of tuples (dim, inner, outer), giving a number of inner and outer repetitions
    # and the dimension to repeat along.
    #
    # An inner repetition repeats values one at a time. For example, an inner repetition
    # of 2 turns [a, b, c] into [a, a, b, b, c, c].
    #
    # An outer repetition repeats the whole array. For example, an outer repetition of
    # 2 turns [a, b, c] into [a, b, c, a, b, c].
    #
    # This is guaranteed to be "complete" and sorted by dimension.
    # By complete we mean that it contains elements (dim, 1, 1) as needed
    # to fill in entries until the maximum value of dim.
    repetitions::Array{Tuple{Int, Int, Int}}

    maxdim::Int
end

Base.eltype(R::RepeatArray{T}) where {T} = T

Base.similar(R::RepeatArray{T}, eltype) where {T} = RepeatArray(similar(R.A, eltype), copy(R.repetitions), R.maxdim)

myrepeat(R::RepeatArray{T}; inner=[], outer=[], dims=nothing) where {T} = begin
    repetitions, maxdim = compute_repetitions(basearray=R, ndims=length(size(R)), inner=inner, outer=outer, dims=dims)
    fused_repetitions = []
    for i in 1:max(R.maxdim, maxdim)
        repspec = get(repetitions, i, (i, 1, 1))
        inner_reps = repspec[2] * get_inner_reps(R, i)
        outer_reps = repspec[3] * get_outer_reps(R, i)
        push!(fused_repetitions, (i, inner_reps, outer_reps))
    end
    RepeatArray{T}(R.A, fused_repetitions, maxdim)
end

myrepeat(A::Array{T, N}; inner=[], outer=[], dims=nothing) where {T, N} = begin
    repetitions, maxdim = compute_repetitions(basearray=A, ndims=N, inner=inner, outer=outer, dims=dims)
    RepeatArray{T}(A, repetitions, maxdim)
end

# Cheat for now until we just override repeat
Base.repeat(R::RepeatArray{T}; inner=[], outer=[]) where {T} = myrepeat(R, inner=inner, outer=outer)

compute_repetitions(; basearray, ndims::Int, inner=[], outer=[], dims=nothing) = begin
    inner = [inner...]
    outer = [outer...]
    if !isnothing(dims)
        length(dims) == length(unique(dims)) || throw(ArgumentError("Repeated dimension specifier"))
        length(inner) == length(outer) || throw(ArgumentError("If dimensions are specified, inner and outer specifications must match in length. Use 1 to indicate no repetition."))
        length(dims) == length(inner) || throw(ArgumentError("Number of dimension specifiers must match number of repetition specifiers"))
    end

    if length(inner) < length(outer)
        append!(inner, ones(length(outer) - length(inner)))
    end
    if length(outer) < length(inner)
        append!(outer, ones(length(inner) - length(outer)))
    end

    dims = if isnothing(dims)
        1:max(length(inner), length(outer))
    else
        dims
    end

    maxdim = max(maximum(dims), ndims)
    repetitions = Tuple{Int, Int, Int}[]

    for d in 1:maxdim
        if d in dims
            inner_rep = inner[findfirst(x -> x == d, dims)]
            outer_rep = outer[findfirst(x -> x == d, dims)]

            # Make sure we treat "inner" repetitions of a singleton dimension as outer,
            # in order to take better advantage of repetition structure during broadcasting.
            if size(basearray, d) == 1
                outer_rep *= inner_rep
                inner_rep = 1
            end
            push!(repetitions, (d, inner_rep, outer_rep))
        else
            push!(repetitions, (d, 1, 1))
        end
    end

    [repetitions, maxdim]
end

size(R::RepeatArray{T}) where {T} = begin
    result = Int[]
    for (dim, inner, outer) in R.repetitions
        if dim <= length(size(R.A))
            push!(result, size(R.A)[dim] * inner * outer)
        else
            push!(result, inner * outer)
        end
    end
    tuple(result...)
end

size(R::RepeatArray{T}, i::Int) where {T} = begin
    s = size(R)
    if i > length(s)
        return 1
    else
        return s[i]
    end
end

get_outer_reps(R::RepeatArray{T}, dim::Int; broadcasting_to=nothing) where {T} = begin
    if dim > length(R.repetitions)
        if isnothing(broadcasting_to)
            return 1
        else
            return broadcasting_to
        end
    else
        return R.repetitions[dim][3]
    end
end

get_inner_reps(R::RepeatArray{T}, dim::Int) where {T} = begin
    if dim > length(R.repetitions)
        return 1
    else
        return R.repetitions[dim][2]
    end
end

Base.collect(R::RepeatArray{T}) where {T} = begin
    if size(R) == size(R.A)
        return R.A
    else
        return collect_by_repeating(R, R.A)
    end
end

collect_by_repeating(R::RepeatArray{T}, A::Array{T, N}) where {T, N} = begin
    missing_dims = length(size(R)) - length(size(A))
    dims = [size(A)...]
    append!(dims, Array{Int}(ones(missing_dims)))

    reshaped = reshape(A, dims...)

    Base.repeat(reshaped, inner=[rep[2] for rep in R.repetitions], outer=[rep[3] for rep in R.repetitions])
end

###########
# Broadcasting
###########

struct RepeatStyle <: Broadcast.BroadcastStyle end

Base.BroadcastStyle(::Type{<:RepeatArray}) = RepeatStyle()

Base.Broadcast.broadcastable(R::RepeatArray) = R

Base.BroadcastStyle(::RepeatStyle, ::Broadcast.DefaultArrayStyle{N}) where {N} = RepeatStyle()

#########
# Unary operators
#########

Base.Broadcast.broadcasted(::RepeatStyle, ::typeof(*), x::Number, R::RepeatArray{T}) where {T} = RepeatArray(x .* R.A, R.repetitions, R.maxdim)
Base.Broadcast.broadcasted(::RepeatStyle, ::typeof(*), R::RepeatArray{T}, x::Number) where {T} = RepeatArray(R.A .* x, R.repetitions, R.maxdim)

Base.Broadcast.broadcasted(::RepeatStyle, ::typeof(/), x::Number, R::RepeatArray{T}) where {T} = RepeatArray(R.A ./ x, R.repetitions, R.maxdim)
Base.Broadcast.broadcasted(::RepeatStyle, ::typeof(/), R::RepeatArray{T}, x::Number) where {T} = RepeatArray(x ./ R.A, R.repetitions, R.maxdim)

Base.Broadcast.broadcasted(::RepeatStyle, ::typeof(Base.literal_pow), ::Base.RefValue{typeof(^)}, R::RepeatArray{T}, v::Base.RefValue{Val{x}}) where {T, InnerType, x} =
    RepeatArray(Base.literal_pow.(^, R.A, v), R.repetitions, R.maxdim)
Base.Broadcast.broadcasted(::RepeatStyle, ::typeof(^), R::RepeatArray{T}, x::Number) where {T} = RepeatArray(R.A .^ x, R.repetitions, R.maxdim)

binaryBroadcast(f::Any, R::RepeatArray{T1}, A::AbstractArray{T2}) where {T1, T2} = begin
    output_size = Base.Broadcast._bcs(size(R), size(A))
    plan = planAbstractArrayBroadcast(R, A)

    intermediate = RepeatArray(R.A, plan.pre_application_reps, length(output_size))

    applied = broadcast(f, collect(intermediate), A)
    return RepeatArray(applied, plan.post_application_reps, length(output_size))
end

binaryBroadcast(f::Any, A::AbstractArray{T1}, R::RepeatArray{T2}) where {T1, T2} = binaryBroadcast((r, a) -> f(a, r), R, A)

Base.Broadcast.broadcasted(::RepeatStyle, ::typeof(+), R::RepeatArray{T1}, A::AbstractArray{T2}) where {T1, T2} = binaryBroadcast(+, R, A)
Base.Broadcast.broadcasted(::RepeatStyle, ::typeof(-), R::RepeatArray{T1}, A::AbstractArray{T2}) where {T1, T2} = binaryBroadcast(-, R, A)
Base.Broadcast.broadcasted(::RepeatStyle, ::typeof(*), R::RepeatArray{T1}, A::AbstractArray{T2}) where {T1, T2} = binaryBroadcast(*, R, A)
Base.Broadcast.broadcasted(::RepeatStyle, ::typeof(/), R::RepeatArray{T1}, A::AbstractArray{T2}) where {T1, T2} = binaryBroadcast(/, R, A)
Base.Broadcast.broadcasted(::RepeatStyle, ::typeof(+), A::AbstractArray{T1}, R::RepeatArray{T2}) where {T1, T2} = binaryBroadcast(+, R, A)
Base.Broadcast.broadcasted(::RepeatStyle, ::typeof(-), A::AbstractArray{T1}, R::RepeatArray{T2}) where {T1, T2} = binaryBroadcast(-, R, A)
Base.Broadcast.broadcasted(::RepeatStyle, ::typeof(*), A::AbstractArray{T1}, R::RepeatArray{T2}) where {T1, T2} = binaryBroadcast(*, R, A)
Base.Broadcast.broadcasted(::RepeatStyle, ::typeof(/), A::AbstractArray{T1}, R::RepeatArray{T2}) where {T1, T2} = binaryBroadcast(/, R, A)

struct AbstractArrayBroadcastPlan
    pre_application_reps::Array{Tuple{Int, Int, Int}}
    post_application_reps::Array{Tuple{Int, Int, Int}}
end

planAbstractArrayBroadcast(R::RepeatArray{T1}, A::AbstractArray{T2}) where {T1, T2} = begin
    output_size = Base.Broadcast._bcs(size(R), size(A))

    pre_application_reps = Tuple{Int, Int, Int}[]
    post_application_reps = Tuple{Int, Int, Int}[]

    for i in 1:length(output_size)
        if size(R, i) == 1 && size(A, i) == 1
            push!(pre_application_reps, (i, 1, 1))
            push!(post_application_reps, (i, 1, 1))
        end
        # Need to repeat this dimension before application of f
        if size(R, i) == 1 && size(A, i) > 1
            push!(pre_application_reps, (i, 1, size(A, i)))
            push!(post_application_reps, (i, 1, 1))
        end
        # Can delay repetition of this dim until after application of f
        if size(R, i) > 1 && size(A, i) == 1
            push!(pre_application_reps, (i, 1, 1))
            # There's no risk of this array access failing because R must have a repetition
            # specification for every dimension wider than a singleton
            push!(post_application_reps, R.repetitions[i])
        end
        if size(R, i) > 1 && size(A, i) > 1
            push!(pre_application_reps, R.repetitions[i])
            push!(post_application_reps, (i, 1, 1))
        end
    end

    return AbstractArrayBroadcastPlan(pre_application_reps, post_application_reps)
end

#########
# Broadcasting between RepeatArrays
#########

struct BinaryBroadcastPlan
    outer_reps::Array{Int}
    inner_reps::Array{Int}
    arg1Expansions::Array{Int}
    arg1Elongations::Array{Int}
    arg2Expansions::Array{Int}
    arg2Elongations::Array{Int}
end

planBinaryBroadcast(R::RepeatArray{T1}, S::RepeatArray{T2}) where {T1, T2} = begin
    output_size = Base.Broadcast._bcs(size(R), size(S))

    outer_reps = Int[]
    for i in 1:length(output_size)
        push!(outer_reps, gcd(get_outer_reps(R, i, broadcasting_to=output_size[i]),
                              get_outer_reps(S, i, broadcasting_to=output_size[i])))
    end

    inner_reps = Int[]

    R_expansions = Int[]
    R_elongations = Int[]
    S_expansions = Int[]
    S_elongations = Int[]

    for i in 1:length(output_size)
        chunk_size = Int(output_size[i] / outer_reps[i])

        inner_reps1 = get_inner_reps(R, i)
        inner_reps2 = get_inner_reps(S, i)

        inner_rep = gcd(inner_reps1, inner_reps2)
        base_dim = Int(chunk_size / inner_rep)

        push!(inner_reps, inner_rep)

        push!(R_expansions, Int(inner_reps1 / inner_rep))
        push!(R_elongations, Int(chunk_size / (size(R.A, i) * inner_reps1)))
        push!(S_expansions, Int(inner_reps2 / inner_rep))
        push!(S_elongations, Int(chunk_size / (size(S.A, i) * inner_reps2)))
    end

    BinaryBroadcastPlan(outer_reps, inner_reps, R_expansions, R_elongations, S_expansions, S_elongations)
end

doBinaryBroadcast(f::Any, R::RepeatArray{T1}, S::RepeatArray{T2}) where {T1, T2} = begin
    plan = planBinaryBroadcast(R, S)
    arg1 = repeat(R.A, inner=plan.arg1Expansions, outer=plan.arg1Elongations)
    arg2 = repeat(S.A, inner=plan.arg2Expansions, outer=plan.arg2Elongations)
    A = broadcast(f, arg1, arg2)
    myrepeat(A, inner=plan.inner_reps, outer=plan.outer_reps)
end

Base.Broadcast.broadcasted(::RepeatStyle, ::typeof(+), R::RepeatArray{T1}, S::RepeatArray{T2}) where {T1, T2} = doBinaryBroadcast(+, R, S)
Base.Broadcast.broadcasted(::RepeatStyle, ::typeof(-), R::RepeatArray{T1}, S::RepeatArray{T2}) where {T1, T2} = doBinaryBroadcast(-, R, S)
Base.Broadcast.broadcasted(::RepeatStyle, ::typeof(*), R::RepeatArray{T1}, S::RepeatArray{T2}) where {T1, T2} = doBinaryBroadcast(*, R, S)
Base.Broadcast.broadcasted(::RepeatStyle, ::typeof(/), R::RepeatArray{T1}, S::RepeatArray{T2}) where {T1, T2} = doBinaryBroadcast(/, R, S)

###########
# Linear Algebra
###########

Base.:*(R::RepeatArray{T1}, S::AbstractArray{T2, N}) where {T1, T2, N} = collect(R) * S
Base.:*(S::AbstractArray{T2, N}, R::RepeatArray{T1}) where {T1, T2, N} = S * collect(R)

end


text:
@testset "Molecule" begin

    molstringA = """
    C       -2.1315511243      2.2861688237      0.0000000000
    H       -1.0615511243      2.2861688237      0.0000000000
    H       -2.4882139062      1.4081046164      0.4966839113
    H       -2.4882187621      2.2950594327     -1.0087661539
    H       -2.4882200570      3.1553408443      0.5120813130"""

    mol = Fermi.Molecule(molstring=molstringA, unit=:angstrom, charge=0, multiplicity=1)
    @test begin
        length(mol.atoms) == 5 &&
        mol.charge == 0 &&
        mol.multiplicity == 1 &&
        mol.Nα == mol.Nβ == 5 &&
        mol.Vnuc ≈ 13.6865180253
    end

    molstringB = """
    C       -6.7482761629      2.9988214333      0.0000000000
    H       -5.0743010691      4.2609493569      0.0000000000
    H       -8.4961890559      4.1563877616      0.0000000000
    """
    mol = Fermi.Molecule(molstring=molstringB, unit=:bohr, charge=0, multiplicity=3)
    @test begin
        length(mol.atoms) == 3 &&
        mol.charge == 0 &&
        mol.multiplicity == 3 &&
        mol.Nα == 5 &&
        mol.Nβ == 3 &&
        mol.Vnuc ≈  6.0160249492
    end

    # Test wrong charge, mult and unit
    @test_throws Fermi.Options.FermiException Fermi.Molecule(molstring=molstringA, charge=10)
    @test_throws Fermi.Options.FermiException Fermi.Molecule(molstring=molstringA, charge=0, multiplicity=2)

    # Test print molecule
    mol = Fermi.Molecule(molstring=molstringB, unit=:bohr, charge=0, multiplicity=3)

    x = Fermi.string_repr(mol)

    @test begin
        occursin("Molecule:", x)          &&
        occursin(r"C\s+?[+-]{0,1}\d+?\.\d+?\s+?[+-]{0,1}\d+?\.\d+?\s+[+-]{0,1}\d+?\.\d+?",x) &&
        occursin(r"H\s+?[+-]{0,1}\d+?\.\d+?\s+?[+-]{0,1}\d+?\.\d+?\s+[+-]{0,1}\d+?\.\d+?",x) &&
        occursin(r"Charge:\s+?0", x)      &&
        occursin(r"Multiplicity:\s+?", x) &&
        occursin(r"Nuclear\s+?repulsion:\s+?\d+?\.\d+?", x)
    end
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
#= This is the assembly language parser
Created by: Kino Roy for CMPT 276, Assignment 1. =#
module AssemblyParser

function Parse(arguments)
  #=-----------------------------------------
  Initializing data structures for holding data
  ------------------------------------------=#
  instructions = Array{AbstractString}(2048,4)
  labelDict = Dict() #Creates the dictionary for labels

  global pathToInstructions = "instructions.txt"
if(size(arguments)[1]>0)
  pathToInstructions = arguments[1]
end

  fill!(instructions, "")

  lines = readFile()

  counter = 1 #Start the counter

  for l in lines
    if(l[1] == ";") #Allows starting the file with comment only lines
      continue
    end
    l = split(l,r";")[1] #eliminate the comments (garbage)
    #Determines type of arguments in the line

    splitLine = split(l,r":| |,",keep=false) #Filters text

    for i in 1:length(splitLine) #Brings all instructions to uppercase
      splitLine[i] = uppercase(splitLine[i])
    end

  #---- Begin Parsing ----#

    #CASE: "INST"
    if size(splitLine)[1] == 1

      instructions[counter,2] = chomp(splitLine[1])

    #CASE: "LABEL + INST"
    elseif size(splitLine)[1] == 2 && contains(l,":")
      instructions[counter,1] = splitLine[1]
      instructions[counter,2] = chomp(splitLine[2])
      #ADDS LABEL TO DICTIONARY
      labelDict["$(splitLine[1])"] = counter

    #CASE: "INST + OP1"
    elseif size(splitLine)[1] == 2
      instructions[counter,2] = splitLine[1]
      instructions[counter,3] = chomp(splitLine[2])
      #ADDS LABEL TO DICTIONARY
      labelDict["$(splitLine[1])"] = counter

    #CASE: "LABEL + INST + OP1"
    elseif size(splitLine)[1] == 3 && contains(l,":")
      instructions[counter,1] = splitLine[1]
      instructions[counter,2] = splitLine[2]
      instructions[counter,3] = chomp(splitLine[3])
      #ADDS LABEL TO DICTIONARY
      labelDict["$(splitLine[1])"] = counter

    #CASE: "INST + OP1 + OP2"
    elseif size(splitLine)[1] == 3 && contains(l,",")
      instructions[counter,2] = splitLine[1]
      instructions[counter,3] = splitLine[2]
      instructions[counter,4] = chomp(splitLine[3])

    #CASE: "LABEL + INST + OP1 + OP2"
    elseif size(splitLine)[1] == 4
      instructions[counter,1] = splitLine[1]
      instructions[counter,2] = splitLine[2]
      instructions[counter,3] = splitLine[3]
      instructions[counter,4] = chomp(splitLine[4])
      #ADDS LABEL TO DICTIONARY
      labelDict["$(splitLine[1])"] = counter
    end

     counter += 1
  end #endfor

  return instructions,labelDict
end

function readFile()
  global pathToInstructions
  stream = open(pathToInstructions) #Open the input file and set the input stream
  lines = readlines(stream) # Create an array with each element, an instruction
  close(stream)

  return lines
end

end
#=------------
END ASSEMBLY PARSER
-------------=#

max_stars_repo_name: kinoroy/KCPSM-Simulator

max_issues_repo_name: kinoroy/KCPSM-Simulator

max_forks_repo_name: kinoroy/KCPSM-Simulator


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -13.266290625856113 <= q <= 11.094199160852137)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.2977454197977769 + 0.7722469146312885*1/(1 + exp(-(0.9100824262656433 + 0.0704650714228956*1/(1 + exp(-(-0.7110372864748489 + -0.3189000958940382*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.039301419137926175*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + -0.19790129907631693*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + 0.8966302258455006*1/(1 + exp(-(0.7299768289984225 + 0.7045346620386517*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.8068588624409387*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + 0.8888276792068996*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + 0.7890296052704904*1/(1 + exp(-(-0.5099367345236407 + -0.3810722083004521*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.4850891682061631*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + -0.158929561309892*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4])))))))))) + 0.5815923034237733*1/(1 + exp(-(-0.15598180224250147 + 0.7627892034632455*1/(1 + exp(-(-0.7110372864748489 + -0.3189000958940382*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.039301419137926175*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + -0.19790129907631693*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + 0.6806349200893864*1/(1 + exp(-(0.7299768289984225 + 0.7045346620386517*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.8068588624409387*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + 0.8888276792068996*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + -0.163059274878961*1/(1 + exp(-(-0.5099367345236407 + -0.3810722083004521*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.4850891682061631*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + -0.158929561309892*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4])))))))))) + -0.731486602090937*1/(1 + exp(-(-0.9577813343303001 + 0.660774759467007*1/(1 + exp(-(-0.7110372864748489 + -0.3189000958940382*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.039301419137926175*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + -0.19790129907631693*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + -0.16269598559087406*1/(1 + exp(-(0.7299768289984225 + 0.7045346620386517*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.8068588624409387*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + 0.8888276792068996*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + 0.9308521240523335*1/(1 + exp(-(-0.5099367345236407 + -0.3810722083004521*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.4850891682061631*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + -0.158929561309892*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))))))))) + 1/(1 + exp(-(0.23088458854123672 + 0.9025271253252471*1/(1 + exp(-(0.9100824262656433 + 0.0704650714228956*1/(1 + exp(-(-0.7110372864748489 + -0.3189000958940382*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.039301419137926175*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + -0.19790129907631693*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + 0.8966302258455006*1/(1 + exp(-(0.7299768289984225 + 0.7045346620386517*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.8068588624409387*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + 0.8888276792068996*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + 0.7890296052704904*1/(1 + exp(-(-0.5099367345236407 + -0.3810722083004521*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.4850891682061631*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + -0.158929561309892*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4])))))))))) + 0.9612479731086991*1/(1 + exp(-(-0.15598180224250147 + 0.7627892034632455*1/(1 + exp(-(-0.7110372864748489 + -0.3189000958940382*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.039301419137926175*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + -0.19790129907631693*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + 0.6806349200893864*1/(1 + exp(-(0.7299768289984225 + 0.7045346620386517*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.8068588624409387*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + 0.8888276792068996*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + -0.163059274878961*1/(1 + exp(-(-0.5099367345236407 + -0.3810722083004521*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.4850891682061631*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + -0.158929561309892*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4])))))))))) + 0.9633061422344626*1/(1 + exp(-(-0.9577813343303001 + 0.660774759467007*1/(1 + exp(-(-0.7110372864748489 + -0.3189000958940382*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.039301419137926175*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + -0.19790129907631693*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + -0.16269598559087406*1/(1 + exp(-(0.7299768289984225 + 0.7045346620386517*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.8068588624409387*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + 0.8888276792068996*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + 0.9308521240523335*1/(1 + exp(-(-0.5099367345236407 + -0.3810722083004521*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.4850891682061631*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + -0.158929561309892*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))))))))) + 1/(1 + exp(-(-0.8299198997441475 + 0.7709117851978275*1/(1 + exp(-(0.9100824262656433 + 0.0704650714228956*1/(1 + exp(-(-0.7110372864748489 + -0.3189000958940382*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.039301419137926175*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + -0.19790129907631693*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + 0.8966302258455006*1/(1 + exp(-(0.7299768289984225 + 0.7045346620386517*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.8068588624409387*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + 0.8888276792068996*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + 0.7890296052704904*1/(1 + exp(-(-0.5099367345236407 + -0.3810722083004521*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.4850891682061631*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + -0.158929561309892*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4])))))))))) + 0.5338511667423695*1/(1 + exp(-(-0.15598180224250147 + 0.7627892034632455*1/(1 + exp(-(-0.7110372864748489 + -0.3189000958940382*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.039301419137926175*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + -0.19790129907631693*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + 0.6806349200893864*1/(1 + exp(-(0.7299768289984225 + 0.7045346620386517*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.8068588624409387*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + 0.8888276792068996*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + -0.163059274878961*1/(1 + exp(-(-0.5099367345236407 + -0.3810722083004521*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.4850891682061631*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + -0.158929561309892*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4])))))))))) + 0.3695524647159538*1/(1 + exp(-(-0.9577813343303001 + 0.660774759467007*1/(1 + exp(-(-0.7110372864748489 + -0.3189000958940382*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.039301419137926175*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + -0.19790129907631693*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + -0.16269598559087406*1/(1 + exp(-(0.7299768289984225 + 0.7045346620386517*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.8068588624409387*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + 0.8888276792068996*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + 0.9308521240523335*1/(1 + exp(-(-0.5099367345236407 + -0.3810722083004521*1/(1 + exp(-(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])))) + 0.4850891682061631*1/(1 + exp(-(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])))) + -0.158929561309892*1/(1 + exp(-(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
@reexport using Distributions
mutable struct Poisson_Bernoulli_VCObs{T <: BlasReal, VD, VL}
    # data
    y::Vector{T}
    X::Matrix{T}
    V::Vector{Matrix{T}}
    n::Int
    p::Int          # number of mean parameters in linear regression
    m::Int          # number of variance components
    # working arrays
    ∇β::Vector{T}   # gradient wrt β
    ∇resβ::Matrix{T}# residual gradient matrix d/dβ_p res_ij (each observation has a gradient of residual is px1)
    ∇τ::Vector{T}   # gradient wrt τ
    ∇θ::Vector{T}   # gradient wrt θ2
    Hβ::Matrix{T}   # Hessian wrt β
    Hθ::Matrix{T}   # Hessian wrt variance components θ
    Hτ::Matrix{T}   # Hessian wrt τ
    res::Vector{T}  # residual vector res_i
    t::Vector{T}    # t[k] = tr(V_i[k]) / 2
    q::Vector{T}    # q[k] = res_i' * V_i[k] * res_i / 2
    xtx::Matrix{T}  # Xi'Xi
    storage_n::Vector{T}
    m1::Vector{T}
    m2::Vector{T}
    storage_p1::Vector{T}
    storage_p2::Vector{T}
    storage_np::Matrix{T}
    storage_pp::Matrix{T}
    added_term_numerator::Matrix{T}
    added_term2::Matrix{T}
    η::Vector{T}    # η = Xβ systematic component
    μ::Vector{T}    # μ(β) = ginv(Xβ) # inverse link of the systematic component
    varμ::Vector{T} # v(μ_i) # variance as a function of the mean
    dμ::Vector{T}   # derivative of μ
    # d::D            # distribution()
    vecd::VD
    veclink::VL      # link function ()
    wt::Vector{T}   # weights wt for GLM.jl
    w1::Vector{T}   # working weights in the gradient = dμ/v(μ)
    w2::Vector{T}   # working weights in the information matrix = dμ^2/v(μ)
end

function Poisson_Bernoulli_VCObs(
    y::Vector{T},
    x::Vector{T},
    V::Vector{Matrix{T}},
    vecd::VD,
    veclink::VL) where {T <: BlasReal, VD, VL}
    X = [transpose(x) zeros(size(transpose(x))); zeros(size(transpose(x))) transpose(x)]
    n, p, m = size(X, 1), size(X, 2), length(V)
    @assert length(y) == n "length(y) should be equal to size(X, 1)"
    # working arrays
    ∇β  = Vector{T}(undef, p)
    ∇resβ  = Matrix{T}(undef, n, p)
    ∇τ  = Vector{T}(undef, 1)
    ∇θ  = Vector{T}(undef, m)
    Hβ  = Matrix{T}(undef, p, p)
    Hθ  = Matrix{T}(undef, m, m)
    Hτ  = Matrix{T}(undef, 1, 1)
    res = Vector{T}(undef, n)
    t   = [tr(V[k])/2 for k in 1:m]
    q   = Vector{T}(undef, m)
    xtx = transpose(X) * X
    storage_n = Vector{T}(undef, n)
    m1        = Vector{T}(undef, m)
    m2        = Vector{T}(undef, m)
    storage_p1 = Vector{T}(undef, p)
    storage_p2 = Vector{T}(undef, p)
    storage_np = Matrix{T}(undef, n, p)
    storage_pp = Matrix{T}(undef, p, p)
    added_term_numerator = Matrix{T}(undef, n, p)
    added_term2 = Matrix{T}(undef, p, p)
    η = Vector{T}(undef, n)
    μ = Vector{T}(undef, n)
    varμ = Vector{T}(undef, n)
    dμ = Vector{T}(undef, n)
    wt = Vector{T}(undef, n)
    fill!(wt, one(T))
    w1 = Vector{T}(undef, n)
    w2 = Vector{T}(undef, n)
    # constructor
    Poisson_Bernoulli_VCObs{T, VD, VL}(y, X, V, n, p, m, ∇β, ∇resβ, ∇τ, ∇θ, Hβ, Hθ,
        Hτ, res, t, q, xtx, storage_n, m1, m2, storage_p1, storage_p2, storage_np, storage_pp, added_term_numerator, added_term2, η, μ, varμ, dμ, vecd, veclink, wt, w1, w2)
end

"""
Poisson_Bernoulli_VCModel
Poisson_Bernoulli_VCModel(gcs)
Bivariate Mixed Poisson, Bernoulli variance component model, which contains a vector of
`Poisson_Bernoulli_VCObs` as data, model parameters, and working arrays.
"""
struct Poisson_Bernoulli_VCModel{T <: BlasReal, VD, VL} <: MathProgBase.AbstractNLPEvaluator
    # data
    data::Vector{Poisson_Bernoulli_VCObs{T, VD, VL}}
    Y1total::T
    Y2total::T
    ntotal::Int     # total number of singleton observations
    p::Int          # number of mean parameters in linear regression
    m::Int          # number of variance components
    # parameters
    β::Vector{T}    # p-vector of mean regression coefficients
    τ::Vector{T}    # inverse of linear regression variance parameter
    θ::Vector{T}
    # working arrays
    ∇β::Vector{T}   # gradient from all observations
    ∇τ::Vector{T}
    ∇θ::Vector{T}
    XtX::Matrix{T}  # X'X = sum_i Xi'Xi
    Hβ::Matrix{T}    # Hessian from all observations
    Hθ::Matrix{T}
    Hτ::Matrix{T}
    Ainv::Matrix{T}
    Aevec::Matrix{T}
    M::Matrix{T}
    vcov::Matrix{T}
    ψ::Vector{T}
    TR::Matrix{T}         # n-by-m matrix with tik = tr(Vi[k]) / 2
    QF::Matrix{T}         # n-by-m matrix with qik = res_i' Vi[k] res_i
    storage_n::Vector{T}
    storage_m::Vector{T}
    storage_θ::Vector{T}
    vecd::Vector{VD}
    veclink::Vector{VL}
    penalized::Bool
end

function Poisson_Bernoulli_VCModel(gcs::Vector{Poisson_Bernoulli_VCObs{T, VD, VL}}; penalized::Bool = false) where {T <: BlasReal, VD, VL}
    n, p, m = length(gcs), size(gcs[1].X, 2), length(gcs[1].V)
    β       = Vector{T}(undef, p)
    τ       = [1.0]
    θ       = Vector{T}(undef, m)
    ∇β      = Vector{T}(undef, p)
    ∇τ      = Vector{T}(undef, 1)
    ∇θ      = Vector{T}(undef, m)
    XtX     = zeros(T, p, p) # sum_i xi'xi
    Hβ      = Matrix{T}(undef, p, p)
    Hθ      = Matrix{T}(undef, m, m)
    Hτ      = Matrix{T}(undef, 1, 1)
    Ainv    = zeros(T, p + m, p + m)
    Aevec   = zeros(T, p + m, p + m)
    M       = zeros(T, p + m, p + m)
    vcov    = zeros(T, p + m, p + m)
    ψ       = Vector{T}(undef, p + m)
    TR      = Matrix{T}(undef, n, m) # collect trace terms
    Y1total  = 0.0
    Y2total  = 0.0
    ntotal  = 0.0
    vecd       = Vector{VD}(undef, n)
    veclink    = Vector{VL}(undef, n)
    for i in eachindex(gcs)
        ntotal  += length(gcs[i].y)
        Y1total  += sum(gcs[i].y[1])
        Y2total  += sum(gcs[i].y[2])
        BLAS.axpy!(one(T), gcs[i].xtx, XtX)
        TR[i, :] = gcs[i].t
        vecd[i] = gcs[i].vecd
        veclink[i] = gcs[i].veclink
    end
    QF        = Matrix{T}(undef, n, m)
    storage_n = Vector{T}(undef, n)
    storage_m = Vector{T}(undef, m)
    storage_θ = Vector{T}(undef, m)
    Poisson_Bernoulli_VCModel{T, VD, VL}(gcs, Y1total, Y2total, ntotal, p, m, β, τ, θ,
        ∇β, ∇τ, ∇θ, XtX, Hβ, Hθ, Hτ, Ainv, Aevec, M, vcov, ψ, TR, QF,
        storage_n, storage_m, storage_θ, vecd, veclink, penalized)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ e704feae-919c-11eb-35ea-0932dbb413dd
begin
	using Pkg; Pkg.activate("MLJ_env", shared=true)
	# Pkg.add("StatsBase")
	# Pkg.add("BenchmarkTools");

	using StatsBase
	using BenchmarkTools
	using Test
	using Random
	using PlutoUI
end

# ╔═╡ 45d72214-919c-11eb-2e69-23c0786baacb
md"""
### Finding index of a random minimum element in an array

Src: [Julia 1.0 Cookbook -  Bogumił Kamiński , Przemysław Szufel - 2018](https://www.packtpub.com/product/julia-1-0-programming-cookbook/9781788998369)

**Context**
In many applications, we need to find the index of the minimum element of some array. The built-in argmin function which is designed to perform this task returns the index of the minimum element in a collection. However, if there are multiple minimal elements, only the first one is returned.
There are situations when we need to get all the indices of a minimal element or a single index is chosen uniformly at random from this set.

Let us implement such a function in two ways.
"""

# ╔═╡ e849191c-919c-11eb-0d04-b9ec7a4ac478
function all_argmin(ary::AbstractArray)::Vector{Integer}
	@assert length(ary) > 0
	m = minimum(ary)
	filter(ix -> ary[ix] == m, eachindex(ary))
end

# ╔═╡ 08f24660-919e-11eb-37c9-155375b48cf4
rand_argmin₁(ary::AbstractArray)::Integer = all_argmin(ary) |> rand

# ╔═╡ 61c23c88-919d-11eb-0d37-23783e1fd854
begin
	@test_throws AssertionError all_argmin([])

	x = [3, 2, 1, 4, 1, 2, 3, 4, 4, 3, 2, 1, 2, 1, 3, 4]
	@test all_argmin([10, 10, 1, 2, 3, 5, 10]) == [3]
	@test all_argmin(x) == [3, 5, 12, 14]

	n = 10^6
	f = n ÷ 4
	@test rand_argmin₁([10, 10, 1, 2, 3, 5, 10]) == 3

	## 10^6 samples - should be roughly equally distributed 1/4 as x has f
	## our 1 (min element)
	@test all(map(p -> abs(p[2] - f) ≤ 1200,
			countmap([rand_argmin₁(x) for _ ∈1:10^6]) |> collect))
end

# ╔═╡ 4f47a640-919f-11eb-0072-e9855b338704
md"""
Let us define a more efficient (and generic function).
"""

# ╔═╡ 67982c78-91a0-11eb-29d7-613d4624cf20
function rand_argmin₂(ary::AbstractArray)::Integer
	indices = eachindex(ary)     ## take into account cartesian and linear indexing
	y = iterate(indices)
	y === nothing && throw(ArgumentError("collection must be non-empty"))

	(ix, state) = y              ## first elem iof iteratble, state
	minval, best_ix, best_cnt = ary[ix], ix, 1
	y = iterate(indices, state)
	while y !== nothing
		(ix, state) = y
		curval = ary[ix]

		if isless(curval, minval)    ## isless moree genral than <
			minval, best_ix = curval, ix
			best_cnt = 1             ## reset

		elseif isequal(curval, minval)
			best_cnt += 1
			## updated prob. of choosing this ix:
			rand() * best_cnt < 1. && (best_ix = ix)

		end
		y = iterate(indices, state)
	end
	best_ix
end

# ╔═╡ 56c60dc2-91a3-11eb-301d-1d0ce38e5cfa
md"""
Notes:
  - `y === nothing`       faster than `isnothing(y)`
  - `while y !== nothing` faster than `while !isnothing(y)`
"""

# ╔═╡ 3a4d6e52-91a2-11eb-0b87-bdc9ff551759
Random.seed!(42)

# ╔═╡ 89cfa838-91a1-11eb-2753-1d838efe223c
(countmap([rand_argmin₁(x) for _ ∈ 1:10^6]),
	countmap([rand_argmin₂(x) for _ ∈ 1:10^6]))

# ╔═╡ 4cc66048-91a2-11eb-2bf3-5b39967d17bb
begin
	const y = rand(1:1000, 10_000)

	with_terminal() do
		@btime rand_argmin₂(y)

		@btime rand_argmin₁(y)
	end
end

# ╔═╡ a15a8f10-91a4-11eb-3489-45611ec3c9d1
md"""
**Comments**

if we encounter a new minimum value, we simply store it. The challenge is what to do if we encounter a minimal element more than once.
Assume that we have already seen some element best_cnt times. A variable best_ix, with probability 1. / best_cnt, holds each previously encountered minimal index.
If we encounter a new index ix for which a[ix] is equal to the current minimum minval, then in the code with probability 1., / (bestcnt + 1), we set best_ix to this value. This means that this new value has the correct probability of being stored.
Also, notice that all earlier values have the following probability of being stored:

$$\frac{1.}{bestcnt} × \frac{bestcnt}{(bestcnt + 1)} = \frac{1.}{(bestcnt + 1)}$$

This means that, ultimately, our algorithm properly randomizes the choice of the returned index.

The second challenge is related to the fact that Julia features many different types of collections that can be efficiently accessed using different indexing, for example, linear indexing or Cartesian indexing.

Therefore, in the code, we traverse a collection of indices into `ary`, returned by `eachindex(ary)`, using an iteration protocol. The protocol is based on the `iterate` function. If we pass an iterable to this function (for example, iterate(indices)), it returns a tuple containing the following:
  - The first element of the iterable
  - A state variable that can be used in subsequent calls

Then, we can iterate over the collection by passing two arguments to iterate — a collection and an earlier returned state. Finally, we know that we have reached the end of the collection if iterate returns nothing.

Note that we have iterated over `eachindex(ary)` and not over `ary`, since we needed to work with the indices of collection a rather than the values contained in it.
"""

# ╔═╡ 9d3d8f1c-91a5-11eb-3f37-79ae44a88223
md"""
We use the `isless` and `isequal` functions in tests and not the `<` and `== ` functions. The main reason for this is that `isless` and `isequal` guarantee to always return a `Bool` result while `<` and `==`  can return missing when one of their arguments is missing.

Additionally, there are subtle differences between those functions when comparing `0.0` and `-0.0` and `NaN`:
"""

# ╔═╡ d73c9186-91a5-11eb-2bce-47367fc9522e
(.0 == -.0, -.0 < .0), (isequal(.0, -.0), isless(-.0, .0))

# ╔═╡ eda015e2-91a5-11eb-32af-2778d08859ed
(NaN == NaN, NaN < NaN), (isequal(NaN, NaN), isless(NaN, NaN))

# ╔═╡ Cell order:
# ╟─45d72214-919c-11eb-2e69-23c0786baacb
# ╠═e704feae-919c-11eb-35ea-0932dbb413dd
# ╠═e849191c-919c-11eb-0d04-b9ec7a4ac478
# ╠═08f24660-919e-11eb-37c9-155375b48cf4
# ╠═61c23c88-919d-11eb-0d37-23783e1fd854
# ╟─4f47a640-919f-11eb-0072-e9855b338704
# ╠═67982c78-91a0-11eb-29d7-613d4624cf20
# ╟─56c60dc2-91a3-11eb-301d-1d0ce38e5cfa
# ╠═3a4d6e52-91a2-11eb-0b87-bdc9ff551759
# ╠═89cfa838-91a1-11eb-2753-1d838efe223c
# ╠═4cc66048-91a2-11eb-2bf3-5b39967d17bb
# ╟─a15a8f10-91a4-11eb-3489-45611ec3c9d1
# ╟─9d3d8f1c-91a5-11eb-3f37-79ae44a88223
# ╠═d73c9186-91a5-11eb-2bce-47367fc9522e
# ╠═eda015e2-91a5-11eb-32af-2778d08859ed

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
export eigmax

function eigmax(A::Symmetric{T,<:BandedMatrix{T}}) where T <: Real
    maximum(KrylovKit.eigsolve(A, issymmetric=true)[1])
end


text:
__precompile__()

module RandomCorrelationMatrices

    using Distributions

    export randcormatrix
    """
        randcormatrix(d, η)

    Generate a d-by-d random correlation matrix using the
    method described in:

        Lewandowski, Kurowicka, Joe
          "Generating random correlation matrices based on vines
            and extended onion method"
        Journal of Multivariate Analysis 100 (2009)
        doi:10.1016/j.jmva.2009.04.008
    """
    function randcormatrix(d, η)
        β = η + (d-2)/2
        u = rand(Beta(β,β))
        r₁₂ = 2*u - 1
        r = [   1 r₁₂ ;
              r₁₂   1 ]
        # In published paper this index is n = 2:d-1
        # n is never mentioned again, but a mysterious k is
        # We'll thus use k
        # Similar approach was taken here:
        # http://stats.stackexchange.com/a/125017/58921
        for k in 2:d-1
            β = β - 1/2
            y = rand(Beta(k/2,β))
            u = randn(k)
            u /= norm(u)
            w = sqrt(y)*u
            A = chol(r, Val{:L})
            z = A*w
            r = [ r  z ;
                  z' 1 ]
        end
        return r
    end

    export randcovmatrix
    """
        randcovmatrix(d, η, σ)

    Use randcormatrix with desired vector of standard deviations σ
    for each term to generate a random covariance matrix.
    """
    function randcovmatrix(d, η, σ)
        length(σ) != d && throw(DimensionMismatch("length(σ) doesn't match d"))
        r = randcormatrix(d, η)
        Σ = zeros(d,d)
        @inbounds for i in 1:d, j in 1:d
            Σ[i,j] = r[i,j] * σ[i] * σ[j]
        end
        return Σ
    end

end # module