text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests of renderers.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Renderers - printf" begin
    matrix = Any[BigFloat(pi) float(pi) 10.0f0  Float16(1)
                 0x01         0x001     0x00001 0x000000001
                 true         false     true    false
                 "Teste" "Teste\nTeste" "Teste \"quote\" Teste" "Teste\n\"quote\"\nTeste"]

    header = ["C1" "C2" "C3" "C4"
              "S1" "S2" "S3" "S4"]

    row_names = [1,2,"3",'4']

    expected = """

    #
    Test
    C1
    C2
    C3
    C4

    S1
    S2
    S3
    S4

    1
    1
    3.141592653589793238462643383279502884197169399375105820974944592307816406286198
    3.14159
    10.0
    1.0

    2
    2
    1
    1
    1
    1

    3
    3
    true
    false
    true
    false

    4
    4
    Teste
    Teste
Teste
    Teste "quote" Teste
    Teste
"quote"
Teste

"""

    result = pretty_table(String, matrix, header;
                          backend = :html,
                          linebreaks = true,
                          row_names = row_names,
                          row_name_column_title = "Test",
                          row_number_column_title = "#",
                          show_row_number = true,
                          standalone = false)

    @test expected == result

    expected = """

    #
    Test
    C1
    C2
    C3
    C4

    S1
    S2
    S3
    S4

    1
    1
    3.141592653589793238462643383279502884197169399375105820974944592307816406286198
    3.141592653589793
    10.0
    1.0

    2
    2
    1
    1
    1
    1

    3
    3
    true
    false
    true
    false

    4
    4
    Teste
    Teste
Teste
    Teste "quote" Teste
    Teste
"quote"
Teste

"""

    result = pretty_table(String, matrix, header;
                          backend = :html,
                          compact_printing = false,
                          linebreaks = true,
                          row_names = row_names,
                          row_name_column_title = "Test",
                          row_number_column_title = "#",
                          show_row_number = true,
                          standalone = false)

    @test expected == result
end

@testset "Renderers - show" begin
    matrix = Any[BigFloat(pi) float(pi) 10.0f0  Float16(1)
                 0x01         0x001     0x00001 0x000000001
                 true         false     true    false
                 "Teste" "Teste\nTeste" "Teste \"quote\" Teste" "Teste\n\"quote\"\nTeste"]

    header = ["C1" "C2" "C3" "C4"
              "S1" "S2" "S3" "S4"]

    row_names = [1,2,"3",'4']

    expected = """

    #
    Test
    C1
    C2
    C3
    C4

    S1
    S2
    S3
    S4

    1
    1
    3.14159
    3.14159
    10.0
    1.0

    2
    2
    0x01
    0x0001
    0x00000001
    0x0000000000000001

    3
    3
    true
    false
    true
    false

    4
    '4'
    Teste
    Teste
Teste
    Teste "quote" Teste
    Teste
"quote"
Teste

"""

    result = pretty_table(String, matrix, header;
                          backend = :html,
                          linebreaks = true,
                          renderer = :show,
                          row_names = row_names,
                          row_name_column_title = "Test",
                          row_number_column_title = "#",
                          show_row_number = true,
                          standalone = false)

    @test expected == result

    expected = """

    #
    Test
    C1
    C2
    C3
    C4

    S1
    S2
    S3
    S4

    1
    1
    3.141592653589793238462643383279502884197169399375105820974944592307816406286198
    3.141592653589793
    10.0f0
    Float16(1.0)

    2
    2
    0x01
    0x0001
    0x00000001
    0x0000000000000001

    3
    3
    true
    false
    true
    false

    4
    '4'
    Teste
    Teste
Teste
    Teste "quote" Teste
    Teste
"quote"
Teste

"""

    result = pretty_table(String, matrix, header;
                          backend = :html,
                          compact_printing = false,
                          linebreaks = true,
                          renderer = :show,
                          row_names = row_names,
                          row_name_column_title = "Test",
                          row_number_column_title = "#",
                          standalone = false,
                          show_row_number = true)

    @test expected == result
end


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

name = "HAMdetector_model_binaries"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/HAMdetector/Escape.jl.git", "27c98e5cd34a015d66ab0898e49b50736acd89aa"),
    ArchiveSource("https://github.com/stan-dev/cmdstan/releases/download/v2.27.0/cmdstan-2.27.0.tar.gz", "ff71f4d255cf26c2d8366a8173402656a659399468871305056aa3d56329c1d5")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd cmdstan-2.27.0/

echo 'CXX=c++' >> make/local
echo 'TBB_CXX_TYPE=c++' >> make/local
echo 'STAN_THREADS=true' >> make/local

cd stan/lib/stan_math/lib/tbb_2020.3/build/
cp linux.inc x86_64-linux-gnu.inc
cp linux.gcc.inc x86_64-linux-gnu.c++.inc
cd ../../../../../..

make build
make runtime=cc5.2.0_libc2.12.2_kernel5.11.12177.current ../Escape.jl/models/model_4
make runtime=cc5.2.0_libc2.12.2_kernel5.11.12177.current ../Escape.jl/models/model_3
make runtime=cc5.2.0_libc2.12.2_kernel5.11.12177.current ../Escape.jl/models/model_2
make runtime=cc5.2.0_libc2.12.2_kernel5.11.12177.current ../Escape.jl/models/model_1

cd ../Escape.jl/models/

mkdir ${bindir}
mkdir ${libdir}

cp model_4 ${bindir}/model_4
cp model_3 ${bindir}/model_3
cp model_2 ${bindir}/model_2
cp model_1 ${bindir}/model_1

cd ../../cmdstan-2.27.0/stan/lib/stan_math/lib/tbb

cp libtbb.so.2 ${libdir}/libtbb.so.2
cp libtbbmalloc.so.2 ${libdir}/libtbbmalloc.so.2
cp libtbbmalloc_proxy.so.2 ${libdir}/libtbbmallo_proxy.so.2

install_license ${WORKSPACE}/srcdir/cmdstan-2.27.0/LICENSE

exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc = "glibc")
]

platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("model_2", :model_2),
    ExecutableProduct("model_1", :model_1),
    ExecutableProduct("model_3", :model_3),
    ExecutableProduct("model_4", :model_4)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"5.2.0")

ext: jl


text:
module Agents

using Requires
using Distributed
using DataStructures
using Graphs
using DataFrames
using Random
import ProgressMeter

import Base.length
import LinearAlgebra

# Core structures of Agents.jl
include("core/agents.jl")
include("core/model.jl")
include("core/space_interaction_API.jl")

# Existing spaces
include("spaces/nothing.jl")
include("spaces/graph.jl")
include("spaces/grid.jl")
include("spaces/discrete.jl")
include("spaces/continuous.jl")
include("spaces/openstreetmap.jl")
include("spaces/utilities.jl")

# Stepping and data collection functionality
include("simulations/step.jl")
include("simulations/collect.jl")
include("simulations/paramscan.jl")
include("simulations/sample.jl")
include("simulations/ensemblerun.jl")

# Other features that exist in submodules
include("submodules/pathfinding/all_pathfinders.jl")
include("submodules/schedulers.jl")
include("submodules/io/AgentsIO.jl")
include("deprecations.jl")

# Predefined models
include("models/Models.jl")
export Models

# Update messages:
using Scratch
display_update = true
version_number = "4.6"
update_name = "update_v$(version_number)"

function __init__()
if display_update
    # Get scratch space for this package
    versions_dir = @get_scratch!("versions")
    if !isfile(joinpath(versions_dir, update_name))
        printstyled(
            stdout,
            """
            \nUpdate message: Agents v$(version_number)
            Please see the changelog online. Some key features:

            * Agents.jl now uses [Graphs.jl](https://github.com/JuliaGraphs/Graphs.jl) instead of LightGraphs.jl
            """;
            color = :light_magenta,
        )
        touch(joinpath(versions_dir, update_name))
    end
end
end

end # module


text:
x = rand(100, 100)
y = rand(100, 100)

@inbounds begin
    for i = 1:100
        for j = 1:100
            x[i, j] *= y[i, j]
            y[i, j] += x[i, j]
        end
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# This file is part of Kpax3. License is MIT.

function computeclusteriseqprobs!(x::UInt8,
                                  b::Int,
                                  priorC::AminoAcidPriorCol,
                                  support::MCMCSupport)
  support.wi.z[1, b] = support.wi.w[1, b] +
                       logcondmarglik(x, support.ni[b], support.vi,
                                      priorC.A[1, b], priorC.B[1, b])
  support.wi.z[2, b] = support.wi.w[2, b] +
                       logcondmarglik(x, support.ni[b], support.vi,
                                      priorC.A[2, b], priorC.B[2, b])
  support.wi.z[3, b] = support.wi.w[3, b] +
                       logcondmarglik(x, support.ni[b], support.vi,
                                      priorC.A[3, b], priorC.B[3, b])
  support.wi.z[4, b] = support.wi.w[4, b] +
                       logcondmarglik(x, support.ni[b], support.vi,
                                      priorC.A[4, b], priorC.B[4, b])

  support.tmp[1] = support.wi.z[1, b] - support.wi.c[b]
  support.tmp[2] = support.wi.z[2, b] - support.wi.c[b]
  support.tmp[3] = support.wi.z[3, b] - support.wi.c[b]
  support.tmp[4] = support.wi.z[4, b] - support.wi.c[b]

  M = max(support.tmp[1], support.tmp[2], support.tmp[3], support.tmp[4])

  M + log(exp(support.tmp[1] - M) + exp(support.tmp[2] - M) +
          exp(support.tmp[3] - M) + exp(support.tmp[4] - M))
end

function computeclusterjseqprobs!(x::UInt8,
                                  b::Int,
                                  priorC::AminoAcidPriorCol,
                                  support::MCMCSupport)
  support.wj.z[1, b] = support.wj.w[1, b] +
                       logcondmarglik(x, support.nj[b], support.vj,
                                      priorC.A[1, b], priorC.B[1, b])
  support.wj.z[2, b] = support.wj.w[2, b] +
                       logcondmarglik(x, support.nj[b], support.vj,
                                      priorC.A[2, b], priorC.B[2, b])
  support.wj.z[3, b] = support.wj.w[3, b] +
                       logcondmarglik(x, support.nj[b], support.vj,
                                      priorC.A[3, b], priorC.B[3, b])
  support.wj.z[4, b] = support.wj.w[4, b] +
                       logcondmarglik(x, support.nj[b], support.vj,
                                      priorC.A[4, b], priorC.B[4, b])

  support.tmp[1] = support.wj.z[1, b] - support.wj.c[b]
  support.tmp[2] = support.wj.z[2, b] - support.wj.c[b]
  support.tmp[3] = support.wj.z[3, b] - support.wj.c[b]
  support.tmp[4] = support.wj.z[4, b] - support.wj.c[b]

  M = max(support.tmp[1], support.tmp[2], support.tmp[3], support.tmp[4])

  M + log(exp(support.tmp[1] - M) + exp(support.tmp[2] - M) +
          exp(support.tmp[3] - M) + exp(support.tmp[4] - M))
end

function updateclusteriweights!(x::UInt8,
                                b::Int,
                                support::MCMCSupport)
  support.ni[b] += float(x)

  support.wi.w[1, b] = support.wi.z[1, b]
  support.wi.w[2, b] = support.wi.z[2, b]
  support.wi.w[3, b] = support.wi.z[3, b]
  support.wi.w[4, b] = support.wi.z[4, b]

  M = max(support.wi.w[1, b], support.wi.w[2, b], support.wi.w[3, b],
          support.wi.w[4, b])

  support.wi.c[b] = M + log(exp(support.wi.w[1, b] - M) +
                            exp(support.wi.w[2, b] - M) +
                            exp(support.wi.w[3, b] - M) +
                            exp(support.wi.w[4, b] - M))

  nothing
end

function updateclusterjweights!(x::UInt8,
                                b::Int,
                                support::MCMCSupport)
  support.nj[b] += float(x)

  support.wj.w[1, b] = support.wj.z[1, b]
  support.wj.w[2, b] = support.wj.z[2, b]
  support.wj.w[3, b] = support.wj.z[3, b]
  support.wj.w[4, b] = support.wj.z[4, b]

  M = max(support.wj.w[1, b], support.wj.w[2, b], support.wj.w[3, b],
          support.wj.w[4, b])

  support.wj.c[b] = M + log(exp(support.wj.w[1, b] - M) +
                            exp(support.wj.w[2, b] - M) +
                            exp(support.wj.w[3, b] - M) +
                            exp(support.wj.w[4, b] - M))

  nothing
end

function updateclusteri!(u::Int,
                         data::Matrix{UInt8},
                         support::MCMCSupport)
  support.vi += 1
  support.ui[support.vi] = u

  for b in 1:support.m
    updateclusteriweights!(data[b, u], b, support)
  end

  nothing
end

function updateclusterj!(u::Int,
                         data::Matrix{UInt8},
                         support::MCMCSupport)
  support.vj += 1
  support.uj[support.vj] = u

  for b in 1:support.m
    updateclusterjweights!(data[b, u], b, support)
  end

  nothing
end

function updatelogmargliki!(priorC::AminoAcidPriorCol,
                            support::MCMCSupport)
  for b in 1:support.m
    support.lpi[1, b] = logmarglik(support.ni[b], support.vi, priorC.A[1, b],
                                   priorC.B[1, b])
    support.lpi[2, b] = logmarglik(support.ni[b], support.vi, priorC.A[2, b],
                                   priorC.B[2, b])
    support.lpi[3, b] = logmarglik(support.ni[b], support.vi, priorC.A[3, b],
                                   priorC.B[3, b])
    support.lpi[4, b] = logmarglik(support.ni[b], support.vi, priorC.A[4, b],
                                   priorC.B[4, b])
  end

  nothing
end

function updatelogmargliki!(ni::Vector{Float64},
                            vi::Int,
                            priorC::AminoAcidPriorCol,
                            support::MCMCSupport)
  for b in 1:support.m
    support.lpi[1, b] = logmarglik(ni[b], vi, priorC.A[1, b], priorC.B[1, b])
    support.lpi[2, b] = logmarglik(ni[b], vi, priorC.A[2, b], priorC.B[2, b])
    support.lpi[3, b] = logmarglik(ni[b], vi, priorC.A[3, b], priorC.B[3, b])
    support.lpi[4, b] = logmarglik(ni[b], vi, priorC.A[4, b], priorC.B[4, b])
  end

  nothing
end

function updatelogmarglikj!(priorC::AminoAcidPriorCol,
                            support::MCMCSupport)
  for b in 1:support.m
    support.lpj[1, b] = logmarglik(support.nj[b], support.vj, priorC.A[1, b],
                                   priorC.B[1, b])
    support.lpj[2, b] = logmarglik(support.nj[b], support.vj, priorC.A[2, b],
                                   priorC.B[2, b])
    support.lpj[3, b] = logmarglik(support.nj[b], support.vj, priorC.A[3, b],
                                   priorC.B[3, b])
    support.lpj[4, b] = logmarglik(support.nj[b], support.vj, priorC.A[4, b],
                                   priorC.B[4, b])
  end

  nothing
end

function updatelogmarglikj!(nj::Vector{Float64},
                            vj::Int,
                            priorC::AminoAcidPriorCol,
                            support::MCMCSupport)
  for b in 1:support.m
    support.lpj[1, b] = logmarglik(nj[b], vj, priorC.A[1, b], priorC.B[1, b])
    support.lpj[2, b] = logmarglik(nj[b], vj, priorC.A[2, b], priorC.B[2, b])
    support.lpj[3, b] = logmarglik(nj[b], vj, priorC.A[3, b], priorC.B[3, b])
    support.lpj[4, b] = logmarglik(nj[b], vj, priorC.A[4, b], priorC.B[4, b])
  end

  nothing
end

function initclusteriweights!(x::UInt8,
                              b::Int,
                              k::Int,
                              priorC::AminoAcidPriorCol,
                              support::MCMCSupport)
  support.ni[b] = float(x)

  support.wi.w[1, b] = priorC.logγ[1] +
                       logmarglik(x, 1, priorC.A[1, b], priorC.B[1, b])
  support.wi.w[2, b] = priorC.logγ[2] +
                       logmarglik(x, 1, priorC.A[2, b], priorC.B[2, b])
  support.wi.w[3, b] = priorC.logγ[3] + priorC.logω[k][1] +
                       logmarglik(x, 1, priorC.A[3, b], priorC.B[3, b])
  support.wi.w[4, b] = priorC.logγ[3] + priorC.logω[k][2] +
                       logmarglik(x, 1, priorC.A[4, b], priorC.B[4, b])

  M = max(support.wi.w[1, b], support.wi.w[2, b],
          support.wi.w[3, b], support.wi.w[4, b])

  support.wi.c[b] = M + log(exp(support.wi.w[1, b] - M) +
                            exp(support.wi.w[2, b] - M) +
                            exp(support.wi.w[3, b] - M) +
                            exp(support.wi.w[4, b] - M))

  nothing
end

function initclusterjweights!(x::UInt8,
                              b::Int,
                              k::Int,
                              priorC::AminoAcidPriorCol,
                              support::MCMCSupport)
  support.nj[b] = float(x)

  support.wj.w[1, b] = priorC.logγ[1] +
                       logmarglik(x, 1, priorC.A[1, b], priorC.B[1, b])
  support.wj.w[2, b] = priorC.logγ[2] +
                       logmarglik(x, 1, priorC.A[2, b], priorC.B[2, b])
  support.wj.w[3, b] = priorC.logγ[3] + priorC.logω[k][1] +
                       logmarglik(x, 1, priorC.A[3, b], priorC.B[3, b])
  support.wj.w[4, b] = priorC.logγ[3] + priorC.logω[k][2] +
                       logmarglik(x, 1, priorC.A[4, b], priorC.B[4, b])

  M = max(support.wj.w[1, b], support.wj.w[2, b],
          support.wj.w[3, b], support.wj.w[4, b])

  support.wj.c[b] = M + log(exp(support.wj.w[1, b] - M) +
                            exp(support.wj.w[2, b] - M) +
                            exp(support.wj.w[3, b] - M) +
                            exp(support.wj.w[4, b] - M))

  nothing
end


text:
module Abm

# TODO: make this cleaner, see https://github.com/sbromberger/LightGraphs.jl/blob/master/src/LightGraphs.jl
using Agents:
    @agent, ABM, AbstractAgent, DiscreteSpace, GridSpace, add_agent!,
    add_agent_pos!, allids, dummystep, genocide!, ids_in_position, kill_agent!,
    move_agent!, nagents, nearby_ids, nearby_positions, positions, random_position, run!

using InteractiveDynamics: abmexploration
using GLMakie
using Random: Xoshiro, rand

include("impl.jl")
include("parameters.jl")
include("step.jl")
include("utils.jl")

function initialize(; M=33, seed=125)
    space = GridSpace((M, M))
    properties = Parameters()
    rng = Xoshiro(seed)

    model = ABM(
        Organism, space;
        properties, rng
    )

    for n ∈ 1:properties.bacteria_count
        agent = Organism(n, (1, 1), rand(model.rng) < 0.5 ? :a : :b, Vector{Int}(), :nothing, :nothing, -1, :bacterium)
        add_bacterium_single!(agent, model)
    end
    for _ ∈ 1:properties.phages_count
        roll = rand(model.rng)
        if roll < 0.45
            kind = :temperate
        elseif roll < 0.85
            kind = :virulent
        else
            kind = :deficient
        end
        add_agent!(Organism, model, roll < 0.5 ? :a : :b, Vector{Int}(), kind, :free, 0, :phage)
    end

    return model
end

function plot(model::ABM)
    params = Dict(
        :a => 0.0:0.1:1.0,
        :b => 0.0:0.1:1.0,
        :m => 0.0:1.0:40.0,
        :α => 0.0:1.0:100.0,
        :κ => 0.01:0.01:0.5,
        :moi_proxy_radius => 1:3,
        :infection_distance => 1:3,
        :latent_period => 1:10,
        :burst_size => 3:8,
        :growth_rate => 0.2:0.1:1.0,
        :decay_factor => 0.01:0.01:0.2,
        :p_burst => 0.4:0.1:0.9
    )

    ac(a::A) where {A<:AbstractAgent} = a.type === :bacterium ? :cyan : :magenta
    as = 10
    am(a::A) where {A<:AbstractAgent} = a.type === :bacterium ? :circle : :diamond
    f(model::ABM{<:DiscreteSpace,A}) where {A<:AbstractAgent} = map(v -> length(v), model.space.s)

    heatarray = f
    heatkwargs = (colorrange=(0, 7), colormap=:thermal)
    plotkwargs = (;
        ac, as, am,
        scatterkwargs=(strokewidth=1.0,),
        heatarray, heatkwargs
    )

    mdata = [:bacteria_count, :phages_count]
    fig, p = abmexploration(model; (model_step!)=complex_step!, params, plotkwargs..., mdata, figure = (resolution = (600,600),))
    return fig
end

function run(model, n)
    run!(model, dummystep, complex_step!, n;
        mdata=[:bacteria_count, :phages_count])
end

end


text:
#=
Takes one step in the Markov chain according to the given state's transition
vector.
Inputs:
- transition_vector: An array of length p where the ith element denotes the
                     probability of transitioning into state i from the current
                     state (which is implicitly given).
Outputs:
- An integer in {1,2,..., p} denoting the index of the next state which occurs
with the probability specified in transition_vector.
=#
function one_step(transition_vector)
    cumulative_vector = cumsum(transition_vector)
    state_index = 1
    U = rand()
    while(U > cumulative_vector[state_index])
        state_index += 1
    end

    return state_index
end

#=
Runs a Markov chain a fixed number of times and returns the final state.
Inputs:
- states: An array of length p where each element represents a state of the
          chain. Assumes the ordering is such that the probability of
          transitioning from state i to state j is given by P[i,j].
          Examples:
          [0,1,2,3,4,5]
          ["red", "green", "blue"]
          [0,"blue",3,2]
- P: A (p x p) matrix denoting the transition matrix of the Markov chain.
     Assumes the rows of P sum to 1 and P[i,j] >= 0 for all i,j.
- S0: The initial state of the Markov chain. Assumed to be an element of the
      state variable.
- steps: A nonnegative integer specifying the number of steps for the Markov
         chain to go through.
Outputs:
- An element of the states variable denoting the state the chain is in after
  the specified number of steps.
=#
function run_fixed_simulation(states, P, S0, steps)

    # Define a dictionary mapping state to their index in P
    state_dict = Dict(states[i] => i for i = 1:length(states))

    # Define the current state of the chain
    Sn = S0

    # Run through the specified number of steps
    for i = 1:steps
        state_index = get(state_dict, Sn, -1)
        next_state_index = one_step(P[state_index,:])
        Sn = states[next_state_index]
    end

    return Sn
end

# State space (χ)
states = [1,2,3,4]

# Transition matrix for gambler's ruin
P = [.1 .1 .4 .4; .1 .2 .3 .4; .4 .3 .15 .15; .4 .4 .15 .05]

# Create histogram of S_100
num_simulations = 100000
state100 = [run_fixed_simulation(states, P, 4, 100) for i = 1:num_simulations]

using Plots
histogram(state100, legend = false, grid = false,
    bg = RGB(247/255, 236/255, 226/255))

cd("/Users/elvis/Documents/DSI/courses/content/stochastic-approximations/julia-code/")
savefig("../images/invariant_measure_hist_2.svg")

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
using LinearAlgebra: Hermitian, Symmetric, UpperTriangular, LowerTriangular
using LinearAlgebra: _chol!, Cholesky, checkpositivedefinite
import LinearAlgebra: cholesky!

function cholesky!(A::Union{Hermitian{NiceNumber,S}, Symmetric{NiceNumber,S}} where S, ::Val{false}=Val(false); check::Bool = true)
    C, info = _chol!(A.data, A.uplo == 'U' ? UpperTriangular : LowerTriangular)
    check && checkpositivedefinite(info)
    return Cholesky(C.data, A.uplo, info)
end


text:
mutable struct BarplotGraphics{R<:Number} <: GraphicsArea
    bars::Vector{R}
    color::ColorType
    char_width::Int
    max_freq::Number
    max_len::Int
    symb::String
    transform

    function BarplotGraphics(
            bars::AbstractVector{R},
            char_width::Int,
            color::UserColorType,
            symb::Union{Char,String},
            transform) where {R}
        length(symb) == 1 || throw(ArgumentError("The symbol to print has to be a single character, got: \"" * symb * "\""))
        char_width = max(char_width, 10)
        max_freq, i = findmax(transform.(bars))
        max_len = length(string(bars[i]))
        char_width = max(char_width, max_len + 7)
        new{R}(bars, crayon_256_color(color), char_width, max_freq, max_len, string(symb), transform)
    end
end

nrows(g::BarplotGraphics) = length(g.bars)
ncols(g::BarplotGraphics) = g.char_width

function BarplotGraphics(
        bars::AbstractVector{R},
        char_width::Int,
        transform = identity;
        color::UserColorType = :green,
        symb = "■") where {R <: Number}
    BarplotGraphics(bars, char_width, crayon_256_color(color), symb, transform)
end

function addrow!(g::BarplotGraphics{R}, bars::AbstractVector{R}) where {R <: Number}
    append!(g.bars, bars)
    g.max_freq, i = findmax(g.transform.(g.bars))
    g.max_len = length(string(g.bars[i]))
    g
end

function addrow!(g::BarplotGraphics{R}, bar::Number) where {R <: Number}
    push!(g.bars, R(bar))
    g.max_freq, i = findmax(g.transform.(g.bars))
    g.max_len = length(string(g.bars[i]))
    g
end

function printrow(io::IO, g::BarplotGraphics, row::Int)
    0 < row <= nrows(g) || throw(ArgumentError("Argument \"row\" out of bounds: $row"))
    bar = g.bars[row]
    max_freq = g.max_freq
    max_bar_width = max(g.char_width - 2 - g.max_len, 1)
    val = g.transform(bar)
    bar_len = max_freq > 0 ? round(Int, max(val, zero(val))/max_freq * max_bar_width, RoundNearestTiesUp) : 0
    bar_str = max_freq > 0 ? repeat(g.symb, bar_len) : ""
    bar_lbl = string(bar)
    print_color(g.color, io, bar_str)
    print_color(:normal, io, " ", bar_lbl)
    pan_len = max(max_bar_width + 1 + g.max_len - bar_len - length(bar_lbl), 0)
    pad = repeat(" ", round(Int, pan_len))
    print(io, pad)
end


text:
include("../../utils.jl")

function pollardfactor(n::T=big(1208925819691594988651321)) where T<:Integer
    for c in T(1):(n - 3)
        G, r, q = ones(T,3)
        y = 2
        m::T = 1900
        ys::T = 0
        x::T = 0
        while G == 1
            x = y
            for i in 1:r
                y = (y^2 + c) % n
            end
            k = T(0)
            G = T(1)
            while k < r && G == 1
                for i in 1:min(r - k, m)
                    ys = y
                    y = (y^2 + c) % n
                    q = (q * abs(x - y)) % n
                end
                G = gcd(q, n)
                k += m
            end
            r *= 2
        end
        G == n && (G = T(1))
        while G == 1
            ys = (ys^2 + c) % n
            G = gcd(abs(x - ys), n)
        end
        if G != n
            return G
        end
    end
end

@gctime pollardfactor()

max_stars_repo_name: vtjnash/GCBenchmarks

max_issues_repo_name: vtjnash/GCBenchmarks

max_forks_repo_name: vtjnash/GCBenchmarks


text:
#=
Jacob-Leal model
- Julia version: 1.7.1
- Authors: Ivan Jericevich, Patrick Chang, Tim Gebbie
=#

#----- Structures -----#
mutable struct JacobLealParameters <: Parameters
	# Fixed parameters
    Nᴸ::Int64 # Number of low-frequency agents
    Nᴴ::Int64 # Number of high-frequency agents
    θ::Int64 # Mean trading frequency for LF trader (minutes)
	θ⁺::Int64 # Upper-bound trading frequency for LF traders (minutes)
	θ⁻::Int64 # Lower-bound trading frequency for LF traders (minutes)
	σ_chartist::Float64 # Chartists' random normal order size standard deviation
	α_chartist::Float64 # Chartists' order size scaling factor (0 < α < 1)
	σ_fundamentalist::Float64 # Fundamentlists' random normal order size standard deviation
	α_fundamentalist::Float64 # Fundamentalists' order size scaling factor
	σʸ::Float64 # Fundamental value standard deviation
	σᶻ::Float64 # LF traders' limit-price tick standard deviation
	δ::Float64 # Price drift for LF traders (Constrained to be greater than 0)
	γᴸ::Int64 # LF traders' unexectuted order resting period until expiry
	γᴴ::Int64 # HF traders' unexectuted order resting period until expiry
	κ⁺::Float64 # Uniform distribution upper-bound for HF traders' order-price support
	κ⁻::Float64 # Uniform distribution lower-bound for HF traders' order-price support
	η⁺::Float64 # Uniform distribution upper-bound for HF traders' activation
	η⁻::Float64 # Uniform distribution lower-bound for HF traders' activation
	λ::Float64 # Market volumes weight in HF traders' order size distribution
	ζ::Float64 # Intensity of switching parameter for chartist and fundamentalist strategies (Constrained to be greater than 0)
	# Time dependent parameters
	ϵₜᶜ::Float64 # Chartist order size noise
	ϵₜᶠ::Float64 # Fundamentalist order size noise
	function JacobLealParameters(; Nᴸ::Int64 = 10000, Nᴴ::Int64 = 100, θ::Int64 = 20, θ⁺::Int64 = 40,
		θ⁻::Int64 = 10, σ_chartist = 0.05::Float64, α_chartist::Float64 = 0.04, σ_fundamentalist::Float64 = 0.01,
		α_fundamentalist::Float64 = 0.04, σʸ::Float64 = 0.01, σᶻ::Float64 = 0.01, δ::Float64 = 0.0001,
		γᴸ::Int64 = 20, γᴴ::Int64 = 1, κ⁺::Float64 = 0.01, κ⁻::Float64 = 0.0, η⁺::Float64 = 0.2, η⁻::Float64 = 0.0,
		λ = 0.625, ζ = 1) # Constructor with default parameter values from Jacob-Leal
		if θ⁻ > θ⁺ # Lower bound must be less than upper bound
			error("Incorrect parameter values. Ensure θ⁻ < θ⁺")
		end
		if κ⁻ > κ⁺ # Lower bound must be less than upper bound
			error("Incorrect parameter values. Ensure κ⁻ < κ⁺")
		end
		if η⁻ > η⁺ # Lower bound must be less than upper bound
			error("Incorrect parameter values. Ensure η⁻ < η⁺")
		end
		if α_chartist > 1 || α_fundamentalist > 1 # Constrained between 0 and 1
			error("Incorrect parameter values. Ensure 0 < α < 1")
		end
		if λ > 1 # Constrained between 0 and 1
			error("Incorrect parameter values. Ensure 0 < λ < 1")
		end
		if γᴸ < γᴴ # HF orders expire faster than LF orders
			error("Incorrect parameter values. Ensure γᴸ > γᴴ")
		end
		new(abs(Nᴸ), abs(Nᴴ), abs(θ), abs(θ⁺), abs(θ⁻), abs(σ_chartist), abs(α_chartist), abs(σ_fundamentalist), abs(α_fundamentalist),
		abs(σʸ), abs(σᶻ), abs(δ), abs(γᴸ), abs(γᴴ), abs(κ⁺), abs(κ⁻), abs(η⁺), abs(η⁻), abs(λ), abs(ζ))
	end
end
@enum Strategy begin
	chartist = 1
	fundamentalist = 2
end
mutable struct LowFrequency <: Agent
	type::Strategy # Chartist or Fundamentalist
	activated::Bool # Activated or not for next trading sesion
	θ::Int64 # Agent trading frequency randomly sampled from truncated exponential
	πₜᶜ::Float64 # Agent's chartist strategy profit at time t
	πₜᶠ::Float64 # Agent's fundamentalist strategy profit at time t
	ϕₜᶜ::Float64 # Probability of agent adopting chartist strategy at time t
	ϕₜᶠ::Float64 # Probability of agent adopting fundamentalist strategy at time t
	dₜᶜ::Float64 # Limit order size for chartist strategy placed at time t
	dₜᶠ::Float64 # Limit order size for fundamentalist strategy placed at time t
	pₜ::Float64 # Limit order price placed at time t
end
mutable struct HighFrequency <: Agent
	activated::Bool # Activated or not for next trading sesion
	πₜ::Float64 # Agent profit at time t (Not used ??CORRECT??)
	Δx::Float64 # Relative price difference activation threshold
	κ::Float64 # Order depth parameter
	dₜ::Float64 # Limit order size placed at time t
	pₜ::Float64 # Limit order price placed at time t
end
mutable struct LimitOrder
	periods_until_expiry::Int64
	size::Float64
	price::Float64
end
mutable struct LimitOrderBook
	t::Int64
	bₜ::Float64 # Best bid
    aₜ::Float64 # Best ask
    bids::Vector{LimitOrder} # Stores all the active bids
    asks::Vector{LimitOrder} # Stores all the active asks
	fₜ::Float64 # Fundamental value of stock at time t
	fₜ₋₁::Float64 # Fundamental value of stock at time t - 1
	p̄::Vector{Float64} # List of closing prices
end
#---------------------------------------------------------------------------------------------------

#----- Results and tracing -----#
mutable struct Options <: SimulationOptions
	time::Int64 # Simulation time (number of trading sessions in days)
	p₀::Float64 # Initial price
	f₀::Float64 # Initial value
	store_trace::Bool
	show_trace::Bool
	show_every::Int64
	seed::Int64
	function Options(; time::Int64 = 100, p₀::Float64 = 100.0, f₀::Float64 = 100.0, store_trace::Bool = false, show_trace::Bool = false,
		show_every::Int64 = 1, seed::Int64 = 1)
		new(time, p₀, f₀, store_trace, show_trace, show_every, seed)
	end
end
mutable struct State <: SimulationState
	t::Int64
	LOB_state::Dict{String, Float64}
	agents::Dict{String, Int64}
end
const Trace = Vector{State}
struct JacobLealResults <: SimulationResults
	p̄::Vector{Float64}
	trace::Trace
end
function Trace!(trace::Trace, options::Options, t::Int64, agents::Vector{Agent}, LOB::LimitOrderBook)
	ask_volume = isempty(LOB.asks) ? 0.0 : sum(ask.size for ask in LOB.asks)
	bid_volume = isempty(LOB.bids) ? 0.0 : sum(bid.size for bid in LOB.bids)
	fundamentalists, chartists, highfrequency = 0, 0, 0
	for agent in agents
		if agent.activated
			if isa(agent, LowFrequency)
				if agent.type == fundamentalist
					fundamentalists += 1
				else
					chartists += 1
				end
			else
				highfrequency += 1
			end
		end
	end
	agents = Dict{String, Int64}("Fundamentalists" => fundamentalists, "Chartsists" => chartists, "High-Frequency Traders" => highfrequency)
	LOB_state = Dict{String, Float64}("price" => LOB.p̄[end], "value" => LOB.fₜ, "ask volume" => ask_volume, "bid volume" => bid_volume)
	state = State(t, LOB_state, agents)
	if options.store_trace
        push!(trace, state)
    end
	if options.show_trace
        if state.t % options.show_every == 0
			LOB_imbalance = round((bid_volume - ask_volume) / (bid_volume + ask_volume), digits = 3)
			@printf("%-5s    %-6d:%6d                  %-5d   %-6.2f          %-6.2f  %-4.3f\n",
			t, fundamentalists, chartists, highfrequency, round(LOB.p̄[end], digits = 2), round(LOB.fₜ, digits = 2), LOB_imbalance)
        end
    end
end
#---------------------------------------------------------------------------------------------------

#----- Simulation initialisation -----#
function InitialiseSimulation(parameters::JacobLealParameters, options::Options)::Tuple{Vector{Agent}, LimitOrderBook}
	agents = Vector{Agent}() # Initialization
	for i in 1:parameters.Nᴸ # Iterate and populate low-frequency agents
		type = rand() < 0.5 ? chartist : fundamentalist # Choose strategy with 50/50 probability
		θ = convert(Int64, trunc(rand(Truncated(Exponential(parameters.θ), parameters.θ⁻, parameters.θ⁺)))) # Sample trading speed from truncated exponential
		push!(agents, LowFrequency(type, true, θ, 0, 0, 0.5, 0.5, 0, 0, 0)) # All agents activated at start
	end
	for j in 1:parameters.Nᴴ # Iterate and populate high-frequency agents
		Δx = rand(Uniform(parameters.η⁻, parameters.η⁺))
		κ = rand(Uniform(parameters.κ⁻, parameters.κ⁺))
		push!(agents, HighFrequency(true, 0, Δx, κ, 0, 0)) # All agents actiivated at start
	end
	p₋₁ = rand(Normal(options.p₀, parameters.σʸ))
	f₋₁ = rand(Normal(options.f₀, parameters.σʸ))
	LOB = LimitOrderBook(0, options.p₀, options.p₀, Vector{LimitOrder}(), Vector{LimitOrder}(), options.f₀, f₋₁, [p₋₁, options.p₀]) # Initial prices and fundamental values
	return agents, LOB
end
#---------------------------------------------------------------------------------------------------

#----- Order creation -----#
# Set order sizes and prices
# Post orders to LOB
function SubmitOrder!(agent::LowFrequency, LOB::LimitOrderBook, parameters::JacobLealParameters)
	if agent.activated
		zₜ = rand(Normal(0, parameters.σᶻ))
		agent.dₜᶜ = parameters.α_chartist * (LOB.p̄[end] - LOB.p̄[end - 1]) + parameters.ϵₜᶜ # Chartist strategy limit order size
		agent.dₜᶠ = parameters.α_fundamentalist * (LOB.fₜ - LOB.p̄[end]) + parameters.ϵₜᶠ # Fundamentalist strategy limit order size
		dₜ = agent.type == chartist ? agent.dₜᶜ : agent.dₜᶠ # Choose order size corresponding to agent strategy
		agent.pₜ = LOB.p̄[end] * (1 + parameters.δ) * (1 + zₜ) # Limit order price
		order = LimitOrder(parameters.γᴸ, abs(dₜ), agent.pₜ) # Create order
		# Push order to LOB
		if dₜ < 0 # Sell order
			push!(LOB.asks, order)
		elseif dₜ > 0 # Buy order
			push!(LOB.bids, order)
		end
	end
end
function SubmitOrder!(agent::HighFrequency, LOB::LimitOrderBook, parameters::JacobLealParameters)
	if agent.activated
		if rand() < 0.5 # Buy order
			if !isempty(LOB.asks) # Only submit order if there is liquidity on opposite side
				# Sample order size from Poisson with mean depending on mean volume on opposite side tuncated from above by the total volume on the opposite side
				total_ask_volume = sum(order.size for order in LOB.asks)
				λ = total_ask_volume / length(LOB.asks)
				agent.dₜ = abs(rand(truncated(Poisson(λ), upper = total_ask_volume)))
				agent.pₜ = LOB.bₜ *  (1 - agent.κ)
				if agent.dₜ > 0 # Only post if order size > 0
					push!(LOB.bids, LimitOrder(parameters.γᴴ, agent.dₜ, agent.pₜ))
				end
			end
		else # Sell order
			if !isempty(LOB.bids) # Only submit order if there is liquidity on opposite side
				# Sample order size from Poisson with mean depending on mean volume on opposite side tuncated from above by the total volume on the opposite side
				total_bid_volume = sum(order.size for order in LOB.bids)
				λ = total_bid_volume / length(LOB.bids)
				agent.dₜ = abs(rand(truncated(Poisson(λ), upper = total_bid_volume)))
				agent.pₜ = LOB.aₜ *  (1 + agent.κ)
				if agent.dₜ > 0  # Only post if order size > 0
					push!(LOB.asks, LimitOrder(parameters.γᴴ, agent.dₜ, agent.pₜ))
				end
			end
		end
	end
end
#---------------------------------------------------------------------------------------------------

#----- Start of session -----#
# Set fundamental market value
# Update time-dependent parameters
# Update LOB parameters
function PreOpen!(LOB::LimitOrderBook, parameters::JacobLealParameters)
	LOB.t += 1
	parameters.ϵₜᶜ = rand(Normal(0, parameters.σ_chartist))
	parameters.ϵₜᶠ = rand(Normal(0, parameters.σ_fundamentalist))
	yₜ = rand(Normal(0, parameters.σʸ))
	LOB.fₜ₋₁ = LOB.fₜ
	LOB.fₜ = LOB.fₜ₋₁ * (1 + parameters.δ) * (1 + yₜ)
end
#---------------------------------------------------------------------------------------------------

#----- Switch low-frequency trader strategy -----#
# Adapt trader strategies for session based on profit from previous session and strategy switching probability
function AdaptStrategy!(agent::LowFrequency, LOB::LimitOrderBook)
	agent.type = rand() < agent.ϕₜᶜ ? chartist : fundamentalist
	agent.activated = mod(LOB.t, agent.θ) == 0
end
function AdaptStrategy!(agent::HighFrequency, LOB::LimitOrderBook)
	agent.activated = abs((LOB.p̄[end] -  LOB.p̄[end - 1]) / LOB.p̄[end - 1]) > agent.Δx # Activate HF agent if relative price difference exceeds threshold
end
#---------------------------------------------------------------------------------------------------

#----- End of session -----#
# Calculate profits
# Expire LOB orders
function PostClose!(agents::Vector{Agent}, LOB::LimitOrderBook, parameters::JacobLealParameters)
	for agent in agents
		if isa(agent, LowFrequency) # Compute LF agent profitability
			# LF agents must compute profits for both strategies (whether adopted or not) ??CORRECT??
			agent.πₜᶠ = (LOB.p̄[end] - agent.pₜ) * abs(agent.dₜᶠ)
			agent.πₜᶜ = (LOB.p̄[end] - agent.pₜ) * abs(agent.dₜᶜ)
			agent.ϕₜᶜ = exp(agent.πₜᶜ / parameters.ζ) / (exp(agent.πₜᶜ / parameters.ζ) + exp(agent.πₜᶠ / parameters.ζ))
			agent.ϕₜᶠ = 1 - agent.ϕₜᶜ
		else # Compute HF agent profitability
			agent.πₜ = (LOB.p̄[end] - agent.pₜ) * agent.dₜ # (Not used ??CORRECT??)
		end
	end
	expired_orders = Vector{Int64}()
	for (index, order) in enumerate(LOB.asks)
		order.periods_until_expiry -= 1
		if order.periods_until_expiry == 0
			push!(expired_orders, index)
		end
	end
	deleteat!(LOB.asks, expired_orders)
	expired_orders = Vector{Int64}()
	for (index, order) in enumerate(LOB.bids)
		order.periods_until_expiry -= 1
		if order.periods_until_expiry == 0
			push!(expired_orders, index)
		end
	end
	deleteat!(LOB.bids, expired_orders)
end
#---------------------------------------------------------------------------------------------------

#----- Closing Auction -----#
function ClosingAuction!(LOB::LimitOrderBook)
	if !(isempty(LOB.bids) || isempty(LOB.asks))
		# Sort by price
		sort!(LOB.asks, by = x -> x.price)
		sort!(LOB.bids, by = x -> x.price, rev = true)
		LOB.aₜ, LOB.bₜ = LOB.asks[1].price, LOB.bids[1].price # Set best bid and best ask
		closing_price = 0
		while LOB.aₜ <= LOB.bₜ # Match order until none cross the spread
			ask_size, bid_size = LOB.asks[1].size, LOB.bids[1].size # Get order sizes on TOB
			if bid_size > ask_size # Bid is partially filled
				bid_size -= ask_size
				LOB.bids[1].size = bid_size
				popfirst!(LOB.asks) # Remove matched orders from TOB on ask side
			elseif bid_size < ask_size # Ask is partially filled
				ask_size -= bid_size
				LOB.asks[1].size = ask_size
				popfirst!(LOB.bids) # Remove matched orders from TOB on bid side
			else # Both sides are completely filled
				# Remove matched orders from TOB on both sides
				popfirst!(LOB.asks)
				popfirst!(LOB.bids)
			end
			closing_price = (LOB.aₜ + LOB.bₜ) / 2 # Compute mid-price or last traded price
			if !(isempty(LOB.asks) || isempty(LOB.bids))
				LOB.aₜ, LOB.bₜ = LOB.asks[1].price, LOB.bids[1].price # Update best ask and best bid
			else
				if !isempty(LOB.asks)
					LOB.aₜ = minimum(x -> x.price, LOB.asks)
				end
				if !isempty(LOB.bids)
					LOB.bₜ = maximum(x -> x.price, LOB.bids)
				end
				break # Break the while loop; no more crossing
			end
		end
		# Append new closing price
		if closing_price > 0
			push!(LOB.p̄, closing_price)
		else
			push!(LOB.p̄, LOB.p̄[end])
		end
	else # We still update new best ask and best bid because of new orders even though no crossing occurs
		if !isempty(LOB.asks)
			LOB.aₜ = minimum(x -> x.price, LOB.asks)
		end
		if !isempty(LOB.bids)
			LOB.bₜ = maximum(x -> x.price, LOB.bids)
		end
		push!(LOB.p̄, LOB.p̄[end])
	end
end
#---------------------------------------------------------------------------------------------------

#----- Model -----#
function Simulate(parameters::JacobLealParameters = JacobLealParameters(); options::Options = Options())::JacobLealResults
	Random.seed!(options.seed)
	tracing = options.store_trace || options.show_trace
	trace = Trace()
	agents, LOB = InitialiseSimulation(parameters, options) # Initialise agents and parameters
	if options.show_trace # Print header
        println("t        Fundamentalists:Chartists\tNᴴ\tpₜ\t\tfₜ\tLOB Imbalance")
		println("-------- ------------------------------ ------- --------------- ------- -------------")
    end
	Trace!(trace, options, 0, agents, LOB)
	for t in 1:options.time # Iterate through days/sessions
		PreOpen!(LOB, parameters) # Update stock fundamental value and parameters for session (required for first session)
		for agent in agents
			AdaptStrategy!(agent, LOB) # Switch strategy for current session based on previous session profits
			SubmitOrder!(agent, LOB, parameters) # Post order to LOB
		end
		ClosingAuction!(LOB) # Execute order matching for current trading session
		PostClose!(agents, LOB, parameters) # Calculate agent profit and probability of switching strategies
		if tracing
			Trace!(trace, options, t, agents, LOB)
		end
	end
	return JacobLealResults(LOB.p̄, trace)
end
#---------------------------------------------------------------------------------------------------

max_stars_repo_name: IvanJericevich/FinancialModels

max_issues_repo_name: IvanJericevich/FinancialModels

max_forks_repo_name: IvanJericevich/FinancialModels


text:
@testset "interval.jl" begin
    it = Interval(2.0, 5.0)
    @test minimum(it) == 2.0
    @test maximum(it) == 5.0

    @test in(3.0, it) == true
    @test in(10.0, it) == false
    @test in(5.0, it) == true

    @test isempty(it) == false
    @test isempty(Interval(2.0,1.0)) == true
    @test isempty(Interval(0.0, 0.0)) == true

    @test intersect(Interval(2.0, 5.0), Interval(3.0, 6.0)) == Interval(3.0, 5.0)
    @test intersect(Interval(2.0, 5.0), Interval(6.0, 9.0)) == Interval(0.0, 0.0)

    @test intersect(Interval(2.0, 5.0), Interval(0.0, 0.0)) == Interval(0.0, 0.0)
    @test intersect(Interval(5.0, 5.0), Interval(2.0, 5.0)) == Interval(0.0, 0.0)
end


text:
# Simplistic kernel calculation
function kernelize!(out::AbstractArray, x::AbstractArray, y::AbstractArray, kernel::Function;
		    center::Bool=false, symmetric::Bool=false)
	m = size(x, 2)
	n = size(y, 2)

	if !symmetric   # generic case, x and y are different
		@inbounds for j in 1:n
        		yj = view(y,:,j)
        		@simd for i in 1:m
				xi = view(x,:,i)
				out[i,j] = kernel(xi, yj)
        		end
		end
	else
    		@inbounds for j = 1:n
        		yj = view(y,:,j)
        		@simd for i in 1:j
				xi = view(x,:,i)
            			out[i,j] = kernel(xi, yj)
				out[j,i] = out[i,j]
        		end
    		end
	end

	if center
		μi = mean(out, 2)	# kernel mean of individual rows
		μj = mean(out, 1)	# kernel mean of individual columns
		μ = mean(out)		# overall mean
		@inbounds @simd for i in 1:m
			for j in 1:n
				out[i,j] = out[i,j] - μi[i] - μj[j] + μ
			end
		end
	end
	return out
end

kernelize!(out::AbstractArray, x::AbstractArray, kernel::Function; center::Bool=false, symmetric::Bool=true) =
	kernelize!(out, x, x, kernel; center=center, symmetric=symmetric)

function kernelize(x::AbstractArray, y::AbstractArray, kernel::Function; center::Bool=false, symmetric::Bool=false)
	m = size(x, 2)
	n = size(y, 2)
	out = similar(x, m, n)
	kernelize!(out, x, y, kernel; center=center, symmetric=symmetric)
end

kernelize(x::AbstractArray, kernel::Function; center::Bool=false, symmetric::Bool=true) =
	kernelize(x, x, kernel; center=center, symmetric=symmetric)

##########################
# FunctionCell Interface #
##########################
"""
	kernel([f];[center=false, symmetric=false])

Constructs an untrained function cell using the kernel function `f` that accepts two vector arguments `x` and `y`,
and returns a scalar value (default `(x,y)->x'y`). When piped data into, the untrained function cell returns a
trained function cell, retaining the training data.

# Keyword arguments
  * `symmetric::Bool` specifies whether the kernel should be symmetric i.e. calculated using the same observations (default `false`)
  * `center::Bool` specifies whether the center the kernel (default `false`)

A list of common used kernels can be found below (definitions from the documentation of to MLKernels.jl):
# Kernels
  * `(x,y)-> a⋅xᵀy + c		Linear, a>0, c≥0`
  * `(x,y)-> (a⋅xᵀy+c)ᵈ		Polynomial, a>0, c≥0, d∈Ζ₊`
  * `(x,y)-> tanh(a⋅xᵀy+c)		Sigmoid, a>0, c≥0`
  * `(x,y)-> exp(a⋅xᵀy)		Exponentiated, a>0`
  * `(x,y)-> exp(-α*‖x-y‖)		Exponential, α>0`
  * `(x,y)-> exp(-α⋅‖x-y‖²)		Squared exponential, α>0`
  * `(x,y)-> exp(-α⋅‖x-y‖ᵞ)		Gamma exponential, α>0, 0<γ<1`
  * `(x,y)-> (1 + α⋅‖x-y‖²)⁻ᵝ	Rational quadratic, α>0, β>0`
  * `(x,y)-> (1 + α⋅‖x-y‖ᵞ)⁻ᵝ	Gamma-rational, α>0, β>0, 0<γ,1`
  * `(x,y)-> exp{-α⋅Σᵢsin²[p⋅(xᵢ-yᵢ)]}	Periodic kernel, p>0, α>0`
  * `(x,y)-> ‖x-y‖²ᵞ		Power, 0<γ≤1 (negative definite)`
  * `(x,y)-> log(1+α⋅‖x-y‖²ᵞ)	Log, α>0, 0<γ≤1 (negative definite)`

# Examples
```
julia> k(a,c) = (x,y)->a*x'y+c # linear kernel
k (generic function with 1 method)

julia> K = kernel(k(1.0,2.0))
Kernel: symmetric=false, varying I/O dimensions, untrained

julia> x=[1,2,3]; y=[4,5,6,7];

julia> y |> K(x)
3×4 Array{Int64,2}:
  6   7   8   9
 10  12  14  16
 14  17  20  23
```
"""
kernel(f::Function=(x,y)->x'y; center::Bool=false, symmetric::Bool=false) =
	FunctionCell(kernel, (f,), ModelProperties(), "Kernel: center=$center, symmetric=$symmetric";
	      center=center, symmetric=symmetric)

############################
# DataCell/Array Interface #
############################
"""
	kernel(x [,f] [;symetric=false])

Trains the function cell using the kernel function `f` and data `x` which will be used
for future kernel calculations.
"""
# Training
kernel(x::T where T<:CellData, f::Function=(x,y)->x'y; center::Bool=false, symmetric::Bool=false) =
	kernel(getx!(x), f; center=center, symmetric=symmetric)
kernel(x::Tuple{T,S} where T<:AbstractArray where S<:AbstractArray, f::Function=(x,y)->x'y;
       center::Bool=false, symmetric::Bool=false) = kernel(x[1], f; center=center, symmetric=symmetric)
kernel(x::T where T<:AbstractVector, f::Function=(x,y)->x'y; center::Bool=false, symmetric::Bool=false) =
	kernel(mat(x, LearnBase.ObsDim.Constant{2}()), f; center=center, symmetric=symmetric)
kernel(x::T where T<:AbstractMatrix, f::Function=(x,y)->x'y; center::Bool=false, symmetric::Bool=false) = begin

	# Build model properties
	modelprops = ModelProperties(nvars(x),nobs(x))

	# Returned trained cell
	FunctionCell(kernel, Model((f, getobs(x), center, symmetric), modelprops),
	      "Kernel: center=$center, symmetric=$symmetric")
end

# Execution
kernel(x::T where T<:CellData, model::Model{<:Tuple{<:Function,Matrix,Bool,Bool}}) =
	datacell(kernel(getx!(x), model), gety(x))
kernel(x::T where T<:AbstractVector, model::Model{<:Tuple{<:Function,Matrix,Bool,Bool}}) =
	kernel(mat(x, LearnBase.ObsDim.Constant{2}()), model)
kernel(x::T where T<:AbstractMatrix, model::Model{<:Tuple{<:Function,Matrix,Bool,Bool}}) =
	kernelize(model.data[2], getobs(x), model.data[1] ;center=model.data[3], symmetric=model.data[4])


text:
## Exercise 4-1

## 🐢 = Turtle()
## @svg begin
##     forward(🐢, 100)
##     turn(🐢, -90)
##     forward(🐢, 100)
## end

## Now modify the macro to draw a square. Don’t go on until you’ve got it working!
println("Ans: ")

using ThinkJulia

turtle = Turtle()

@svg begin
    forward(turtle, 100)
    turn(turtle, -90)
    forward(turtle, 100)
    turn(turtle, -90)
    forward(turtle, 100)
    turn(turtle, -90)
    forward(turtle, 100)
end

println("End.")


text:
function trim(polymer)
    chars = Array{Char, 1}()
    for char in polymer
        if char == '\n'
            continue
        elseif !(isempty(chars)) && lowercase(char) == lowercase(chars[end]) && chars[end] != char
            pop!(chars)
        else
            push!(chars, char)
        end
    end
    return join(chars)
end

function purge(polymer, char)
    return join(x for x in polymer if lowercase(x) != char)
end

function foo()
    polymer = trim(convert.(Char, read(open("data/5.txt", "r"))))
    min_length = 999999
    for char in 'a':'z'
        len = length(trim(purge(polymer, char)))
        if min_length > len
            min_length = len
        end
    end
    return min_length
end

println(foo())


text:
function print_header(params::Params, lpip_pb::LPIPLinearProblem{T}) where T
    if params.verbose
        println("                                                   ")
        println(" LPIP.jl v0.1.0  (c) Guilherme Bodin, Iago Chávarry 2019        ")
        println("                                                   ")
        println("Problem with: ", lpip_pb.n - lpip_pb.l, " variables,")
        println("              ", lpip_pb.f, " equality constraints,")
        println("              ", lpip_pb.l, " nonnegative constraints.")
        println("                                                   ")
        println("Parameters:  rho = ", params.rho, " alpha = ", params.alpha, " tol = ", params.tol)
        println("             max_iter = ", params.max_iter, " time limit (s) = ", params.time_limit)
        println("                                                   ")
        println("||     obj    |    gap    | iter |   time (s) ||")
    end
end

function print_evolution(params::Params, lpip_pb::LPIPLinearProblem{T}, iter::Int, t0::Float64) where T
    if params.verbose
        t1 = time() - t0
        obj = dot(lpip_pb.c, lpip_pb.variables.x) + lpip_pb.d
        gap = abs(dot(lpip_pb.variables.x, lpip_pb.variables.s))

        str_obj = " " * @sprintf("%.6f", obj) * "  |"
        str_gap   = " " * @sprintf("%.6f", gap) * "  |"
        str_iter = @sprintf("%d", iter) * " |"
        str_time   = "   " * @sprintf("%.4f", t1) * "   ||"

        str = "||"
        str *= " "^max(0, 11 - length(str_obj)) * str_obj
        str *= " "^max(0, 11 - length(str_gap)) * str_gap
        str *= " "^max(0, 7 - length(str_iter)) * str_iter
        str *= " "^max(0, 12 - length(str_time)) * str_time

        println(str)
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 76b8f2b0-2719-11eb-37f6-01188db45771
using TikzNeuralNetworks

# ╔═╡ 36fce9f0-27c4-11eb-2ee9-bb76b493374e
using PyPlot

# ╔═╡ 643b6720-27c4-11eb-1c36-6f9262dc4995
using Seaborn

# ╔═╡ b5d7170e-2719-11eb-02ca-172fa61417d1
md"""
## Autoencoder
"""

# ╔═╡ 8263b410-2719-11eb-03ad-4fe38be29532
begin
	tikznn = TikzNeuralNetwork(
		input_size=5,
		input_arrows=false,
		input_label=i->i == 3 ? "\$\\mathbf{x}\$" : "",
		# activation_functions=["g","g","g"],0
		hidden_layer_sizes=[4,2,4],
		hidden_layer_labels=["", "\$\\tilde{\\mathbf{x}}\$", ""],
		output_size=5,
		output_arrows=false,
		output_label=i->i == 3 ?  "\$\\mathbf{x}'\$" : "")
	tikznn.tikz.width="12cm"
	tikznn
end

# ╔═╡ 45ff5610-27a9-11eb-1592-69eb568a2e55
save(TEX("../tex/autoencoder-nn.tex"), tikznn)

# ╔═╡ be9ee750-271a-11eb-08aa-9377c292341d
md"""
## Adversary

$$\mathcal{J}(\mathbf{\hat{y}}, \mathbf{y}) = \frac{1}{m}\sum_{i=1}^m \mathcal{L}(\hat{y}_i, y_i)$$

$$\mathcal{L}(\hat{y}, y) = - y \log(\hat y) - (1 - y)\log(1-\hat y)$$

$$\mathcal{L}(\hat{y}, y) = - y \log(\hat y) \Omega - (1 - y)\log(1-\hat y)$$

"""

# ╔═╡ c15bbd60-271a-11eb-08f7-2349fae8fd8f
begin
	tikzadv = TikzNeuralNetwork(
		input_size=2,
		input_arrows=false,
		input_label=i->"\$x_{$i}\$",
		hidden_layer_sizes=[4,2,1],
		hidden_layer_labels=["ReLU", "ReLU", L"\sigma"],
		output_size=1,
		output_arrows=false,
		output_label=i->"\$\\hat{y}\$")
	tikzadv.tikz.width="12cm"
	tikzadv
end

# ╔═╡ 45bb62f0-271f-11eb-2462-9bbd423021d4
md"""
## System Under Test
"""

# ╔═╡ 49b91d70-271f-11eb-3c2e-81f1d12885c5
begin
	tikzsut = TikzNeuralNetwork(
		input_size=4,
		input_arrows=false,
		input_label=i->"\$x_{$i}\$",
		hidden_layer_sizes=[2],
		hidden_layer_labels=["ReLU"],
		output_size=1,
		output_arrows=false,
		output_label=i->"\$\\mathbf{\\hat{y}}\$")
	tikzsut.tikz.width="12cm"
	tikzsut
end

# ╔═╡ 3893b690-27c4-11eb-084f-1fe93bb22697
M = [10 0 0 0 0;
     0 10 0 0 0;
     0 0 10 0 0;
     0 0 0 10 0;
     0 0 0 0 10]

# ╔═╡ 99c0e13e-27c4-11eb-0ce5-3521508cfbe4
close("all"); heatmap(M, annot=true, cmap="viridis", xticklabels=1:5, yticklabels=1:5); xlabel(L"x"); ylabel(L"y"); gcf()

# ╔═╡ fbf23580-27c4-11eb-2f92-0bd160f6cdc4
savefig("confusion.pgf")

# ╔═╡ Cell order:
# ╠═76b8f2b0-2719-11eb-37f6-01188db45771
# ╟─b5d7170e-2719-11eb-02ca-172fa61417d1
# ╠═8263b410-2719-11eb-03ad-4fe38be29532
# ╠═45ff5610-27a9-11eb-1592-69eb568a2e55
# ╠═be9ee750-271a-11eb-08aa-9377c292341d
# ╠═c15bbd60-271a-11eb-08f7-2349fae8fd8f
# ╟─45bb62f0-271f-11eb-2462-9bbd423021d4
# ╠═49b91d70-271f-11eb-3c2e-81f1d12885c5
# ╠═36fce9f0-27c4-11eb-2ee9-bb76b493374e
# ╠═643b6720-27c4-11eb-1c36-6f9262dc4995
# ╠═3893b690-27c4-11eb-084f-1fe93bb22697
# ╠═99c0e13e-27c4-11eb-0ce5-3521508cfbe4
# ╠═fbf23580-27c4-11eb-2f92-0bd160f6cdc4


text:
function pcontrib(p::Int64, a::Int64)
    n = one(p)
    pcon = one(p)
    for i in 1:a
        n *= p
        pcon += n
    end
    return pcon
end

function divisorsum(n::Int64)
    dsum = one(n)
    for (p, a) in factor(n)
        dsum *= pcontrib(p, a)
    end
    dsum -= n
end


text:
export TreeNode;
export EMPTY_LABEL, EMPTY_NODE, ROOT_TOKEN, NO_HEAD

const ROOT_TOKEN = "-ROOT-"
const NO_HEAD = "-NO_HEAD-"
const EMPTY_LABEL = "EMPTY_LABEL"
const EMPTY_NODE = -1

mutable struct TreeNode
  word_id::Integer
  label::String
  head_id::Integer
  token::String

  TreeNode(word_id::Integer) = new(word_id, EMPTY_LABEL, EMPTY_NODE, ROOT_TOKEN)
  TreeNode(word_id::Integer, token::String) = new(word_id, EMPTY_LABEL, EMPTY_NODE, token)
  TreeNode(word_id::Integer, label::String, head_id::Integer) = new(word_id, label, head_id)
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
include("abstract_learner.jl")
include("temporal_difference_learner.jl")
include("basic_dqn_learner.jl")
include("dqn_learner.jl")
include("prioritized_dqn_learner.jl")
include("monte_carlo_learner.jl")
include("gradient_bandit_learner.jl")
include("rainbow.jl")


text:
# surely, there is a more straightforward way
# These are meant to collect and convert arguments for use with PermutationsA

collect1(a::Vector{<:String}) =  map(BigInt,a)
collect1(a::Vector{<:Real}) = a

function collect1(::Type{T}, a::Vector{V}) where {T, V}
    return T == V ? a : collect(T, a)  # note that it may return a copy!
end

function tupcollect(::Type{T}, a::Tuple) where T
    aout = Array{Array{T,1}}(undef, 0)
    for x in a
        a1 = Array{T}(undef, 0)
        for x1 in x
            push!(a1,convert(T,x1))
        end
        push!(aout,a1)
    end
    return aout
end

for (f,t) in ((:int32, Int32), (:int64, Int64), (:BigInt, BigInt))
    @eval begin
        function tupcollect(::Type{$t}, a::Tuple)
            aout = Array{Array{$t,1}}(undef, 0)
            for x in a
                a1 = Array{$t}(undef, 0)
                for x1 in x
                    push!(a1, ($f)(x1))
                end
                push!(aout, a1)
            end
            return aout
        end
    end
end

function tupcollect(a::Tuple)
    tupcollect(Int,a)
end

function collect2(a::AbstractVector...)
    collect(map(collect1,a))
end

function collect2(::Type{T}, a::AbstractVector...) where T
    collect(map((x)->collect1(T,x) ,a))
end


text:
include("setup_ui_test_module.jl")

clear_breakpoints!(); clear_uninstrumenteds!()
#@testset "No breakpoints" begin
    make_readline_patch([])
    record = make_recording_breakpoint_hit_patch()

    @run eg1()
    @test record == []
#end

#@testset "1 breakpoint" begin
    make_readline_patch(["CC"])
    record = make_recording_breakpoint_hit_patch()
    set_breakpoint!(eg2)

    @show list_breakpoints()
    @run eg1()
    @test first(record).f == eg2
#end

clear_breakpoints!(); clear_uninstrumenteds!()
#@testset "Two breakpoints" begin
    make_readline_patch(["CC", "CC"])
    record = make_recording_breakpoint_hit_patch()
    set_breakpoint!(eg2)
    set_breakpoint!(eg3)
    @run eg1()
    @test first.(record) == [eg2, eg3]
#end

###################################################
clear_breakpoints!(); clear_uninstrumenteds!()
#@testset "No breakpoints, With no instrumenting of Base" begin
    make_readline_patch([])
    record = make_recording_breakpoint_hit_patch()
    set_uninstrumented!(Base)
    @run eg1()
    @test record == []
#end
###########################

clear_breakpoints!(); clear_uninstrumenteds!()
#@testset "breakpoint by file and linenum" begin
    make_readline_patch(["CC"])
    record = make_recording_breakpoint_hit_patch()

    set_breakpoint!("demo.jl", 9)
    @test !isempty(list_breakpoints())
    @run eg1()
    @test first.(record) == [eg2]
#end

#########################################################
 # Stepping Mode

clear_breakpoints!(); clear_uninstrumenteds!()
#@testset "step in" begin
    make_readline_patch(["SI", "CC"])
    record = make_recording_breakpoint_hit_patch()

    set_breakpoint!(eg2)
    @run eg1()
    @test first.(record) == [eg2, eg21]
#end

clear_breakpoints!(); clear_uninstrumenteds!()
#@testset "step out" begin
    set_breakpoint!(eg2)
    make_readline_patch(["SO", "CC"])
    record = make_recording_breakpoint_hit_patch()

    @run eg1()
    @test first.(record) == [eg2, eg1]
#end

clear_breakpoints!(); clear_uninstrumenteds!()
#@testset "step next" begin
    set_breakpoint!(eg2)
    make_readline_patch(["SN", "SN", "CC"])
    record = make_recording_breakpoint_hit_patch()

    @run eg1()
    @test first.(record) == [eg2, eg2, eg2]
#end

###############################################

clear_breakpoints!(); clear_uninstrumenteds!()
#@testset "@enter" begin
    make_readline_patch(["CC"])
    record = make_recording_breakpoint_hit_patch()

    @enter eg1()
    @test first.(record) == [eg1]
    @test isempty(list_breakpoints())
#end

clear_breakpoints!(); clear_uninstrumenteds!()
#@testset "@enter, too complex" begin
    make_readline_patch(["CC"])
    record = make_recording_breakpoint_hit_patch()

    @test_throws ErrorException (@enter (()->eg1())())
    @test isempty(list_breakpoints())
#end

###############################################
clear_breakpoints!(); clear_uninstrumenteds!()
#@testset "Influence calling enviroment" begin
    global zzz = 10
    make_readline_patch(["zzz = 20", "CC"])

    set_breakpoint!(eg2)
    @run eg1()
    @test zzz==20
#end

clear_breakpoints!(); clear_uninstrumenteds!()
#@testset "Abort" begin
    make_readline_patch(["XX"])

    set_breakpoint!(eg2)
    @test nothing==@run eg1()
#end

########################################
# Variables
function var_demo1(x)
    local y
    z=2x
    y=z
    z=1
    return 0
end

clear_breakpoints!(); clear_uninstrumenteds!()
#@testset "Variables stepping" begin
    set_breakpoint!(var_demo1)
    make_readline_patch(["SN", "SN", "SN", "CC"])
    record = make_recording_breakpoint_hit_patch()
    @run var_demo1(5)

    @test record[1].variables == LittleDict(:x=>5)
    @test record[2].variables == LittleDict(:x=>5, :z=>10)
    @test record[3].variables == LittleDict(:x=>5, :z=>10, :y=>10)
    @test record[4].variables == LittleDict(:x=>5, :y=>10, :z=>1)
#end

#######
# Done
reset_patched_functions!()
println("\tUser interaction tests complete")


text:
cd(joinpath(ENV["HOME"],"CODE/DCAlign"))
using Pkg,DelimitedFiles,BenchmarkTools, Profile
Pkg.activate(".")
using Revise,DCAlign
m = DCAlign

q=21; N=87; gap=0;
ctype=Symbol("amino")
typel=Symbol("bm")

J,H =  m.read_parameters("../AlignPotts/test/PF00677/Parameters_bm_PF00677potts.dat",q,N,gap=0)
Lambda_all = readdlm("test/PF00677/Lambda_PF00677.dat")
lambda_o = Lambda_all[:,1]
lambda_e = Lambda_all[:,2]

delta = 50;
al = m.enveloptoalign( "test/PF00677/PF00677.full", "test/PF00677/PF00677.align", "test/PF00677/PF00677.ins", delta = delta, ctype = ctype);
M = length(al)

idx0 = 1
(aux,garb) = split(al[idx0][1], "/")
seqhmm = al[idx0][3];
seqins = al[idx0][4];
seq = m.Seq(seqhmm,al[idx0][2],ctype)
N = length(al[idx0][2])
display(al[idx0][1])
display(al[idx0][2])
display(seqhmm)
display(seqins)

μext = 0.00;
μint = 0.00;
β=1; T0=false; seed=33; damp=0.0;
@benchmark resalign=m.palign($seq, $J, $H, $lambda_o, $lambda_e,
    ctype = ctype ,mindec = 50, nprint=1, maxiter=2, T0 = T0, μext=μext, μint =  μint,seed=seed,
    damp=damp, epsconv=1e-10 ,verbose=false)

@benchmark resalign=m.palign($seq, $J, $H, $lambda_o, $lambda_e,
        ctype = ctype ,mindec = 50, nprint=1, maxiter=2, T0 = true, μext=μext, μint =  μint,seed=seed,
        damp=damp, epsconv=1e-10 ,verbose=false)

T0=false; seed=33; damp=0.0;
Profile.clear()
@profiler resalign=m.palign(seq, J, H, lambda_o, lambda_e,
        ctype = ctype ,mindec = 50, nprint=1, maxiter=20,
        T0 = T0, μext=μext, μint =  μint,seed=seed,
        damp=damp, epsconv=1e-10 ,verbose=false) combine=true

resalign=m.palign(seq, J, H, lambda_o, lambda_e,
                ctype = ctype ,mindec = 50, nprint=1, maxiter=20000,
                T0 = false, μext=μext, μint =  μint,seed=seed,
                damp=damp, epsconv=1e-10 ,verbose=true)


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

portfolio rebalancing problem
=#

struct PortfolioJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    num_stocks::Int
    epinormeucl_constr::Bool # add an L2 ball constraint
    epinorminf_constrs::Bool # add L1 and Linfty ball constraints
end

function build(inst::PortfolioJuMP{T}) where {T <: Float64}
    @assert xor(inst.epinormeucl_constr, inst.epinorminf_constrs)
    num_stocks = inst.num_stocks

    returns = rand(num_stocks)
    sigma_half = randn(num_stocks, num_stocks)
    x = randn(num_stocks)
    gamma = sum(abs, sigma_half * x) / norm(x)
    A = randn(div(num_stocks, 2), num_stocks)

    model = JuMP.Model()
    JuMP.@variable(model, invest[1:num_stocks])
    JuMP.@objective(model, Max, dot(returns, invest))
    JuMP.@constraint(model, sum(invest) == 0)
    JuMP.@constraint(model, A * invest .== 0)

    if inst.epinormeucl_constr
        JuMP.@constraint(model, vcat(gamma / sqrt(num_stocks),
            sigma_half * invest) in JuMP.SecondOrderCone())
    end
    if inst.epinorminf_constrs
        JuMP.@constraint(model, vcat(1, invest) in
            MOI.NormInfinityCone(1 + num_stocks))
        JuMP.@constraint(model, vcat(gamma, sigma_half * invest) in
            MOI.NormOneCone(1 + num_stocks))
    end

    return model
end


text:
function get_absolute_accuracy(obj::BaseUnivariateSolver)
    return jcall(obj, "getAbsoluteAccuracy", jdouble, ())
end

function get_evaluations(obj::BaseUnivariateSolver)
    return jcall(obj, "getEvaluations", jint, ())
end

function get_function_value_accuracy(obj::BaseUnivariateSolver)
    return jcall(obj, "getFunctionValueAccuracy", jdouble, ())
end

function get_max_evaluations(obj::BaseUnivariateSolver)
    return jcall(obj, "getMaxEvaluations", jint, ())
end

function get_relative_accuracy(obj::BaseUnivariateSolver)
    return jcall(obj, "getRelativeAccuracy", jdouble, ())
end

function solve(obj::BaseUnivariateSolver, arg0::jint, arg1::UnivariateFunction, arg2::jdouble)
    return jcall(obj, "solve", jdouble, (jint, UnivariateFunction, jdouble), arg0, arg1, arg2)
end

function solve(obj::BaseUnivariateSolver, arg0::jint, arg1::UnivariateFunction, arg2::jdouble, arg3::jdouble)
    return jcall(obj, "solve", jdouble, (jint, UnivariateFunction, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function solve(obj::BaseUnivariateSolver, arg0::jint, arg1::UnivariateFunction, arg2::jdouble, arg3::jdouble, arg4::jdouble)
    return jcall(obj, "solve", jdouble, (jint, UnivariateFunction, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4)
end


text:
using Test
using VariantsNMF

printstyled("Julia version: $VERSION\n", color=:blue)
printstyled("Running tests: \n", color=:blue)

const tests = [
    "convergence",
    "clustering",
    "divergences",
    "nmf_FI",
    "rank_by_stability",
    "rank_by_consensus",
]

# test utilities to avoid redundancies
include("testutils.jl")

for t in tests
    @testset "Test $t" begin
        include("$t.jl")
    end
end


text:
"""
    latin_hypercube_sampling(mins, maxs, n)

Randomly sample `n` vectors from the parallelogram defined
by `mins` and `maxs` using the Latin hypercube algorithm.

Returns `dims`×`n` matrix.
"""
function latin_hypercube_sampling(mins::AbstractVector{T},
                                  maxs::AbstractVector{T},
                                  n::Integer) where T<:Real
    length(mins) == length(maxs) ||
        throw(DimensionMismatch("mins and maxs should have the same length"))
    all(xy -> xy[1] <= xy[2], zip(mins, maxs)) ||
        throw(ArgumentError("mins[i] should not exceed maxs[i]"))
    dims = length(mins)
    result = Matrix{Float64}(undef, dims, n)
    cubedim = Vector{Float64}(undef, n)
    @inbounds for i in 1:dims
        imin = mins[i]
        dimstep = (maxs[i] - imin) / n
        for j in 1:n
            cubedim[j] = imin + dimstep * (j - 1 + rand())
        end
        result[i, :] .= shuffle!(cubedim)
    end
    return result
end


text:
#!/usr/bin/env julia

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "-M"
    help = "Print out the list of files this program will create."
    action = :store_true
    "config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using DiskTracer.model
using DiskTracer.constants

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]
pars = convert_dict(config["parameters"], config["model"])

mol = mols[species*transition]

using LaTeXStrings
using Plots
pyplot()

# Make a 2D plot of the source function and Upsilon
# density structure
function plot_S(pars::AbstractParameters)

    nr = 128
    nz = 128
    rs = linspace(0.1, 600 * AU, nr)
    zs = linspace(0.1, 600 * AU, nz)

    rr = rs./AU
    zz = zs./AU

    Ss = Array{Float64}(nz, nr)

    for i=1:nz
        for j=1:nr
            Ss[i,j] = DiskTracer.model.S_nu(rs[j], zs[i], mol.nu_0, pars)
        end
    end

    p = heatmap(log10.(Ss), xlabel=L"$r$ [AU]", ylabel=L"$z$ [AU]") #, size=(800,800))
    savefig(p, "source_function.png")
end

function plot_Upsilon(pars::AbstractParameters)

    nr = 128
    nz = 128
    rs = linspace(0.1, 200 * AU, nr)
    zs = linspace(0.1, 200 * AU, nz)

    rr = rs./AU
    zz = zs./AU

    Us = Array{Float64}(nz, nr)

    for i=1:nz
        for j=1:nr
            Us[i,j] = DiskTracer.model.Upsilon_nu(rs[j], zs[i], mol.nu_0, pars, mol)
        end
    end

    # show 5 decades of scale
    p = heatmap(log10.(Us+ 1e-5 * maximum(Us)), xlabel=L"$r$ [AU]", ylabel=L"$z$ [AU]")
    savefig(p, "Upsilon.png")

end

plot_S(pars)
plot_Upsilon(pars)


text:
@testset "119.pascals-triangle-ii.jl" begin
    @test get_pascal_row(3) == [1, 3, 3, 1]
    @test get_pascal_row(5) == [1, 5, 10, 10, 5, 1]
    @test get_pascal_row(6) == [1, 6, 15, 20, 15, 6, 1]
    @test get_pascal_row(23) ==
          [1, 23, 253, 1771, 8855, 33649, 100947, 245157, 490314, 817190, 1144066, 1352078,
           1352078, 1144066, 817190, 490314, 245157, 100947, 33649, 8855, 1771, 253, 23, 1]
end


text:
include("solver_preamble.jl")
using MosekTools
# It cannot be used in Direct mode yet as it gets:
# MathOptInterface.DeleteNotAllowed{MathOptInterface.ConstraintIndex{MathOptInterface.VectorOfVariables,MathOptInterface.RotatedSecondOrderCone}}
# for `sdsos_bivariate_quadratic`
factory = with_optimizer(Mosek.Optimizer, QUIET=true)
config = MOI.Test.TestConfig(atol=1e-5, rtol=1e-5)
@testset "Linear" begin
    Tests.linear_test(factory, config)
end
@testset "SOC" begin
    Tests.soc_test(factory, config, String[
        # FIXME MethodError: Cannot `convert` an object of type Nothing to an object of type MathOptInterface.VariableIndex
        "sdsos_horn"
    ])
end
@testset "SDP" begin
    Tests.sd_test(factory, config, [
        # Expression: termination_status(model) == MOI.OPTIMAL
        #  Evaluated: MathOptInterface.SLOW_PROGRESS == MathOptInterface.OPTIMAL
        # Expression: ≈(objective_value(model), α_value, atol=atol, rtol=rtol)
        #  Evaluated: 0.6847213936286884 ≈ 10.0 (atol=1.0e-5, rtol=1.0e-5)
        "BPT12e399_maxdegree",
        # FIXME AssertionError: m.x_type[ref2id(vi)] == Deleted
        "maxcut", "sos_horn", "motzkin",
        # Expression: JuMP.termination_status(model) == MOI.INFEASIBLE
        # Evaluated: MathOptInterface.SLOW_PROGRESS == MathOptInterface.INFEASIBLE
        # Contacted Mosek and they replied that there is nothing wrong as the
        # PrimalStatus and DualStatus are correct
        "quartic_infeasible_lyapunov_switched_system", "quartic_infeasible_scaled_lyapunov_switched_system"
    ])
end


text:
println("foo!")
println(STDERR, "bar!")
function calculate100(f, a, b)
    println("f: $f")
    println("a: $a, b: $b")
    return f.(linspace(a, b, 100))
end


text:
# todo change
__precompile__(false)

module Bayesian_EMR_MSM

using CmdStan
using LinearAlgebra
using Distributions
using Random
using DataFrames
using CSV

include("timeseries/timeseries.jl")
include("estimations/estimations.jl")
include("predictions/predictions.jl")

end # module


text:
module ColorVectorSpaceTests

using FactCheck, Base.Test, ColorVectorSpace, ColorTypes, FixedPointNumbers, Compat

macro test_colortype_approx_eq(a, b)
    :(test_colortype_approx_eq($(esc(a)), $(esc(b)), $(string(a)), $(string(b))))
end

u8sum(x,y) = Float64(UFixed8(x)) + Float64(UFixed8(y))

FactCheck.roughly(x::Gray) = (y::Gray) -> isapprox(y, x)
FactCheck.roughly(x::RGB) = (y::RGB) -> isapprox(y, x)

facts("Colortypes") do
    function test_colortype_approx_eq(a::Colorant, b::Colorant, astr, bstr)
        @fact typeof(a) --> typeof(b)
        n = length(fieldnames(typeof(a)))
        for i = 1:n
            @fact getfield(a, i) --> roughly(getfield(b,i))
        end
    end

    context("nan") do
        function make_checked_nan{T}(::Type{T})
            x = nan(T)
            isa(x, T) && isnan(x)
        end
        for S in (Float32, Float64)
            @fact make_checked_nan(S) --> true
            @fact make_checked_nan(Gray{S}) --> true
            @fact make_checked_nan(AGray{S}) --> true
            @fact make_checked_nan(GrayA{S}) --> true
            @fact make_checked_nan(RGB{S}) --> true
            @fact make_checked_nan(ARGB{S}) --> true
            @fact make_checked_nan(ARGB{S}) --> true
        end
    end

    context("Arithmetic with Gray") do
        cf = Gray{Float32}(0.1)
        ccmp = Gray{Float32}(0.2)
        @fact 2*cf --> ccmp
        @fact cf*2 --> ccmp
        @fact ccmp/2 --> cf
        @fact 2.0f0*cf --> ccmp
        @test_colortype_approx_eq cf*cf Gray{Float32}(0.01)
        @test_colortype_approx_eq cf^2 Gray{Float32}(0.01)
        @test_colortype_approx_eq cf^3.0f0 Gray{Float32}(0.001)
        @fact eltype(2.0*cf) --> Float64
        @fact abs2(ccmp) --> 0.2f0^2
        @fact norm(cf) --> 0.1f0
        @fact sumabs2(ccmp) --> 0.2f0^2
        cu = Gray{U8}(0.1)
        @fact 2*cu --> Gray(2*cu.val)
        @fact 2.0f0*cu --> Gray(2.0f0*cu.val)
        f = U8(0.5)
        @fact (f*cu).val --> roughly(f*cu.val)
        @fact 2.*cf --> ccmp
        @fact cf.*2 --> ccmp
        @fact cf/2.0f0 --> Gray{Float32}(0.05)
        @fact cu/2 --> Gray(cu.val/2)
        @fact cu/0.5f0 --> Gray(cu.val/0.5f0)
        @fact cf+cf --> ccmp
        @fact cf --> isfinite
        @fact cf --> not(isinf)
        @fact cf --> not(isnan)
        @fact Gray(NaN) --> not(isfinite)
        @fact Gray(NaN) --> not(isinf)
        @fact Gray(NaN) --> isnan
        @fact Gray(Inf) --> not(isfinite)
        @fact Gray(Inf) --> isinf
        @fact Gray(Inf) --> not(isnan)
        @fact abs(Gray(0.1)) --> roughly(0.1)
        @fact eps(Gray{U8}) --> Gray(eps(U8))  # #282

        acu = Gray{U8}[cu]
        acf = Gray{Float32}[cf]
        @fact typeof(acu+acf) --> Vector{Gray{Float32}}
        @fact typeof(acu-acf) --> Vector{Gray{Float32}}
        @fact typeof(acu.+acf) --> Vector{Gray{Float32}}
        @fact typeof(acu.-acf) --> Vector{Gray{Float32}}
        @fact typeof(acu+cf) --> Vector{Gray{U8}}
        @fact typeof(acu-cf) --> Vector{Gray{U8}}
        @fact typeof(acu.+cf) --> Vector{Gray{U8}}
        @fact typeof(acu.-cf) --> Vector{Gray{U8}}
        @fact typeof(2*acf) --> Vector{Gray{Float32}}
        @fact typeof(2.*acf) --> Vector{Gray{Float32}}
        @fact typeof(0x02*acu) --> Vector{Gray{Float32}}
        @fact typeof(acu/2) --> Vector{Gray{typeof(U8(0.5)/2)}}
        @fact typeof(acf.^2) --> Vector{Gray{Float32}}
        @fact (acu/Gray{U8}(0.5))[1] --> gray(acu[1])/U8(0.5)
        @fact (acf/Gray{Float32}(2))[1] --> roughly(0.05f0)
        @fact (acu/2)[1] --> Gray(gray(acu[1])/2)
        @fact (acf/2)[1] --> roughly(Gray{Float32}(0.05f0))
        @fact sumabs2([cf,ccmp]) --> roughly(0.05f0)

        @fact gray(0.8) --> 0.8

        a = Gray{U8}[0.8,0.7]
        @fact sum(a) --> Gray(u8sum(0.8,0.7))
        @fact abs( var(a) - (a[1]-a[2])^2 / 2 ) --> less_than(0.001)
        @fact isapprox(a, a) --> true
        @fact real(Gray{Float32}) <: Real --> true
        @fact zero(ColorTypes.Gray)-->0
        @fact one(ColorTypes.Gray)-->1
        a = Gray{U8}[0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
        @fact histrange(a,10)-->0.1:0.1:1

    end

    context("Comparisons with Gray") do
        g1 = Gray{U8}(0.2)
        g2 = Gray{U8}(0.3)
        @fact isless(g1, g2) --> true
        @fact isless(g2, g1) --> false
        @fact g1 < g2 --> true
        @fact g2 < g1 --> false
        @fact isless(g1, 0.5) --> true
        @fact isless(0.5, g1) --> false
        @fact g1 < 0.5 --> true
        @fact 0.5 < g1 --> false
        @fact @inferred(max(g1, g2)) --> g2
        @fact @inferred(max(g1, 0.1)) --> 0.2
        @fact @inferred(min(g1, g2)) --> g1
        @fact @inferred(min(g1, 0.1)) --> 0.1
        a = Gray{Float64}(0.9999999999999999)
        b = Gray{Float64}(1.0)

        @fact isapprox(a, b) --> true
        a = Gray{Float64}(0.99)
        @fact isapprox(a, b, rtol = 0.01) --> false
        @fact isapprox(a, b, rtol = 0.1) --> true
    end

    context("Unary operations with Gray") do
        for g in (Gray(0.4), Gray{U8}(0.4))
            for op in ColorVectorSpace.unaryOps
                try
                    v = @eval $op(gray(g))  # if this fails, don't bother
                    @fact $op(g) --> v
                end
            end
        end
        u = U8(0.4)
        @fact ~Gray(u) --> Gray(~u)
        @fact -Gray(u) --> Gray(-u)
    end

    context("Arithmetic with GrayA") do
        p1 = GrayA{Float32}(Gray(0.8), 0.2)
        p2 = GrayA{Float32}(Gray(0.6), 0.3)
        @test_colortype_approx_eq p1+p2 GrayA{Float32}(Gray(1.4),0.5)
        @test_colortype_approx_eq (p1+p2)/2 GrayA{Float32}(Gray(0.7),0.25)
        @test_colortype_approx_eq 0.4f0*p1+0.6f0*p2 GrayA{Float32}(Gray(0.68),0.26)
        @test_colortype_approx_eq ([p1]+[p2])[1] GrayA{Float32}(Gray(1.4),0.5)
        @test_colortype_approx_eq ([p1].+[p2])[1] GrayA{Float32}(Gray(1.4),0.5)
        @test_colortype_approx_eq ([p1]+p2)[1] GrayA{Float32}(Gray(1.4),0.5)
        @test_colortype_approx_eq ([p1].+p2)[1] GrayA{Float32}(Gray(1.4),0.5)
        @test_colortype_approx_eq ([p1]-[p2])[1] GrayA{Float32}(Gray(0.2),-0.1)
        @test_colortype_approx_eq ([p1].-[p2])[1] GrayA{Float32}(Gray(0.2),-0.1)
        @test_colortype_approx_eq ([p1]-p2)[1] GrayA{Float32}(Gray(0.2),-0.1)
        @test_colortype_approx_eq ([p1].-p2)[1] GrayA{Float32}(Gray(0.2),-0.1)
        @test_colortype_approx_eq ([p1]/2)[1] GrayA{Float32}(Gray(0.4),0.1)
        @test_colortype_approx_eq (0.4f0*[p1]+0.6f0*[p2])[1] GrayA{Float32}(Gray(0.68),0.26)

        a = GrayA{U8}[GrayA(0.8,0.7), GrayA(0.5,0.2)]
        @fact sum(a) --> GrayA(u8sum(0.8,0.5), u8sum(0.7,0.2))
        @fact isapprox(a, a) --> true
        a = AGray{Float64}(1.0, 0.9999999999999999)
        b = AGray{Float64}(1.0, 1.0)

        @fact isapprox(a, b) --> true
        a = AGray{Float64}(1.0, 0.99)
        @fact isapprox(a, b, rtol = 0.01) --> false
        @fact isapprox(a, b, rtol = 0.1) --> true

    end

    context("Arithemtic with RGB") do
        cf = RGB{Float32}(0.1,0.2,0.3)
        ccmp = RGB{Float32}(0.2,0.4,0.6)
        @fact 2*cf --> ccmp
        @fact cf*2 --> ccmp
        @fact ccmp/2 --> cf
        @fact 2.0f0*cf --> ccmp
        @fact eltype(2.0*cf) --> Float64
        cu = RGB{U8}(0.1,0.2,0.3)
        @test_colortype_approx_eq 2*cu RGB(2*cu.r, 2*cu.g, 2*cu.b)
        @test_colortype_approx_eq 2.0f0*cu RGB(2.0f0*cu.r, 2.0f0*cu.g, 2.0f0*cu.b)
        f = U8(0.5)
        @fact (f*cu).r --> roughly(f*cu.r)
        @fact 2.*cf --> ccmp
        @fact cf.*2 --> ccmp
        @fact cf/2.0f0 --> RGB{Float32}(0.05,0.1,0.15)
        @fact cu/2 --> RGB(cu.r/2,cu.g/2,cu.b/2)
        @fact cu/0.5f0 --> RGB(cu.r/0.5f0, cu.g/0.5f0, cu.b/0.5f0)
        @fact cf+cf --> ccmp
        @fact cu * 1//2 --> roughly(RGB{Float64}(U8(0.1)/2, U8(0.2)/2, U8(0.3)/2))
        @test_colortype_approx_eq (cf*[0.8f0])[1] RGB{Float32}(0.8*0.1,0.8*0.2,0.8*0.3)
        @test_colortype_approx_eq ([0.8f0]*cf)[1] RGB{Float32}(0.8*0.1,0.8*0.2,0.8*0.3)
        @test_colortype_approx_eq (cf.*[0.8f0])[1] RGB{Float32}(0.8*0.1,0.8*0.2,0.8*0.3)
        @test_colortype_approx_eq ([0.8f0].*cf)[1] RGB{Float32}(0.8*0.1,0.8*0.2,0.8*0.3)
        @fact cf --> isfinite
        @fact cf --> not(isinf)
        @fact cf --> not(isnan)
        @fact RGB(NaN, 1, 0.5) --> not(isfinite)
        @fact RGB(NaN, 1, 0.5) --> not(isinf)
        @fact RGB(NaN, 1, 0.5) --> isnan
        @fact RGB(1, Inf, 0.5) --> not(isfinite)
        @fact RGB(1, Inf, 0.5) --> isinf
        @fact RGB(1, Inf, 0.5) --> not(isnan)
        @fact abs(RGB(0.1,0.2,0.3)) --> roughly(0.6)
        @fact sumabs2(RGB(0.1,0.2,0.3)) --> roughly(0.14)
        @fact norm(RGB(0.1,0.2,0.3)) --> roughly(sqrt(0.14))

        acu = RGB{U8}[cu]
        acf = RGB{Float32}[cf]
        @fact typeof(acu+acf) --> Vector{RGB{Float32}}
        @fact typeof(acu-acf) --> Vector{RGB{Float32}}
        @fact typeof(acu.+acf) --> Vector{RGB{Float32}}
        @fact typeof(acu.-acf) --> Vector{RGB{Float32}}
        @fact typeof(acu+cf) --> Vector{RGB{U8}}
        @fact typeof(acu-cf) --> Vector{RGB{U8}}
        @fact typeof(acu.+cf) --> Vector{RGB{U8}}
        @fact typeof(acu.-cf) --> Vector{RGB{U8}}
        @fact typeof(2*acf) --> Vector{RGB{Float32}}
        @fact typeof(convert(UInt8, 2)*acu) --> Vector{RGB{Float32}}
        @fact typeof(acu/2) --> Vector{RGB{typeof(U8(0.5)/2)}}

        a = RGB{U8}[RGB(1,0,0), RGB(1,0.8,0)]
        @fact sum(a) --> RGB(2.0,0.8,0)
        @fact isapprox(a, a) --> true
        a = RGB{Float64}(1.0, 1.0, 0.9999999999999999)
        b = RGB{Float64}(1.0, 1.0, 1.0)

        @fact isapprox(a, b) --> true
        a = RGB{Float64}(1.0, 1.0, 0.99)
        @fact isapprox(a, b, rtol = 0.01) --> false
        @fact isapprox(a, b, rtol = 0.1) --> true
    end

    context("Arithemtic with RGBA") do
        cf = RGBA{Float32}(0.1,0.2,0.3,0.4)
        ccmp = RGBA{Float32}(0.2,0.4,0.6,0.8)
        @fact 2*cf --> ccmp
        @fact cf*2 --> ccmp
        @fact ccmp/2 --> cf
        @fact 2.0f0*cf --> ccmp
        @fact eltype(2.0*cf) --> Float64
        cu = RGBA{U8}(0.1,0.2,0.3,0.4)
        @test_colortype_approx_eq 2*cu RGBA(2*cu.r, 2*cu.g, 2*cu.b, 2*cu.alpha)
        @test_colortype_approx_eq 2.0f0*cu RGBA(2.0f0*cu.r, 2.0f0*cu.g, 2.0f0*cu.b, 2.0f0*cu.alpha)
        f = U8(0.5)
        @fact (f*cu).r --> roughly(f*cu.r)
        @fact 2.*cf --> ccmp
        @fact cf.*2 --> ccmp
        @fact cf/2.0f0 --> RGBA{Float32}(0.05,0.1,0.15,0.2)
        @fact cu/2 --> RGBA(cu.r/2,cu.g/2,cu.b/2,cu.alpha/2)
        @fact cu/0.5f0 --> RGBA(cu.r/0.5f0, cu.g/0.5f0, cu.b/0.5f0, cu.alpha/0.5f0)
        @fact cf+cf --> ccmp
        @test_colortype_approx_eq (cf*[0.8f0])[1] RGBA{Float32}(0.8*0.1,0.8*0.2,0.8*0.3,0.8*0.4)
        @test_colortype_approx_eq ([0.8f0]*cf)[1] RGBA{Float32}(0.8*0.1,0.8*0.2,0.8*0.3,0.8*0.4)
        @test_colortype_approx_eq (cf.*[0.8f0])[1] RGBA{Float32}(0.8*0.1,0.8*0.2,0.8*0.3,0.8*0.4)
        @test_colortype_approx_eq ([0.8f0].*cf)[1] RGBA{Float32}(0.8*0.1,0.8*0.2,0.8*0.3,0.8*0.4)
        @fact cf --> isfinite
        @fact cf --> not(isinf)
        @fact cf --> not(isnan)
        @fact RGBA(NaN, 1, 0.5, 0.8) --> not(isfinite)
        @fact RGBA(NaN, 1, 0.5) --> not(isinf)
        @fact RGBA(NaN, 1, 0.5) --> isnan
        @fact RGBA(1, Inf, 0.5) --> not(isfinite)
        @fact RGBA(1, Inf, 0.5) --> isinf
        @fact RGBA(1, Inf, 0.5) --> not(isnan)
        @fact RGBA(0.2, 1, 0.5, NaN) --> not(isfinite)
        @fact RGBA(0.2, 1, 0.5, NaN) --> not(isinf)
        @fact RGBA(0.2, 1, 0.5, NaN) --> isnan
        @fact RGBA(0.2, 1, 0.5, Inf) --> not(isfinite)
        @fact RGBA(0.2, 1, 0.5, Inf) --> isinf
        @fact RGBA(0.2, 1, 0.5, Inf) --> not(isnan)
        @fact abs(RGBA(0.1,0.2,0.3,0.2)) --> roughly(0.8)

        acu = RGBA{U8}[cu]
        acf = RGBA{Float32}[cf]
        @fact typeof(acu+acf) --> Vector{RGBA{Float32}}
        @fact typeof(acu-acf) --> Vector{RGBA{Float32}}
        @fact typeof(acu.+acf) --> Vector{RGBA{Float32}}
        @fact typeof(acu.-acf) --> Vector{RGBA{Float32}}
        @fact typeof(acu+cf) --> Vector{RGBA{U8}}
        @fact typeof(acu-cf) --> Vector{RGBA{U8}}
        @fact typeof(acu.+cf) --> Vector{RGBA{U8}}
        @fact typeof(acu.-cf) --> Vector{RGBA{U8}}
        @fact typeof(2*acf) --> Vector{RGBA{Float32}}
        @fact typeof(convert(UInt8, 2)*acu) --> Vector{RGBA{Float32}}
        @fact typeof(acu/2) --> Vector{RGBA{typeof(U8(0.5)/2)}}

        a = RGBA{U8}[RGBA(1,0,0,0.8), RGBA(0.7,0.8,0,0.9)]
        @fact sum(a) --> RGBA(u8sum(1,0.7),0.8,0,u8sum(0.8,0.9))
        @fact isapprox(a, a) --> true
        a = ARGB{Float64}(1.0, 1.0, 1.0, 0.9999999999999999)
        b = ARGB{Float64}(1.0, 1.0, 1.0, 1.0)

        @fact isapprox(a, b) --> true
        a = ARGB{Float64}(1.0, 1.0, 1.0, 0.99)
        @fact isapprox(a, b, rtol = 0.01) --> false
        @fact isapprox(a, b, rtol = 0.1) --> true
    end
end

isinteractive() || FactCheck.exitstatus()

end


text:
mutable struct UMVE{T<:Real}
    A::Matrix{T}
    B::Matrix{T}
    C::OutputMatrix{T}
    G::Matrix{T}
    W::Matrix{T}
    V::Matrix{T}
    Ebar::Matrix{T}
    z::Vector{T}
    xi::Vector{T}
    Pi::Matrix{T} # state estimation error covariance
    Delta::Union{Matrix{T}, Missing} # input estimation error covariance

    function UMVE(A::S,
        B::S,
        C::OutputMatrix{T},
        E::S,
        W::AbstractMatrix{T},
        V::AbstractMatrix{T},
        decompose::Bool = true) where {T<:Real, S<:Matrix{T}}
        if isposdef(-W) #this does not exist
            throw(DomainError(W, "W must be positive definite"))
        end
        if !isposdef(V)
            throw(DomainError(V, "V must be positive definite"))
        end
        W = Matrix(W)
        V = Matrix(V)
        feasible, G, Ebar = checkUMVEFeasibility(C,E, decompose)
        nx = size(A,1)
        nxi = size(G,2)
        if feasible
            Pk = W
            Rk = R(A, C, Pk, W, V)
            Mk = M(G, C, Rk)
            new{T}(
                A,
                B,
                C,
                G,
                W,
                V,
                Ebar,
                zeros(T, nx),
                zeros(T, nxi),
                W,
                missing
            )
        else
            error("unfeasible")
        end
    end
end

function checkUMVEFeasibility(C::OutputMatrix, E::Matrix, decompose::Bool)
    rE = rank(E)
    feasible = rank(C*E) == rE
    if decompose
        SVD = svd(E)
        G = SVD.U[:,1:rE]
        Ebar = (Diagonal(SVD.S) * SVD.Vt)[1:rE,:]
    else
        nz = size(E,1)
        G = isempty(E) ? zeros(nz, 0) : Matrix{eltype(E)}(I, nz, nz)
        Ebar = E
    end
    feasible, G, Ebar
end

R(A::T, C::OutputMatrix{S}, Pk::T, W::T, V::T) where {S<:Real, T<:Matrix{S}} = C*(A*Pk*A' + W)*C' + V
M(G::T, C::OutputMatrix{S}, Rk::T) where {S<:Real, T<:Matrix{S}} = inv( G'C'*inv(Rk)*C*G ) * G'C'inv(Rk)
K(A::T, C::OutputMatrix{S}, Pk::T, W::T, Rk::T) where {S<:Real, T<:Matrix{S}} = C*(A*Pk*A' + W)*C'inv(Rk)
barA(Kk::T, C::OutputMatrix{S}, G::T, Mk::T) where {S<:Real, T<:Matrix{S}} = (I - Kk*C)*(I - G*Mk*C)
barL(Kk::T, C::OutputMatrix{S}, G::T, Mk::T) where {S<:Real, T<:Matrix{S}} = Kk + (I - Kk*C)*G*Mk

function updateState(o::UMVE{T}, u::Vector{T}, y::Vector{T}) where {T}
    Rk = R(o.A, o.C, o.Pi, o.W, o.V)
    # setfield!(o, :R, Rk)
    Mk = M(o.G, o.C, Rk)
    setfield!(o, :Delta, Mk*Rk*Mk')
    Kk = K(o.A, o.C, o.Pi, o.W, Rk)
    Abar = barA(Kk, o.C, o.G, Mk)
    Lbar = barL(Kk, o.C, o.G, Mk)
    biased_z = o.A*o.z + o.B*u
    setfield!(o, :z, Abar*biased_z + Lbar*y)
    setfield!(o, :xi, Mk*(y - o.C*biased_z))
    setfield!(o, :Pi, Abar*o.A*o.Pi*o.A'Abar' + Abar*o.W*Abar' + Lbar*o.V*Lbar')
end
updateState(o::UMVE{T}, u::T, y::Vector{T}) where {T} = updateState(o, [u], y)

outputMap(o::UMVE) = o.C*o.z

function Base.getproperty(o::UMVE, s::Symbol)
    if s === :nz
        return length(o.z)
    elseif s === :nxi
        return length(o.xi)
    else
        return getfield(o, s)
    end
end


text:
const reweightconf = Dict{Symbol,Any}(
:type=>:beammodel,
:method=>:std,
:raymethod=>:bresenham,
#:max_range_meters=>20,
:max_range_meters=>5.6,
:exclude_max_range=>true,
:inv_squash_factor=>2,
:a_short_exp=>0.4,
:a_short_lin=>0.2,
:a_max=>2,
:a_rand=>:0.7,
:a_hit=>0.75,
:sigma_hit=>0.07,
:lambda_short=>1/0.6,
:nrays=>50)

const predictconf = Dict{Symbol,Any}(
:type=>:ackermann,
:sigma_x=>0.01,
:sigma_y=>0.01,
:sigma_theta=>0.05,
:init_sigma_x=>0.3,
:init_sigma_y=>0.3,
:init_sigma_theta=>0.05,
:vel_sigma=>0.000001,
:delta_sigma=>0.000001,
:car_length=>0.33,
:stochastic=>true)

const pfconf=Dict{Symbol,Any}(:max_particles=>1000, :dtype=>Float64, :rng_seed=>1234, :resampler=>:naive, :predictconf=>predictconf, :reweightconf=>reweightconf)

const vescconf = Dict{Symbol, Any}(
:speed2erpm_offset=>0.0,
:speed2erpm_gain=>-1900.0,
#:speed2erpm_gain=>4614.0,
#:steering2servo_offset=>0.5304,
:steering2servo_offset=>0.42,
:steering2servo_gain=>-1.2135)

const rosconf = Dict{Symbol,Any}(
:max_viz_particles=>50,
:static_map_service=>"/static_map",
:ackermann_topic=>"/mux/ackermann_cmd_mux/output",
:servo_pos_topic=>"/vesc/sensors/servo_position_command",
:vesc_topic=>"/vesc/sensors/core",
:pose_topic=>"/pf/inferred_pose",
:particle_topic=>"/pf/viz/particles",
:scan_topic=>"/scan",
:initialpose_topic=>"/initialpose",
:odom_topic=>"/odom",
:vescconf=>vescconf,
:laserpose_topic=>"/pf/viz/laserpose")

const conf = Dict{Symbol,Any}(:rosconf=>rosconf, :pfconf=>pfconf)


text:
"""
    College(j, f, t, ft, omf)

Contains a college's index `j`, admissions probability `f`, utility value `t`,
their product `ft`, and `1 - f = omf`. Used only by `applicationorder()`.
"""
struct College
    j::Int
    f::Float64
    t::Float64
    ft::Float64
    omf::Float64
end

# Overload isless() so that the heap is ordered by expected utility.
isless(c1::College, c2::College) = isless(c1.ft, c2.ft)

"""
    applicationorder_list(mkt::SameCostsMarket) -> X, V

Produce the optimal application order `X` and associated valuations `V`
for the [`SameCostsMarket`](@ref) defined by `mkt`. Uses a list data structure;
typically faster than the equivalent [`applicationorder_heap`](@ref).

All `SameCostsMarket`s satisfy a nestedness property, meaning that the optimal
portfolio when `mkt.h = h` is given by the first `h` entries of the optimal portfolio
when `mkt.h = mkt.m`. For example:

```julia-repl
julia> mkt = SameCostsMarket([0.2, 0.5, 0.1, 0.6, 0.1], [1, 4, 9, 1, 8], 5);

julia> x, v = applicationorder_list(mkt)
([2, 3, 5, 4, 1], [2.0, 2.7, 3.24, 3.483, 3.5154])

julia> x[1:4], v[4] # optimal portfolio and valuation for h = 4
([2, 3, 5, 4], 3.483)
```
"""
function applicationorder_list(
    mkt::SameCostsMarket;
    verbose::Bool=false
)::Tuple{Vector{Int},Vector{Float64}}
    apporder = zeros(Int, mkt.h)
    v = zeros(mkt.h)
    mkt_list = College[College(j, mkt.f[j], mkt.t[j], mkt.ft[j], mkt.omf[j]) for j in 1:mkt.m]

    c_best::College, idx_best::Int = findmax(mkt_list)
    @inbounds for j in 1:mkt.h
        if verbose
            println("Iteration $j")
            println("  ft: ", [mkt_list[i].ft for i in 1:length(mkt_list)])
            println("  Add school $(c_best.j)")
        end

        v[j] = get(v, j - 1, 0) + c_best.ft
        apporder[j] = c_best.j

        next_c_best::College = College(0, 1.0, -1.0, -1.0, 0.0)
        next_idx_best::Int = 0

        for i in 1:idx_best-1
            mkt_list[i] =
                College(
                    mkt_list[i].j,
                    mkt_list[i].f,
                    mkt_list[i].t * c_best.omf,
                    mkt_list[i].ft * c_best.omf,
                    mkt_list[i].omf
                )

            if isless(next_c_best, mkt_list[i])
                next_c_best = mkt_list[i]
                next_idx_best = i
            end
        end

        for i in idx_best+1:length(mkt_list)
            mkt_list[i-1] =
                College(
                    mkt_list[i].j,
                    mkt_list[i].f,
                    mkt_list[i].t - c_best.ft,
                    mkt_list[i].ft - mkt_list[i].f * c_best.ft,
                    mkt_list[i].omf
                )

            if isless(next_c_best, mkt_list[i-1])
                next_c_best = mkt_list[i-1]
                next_idx_best = i - 1
            end
        end

        c_best, idx_best = next_c_best, next_idx_best

        pop!(mkt_list)
    end

    return mkt.perm[apporder], v
end

"""
    applicationorder_list(mkt::SameCostsMarket) -> X, V

Produce the optimal application order `X` and associated valuations `V`
for the [`SameCostsMarket`](@ref) defined by `mkt`. Uses a heap data structure;
typically the equivalent [`applicationorder_list`](@ref) is faster.

All `SameCostsMarket`s satisfy a nestedness property, meaning that the optimal
portfolio when `mkt.h = h` is given by the first `h` entries of the optimal portfolio
when `mkt.h = mkt.m`. For example:

```julia-repl
julia> mkt = SameCostsMarket([0.2, 0.5, 0.1, 0.6, 0.1], [1, 4, 9, 1, 8], 5);

julia> x, v = applicationorder_heap(mkt)
([2, 3, 5, 4, 1], [2.0, 2.7, 3.24, 3.483, 3.5154])

julia> x[1:4], v[4] # optimal portfolio and valuation for h = 4
([2, 3, 5, 4], 3.483)
```
"""
function applicationorder_heap(mkt::SameCostsMarket)::Tuple{Vector{Int},Vector{Float64}}
    apporder = zeros(Int, mkt.h)
    v = zeros(mkt.h)

    mkt_heap = BinaryMaxHeap{College}(collect(College(j, mkt.f[j], mkt.t[j], mkt.ft[j], mkt.omf[j]) for j in 1:mkt.m))

    @inbounds for j in 1:mkt.h
        c_k = first(mkt_heap)
        v[j] = get(v, j - 1, 0) + c_k.ft
        apporder[j] = c_k.j

        # These two implementations perform almost identically
        mkt_heap = BinaryMaxHeap{College}(
            # College{T}[
            #     College(
            #         c.j,
            #         c.f,
            #         c.j < c_k.j ? c.t * c_k.omf : c.t - c_k.ft,
            #         c.j < c_k.j ? c.ft * c_k.omf : c.ft - c.f * c_k.ft,
            #         c.omf
            #     )
            #     for c in mkt_heap.valtree if c.j != c_k.j
            # ]
            map(filter(c -> c.j != c_k.j, mkt_heap.valtree)) do c
                College(
                    c.j,
                    c.f,
                    c.j < c_k.j ? c.t * c_k.omf : c.t - c_k.ft,
                    c.j < c_k.j ? c.ft * c_k.omf : c.ft - c.f * c_k.ft,
                    c.omf
                )
            end
        )
    end

    return mkt.perm[apporder], v
end

# Used by dynamic program below
@inbounds function V_recursor!(
    V_dict::Dict{Tuple{Int,Int},Float64},
    j::Int,
    h::Int,
    mkt::VariedCostsMarket
)::Float64
    return get(V_dict, (j, h)) do
        if j == 0 || h == 0
            return 0.0
        elseif h < mkt.g[j]
            push!(V_dict, (j, h) => V_recursor!(V_dict, j - 1, h, mkt))
            return V_dict[(j, h)]
        else
            jmo = j - 1
            push!(V_dict, (j, h) => max(
                V_recursor!(V_dict, jmo, h, mkt),
                mkt.omf[j] * V_recursor!(V_dict, jmo, h - mkt.g[j], mkt) + mkt.ft[j]
            ))
            return V_dict[(j, h)]
        end
    end
end

"""
    optimalportfolio_dynamicprogram(mkt::VariedCostsMarket) -> X, v

Use the dynamic program on application costs to produce the optimal portfolio `X` and associated
value `v` for the [`VariedCostsMarket`](@ref) defined by `mkt`.

```julia-repl
julia> mkt = VariedCostsMarket([0.2, 0.5, 0.1, 0.6, 0.1], [1, 4, 9, 1, 8], [2, 4, 2, 5, 1], 8);

julia> optimalportfolio_dynamicprogram(mkt)
([3, 5, 2], 3.24)
```
"""
function optimalportfolio_dynamicprogram(
    mkt::VariedCostsMarket;
    verbose::Bool=false
)::Tuple{Vector{Int},Float64}
    V_dict = Dict{Tuple{Int,Int},Float64}()
    sizehint!(V_dict, mkt.m * mkt.m ÷ 2)

    v = V_recursor!(V_dict, mkt.m, mkt.H, mkt)

    h = mkt.H
    X = Int[]

    for j in reverse(1:mkt.m)
        if V_recursor!(V_dict, j - 1, h, mkt) < V_recursor!(V_dict, j, h, mkt)
            push!(X, j)
            h -= mkt.g[j]
        end
    end

    if verbose
        V_table = Array{Union{Missing,Float64}}(missing, mkt.m, mkt.H)
        for (j, h) in keys(V_dict)
            if 0 < j ≤ mkt.m && 0 < h ≤ mkt.H
                V_table[j, h] = V_dict[(j, h)]
            end
        end

        display(V_table)
    end

    return mkt.perm[X], v
end

max_stars_repo_name: maxkapur/OptimalApplication

max_issues_repo_name: maxkapur/OptimalApplication

max_forks_repo_name: maxkapur/OptimalApplication


text:
import MD5
using BenchmarkTools

data = randstring(10^4)

if Pkg.installed("Nettle") != nothing
import Nettle
@show Nettle.hexdigest("md5", data)
display(@benchmark Nettle.hexdigest("md5", data))
    println()
end

@show bytes2hex(MD5.md5(data))
display(@benchmark MD5.md5(data))
println()


text:
using Clp,Cbc

JuDGE_SP_Solver = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0, "ratioGap" => 0.0, "presolve" => "off", "preprocess" => "off")

JuDGE_MP_Solver = (optimizer_with_attributes(Clp.Optimizer, "LogLevel" => 0, "Algorithm" => 4),
				   JuDGE_SP_Solver)

JuDGE_DE_Solver = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 1, "ratioGap" => 0.0)


text:
struct AuAu <: JuLIP.SitePotential end

JuLIP.cutoff(::AuAu) = JuLIP.lattice_constant(:Au)

function JuLIP.evaluate!(tmp, ::AuAu, Rs, Zs, z0)
    Es = 0
    if z0 == 79 # Ensure current site is gold atom
        slice = (Zs .== 79) # Identify gold atoms
        AuRs = @view Rs[slice] # Slice gold atoms
        I = get_nearest_neighbours(AuRs)
        for R in @view AuRs[I]
            i = select_orientation_index(R)
            R -= direction_vectors[i]
            A = force_tensors[i]
            Es += R' * A * R
        end
    end
    return Es / 4
end

function JuLIP.evaluate_d!(dEs, tmp, ::AuAu, Rs, Zs, z0)
    if z0 == 79 # Ensure current site is gold atom
        slice = (Zs .== 79) # Identify gold atoms
        AuRs = @view Rs[slice] # Slice gold atoms
        I = get_nearest_neighbours(AuRs)
        for i in 1:length(dEs)
            dEs[i] = zero(dEs[i])
        end
        for (j, R) in zip(I, @view AuRs[I])
            i = select_orientation_index(R)
            R -= direction_vectors[i]
            A = force_tensors[i]
            dEs[j] = A*R / 2
        end
    end
    return dEs
end

function get_nearest_neighbours(AuRs)
    if length(AuRs) >= 12
        I = partialsortperm(AuRs, 1:12, by=x->norm(x))
    else
        I = partialsortperm(AuRs, 1:9, by=x->norm(x))
    end
    return I
end

"""
Identifies which force tensor is to be used based on symmetry considerations.

We calculate the dot product between the current vector and all possible basis vectors.
The largest dot product identifies the correct tensor.
"""
function select_orientation_index(vec)
    _, index = findmax(dot.(Ref(vec), direction_vectors))
    return index
end