text:
module Points

export Point

using ..Elements
import ..Motions: induce_velocity, mutually_induce_velocity!, self_induce_velocity!, advect,
                  allocate_velocity

#== Type definition ==#

"""
    Point{T,R} <: Elements.Element

An immutable structure representing a point source/vortex

## Fields

- `z::ComplexF64`: position
- `S::T`: strength/circulation
"""
struct Point{S <: Number,R <: Real} <: Element
    z::Complex{R}
    S::Union{R,Complex{R}}
    Point{S}(z::Complex{R}, s::T)           where {S <: Complex, T <: Real, R <: Real} = (U = promote_type(R,T); new{Complex{U},U}(z, im*convert(U,s)))
    Point{S}(z::Complex{R}, s::Complex{T})  where {S <: Complex, T <: Real, R <: Real} = (U = promote_type(R,T); new{Complex{U},U}(z, convert(Complex{U},s)))
    Point{S}(z::Complex{R}, s::T)           where {S <: Real, T <: Real, R <: Real}    = (U = promote_type(R,T); new{U,U}(z, convert(U,s)))
end
# struct Point{T <: Number, R <: Number} <: Element
#     z::Complex{R}
#     S::T
#     Point{T,R}(z, s::Real) where {T <: Complex,R <: Number} = new(z, im*s)
#     Point{T,R}(z, s)       where {T <: Complex,R <: Number} = new(z, s)
#     Point{T,R}(z, s::Real) where {T <: Real, R <: Number} = new(z, s)
# end

#Point(z::Complex{R},S::T) where {T,R} = Point{T,R}(z,S)

Elements.kind(::Point) = Singleton
Elements.kind(::Type{Point{T,R}}) where {T,R} = Singleton

Elements.property_type(::Point{T,R}) where {T,R} = R
Elements.property_type(::Type{Point{T,R}}) where {T,R} = R

#Elements.promote_property_type(::Point{T,R}) where {T,R} = promote_type(T,R)
#Elements.promote_property_type(::Type{Point{T,R}}) where {T,R} = promote_type(T,R)

#== Methods to be extended ==#

Elements.position(p::Point) = p.z
Elements.blobradius(p::Point) = 0.0

Elements.streamfunction(z::Complex{T}, p::Point) where {T} = real(-0.5p.S*log(z - p.z)/π)

Elements.complexpotential(z::Complex{T}, p::Point) where {T} = -0.5im*p.S*log(z - p.z)/π

cauchy_kernel(z) = z != zero(z) ? 0.5im/(π*conj(z)) : zero(z)

# ensures that the velocity is of same type as position
allocate_velocity(v::Vector{Point{T,R}}) where {T,R} = zeros(Complex{R},length(v))

function induce_velocity(z::Complex{T}, p::Point, t) where {T}
    p.S'*cauchy_kernel(z - p.z)
end

function mutually_induce_velocity!(ws₁, ws₂,
                                   points₁::Vector{Point{T₁}},
                                   points₂::Vector{Point{T₂}}, t) where {T₁, T₂}
    for (s, source) in enumerate(points₁)
        for (t, target) in enumerate(points₂)
            K = cauchy_kernel(target.z - source.z)
            ws₂[t] += source.S'*K
            ws₁[s] -= target.S'*K
        end
    end
    nothing
end

function self_induce_velocity!(ws, points::Vector{Point{T,R}}, t) where {T,R}
    N = length(points)

    for s in 1:N, t in s+1:N
        K = cauchy_kernel(points[t].z - points[s].z)
        ws[t] += points[s].S'*K
        ws[s] -= points[t].S'*K
    end
    ws
end

function advect(p::Point{T}, w::ComplexF64, Δt::Float64) where {T}
    Point{T}(p.z + w*Δt, p.S)
end

#function Base.show(io::IO, p::Point)
#    print(io, "Vortex.Point(z = $(round(p.z, digits=3)), Γ = $(round(p.Γ, digits=3)))")
#end

end


text:
# Problems taken from
# Deb et. al - "Scalable Test Problems for Evolutionary Multiobjective Optimization"

# DTLZ8 & DTLZ9 not supported yet due to constraints...
export DTLZ1, DTLZ2, DTLZ3, DTLZ4, DTLZ5, DTLZ6, DTLZ7, DTLZ;
export get_ideal_point;

dtlz_ref_string = """Deb et. al - "Scalable Test Problems for Evolutionary Multiobjective Optimization"."""

#%% -------------- Problem definitions ----------------------------------#
abstract type DTLZ <: MOP end;

dtlz_descriptor( i, n_vars, M) = """
    struct DTLZ$i <: DTLZ
        n_vars :: Int = $n_vars;
        M :: Int = $M ;
        k :: Int = n_vars + 1 - M;
    end

Test problem defined in: $dtlz_ref_string\n
Has `n_vars` decision variables and `M` objectives.
The decision vector `x` is split into two parts: `x[1:M-1]` and `x[M:end]`.
The latter has length `k`.
"""

pareto_front_descriptions = Dict(
    1 => "DTLZ1 has an affine-linear global Pareto-optimal front with Σfᵢ=0.5 and
    **many** local minima because a function `g` akin to Rastrigin's is used.",
    2 => "DTLZ2 has spherical global Pareto-optimal front with Σfᵢ²=1.",
    3 => "DTLZ3 has spherical global Pareto-optimal front with Σfᵢ²=1 and
    **many** local minima because a function `g` akin to Rastrigin's is used.",
    4 => "DTLZ4 has spherical global Pareto-optimal front with Σfᵢ²=1 and
    solutions tend to lie near the ``f_M-f₁`` plane.",
    5 => "DTLZ5 has spherical Pareto-optimal front that is a curve.",
    6 => "DTLZ6 has spherical Pareto-optimal front that is a curve and hard to reach.",
    7 => "DTLZ7 has ``2^{M-1}`` disconnected Pareto-optimal regions in search space.
    It is similar in construction to the ZDT problems.",
);

# dtlz index => (n_vars, M)
dtlz_standard_settings = Dict(
    1 => (10,6),
    2 => (15,6),
    3 => (15,6),
    4 => (15,6),
    5 => (15,6),
    6 => (15,6),
    7 => (22,3)
)

for i = 1 : 7

    n_vars, M = dtlz_standard_settings[i];
    help_text = """
    $(dtlz_descriptor(1, n_vars, M))

    $(pareto_front_descriptions[i])
    """

    typename = Symbol("DTLZ", i);

    @eval begin
        @doc $help_text
        @with_kw struct $typename <: DTLZ
            n_vars :: Int = $(n_vars);
            M :: Int = $(M) ;
            k :: Int = n_vars + 1 - M;
            @assert 1 <= k;
            @assert 1 <= M <= n_vars;
            @assert n_vars == M - 1 + k;
        end
    end
end

#%% --------------------------------------------------------------------#
# some general properties

num_vars( mop :: D where D<:DTLZ) = mop.n_vars;
num_objectives( mop :: D where D<:DTLZ) = mop.M;

for i=1:6
    typename = Symbol("DTLZ", i)
    @eval get_ideal_point( mop :: $typename ) = zeros(mop.M);
end

function constraints( mop :: D where D<:DTLZ)
    Box(
        zeros( num_vars( mop ) ),
        ones( num_vars( mop ) )
    )
end

function get_pareto_set( mop :: D where D<:DTLZ)
    return SamplingFunction( mop, :ParetoSet, [:random] );
end

function get_pareto_front(  mop :: D where D<:DTLZ)
    return SamplingFunction( mop, :ParetoFront, [:random] );
end

function get_points( sample_func :: SamplingFunction, mop :: Union{DTLZ1,DTLZ6,DTLZ7},
    ::Val{:ParetoSet}, num_points :: Int; method = :random )
return [ [ rand(mop.M-1); zeros(mop.k) ] for i = 1 : num_points ]
end

function get_points( sample_func :: SamplingFunction, mop :: Union{DTLZ2,DTLZ3,DTLZ4,DTLZ5},
    ::Val{:ParetoSet}, num_points :: Int; method = :random )
return [ [ rand(mop.M-1); fill(.5, mop.k) ] for i = 1 : num_points ]
end

#%% ----------------------- Objectives --------------------------#
########## DTLZ 1 ##############

# concerning `g(x)` there is a inconsistency in the source:
# they use other factors but then the Pareto Front
# is no longer the hyperplane with Σ fⱼ = .5
function g(x, mop::Union{DTLZ1, DTLZ3})
    χ = x[mop.M : end];
    return (
        .75 * mop.k +
        sum( (χ .- .5 ).^2 ) -
        sum( cos.( (20 * pi) .* ( χ .- .5) ) )
    );
end

function get_objectives( mop :: DTLZ1 )
    objectives = Vector{Function}( undef, mop.M )
    for j = 1 : mop.M
        fⱼ = function( x :: Vector{R} where R <: Real)
            return (
                .5 * prod( x[1 : mop.M-j] ) *
                (1 + g( x , mop )) *
                (1 - x[mop.M-j+1])^(j>1)
            )
        end
        objectives[j] = fⱼ
    end
    return objectives
end

########## DTLZ 2 ##############

θ(x, :: Union{DTLZ2, DTLZ3} ) = (pi/2) .* x;
function g(x, mop:: Union{DTLZ2,DTLZ4,DTLZ5} )
    χ = x[mop.M : end];
    return sum( (χ .- 0.5 ).^2 );
end

function get_objectives( mop :: Union{DTLZ2,DTLZ3,DTLZ4,DTLZ5,DTLZ6} )
    objectives = Vector{Function}( undef, mop.M )
    for j = 1 : mop.M
        fⱼ = function( x :: Vector{R} where R <: Real)
            return (
                (1 + g( x, mop ) ) *
                prod( cos.( θ(x,mop)[1 : mop.M - j] ) ) *
                sin( θ(x,mop)[mop.M-j+1] )^(j>1)
            )
        end
        objectives[j] = fⱼ
    end
    return objectives
end
########## DTLZ 3 ##############

# `θ` as for DTLZ2
# `g` as for DTLZ1

# `get_objectives` defined above

########## DTLZ 4 ##############

θ(x, :: DTLZ4) = (pi/2) .* x.^100;
# `g` same as for DTLZ2 & DTLZ5

# `get_objectives` defined above

########## DTLZ 5 ##############

# `g` same as for DTLZ2 & DTLZ4

function θ(x, mop :: Union{DTLZ5, DTLZ6} )
    G = g( x, mop );
    return [
        (pi/2) * x[1];
        ( pi / ( 4 * (1 + G ) ) ) .* ( 1 .+ (2*G) .* x[2: mop.M - 1] ) ;
        x[mop.M : end]
    ];
end

# `get_objectives` defined above

########## DTLZ 6 ##############

function g(x, mop::DTLZ6)
    χ = x[mop.M : end];
    return sum( χ.^0.1 )
end
# `θ` same as for DTLZ5

# `get_objectives` defined above

########## DTLZ 7 ##############

function g(x, mop::DTLZ7)
    χ = x[mop.M : end];
    return 1 + (9/mop.k) * sum( χ )
end

function get_objectives( mop::DTLZ7 )
    objectives = Vector{Function}( undef, mop.M )

    h(F, g) = mop.M - sum( F[1 : mop.M-1] ./ (1 + g) .* ( 1 .+ sin.( (3*pi) .* F[ 1:mop.M-1 ] ) ) );

    for j=1:mop.M-1
        fⱼ = function( x :: Vector{R} where R<:Real )
            x[j]
        end
        objectives[j] = fⱼ;
    end

    objectives[mop.M] = function( x :: Vector{R} where R<:Real )
        F = [ objectives[j](x) for j=1:mop.M-1 ];
        G = g( x, mop );
        return (1+G) * h(F,G)
    end
    return objectives
end

# TODO check this
function get_ideal_point( mop::DTLZ7)
    [zeros( mop.M - 1); 1.0];
end
##########################################################################

ext: jl


text:
# https://adventofcode.com/2021/day/12
using AdventOfCode
using Graphs
using Unicode
using Combinatorics
using SimpleWeightedGraphs
input = '"'.*replace.(readlines("2021/data/day_12.txt"),"-"=>"""\"=>\"""").*'"' .|>parse
Unicode.isuppercase(str::AbstractString)=isuppercase(str[1])
strdict=Dict(unique([first.(input);last.(input)]).=> 1:length(unique([first.(input);last.(input)])))
uppers=filter(isuppercase,keys(strdict))
input=(edges=input,strdict=strdict,uppers=uppers)
Graphs.add_edge!(g::Graphs.SimpleGraphs.AbstractSimpleGraph, a::AbstractVector{<:Integer})=add_edge!(g, a[begin],a[begin+1])
function makegraph(input)
    strdict=input.strdict
    g=SimpleWeightedGraph(SimpleGraph(strdict|>length), 1)
    for edge in input.edges
        add_edge!(g,strdict[first(edge)],strdict[last(edge)],1)
    end
    for upvex in getindex.(Ref(strdict),uppers) # for each upper case vertice
        nb=neighbors(g,upvex)|>collect
        for edge in combinations(nb,2)
            add_edge!(g,edge[begin],edge[end],g.weights[edge[begin],edge[end]]+1) # connect among neighbors of upvex
        end
        for edge in zip(nb,nb)
            add_edge!(g,edge[begin],edge[end],g.weights[edge[begin],edge[end]]+1) # add self edges for neighbors
        end
        nb=neighbors(g,upvex)|>collect
        rem_edge!.(Ref(g),upvex,nb) # remove edges from upvex
    end
    g
end
function part_1(input)
    strdict=input.strdict
    g=makegraph(input)
    routes=combinations(setdiff(values(strdict),[strdict["start"];strdict["end"];getindex.(Ref(strdict),uppers)])).|>permutations
    npath = x->prod(getindex.(Ref(g.weights),[strdict["start"];x],[x;strdict["end"]]))
    sum(sum(npath, route) for route in routes)
end
onlyonetwicefilter(vectors)=filter(x->length(x)-length(unique(x))<2,vectors)
function part_2(input)
    strdict=input.strdict
    g=makegraph(input)
    caves=setdiff(values(strdict),[strdict["start"];strdict["end"];getindex.(Ref(strdict),uppers)])
    doublecaves=repeat(caves,2)
    routes=combinations(doublecaves).|>sort|>onlyonetwicefilter|>unique.|>permutations
    npath = x->prod(getindex.(Ref(g.weights),[strdict["start"];x],[x;strdict["end"]]))
    sum(sum(npath, unique(route)) for route in routes)+npath(Int[])
end


text:
module TextTiling
include("Utls.jl")
using Statistics
using .Utls

"""
    TextTiling.SegmentObject(window_size, smooth_window_size, tokenizer)
TextTiling is a method for finding segment boundaries based on lexical cohesion and similarity between adjacent blocks.
# Arguments
- `window_size`: Sliding window size.
- `smooth_window_size`: Window size for smoothing depth scores.
- `tokenizer`: Tokenizer for word segmentation.
"""
mutable struct SegmentObject
    window_size::Int
    smooth_window_size::Int
    tokenizer::Any
end

function preprocessing(seg::SegmentObject, document)
    n = length(document)
    @assert n > 0 && length([d for d in document if typeof(d) != String]) == 0
    seg.window_size = maximum([minimum([seg.window_size, n / 3]), 1])
    return [Utls.count_elements(seg.tokenizer(document[i])) for i = 1:n]
end

function calculate_gap_score(seg::SegmentObject, preprocessed_document)
    n = length(preprocessed_document)
    gap_score = [0.0 for _ = 1:n]

    for i = 1:n
        sz = minimum([minimum([i, n - i]), seg.window_size])
        left_side, right_side = Utls.SentenceElements(Dict{String,Int}()),
        Utls.SentenceElements(Dict{String,Int}())

        for j = Int(i - sz + 1):Int(i)
            Utls.merge_elements(left_side, preprocessed_document[j])
        end

        for j = Int(i + 1):Int(i + sz)
            Utls.merge_elements(right_side, preprocessed_document[j])
        end

        gap_score[i] =
            Utls.calculate_cosin_similarity(left_side.elements_dct, right_side.elements_dct)
    end
    return gap_score
end

function calculate_depth_score(seg::SegmentObject, gap_score)
    n = length(gap_score)
    depth_score = [0.0 for _ = 1:n]

    for i = 1:n
        if i <= seg.window_size || i + seg.window_size > n
            continue
        end

        ptr = i - 1
        while ptr >= 1 && gap_score[ptr] >= gap_score[ptr+1]
            ptr -= 1
        end
        lval = gap_score[ptr+1]

        ptr = i + 1
        while ptr < n && gap_score[ptr] >= gap_score[ptr-1]
            ptr += 1
        end
        rval = gap_score[ptr-1]

        depth_score[i] = 1 / 2 * (lval - gap_score[i] + rval - gap_score[i])
    end
    return depth_score
end

function smoothing(seg::SegmentObject, depth_score)
    n = length(depth_score)
    smooth_depth_score = [0.0 for _ = 1:n]

    for i = 1:n
        if i - seg.smooth_window_size < 1 || i + seg.smooth_window_size >= n
            smooth_depth_score[i] = depth_score[i]
        else
            smooth_depth_score[i] =
                mean(depth_score[(i-seg.smooth_window_size):(i+seg.smooth_window_size)])
        end
    end
    return smooth_depth_score
end

function determine_boundaries(seg::SegmentObject, smooth_depth_score, num_topics)
    n = length(smooth_depth_score)
    boundaries = ["0" for _ = 1:n]
    cutoff_threshold = mean(smooth_depth_score) - std(smooth_depth_score) / 2.0
    depth_tuples = [(x, y) for (x, y) in zip(smooth_depth_score, [i for i = 1:n])]
    depth_tuples = reverse(sort(depth_tuples))

    if num_topics == Nothing
        for x in depth_tuples
            if x[1] > cutoff_threshold
                boundaries[x[2]] = "1"
            end
        end
    else
        for (i, x) in enumerate(depth_tuples)
            if x[1] > cutoff_threshold && i <= num_topics - 1
                boundaries[x[2]] = "1"
            end
        end
    end
    return join(boundaries[1:end-1])
end

"""
    TextTiling.segment(seg, document, [num_topics]) -> String
Performs the splitting of the document entered in the `document` argument.
# Arguments
- `seg`: Segment object.
- `document`: The document to be text segmented.
- `num_topics`: num_topics is the number of topics in the document. If this value is specified, segment boundaries are determined by the number of num_topics, starting with the highest depth score.

# Examples
```jldoctest
using TextSegmentation

window_size = 2
smooth_window_size = 1
num_topics = 3

tt = TextTiling.SegmentObject(window_size, smooth_window_size, Utls.tokenize)
result = TextTiling.segment(tt, document, num_topics)
println(result)
000100010000
```
"""
function segment(seg, document, num_topics=Nothing)
    preprocessed_document = preprocessing(seg, document)
    gap_score = calculate_gap_score(seg, preprocessed_document)
    depth_score = calculate_depth_score(seg, gap_score)
    smooth_depth_score = smoothing(seg, depth_score)
    return determine_boundaries(seg, smooth_depth_score, num_topics)
end
end


text:
using LinearAlgebra
using UnicodePlots
include("gemm.jl")

function plot_data(ns, data)
  scaling = 2 * ns.^3

  minv = floatmax()
  maxv = floatmin()
  for (s, v) = pairs(data)
    # if s != :ref
    minv = min(minv, minimum(scaling ./ v))
    maxv = max(maxv, maximum(scaling ./ v))
  # end
  end

  plt = lineplot(ns, scaling ./ data.ref, xlim = (ns[1], ns[end]),
                 ylim = (minv, maxv), title="max GFLOPS", xlabel="n",
                 ylabel="time", name = "ref")

  for (s, v) = pairs(data)
    if s != :ref
      lineplot!(plt, ns, scaling ./ v, name = string(s))
    end
  end
  plt
end

using BenchmarkTools

function run(ns, nsamples = 10)
  ref = zeros(length(ns))
  ji_reg = zeros(length(ns))
  ji = zeros(length(ns))

  for (i, n) = enumerate(ns)
    println()
    @show i,n
    m = k = n
    C = rand(m, n)
    A = rand(m, k)
    B = rand(k, n)
    C2 = copy(C)
    C1 = copy(C)

    t = @elapsed for n = 1:nsamples
      mul!(C2, A, B, 1, 1)
    end
    ref[i] = t / nsamples
    @show ref[i], 2*m^3 / ref[i]

    C .= C1
    t = @elapsed for n = 1:nsamples
      mygemm_NRxMR!(C, m, A, m, B, k, Val(8), Val(4))
    end
    ji_reg[i] = t / nsamples
    @show ji_reg[i], 2m^3/ji_reg[i]
    @show extrema(C - C2)

    C .= C1
    t = @elapsed for n = 1:nsamples
      mygemm!(C, A, B, Val(48), Val(256), Val(256), Val(8), Val(4))
    end
    ji[i] = t / nsamples
    @show ji[i], 2m^3/ji[i]
    @show extrema(C - C2)
  end

  data = (ref = ref,
          ji_reg = ji_reg,
          ji = ji,
         )

  return data
end

ns = 48:48:1000
data = run(ns)
plt = plot_data(ns, data)

max_stars_repo_name: jkozdon/Julia-LAFF-On-PfHP

max_issues_repo_name: jkozdon/Julia-LAFF-On-PfHP

max_forks_repo_name: jkozdon/Julia-LAFF-On-PfHP


text:
# ---
# title: 925. Long Pressed Name
# id: problem925
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Two Pointers, String
# link: <https://leetcode.com/problems/long-pressed-name/description/>
# hidden: true
# ---
#
# Your friend is typing his `name` into a keyboard.  Sometimes, when typing a
# character `c`, the key might get _long pressed_ , and the character will be
# typed 1 or more times.
#
# You examine the `typed` characters of the keyboard.  Return `True` if it is
# possible that it was your friends name, with some characters (possibly none)
# being long pressed.
#
#
#
# **Example 1:**
#
#
#
#     Input: name = "alex", typed = "aaleex"
#     Output: true
#     Explanation: 'a' and 'e' in 'alex' were long pressed.
#
#
# **Example 2:**
#
#
#
#     Input: name = "saeed", typed = "ssaaedd"
#     Output: false
#     Explanation: 'e' must have been pressed twice, but it wasn't in the typed output.
#
#
# **Example 3:**
#
#
#
#     Input: name = "leelee", typed = "lleeelee"
#     Output: true
#
#
# **Example 4:**
#
#
#
#     Input: name = "laiden", typed = "laiden"
#     Output: true
#     Explanation: It's not necessary to long press any character.
#
#
#
#
# **Constraints:**
#
#   * `1 <= name.length <= 1000`
#   * `1 <= typed.length <= 1000`
#   * The characters of `name` and `typed` are lowercase letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#Script to plot Rydberg radial wave functions
#23/07/2017

using Plots, JLD, LaTeXStrings
pyplot()
include("functions.jl")
PyPlot.close("all")

#Input information
atom = "87Rb"
nn = 50
ll = 0
jj = 0.5

#Calculate wave function
normY_sol, rr = numerovfunc(atom,nn,ll,jj)

#Rescale for plotting
plotscale = sqrt(rr)
probamp = (normY_sol.*plotscale).^2

alpha_c = getalpha(atom)

PyPlot.figure()
if nn>20
  plot(plotscale[plotscale .> sqrt(alpha_c^(1/3))],normY_sol[plotscale .> sqrt(alpha_c^(1/3))],linewidth=2)
else
  plot(plotscale,normY_sol,linewidth=2)
end
plot!(xlabel=L"(r/a_0)^{1/2}")
plot!(ylabel=L"r^{1/2}R(r) \,(a_0^{-1})")
plot!(title=string(atom," radial wavefunction |n,l,j⟩ = |",string(nn),",",string(ll),",",string(jj),"⟩"))
plot!(leg=false)
gui()

PyPlot.figure()
if nn>20
  plot(rr[rr .> alpha_c^(1/3)],probamp[rr .> alpha_c^(1/3)],linewidth=2)
else
  plot(rr,probamp,linewidth=2)
end
plot!(xlabel=L"r/a_0")
plot!(ylabel=L"|rR(r)|^2 \,(a_0^{-1})")
plot!(title=string(atom," radial probability density |n,l,j⟩ = |",string(nn),",",string(ll),",",string(jj),"⟩"))
plot!(leg=false)
gui()

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
##  SCRIPT NAME:   solveValueFunction.jl
##  DATE:          January 2015
#
#   This script solves the household problem by value function iteration.
#   The function obtains the value function for each time period and the policy
#   function (i.e. optimal next-period asset choice). From there I can work out
#   the optimal consumption level. The solution approach taken is by backwards
#   recursion. The optimization each period is carried out using NLopt.
#
#
##  ALEXANDROS THELOUDIS, UCL

##  -------------------------------------------------------------------------------------------
cd("C:\\Dropbox\\JuliaMeeting3")

#   Load packages:
using Grid
using NLopt

#   Load functions:
include("objectiveFunction.jl")
include("utilityH.jl")
include("utilityW.jl")
include("equalityConstraints1.jl")
include("mUtilH_cH.jl")
include("mUtilW_cW.jl")
include("mUtilH_k.jl")
include("mUtilW_k.jl")

##  ---------------------------------------------------------------------------------------------
#   Economy constants:
const r                  = 0.01            # Deterministic interest rate
const beta               = 0.98            # Discount factor
const cminH              = .1              # Min consumption husband (normalised /1000)
const cminW              = .1              # Min consumption wife (normalised /1000)
const kmin               = .1              # Min (inputs to) public consumption (normalised /1000)
const maxAnnWorkingHrsH  = 2.0             # 40 hrs/week times 50 weeks/year (normalised /1000)
const maxAnnWorkingHrsW  = 2.0             # 40 hrs/week times 50 weeks/year (normalised /1000)

#   Pareto weights related constants:
const paretoH            = 0.50            # Pareto weight husband
const paretoW            = 1 - paretoH     # Pareto weight wife

#   Preference-related constants (husband):
const gammaCH            = 1.5             # CRRA on private consumption
const gammaKH            = 1.5             # CRRA on public consumption
const alphaCH            = 0.5             # Weight of private consumption

#   Preference-related constants (wife):
const gammaCW            = 1.5             # CRRA on private consumption
const gammaKW            = 1.5             # CRRA on public consumption
const alphaCW            = 0.5             # Weight of private consumption

#   Solution constant:
const lengthA1           = 50              # Length of asset grid tomorrow

##  ---------------------------------------------------------------------------------------------
#   Current and future state:
A         = 10          # current assets
wH        = 5           # husband's wage today
wW        = 5           # wife's wage today
priceC    = 1.5         # price of private consumption

#   Calculate min and max assets tmr:
# ..lowest assets: cannot die with debts
lbA1      = 0
# ..highest assets: everything saved
ubA1      = (A + wH*maxAnnWorkingHrsH + wW*maxAnnWorkingHrsW - priceC*(cminH + cminW) -kmin)*(1+r)
# ..generate an asset grid for tmr assets (to be used in the EVal interpolation)
xcoord   = range(lbA1, (ubA1 - lbA1)/(lengthA1-1), lengthA1)

#   Expected Value function tomorrow (tmr is terminal period):
EVal1     = zeros(lengthA1)

#  ------------------------------------------------------------------------------------
#  Declare the relationship (linear) between next period's assets and Eval1:
VA!i     = CoordInterpGrid(xcoord, EVal1, BCnil, InterpLinear)

#  ------------------------------------------------------------------------------------
#  Define function to be optimized:

# Set properties of NL optimization problem:
valueFunctionOpt = Opt(:LD_SLSQP, 4)          # algorithm and dimensionality (:LN_COBYLA  :LD_SLSQP)
xtol_rel!(valueFunctionOpt,1e-6)   # tolerance

# Set bounds of optimization:
low_bound = [lbA1 ; cminH ; cminW ; kmin]     # lower bounds
up_bound  = [ubA1 ; Inf   ; Inf   ; Inf ]     # upper bounds
lower_bounds!(valueFunctionOpt,low_bound)     # lower bounds
upper_bounds!(valueFunctionOpt,up_bound)      # upper bounds

# Set objective function, min/max mode, and constraints:
max_objective!(valueFunctionOpt, (x,g) -> objectiveFunction(x, g, A, wH, wW, VA!i))
equality_constraint!(valueFunctionOpt, (x,g) -> equalityConstraints1(x, g, A, wH, wW, priceC), 1e-6)

# Set starting value and initialize optimization:
start_point = [ (ubA1 + lbA1)/2 ;
                (ubA1 - lbA1)/3 ;
                (ubA1 - lbA1)/3 ;
                (ubA1 - lbA1)/3 ]
startpt_valueFunctionOpt     = start_point             # starting point
(fval,x_optimal,return_flag) = optimize(valueFunctionOpt,[startpt_valueFunctionOpt])

# Store results of optimization:
policyA1 = maximum(x_optimal[1])    # policy for next period's assets
policyCH = maximum(x_optimal[2])    # policy for husband's consumption
policyCW = maximum(x_optimal[3])    # policy for wife's consumption
policyK  = maximum(x_optimal[4])    # policy for public consumption
Val      = fval                     # value function today

# Print results:
println("NLopt finished solving. It found an fval=",Val," with CH=",policyCH,", CW=",policyCW,", K=",policyK," and next period's assets are ",policyA1)

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
### Constraints
#
# Constraints are specified by the user as
#    lx_i ≤   x[i]  ≤ ux_i  # variable (box) constraints
#    lc_i ≤ c(x)[i] ≤ uc_i  # linear/nonlinear constraints
# and become equality constraints with l_i = u_i. ±∞ are allowed for l
# and u, in which case the relevant side(s) are unbounded.
#
# The user supplies functions to calculate c(x) and its derivatives.
#
# Of course we could unify the box-constraints into the
# linear/nonlinear constraints, but that would force the user to
# provide the variable-derivatives manually, which would be silly.
#
# This parametrization of the constraints gets "parsed" into a form
# that speeds and simplifies the IPNewton algorithm, at the cost of many
# additional variables. See `parse_constraints` for details.

struct ConstraintBounds{T}
    nc::Int            # Number of linear/nonlinear constraints supplied by user
    # Box-constraints on variables (i.e., directly on x)
    eqx::Vector{Int}   # index-vector of equality-constrained x (not actually variable...)
    valx::Vector{T}    # value of equality-constrained x
    ineqx::Vector{Int} # index-vector of other inequality-constrained variables
    σx::Vector{Int8}   # ±1, in constraints σ(v-b) ≥ 0 (sign depends on whether v>b or v<b)
    bx::Vector{T}      # bound (upper or lower) on variable
    # Linear/nonlinear constraint functions and bounds
    eqc::Vector{Int}   # index-vector equality-constrained entries in c
    valc::Vector{T}    # value of the equality-constraint
    ineqc::Vector{Int} # index-vector of inequality-constraints
    σc::Vector{Int8}   # same as σx, bx except for the nonlinear constraints
    bc::Vector{T}
end
function ConstraintBounds(lx, ux, lc, uc)
    _cb(symmetrize(lx, ux)..., symmetrize(lc, uc)...)
end
function _cb(lx::AbstractArray{Tx}, ux::AbstractArray{Tx}, lc::AbstractVector{Tc}, uc::AbstractVector{Tc}) where {Tx,Tc}
    T = promote_type(Tx, Tc)
    ConstraintBounds{T}(length(lc), parse_constraints(T, lx, ux)..., parse_constraints(T, lc, uc)...)
end

Base.eltype(::Type{ConstraintBounds{T}}) where T = T
Base.eltype(cb::ConstraintBounds) = eltype(typeof(cb))

Base.convert(::Type{ConstraintBounds{T}}, cb::ConstraintBounds{S}) where {T,S} =
    ConstraintBounds(cb.nc, cb.eqx, convert(Vector{T}, cb.valx),
                     cb.ineqx, cb.σx, convert(Vector{T}, cb.bx),
                     cb.eqc, convert(Vector{T}, cb.valc), cb.ineqc,
                     cb.σc, convert(Vector{T}, cb.bc))

"""
    nconstraints(bounds) -> nc

The number of linear/nonlinear constraint functions supplied by the
user. This does not include bounds-constraints on variables.

See also: nconstraints_x.
"""
nconstraints(cb::ConstraintBounds) = cb.nc

"""
    nconstraints_x(bounds) -> nx

The number of "meaningful" constraints (not `±Inf`) on the x coordinates.

See also: nconstraints.
"""
function nconstraints_x(cb::ConstraintBounds)
    mi = isempty(cb.ineqx) ? 0 : maximum(cb.ineqx)
    me = isempty(cb.eqx) ? 0 : maximum(cb.eqx)
    nmax = max(mi, me)
    hasconstraint = falses(nmax)
    hasconstraint[cb.ineqx] .= true
    hasconstraint[cb.eqx] .= true
    sum(hasconstraint)
end

function Base.show(io::IO, cb::ConstraintBounds)
    indent = "    "
    print(io, "ConstraintBounds:")
    print(io, "\n  Variables:")
    showeq(io, indent, cb.eqx, cb.valx, 'x', :bracket)
    showineq(io, indent, cb.ineqx, cb.σx, cb.bx, 'x', :bracket)
    print(io, "\n  Linear/nonlinear constraints:")
    showeq(io, indent, cb.eqc, cb.valc, 'c', :subscript)
    showineq(io, indent, cb.ineqc, cb.σc, cb.bc, 'c', :subscript)
    nothing
end

# Synonym constructor for ConstraintBounds with no c(x)
BoxConstraints(lx, ux) = ConstraintBounds(lx, ux, [], [])

abstract type AbstractConstraints end

nconstraints(constraints::AbstractConstraints) = nconstraints(constraints.bounds)

struct OnceDifferentiableConstraints{F,J,T} <: AbstractConstraints
    c!::F         # c!(storage, x) stores the value of the constraint-functions at x
    jacobian!::J  # jacobian!(storage, x) stores the Jacobian of the constraint-functions
    bounds::ConstraintBounds{T}
end

function OnceDifferentiableConstraints(c!, jacobian!,
                                       lx::AbstractArray, ux::AbstractArray,
                                       lc::AbstractArray, uc::AbstractArray)
    b = ConstraintBounds(lx, ux, lc, uc)
    OnceDifferentiableConstraints(c!, jacobian!, b)
end

function OnceDifferentiableConstraints(lx::AbstractArray, ux::AbstractArray)
    bounds = ConstraintBounds(lx, ux, [], [])
    OnceDifferentiableConstraints(bounds)
end

function OnceDifferentiableConstraints(bounds::ConstraintBounds)
    c! = (c, x)->nothing
    J! = (J, x)->nothing
    OnceDifferentiableConstraints(c!, J!, bounds)
end

function checked_chunk(lx)
    if isempty(lx)
        throw(ArgumentError("autodiff on constraints require the full lower bound vector `lx`."))
    end
    ForwardDiff.Chunk(lx)
end

function OnceDifferentiableConstraints(c!, lx::AbstractVector, ux::AbstractVector,
                                       lc::AbstractVector, uc::AbstractVector,
                                       autodiff::Symbol = :central,
                                       chunk::ForwardDiff.Chunk = checked_chunk(lx))

    bounds = ConstraintBounds(lx, ux, lc, uc)
    T = eltype(bounds)
    sizex = size(lx)
    sizec = size(lc)

    xcache = zeros(T, sizex)
    ccache = zeros(T, sizec)

    if is_finitediff(autodiff)
        ccache2 = similar(ccache)
        fdtype = finitediff_fdtype(autodiff)
        jacobian_cache = FiniteDiff.JacobianCache(xcache, ccache,ccache2,fdtype)
        function jfinite!(J, x)
            FiniteDiff.finite_difference_jacobian!(J, c!, x, jacobian_cache)
            J
        end
        return OnceDifferentiableConstraints(c!, jfinite!, bounds)
    elseif is_forwarddiff(autodiff)
        jac_cfg = ForwardDiff.JacobianConfig(c!, ccache, xcache, chunk)
        ForwardDiff.checktag(jac_cfg, c!, xcache)

        function jforward!(J, x)
            ForwardDiff.jacobian!(J, c!, ccache, x, jac_cfg, Val{false}())
            J
        end
        return OnceDifferentiableConstraints(c!, jforward!, bounds)
    else
        error("The autodiff value $autodiff is not support. Use :finite or :forward.")
    end
end

struct TwiceDifferentiableConstraints{F,J,H,T} <: AbstractConstraints
    c!::F # c!(storage, x) stores the value of the constraint-functions at x
    jacobian!::J # jacobian!(storage, x) stores the Jacobian of the constraint-functions
    h!::H   # h!(storage, x) stores the hessian of the constraint functions
    bounds::ConstraintBounds{T}
end

function TwiceDifferentiableConstraints(c!, jacobian!, h!, lx, ux, lc, uc)
    b = ConstraintBounds(lx, ux, lc, uc)
    TwiceDifferentiableConstraints(c!, jacobian!, h!, b)
end

function TwiceDifferentiableConstraints(c!, lx::AbstractVector, ux::AbstractVector,
    lc::AbstractVector, uc::AbstractVector,
    autodiff::Symbol = :central,
    chunk::ForwardDiff.Chunk = checked_chunk(lx))
       if is_finitediff(autodiff)
        fdtype = finitediff_fdtype(autodiff)
        return twicediff_constraints_finite(c!,lx,ux,lc,uc,fdtype,nothing)
    elseif is_forwarddiff(autodiff)
        return twicediff_constraints_forward(c!,lx,ux,lc,uc,chunk,nothing)
    else
        error("The autodiff value $autodiff is not support. Use :finite or :forward.")
    end
end

function TwiceDifferentiableConstraints(c!, con_jac!,lx::AbstractVector, ux::AbstractVector,
    lc::AbstractVector, uc::AbstractVector,
    autodiff::Symbol = :central,
    chunk::ForwardDiff.Chunk = checked_chunk(lx))
       if is_finitediff(autodiff)
        fdtype = finitediff_fdtype(autodiff)
        return twicediff_constraints_finite(c!,lx,ux,lc,uc,fdtype,con_jac!)
    elseif is_forwarddiff(autodiff)
        return  twicediff_constraints_forward(c!,lx,ux,lc,uc,chunk,con_jac!)
    else
        error("The autodiff value $autodiff is not support. Use :finite or :forward.")
    end
end

function TwiceDifferentiableConstraints(lx::AbstractArray, ux::AbstractArray)
    bounds = ConstraintBounds(lx, ux, [], [])
    TwiceDifferentiableConstraints(bounds)
end

function twicediff_constraints_forward(c!, lx, ux, lc, uc,chunk,con_jac! = nothing)
    bounds = ConstraintBounds(lx, ux, lc, uc)
    T = eltype(bounds)
    nc = length(lc)
    nx = length(lx)
    ccache = zeros(T, nc)
    xcache = zeros(T, nx)
    cache_check =  Ref{DataType}(Missing) #the datatype Missing, not the singleton
    ref_f= Ref{Any}() #cache for intermediate jacobian used in the hessian
    cxxcache = zeros(T, nx * nc, nx) #output cache for hessian
    h = reshape(cxxcache, (nc, nx, nx)) #reshaped output
    hi = [@view h[i,:,:] for i in 1:nc]
        #ref_f caches the closure function with its caches. other aproaches include using a Dict, but the
        #cost of switching happens just once per optimize call.

    if isnothing(con_jac!) #if the jacobian is not provided, generate one
        jac_cfg = ForwardDiff.JacobianConfig(c!, ccache, xcache, chunk)
        ForwardDiff.checktag(jac_cfg, c!, xcache)

        jac! = (J, x) -> begin
            ForwardDiff.jacobian!(J, c!, ccache, x, jac_cfg, Val{false}())
            J
        end

        con_jac_cached = x -> begin
            exists_cache = (cache_check[] == eltype(x))
            if exists_cache
                f = ref_f[]
                return f(x)
            else
                jcache = zeros(eltype(x), nc)
                out_cache = zeros(eltype(x), nc, nx)
                cfg_cache = ForwardDiff.JacobianConfig(c!,jcache,x)
                f = z->ForwardDiff.jacobian!(out_cache, c!, jcache, z,cfg_cache,Val{false}())
                ref_f[] = f
                cache_check[]= eltype(x)
                return f(x)
            end
        end

    else
        jac! = (J,x) -> con_jac!(J,x)

        #here, the cache should also include a JacobianConfig
         con_jac_cached = x -> begin
            exists_cache = (cache_check[] == eltype(x))
            if exists_cache
                f = ref_f[]
                return f(x)
            else
                out_cache = zeros(eltype(x), nc, nx)
                f = z->jac!(out_cache,x)
                ref_f[] = f
                cache_check[]= eltype(x)
                return f(x)
            end
        end
    end

    hess_config_cache = ForwardDiff.JacobianConfig(typeof(con_jac_cached),lx)
    function con_hess!(hess, x, λ)
        ForwardDiff.jacobian!(cxxcache, con_jac_cached, x,hess_config_cache,Val{false}())
        for i = 1:nc  #hot hessian loop
            hess+=λ[i].*hi[i]
        end
        return hess
    end

    return TwiceDifferentiableConstraints(c!, jac!, con_hess!, bounds)
end

function twicediff_constraints_finite(c!,lx,ux,lc,uc,fdtype,con_jac! = nothing)
    bounds = ConstraintBounds(lx, ux, lc, uc)
    T = eltype(bounds)
    nx = length(lx)
    nc = length(lc)
    xcache = zeros(T, nx)
    ccache = zeros(T, nc)

    if isnothing(con_jac!)
        jac_ccache = similar(ccache)
        jacobian_cache = FiniteDiff.JacobianCache(xcache, ccache,jac_ccache,fdtype)
        function jac!(J, x)
            FiniteDiff.finite_difference_jacobian!(J, c!, x, jacobian_cache)
            J
        end
    else
        jac! = (J,x) -> con_jac!(J,x)
    end
        cxxcache = zeros(T,nc*nx,nx) # to create cached jacobian
        h = reshape(cxxcache, (nc, nx, nx)) #reshaped output
        hi = [@view h[i,:,:] for i in 1:nc]

        function jac_vec!(J,x) #to evaluate the jacobian of a jacobian, FiniteDiff needs a vector version of that
            j_mat = reshape(J,nc,nx)
            return jac!(j_mat,x)
            return J
        end
        hess_xcache =similar(xcache)
        hess_cxcache =zeros(T,nc*nx) #output of jacobian, as a vector
        hess_cxxcache =similar(hess_cxcache)
        hess_config_cache = FiniteDiff.JacobianCache(hess_xcache,hess_cxcache,hess_cxxcache,fdtype)
        function con_hess!(hess, x, λ)
            FiniteDiff.finite_difference_jacobian!(cxxcache, jac_vec!, x,hess_config_cache)
            for i = 1:nc
                hi = @view h[i,:,:]
                hess+=λ[i].*hi
            end
            return hess
        end
        return TwiceDifferentiableConstraints(c!, jac!, con_hess!, bounds)
end

function TwiceDifferentiableConstraints(bounds::ConstraintBounds)
    c! = (x, c)->nothing
    J! = (x, J)->nothing
    h! = (x, λ, h)->nothing
    TwiceDifferentiableConstraints(c!, J!, h!, bounds)
end

## Utilities

function symmetrize(l, u)
    if isempty(l) && !isempty(u)
        l = fill!(similar(u), -Inf)
    elseif !isempty(l) && isempty(u)
        u = fill!(similar(l), Inf)
    end
    # TODO: change to indices?
    size(l) == size(u) || throw(DimensionMismatch("bounds arrays must be consistent, got sizes $(size(l)) and $(size(u))"))
    _symmetrize(l, u)
end
_symmetrize(l::AbstractArray{T,N}, u::AbstractArray{T,N}) where {T,N} = l, u
_symmetrize(l::Vector{Any}, u::Vector{Any}) = _symm(l, u)
_symmetrize(l, u) = _symm(l, u)

# Designed to ensure that bounds supplied as [] don't cause
# unnecessary broadening of the eltype. Note this isn't type-stable; if
# the user cares, it can be avoided by supplying the same concrete
# type for both l and u.
function _symm(l, u)
    if isempty(l) && isempty(u)
        if eltype(l) == Any
            # prevent promotion from returning eltype Any
            l = Array{Union{}}(undef, 0)
        end
        if eltype(u) == Any
            u = Array{Union{}}(undef, 0)
        end
    end
    promote(l, u)
end

"""
    parse_constraints(T, l, u) -> eq, val, ineq, σ, b

From user-supplied constraints of the form

    l_i ≤  v_i  ≤ u_i

(which include both inequality and equality constraints, the latter
when `l_i == u_i`), convert into the following representation:

    - `eq`, a vector of the indices for which `l[eq] == u[eq]`
    - `val = l[eq] = u[eq]`
    - `ineq`, `σ`, and `b` such that the inequality constraints can be written as
             σ[k]*(v[ineq[k]] - b[k]) ≥ 0
       where `σ[k] = ±1`.

Note that since the same `v_i` might have both lower and upper bounds,
`ineq` might have the same index twice (once with `σ`=-1 and once with `σ`=1).

Supplying `±Inf` for elements of `l` and/or `u` implies that `v_i` is
unbounded in the corresponding direction. In such cases there is no
corresponding entry in `ineq`/`σ`/`b`.

T is the element-type of the non-Int outputs
"""
function parse_constraints(::Type{T}, l, u) where T
    size(l) == size(u) || throw(DimensionMismatch("l and u must be the same size, got $(size(l)) and $(size(u))"))
    eq, ineq = Int[], Int[]
    val, b = T[], T[]
    σ = Array{Int8}(undef, 0)
    for i = 1:length(l)
        li, ui = l[i], u[i]
        li <= ui || throw(ArgumentError("l must be smaller than u, got $li, $ui"))
        if li == ui
            push!(eq, i)
            push!(val, ui)
        else
            if isfinite(li)
                push!(ineq, i)
                push!(σ, 1)
                push!(b, li)
            end
            if isfinite(ui)
                push!(ineq, i)
                push!(σ, -1)
                push!(b, ui)
            end
        end
    end
    eq, val, ineq, σ, b
end

### Compact printing of constraints

function showeq(io, indent, eq, val, chr, style)
    if !isempty(eq)
        print(io, '\n', indent)
        if style == :bracket
            eqstrs = map((i, v)->"$chr[$i]=$v", eq, val)
        else
            eqstrs = map((i, v)->"$(chr)_$i=$v", eq, val)
        end
        foreach(s->print(io, s * ", "), eqstrs[1:end - 1])
        print(io, eqstrs[end])
    end
end

function showineq(io, indent, ineqs, σs, bs, chr, style)
    if !isempty(ineqs)
        print(io, '\n', indent)
        if style == :bracket
            ineqstrs = map((i, σ, b)->"$chr[$i]" * ineqstr(σ, b), ineqs, σs, bs)
        else
            ineqstrs = map((i, σ, b)->"$(chr)_$i" * ineqstr(σ, b), ineqs, σs, bs)
        end
        foreach(s->print(io, s * ", "), ineqstrs[1:end - 1])
        print(io, ineqstrs[end])
    end
end
ineqstr(σ, b) = σ > 0 ? "≥$b" : "≤$b"


text:
module Func
using Random, Test, StatsBase, ProgressMeter, StatsPlots, Distributions

include("types.jl")
include("params.jl")
import .Types, .Params

# NOTE Used to have default parameters: drift=Params.drift, marketvol=Params.marketvol. Consider these for the functions in general, so that function arguments are minimal outside of testing.

# TODO: Add description to each function, according to some standard format
# that should probably inclucde input, output and purpose, plus comments on
# assumptions or non-obvious points.

function marketmove!(
    marketvals,
    t::Int,
    drift=Params.drift,
    marketvol=Params.marketvol)

    marketvals[t] = marketvals[t-1]*(1 + drift) + marketvol * randn()

    return marketvals
end

function marketinit!(
    marketvals,
    marketstartval=Params.marketstartval,
    perfwindow=Params.perfwindow,
    drift=Params.drift,
    marketvol=Params.marketvol)

    horizon = perfwindow[end]
    marketvals[1] = marketstartval
    for t in 2:horizon+1
        Func.marketmove!(marketvals, t, drift, marketvol)
    end

    return marketvals
end

function stockmove!(
    stockvals,
    t::Int,
    mktvals,
    betas,
    stockvolas)

    mktreturn = (mktvals[t] - mktvals[t-1]) / mktvals[t-1]
    stockvals[:, t] = ((1 .+ mktreturn .* betas) .* stockvals[:, t-1]) +
        stockvolas .* randn(length(stockvolas)) .* stockvals[:, t-1]

    return stockvals
end

function betainit!(
    betas,
    bigm=Params.bigm,
    betastd=Params.betastd,
    betamean=Params.betamean)

    betas .= randn(bigm) .* betastd .+ betamean

    return betas
end

# Function for a discrete volatility range
function stockvolinit!(
    stockvol,
    volrange=Params.stockvolrange,
    bigm=Params.bigm)

    stockvol .= rand(volrange, bigm)

    return stockvol
end

function stockvalueinit!(
    stocks,
    marketvals,
    stockstartval=Params.stockstartval,
    perfwindow=Params.perfwindow)

    horizon = perfwindow[end]
    stocks.value[:,1] .= stockstartval
    for t in 2:(horizon + 1)
        stocks.value .= Func.stockmove!(
        stocks.value, t, marketvals, stocks.beta, stocks.vol)
    end

    return stocks.value
end

function stockimpactinit!(
    impact,
    impactrange=Params.impactrange)

    impact .= rand(impactrange, length(impact))

    return impact
end

function invhorizoninit!(
    horizons,
    perfwindow=Params.perfwindow)
    horizons .= rand(perfwindow, length(horizons))

    return horizons
end

function invthreshinit!(
    thresholds,
    mean=Params.thresholdmean,
    std=Params.thresholdstd)

    thresholds .= randn(length(thresholds)) * std .+ mean

    return thresholds
end

# TODO: Sort out the capital range to match other items, as in rand(range, n)
function invassetinit!(
    invassets,
    caprange=Params.invcaprange,
    bigk=Params.bigk)

    ninvs = size(invassets, 1)
    capital = rand(caprange[1]:caprange[2], ninvs)
    for inv in 1:bigk
        invassets[inv, inv] = capital[inv]
    end
    for inv in (bigk+1):ninvs
        invassets[inv, rand(1:bigk)] = capital[inv]
    end

    return invassets
end

function fundcapitalinit!(
    fundvals,
    investorassets)

    fundvals[:,1] .= sum(investorassets, dims=1)[1:end-1]

    return fundvals
end

# TODO: Test of this function didn't throw error when ncols was misspecified
# as ncols=size(investorassets,1)- 1, i.e. wrong axis. Check why, write another
function fundstakeinit!(
    fundstakes,
    investorassets)

    ncols = size(investorassets, 2) - 1
    for col in 1:ncols
        fundstakes[col,:] = investorassets[:,col] ./ sum(investorassets[:,col])
    end

    return fundstakes
end

function fundholdinit!(
    holdings,
    initialcapital,
    stockvals,
    portfsizerange=Params.portfsizerange)

    kfunds = size(holdings, 1)
    mstocks = size(holdings, 2)
    nofm = rand(portfsizerange, kfunds) # Number of stocks in each portfolio
    for k in 1:kfunds # Loop over funds
        selection = rand(1:mstocks, nofm[k]) # Randomly select w/ replacement
        for stock in selection # Loop over stocks
            holdings[k, stock] +=
            (initialcapital[k] / length(selection)) / stockvals[stock]
        end
    end

    return holdings
end

function fundvalinit!(
    fundvals,
    holdings,
    stockvals,
    perfwindow=Params.perfwindow)

    horizon = perfwindow[end]
    kfunds = size(fundvals, 1)
    for t in 2:(horizon + 1) # t=1 is initialised by fundcapitalinit!
        for k in 1:kfunds
            fundvals[k, t] = sum(holdings[k, :] .* stockvals[:, t])
        end
    end

    return fundvals
end

function drawreviewers(
    bign=Params.bign,
    reviewprob=Params.reviewprobability)

    reviewers = rand(bign) .< reviewprob

    return findall(reviewers)
end

# QUESTION: Is it sufficient to compare fund values or do we need to compare the
# value of investors' stakes over time?

# TODO: Test that only a single investors.asset is > 0
function perfreview(
    t,
    reviewers,
    investors,
    fundvals)

    divestments = zeros(Int64, 1,2)
    for rev in reviewers
        horizon = investors.horizon[rev]
        fnds = findall(investors.assets[rev, 1:(end-1)] .> 0)
        for fnd in fnds
            valchange =
            (fundvals[fnd,t] - fundvals[fnd, t-horizon]) / fundvals[fnd,horizon]
            if valchange < investors.threshold[rev]
                reviewerfundpair = [convert(Int64, rev) convert(Int64, fnd)]
                if sum(divestments[1,:]) == 0
                    divestments = reviewerfundpair
                else
                    divestments = [divestments; reviewerfundpair]
                end
            end
        end
    end

    return divestments
end

# TODO: Write test that checks that sellorders and buyorders have the same
# number of rows in .values and .investors/.funds
# QUESTION: Should I implement size and value checks for all intermediate outpts
# NOTE: Passes one column of stock values only

function liquidate!(
    holdings,
    stakes,
    stockvals,
    divestments)

    sellorders = Types.SellMarketOrder(
    Array{Float64}(undef, 0, size(holdings, 2)), Array{Float64}(undef, 0))
    for row in 1:size(divestments, 1)
        investor = divestments[row, 1]
        fund = divestments[row, 2]

        # Note the minus
        salevals = -stakes[fund, investor] .* (holdings[fund, :] .* stockvals)'
        sellorders.values = vcat(sellorders.values, salevals)
        sellorders.investors = vcat(sellorders.investors, investor)

        remainder = (1 - stakes[fund, investor])
        holdings[fund, :] = holdings[fund, :] .* remainder

        stakes[fund, investor] = 0
        if sum(stakes[fund,: ]) > 0
            stakes[fund, :] = stakes[fund, :] ./ sum(stakes[fund,: ])
        end
    end

    return sellorders, holdings, stakes
end

# NOTE: Truncate minimum stock value at 0.000000001
# NOTE: Could move volume calculation into marketmake
function marketmake!(
    stocks,
    t,
    ordervals)

    netimpact = sum(ordervals, dims=1)' .* stocks.impact
    stocks.value[:, t] .= vec((1 .+ netimpact) .* stocks.value[:, t])

    nearzeroidx = findall(stocks.value[:, t] .< 0.000000001)
    stocks.value[nearzeroidx, t] .= 0.000000001

    return stocks.value
end

function executeorder!(
    stocks,
    t,
    orders::Types.SellMarketOrder)

    cashout = Array{Float64}(undef, 0, 2)
    oldstockvals = copy(stocks.value)
    stocks.value .= marketmake!(stocks, t, orders.values)
    priceratio = stocks.value[:, t] ./ oldstockvals[:, t]
    for order in 1:size(orders.values, 1)
        investor = orders.investors[order]
        amount = -sum(orders.values[order, :] .* priceratio)
        cashout = vcat(cashout, [investor amount])
    end

    return stocks.value, cashout
end

# NOTE: Small issue here that sell orders happen first, so sellers get a bad
# price, followed by buyers who get a good price, ~midmarket as their order
# neutralises the prior decline. However, the amount they buy is also diminished
function executeorder!(
    stocks,
    t,
    orders::Types.BuyMarketOrder)

    mstocks = size(stocks.value, 1)
    sharesout = Array{Float64}(undef, 0, mstocks + 1)
    stocks.value .= marketmake!(stocks, t, orders.values)
    for order in 1:size(orders.values, 1)
        fund = orders.funds[order]
        shares = (orders.values[order, :] ./ stocks.value[:, t])'
        sharesout = vcat(sharesout, [fund shares])
    end

    return stocks.value, sharesout
end

# QUESTION: Structure of disburse! methods is not symmetric in args... fix?
# Disburse cash to investors following  sell order / divestment
function disburse!(
    invassets,
    divestments,
    cashout)

    for row in 1:size(divestments,1)
        inv = divestments[row, 1]
        fund = divestments[row, 2]
        #@assert(inv == cashout[row, 1]) # Ensure divestment matches cashout
        invassets[inv, fund] = 0
        invassets[inv, end] = cashout[row, 2]
    end

    return invassets
end

# Disburse shares to funds following buy order/investment, update value history
# QUESTION: Why doesn't the value of funds drop as investors divest?
function disburse!(
    funds::Types.EquityFund,
    sharesout,
    stockvals)

    for row in 1:size(sharesout, 1)
        fund = convert(Int64, sharesout[row, 1])
        funds.holdings[fund, :] += sharesout[row, 2:end]
        fundrevalue!(funds, stockvals, [fund])
    end

    return funds
end

function fundrevalue!(
    funds,
    stockvals,
    targets=1:Params.bigk)

    for k in targets
        funds.value[k, :] .= vec(funds.holdings[k, :]' * stockvals)
    end

    return funds.value
end

function respawn!(
    funds,
    investors,
    t,
    stockvals,
    portfsizerange=Params.portfsizerange)

    buyorders = Types.BuyMarketOrder(
    Array{Float64}(undef, 0, size(funds.holdings, 2)), Array{Float64}(undef, 0))
    spawn = findall(vec(sum(funds.holdings, dims=2) .== 0))
    neggs = length(spawn)
    potentialinvs = findall(vec(investors.assets[:, end] .> 0))
    anchorinvs = Random.shuffle(potentialinvs)[1:neggs]
    for idx in 1:neggs
        inv = anchorinvs[idx]
        egg = spawn[idx]
        capital = investors.assets[inv, end]
        # Anchor investor notes investment in the fund
        investors.assets[inv, egg] = capital
        investors.assets[inv, end] = 0
        # Re-use fund holdings init fctn to draw new stocks and generate order
        buyorder = (fundholdinit!(funds.holdings[egg, :]',
        capital, stockvals[:, t], portfsizerange)' .* stockvals[:, t])'
        buyorders.values = vcat(buyorders.values, buyorder)
        buyorders.funds = vcat(buyorders.funds, egg)
        # Set anchor investor's stake
        funds.stakes[egg, inv] = 1
    end

    return buyorders, investors.assets, funds.stakes
end

# NOTE: Bit of an issue: as many reinvestments as dead funds are randomly drawn,
# which reduces the flow-performance relationship, especially if many funds die

function reinvest!(
    investors,
    funds,
    stockvals,
    t,
    selection="best")

    kfunds = size(investors.assets, 2) - 1

    buyorders = Types.BuyMarketOrder(
    Array{Float64}(undef, 0, size(funds.holdings, 2)), Array{Float64}(undef, 0))

    reinvestors = findall(vec(investors.assets[:, end] .> 0))
    for reinv in reinvestors
        injection = investors.assets[reinv, end]
        if selection == "best"
            # Choose best-performing fund
            choice = bestperformer(funds.value, investors.horizon[reinv], t)
        elseif selection == "goodenough"
            # Choose any fund above return threshold
            choice = goodenoughfund(
            funds.value, investors.horizon[reinv], investors.threshold[reinv],t)
        else
            println("Please specify a valid fund selection method in reinvest!")
        end
        # Move cash into the fund
        investors.assets[reinv, choice] = injection
        investors.assets[reinv, end] = 0

        fundval = funds.holdings[choice, :]' * stockvals[:,t]
        # QUESTION: Should assert hold true, i.e. funds.value be unchanged?
        #@assert fundval == funds.value[choice, t]
        stakevals = funds.stakes[choice, :] .* fundval
        stakevals[reinv] = injection
        newcapital = fundval + injection
        funds.stakes[choice, :] .= stakevals ./ newcapital

        portfoliovalues = funds.holdings[choice, :] .* stockvals[:, t]
        buyorder = (portfoliovalues ./ fundval)' .* injection
        buyorders.values = vcat(buyorders.values, buyorder)
        buyorders.funds = vcat(buyorders.funds, choice)
    end

    return investors.assets, funds.stakes, buyorders
end

function bestperformer(
    fundvals,
    horizon,
    t)

    # Return between investor's horizon and now/t
    horizonreturns =
    (fundvals[:, t] - fundvals[:, t-horizon]) ./ fundvals[:, t-horizon]
    # Index of best-performing fund
    _, bestfund = findmax(horizonreturns)

    return bestfund
end

# TODO: Test for goodenoughfund function
# NOTE: Choose best performer if no funds meet performance threshold
function goodenoughfund(
    fundvals,
    horizon,
    threshold,
    t)

    # Return between investor's horizon and now/t
    horizonreturns =
    (fundvals[:, t] - fundvals[:, t-horizon]) ./ fundvals[:, t-horizon]
    # Index of best-performing fund
    candidates = findall(horizonreturns .> threshold)
    if size(candidates,1) == 0
        _, goodenough = findmax(horizonreturns)
    else
        goodenough = rand(1:size(fundvals,1))
    end

    return goodenough
end

function initialise(market, stocks, investors, funds)
    Func.marketinit!(market.value)
    Func.betainit!(stocks.beta)
    Func.stockvolinit!(stocks.vol)
    Func.stockvalueinit!(stocks, market.value)
    Func.stockimpactinit!(stocks.impact)
    Func.invhorizoninit!(investors.horizon)
    Func.invthreshinit!(investors.threshold)
    Func.invassetinit!(investors.assets)
    Func.fundcapitalinit!(funds.value, investors.assets)
    Func.fundstakeinit!(funds.stakes, investors.assets)
    Func.fundholdinit!(funds.holdings, funds.value[:, 1], stocks.value)
    Func.fundvalinit!(funds.value, funds.holdings, stocks.value)
end

function boundstest(market, stocks, investors, funds)
    @testset "Variable Bounds" begin

        @test all(stocks.value .>= 0)
        @test all(stocks.vol .>= 0)
        @test all(stocks.impact .>= 0)

        @test all(funds.holdings .>= 0)
        @test all(funds.stakes .>= 0)
        @test all(funds.value .>= 0)

        @test all(investors.assets .>= 0)
        @test all(investors.horizon .>= 0)
        @test all(1 .>= investors.threshold .>= -1)
    end
end

function modelrun(market, stocks, investors, funds)
    @showprogress for t in (Params.perfwindow[end]+2):Params.bigt

        # Market index and assets move
        Func.marketmove!(market.value, t)
        # Func.boundstest(market, stocks, investors, funds)
        Func.stockmove!(stocks.value, t, market.value, stocks.beta, stocks.vol)
        # Func.boundstest(market, stocks, investors, funds)
        Func.fundrevalue!(funds, stocks.value)
        # Func.boundstest(market, stocks, investors, funds)
        # print("Time: ", t)
        # QUESTION: Consider whether I should pass t instead of stocks.value[:, t]
        # Investors (may_ sell out of underperforming funds
        # TODO: Test sellorders.investors == cashout[1,:], similar for buyorders
        # TODO: No cashout of value 0
        # QUESTION: Enforce integers for fund/investor IDs in sellorders, cashout..?
        # QUESTION: Is a separate trackvolume! function appropriate?

        # Investors' performance review
        reviewers = Func.drawreviewers()
        divestments = Func.perfreview(t, reviewers, investors, funds.value)

        # Run liquidation/reinvestment cycle if there are divestments
        if sum(divestments[1,:]) > 0

            sellorders, _, _ = Func.liquidate!(funds.holdings, funds.stakes,
            stocks.value[:, t], divestments)

            Func.trackvolume!(stocks.volume, sellorders, t)

            _, cashout = Func.executeorder!(stocks, t, sellorders)
            Func.disburse!(investors.assets, divestments, cashout)

            # TODO: Write test that check correct outocme of the whole process,
            # step by step

            # Generate a new fund for each dead one, if there are dead ones
            if any(sum(funds.stakes, dims=2) .== 0)

                buyorders, _, _ = Func.respawn!(funds, investors, t, stocks.value)

                Func.trackvolume!(stocks.volume, sellorders, t)

                _, sharesout = Func.executeorder!(stocks, t, buyorders)

                Func.disburse!(funds, sharesout, stocks.value)

            end # respawn loop

            # Re-value funds prior to reinvestment decisions
            Func.fundrevalue!(funds, stocks.value)

            # If there are still investors with spare cash, they reinvest it
            if any(investors.assets[:, end] .> 0)

                _, _, buyorders = Func.reinvest!(
                investors, funds, stocks.value, t, "goodenough")

                Func.trackvolume!(stocks.volume, sellorders, t)

                _, sharesout = Func.executeorder!(stocks, t, buyorders)

                Func.disburse!(funds, sharesout, stocks.value)
            end # reinvestment loop

        end # Divestment loop
    end # Model run loop
end # runmodel function

# TODO: Write tests for both assetreturns functions and vecreturns
# NOTE: Duplication in assetreturns and demean, i.e. targets for refactoring
function vecreturns(vector::Array{Float64,1})
    returnsvec = (vector[2:end] - vector[1:(end-1)]) ./ vector[1:(end-1)]
    return returnsvec
end

function assetreturns(vector::Array{Float64,1})
    returnsvec = vecreturns(vector)
    return returnsvec
end

function assetreturns(array::Array{Float64,2})
    nassets = size(array, 1)
    nperiods = size(array, 2)
    returnsarray = zeros(nassets, nperiods-1)
    for asset in 1:nassets
        returnsarray[asset,:] .= vecreturns(array[asset,:])
    end
    return returnsarray
end

function subtractmean(returns::Array{Float64,1})
    demeaned = returns .- mean(returns)
    return demeaned
end

# TODO: Test demean functions
function demean(returns::Array{Float64,1})
    demeanedreturns = subtractmean(returns)
    return demeanedreturns
end

function demean(returns::Array{Float64,2})
    nassets = size(returns, 1)
    nperiods = size(returns, 2)
    demeanedreturns = zeros(nassets, nperiods)
    for asset in 1:nassets
        demeanedreturns[asset, :] .= subtractmean(returns[asset, :])
    end
    return demeanedreturns
end
function calckurtoses(returns, startidx=Params.perfwindow[end]+1)
    nassets = size(returns, 1)
    kurtoses = zeros(nassets)
    for asset in 1:nassets
        kurtoses[asset] = kurtosis(returns[asset, startidx:end])
    end
    return kurtoses
end

# TODO: Test for trackvolume!
# NOTE: Stylised facts should only be calculated for period after initialisation
function trackvolume!(volume, order::Types.BuyMarketOrder, t)
    volume[:, t] += sum(order.values, dims=1)'
    return volume
end

function trackvolume!(volume, order::Types.SellMarketOrder, t)
    volume[:, t] -= sum(order.values, dims=1)'
    return volume
end

# TODO: Test for volavolumecorr
function volavolumecorr(volumes, returns, startidx=Params.perfwindow[end]+1)
    nassets = size(volumes, 1)
    correlation = zeros(nassets)
    for asset in 1:nassets
        correlation[asset] =
        StatsBase.cor(volumes[asset, (startidx+1):end],
        broadcast(abs, returns[asset, startidx:end]))
    end
    return correlation
end

# TODO: Test for lossgainratio
function lossgainratio(returns, cutoff, startidx=Params.perfwindow[end]+1)
    nassets = size(returns, 1)
    ratios = zeros(nassets)
    for asset in 1:nassets
        stock = returns[asset, startidx:end]
        absolutevals = broadcast(abs, stock)
        largeamount = percentile(absolutevals, cutoff)
        upmoves = sum(stock .> largeamount)
        downmoves = sum(stock .< -largeamount)
        ratios[asset] = downmoves/upmoves
    end
    return ratios
end

# TODO: move plots to separate file
# TODO: save plots in plots file

function plot_stylisedfacts(
    marketval, stocksval, stocksvolume, returnspctile=5)
    marketreturns = Func.assetreturns(marketval)
    stockreturns = Func.assetreturns(stocksval)
    demeanedmarketreturns = Func.demean(marketreturns)
    demeanedstockreturns = Func.demean(stockreturns)
    demeanedmarketreturnssquared = demeanedmarketreturns.^2
    demeanedstockreturnssquared = demeanedstockreturns.^2
    lags = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]#collect[1:10]
    choice = rand(1:size(stockreturns,1))
    pyplot() # as GR wasn't working, maybe fix that
    plot_pricehistories(stocksval, marketval)
    plot_marketreturnhistogram(demeanedmarketreturns)
    plot_stockreturnhistogram(demeanedstockreturns, choice)
    plot_marketpacf(demeanedmarketreturns, lags)
    plot_stockpacf(demeanedstockreturns, choice, lags)
    plot_marketvolaclustering(demeanedmarketreturnssquared, lags)
    plot_stockvolaclustering(demeanedstockreturnssquared, choice, lags)
    plot_stockkurtoses(demeanedstockreturns)
    plot_lossgainratio(demeanedstockreturns, returnspctile)
    plot_volavolumecorr(stocksvolume, demeanedstockreturns)
end

# TODO: refactor plots - DRY
function plot_pricehistories(stocksval, marketval)
    periods = 1:length(marketval)
    nstocks = 1:100
    StatsPlots.plot(periods, transpose(stocksval)[:, nstocks], legend=:none,
    title="History of stock returns")
    plt = StatsPlots.plot!(
        periods, marketval, lw=3, lc=:black, label="Market index")
    png(joinpath(Params.plotpath, "returnshistory"))
    display(plt)
end

function plot_marketreturnhistogram(asset)
    StatsPlots.histogram(asset, bins=100, title="Market returns vs normal distribution", label="Market returns", legend=:topright)
    x = findmin(asset)[1]:0.001:findmax(asset)[1]
    d = fit(Normal, asset)
    plt = StatsPlots.plot!(x, pdf.(d, x), lc=:red, label="Normal distn")
    png(joinpath(Params.plotpath, "marketreturnsvsnormal"))
    display(plt)
end

function plot_stockreturnhistogram(demeanedstockreturns, choice)
    asset = demeanedstockreturns[choice, :]
    StatsPlots.histogram(asset, bins=100, title="Example of stock returns", label="Random stock's returns", legend=:topright)
    x = findmin(asset)[1]:0.001:findmax(asset)[1]
    d = fit(Normal, asset)
    plt = StatsPlots.plot!(x, pdf.(d, x), lc=:red, label="Normal distn")
    png(joinpath(Params.plotpath, "stockreturnsvsnormal"))
    display(plt)
end

function plot_marketpacf(demeanedmarketreturns, lags)
    asset = demeanedmarketreturns
    pacfcoeffs = pacf(demeanedmarketreturns, lags)
    StatsPlots.bar(pacfcoeffs, title="Market autocorrelation", label="Partial correlation coefficients")
    plt = StatsPlots.plot!(0:0.01:10, [ones(length(0:0.01:10)) .* 1.96 / sqrt(size(asset, 1)), -ones(length(0:0.01:10)) .* 1.96 / sqrt(size(asset, 1)),], lc=:red, label="critical values", legend=:bottomright)
    png(joinpath(Params.plotpath, "marketautocorrelation"))
    display(plt)
end

function plot_stockpacf(demeanedstockreturns, choice, lags)
    asset = demeanedstockreturns[choice,: ]
    pacfcoeffs = pacf(asset, lags)
    StatsPlots.bar(pacfcoeffs, title="Stock autocorrelation (random stock example)", label="Partial autocorrelation coefficients, demeaned returns", legend=:bottomright)
    plt = StatsPlots.plot!(0:0.01:10, [ones(length(0:0.01:10)) .* 1.96 / sqrt(size(asset, 1)), -ones(length(0:0.01:10)) .* 1.96 / sqrt(size(asset, 1)),], lc=:red, label="Critical values")
    png(joinpath(Params.plotpath, "stockautocorrelation"))
    display(plt)
end

function plot_marketvolaclustering(demeanedmarketreturnssquared, lags)
    asset = demeanedmarketreturnssquared
    pacfcoeffs = pacf(asset, lags)
    StatsPlots.bar(pacfcoeffs, title="Volatility clustering, market", label="Partial autocorrelation coefficients, squared demeaned returns", legend=:bottomright)
    plt = StatsPlots.plot!(0:0.01:10, [ones(length(0:0.01:10)) .* 1.96 / sqrt(size(asset, 1)), -ones(length(0:0.01:10)) .* 1.96 / sqrt(size(asset, 1)),], lc=:red, label="Critical values")
    png(joinpath(Params.plotpath, "marketvolaclustering"))
    display(plt)
end

function plot_stockvolaclustering(demeanedstockreturnssquared, choice, lags)
    asset = demeanedstockreturnssquared[choice,:]
    pacfcoeffs = pacf(asset, lags)
    StatsPlots.bar(pacfcoeffs, title="Volatility clustering, random stock", label="partial autocorrelation coefficients, squared demeaned returns", legend=:bottomright)
    plt = StatsPlots.plot!(0:0.01:10, [ones(length(0:0.01:10)) .* 1.96 / sqrt(size(asset, 1)), -ones(length(0:0.01:10)) .* 1.96 / sqrt(size(asset, 1)),], lc=:red, label="Critical values")
    png(joinpath(Params.plotpath, "stockvolaclustering"))
    display(plt)
end

function plot_stockkurtoses(demeanedstockreturns)
    stockkurtoses = Func.calckurtoses(demeanedstockreturns)
    plt = StatsPlots.bar(sort!(stockkurtoses, rev=true), title="Kurtoses of stock returns", label="Kurtoses")
    png(joinpath(Params.plotpath, "kurtoses"))
    display(plt)
end

function plot_lossgainratio(demeanedstockreturns, returnspercentile)
    ratios = Func.lossgainratio(demeanedstockreturns, returnspercentile)
    plot(ratios, title="Gain-loss asymmetry in stock returns", label="# large losses / # large gains")
    plt = plot!(0:0.01:200, ones(length(0:0.01:200)) * mean(ratios), label="Mean ratio")
    png(joinpath(Params.plotpath, "lossgainratio"))
    display(plt)
end

function plot_volavolumecorr(volume, demeanedstockreturns)
    corrs = Func.volavolumecorr(volume, demeanedstockreturns)
    plt = bar(sort(corrs, rev=true), title="Volume-volatility correlation of stocks", label="Volume - absolute return correlation coefficients")
    png(joinpath(Params.plotpath, "volavolumecorr"))
    display(plt)
end

end  # module Functions


text:
using AxisArrays
using JuMP
using MathOptInterface
using Cbc


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# Illustration of the corner (circled), edge (braces), and face index numbering convention
# used in these functions.
#
#                    ⑧────────────────────────{7}────────────────────────⑦
#                   ╱│                                                   ╱│
#                  ╱ │                                                  ╱ │
#                 ╱  │                                                 ╱  │
#                ╱   │                5 (+z)                          ╱   │
#               ╱    │                                               ╱    │
#              ╱     │                                              ╱     │
#            {12}    │                                            {11}    │
#            ╱       │                                            ╱       │
#           ╱        │                                           ╱        │
#          ╱         │                    2 (+y)                ╱         │
#         ╱          │                                         ╱          │
#        ╱          {8}                                       ╱          {6}
#       ╱            │                                       ╱            │
#      ⑤─────────────────────────{3}───────────────────────⑥    4 (+x)   │
#      │             │                                      │             │
#      │             │                                      │             │
#      │             │                                      │             │
#      │    6 (-x)   │                                      │             │
#      │             │                                      │             │
#      │             │                                      │             │
#      │             │                                      │             │
#      │             │                                      │             │
#      │             ④────────────────────────{5}──────────│─────────────③
#      │            ╱                                       │            ╱
#      │           ╱              1 (-y)                    │           ╱
#     {4}         ╱                                        {2}         ╱
#      │         ╱                                          │         ╱
#      │        ╱                                           │        ╱
#      │      {9}                                           │      {10}
#      │      ╱                                             │      ╱
#      │     ╱                                              │     ╱  Global coordinates:
#      │    ╱                                               │    ╱       z
#      │   ╱                             3 (-z)             │   ╱        ↑   y
#      │  ╱                                                 │  ╱         │  ╱
#      │ ╱                                                  │ ╱          │ ╱
#      │╱                                                   │╱           └─────> x
#      ①───────────────────────{1}─────────────────────────②

# Transfinite mapping formula from a point (xi, eta, zeta) in reference space [-1,1]^3 to a
# physical coordinate (x, y, z) for a hexahedral element with straight sides
function straight_side_hex_map(xi, eta, zeta, corner_points)

  coordinate = zeros(eltype(xi), 3)
  for j in 1:3
    coordinate[j] += (0.125 * ( corner_points[j, 1] * (1 - xi) * (1 - eta) * (1 - zeta)
                              + corner_points[j, 2] * (1 + xi) * (1 - eta) * (1 - zeta)
                              + corner_points[j, 3] * (1 + xi) * (1 + eta) * (1 - zeta)
                              + corner_points[j, 4] * (1 - xi) * (1 + eta) * (1 - zeta)
                              + corner_points[j, 5] * (1 - xi) * (1 - eta) * (1 + zeta)
                              + corner_points[j, 6] * (1 + xi) * (1 - eta) * (1 + zeta)
                              + corner_points[j, 7] * (1 + xi) * (1 + eta) * (1 + zeta)
                              + corner_points[j, 8] * (1 - xi) * (1 + eta) * (1 + zeta) ) )
  end

  return coordinate
end

# Construct the (x, y, z) node coordinates in the volume of a straight sided hexahedral element
function calc_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, element, nodes, corners)

  for k in eachindex(nodes), j in eachindex(nodes), i in eachindex(nodes)
    node_coordinates[:, i, j, k, element] .= straight_side_hex_map(nodes[i], nodes[j], nodes[k], corners)
  end

  return node_coordinates
end

# Transfinite mapping formula from a point (xi, eta, zeta) in reference space [-1,1]^3 to a point
# (x,y,z) in physical coordinate space for a hexahedral element with general curved sides
# See Section 4.3
# - Andrew R. Winters (2014)
#   Discontinous Galerkin spectral element approximations for the reflection and
#   transmission of waves from moving material interfaces
#   [PhD thesis, Florida State University](https://diginole.lib.fsu.edu/islandora/object/fsu%3A185342)
function transfinite_hex_map(xi, eta, zeta, face_curves::AbstractVector{<:CurvedFace})

  coordinate = zeros(eltype(xi), 3)
  face_values = zeros(eltype(xi), (3, 6))
  edge_values = zeros(eltype(xi), (3, 12))
  corners = zeros(eltype(xi), (3, 8))

  # Compute values along the face edges
  edge_values[:, 1] .= evaluate_at(SVector(xi,   -1), face_curves[1])
  edge_values[:, 2] .= evaluate_at(SVector( 1, zeta), face_curves[1])
  edge_values[:, 3] .= evaluate_at(SVector(xi,    1), face_curves[1])
  edge_values[:, 4] .= evaluate_at(SVector(-1, zeta), face_curves[1])

  edge_values[:, 5] .= evaluate_at(SVector(xi,   -1), face_curves[2])
  edge_values[:, 6] .= evaluate_at(SVector( 1, zeta), face_curves[2])
  edge_values[:, 7] .= evaluate_at(SVector(xi,    1), face_curves[2])
  edge_values[:, 8] .= evaluate_at(SVector(-1, zeta), face_curves[2])

  edge_values[:, 9]  .= evaluate_at(SVector(eta, -1), face_curves[6])
  edge_values[:, 10] .= evaluate_at(SVector(eta, -1), face_curves[4])
  edge_values[:, 11] .= evaluate_at(SVector(eta,  1), face_curves[4])
  edge_values[:, 12] .= evaluate_at(SVector(eta,  1), face_curves[6])

  # Compute values on the face
  face_values[:, 1] .= evaluate_at(SVector( xi, zeta), face_curves[1])
  face_values[:, 2] .= evaluate_at(SVector( xi, zeta), face_curves[2])
  face_values[:, 3] .= evaluate_at(SVector( xi,  eta), face_curves[3])
  face_values[:, 4] .= evaluate_at(SVector(eta, zeta), face_curves[4])
  face_values[:, 5] .= evaluate_at(SVector( xi,  eta), face_curves[5])
  face_values[:, 6] .= evaluate_at(SVector(eta, zeta), face_curves[6])

  # Pull the eight corner values and compute the straight sided hex mapping
  corners[:,1] .= face_curves[1].coordinates[:, 1,   1]
  corners[:,2] .= face_curves[1].coordinates[:, end, 1]
  corners[:,3] .= face_curves[2].coordinates[:, end, 1]
  corners[:,4] .= face_curves[2].coordinates[:, 1,   1]
  corners[:,5] .= face_curves[1].coordinates[:, 1,   end]
  corners[:,6] .= face_curves[1].coordinates[:, end, end]
  corners[:,7] .= face_curves[2].coordinates[:, end, end]
  corners[:,8] .= face_curves[2].coordinates[:, 1,   end]

  coordinate_straight = straight_side_hex_map(xi, eta, zeta, corners)

  # Compute the transfinite mapping
  for j in 1:3
    # Linear interpolation between opposite faces
    coordinate[j] = ( 0.5 * ( face_values[j, 6] * (1 - xi  ) + face_values[j, 4] * (1 + xi  )
                            + face_values[j, 1] * (1 - eta ) + face_values[j, 2] * (1 + eta )
                            + face_values[j, 3] * (1 - zeta) + face_values[j, 5] * (1 + zeta) ) )

    # Edge corrections to ensure faces match
    coordinate[j] -= ( 0.25 * ( edge_values[j, 1 ] * (1 - eta) * (1 - zeta)
                              + edge_values[j, 2 ] * (1 + xi ) * (1 - eta )
                              + edge_values[j, 3 ] * (1 - eta) * (1 + zeta)
                              + edge_values[j, 4 ] * (1 - xi ) * (1 - eta )
                              + edge_values[j, 5 ] * (1 + eta) * (1 - zeta)
                              + edge_values[j, 6 ] * (1 + xi ) * (1 + eta )
                              + edge_values[j, 7 ] * (1 + eta) * (1 + zeta)
                              + edge_values[j, 8 ] * (1 - xi ) * (1 + eta )
                              + edge_values[j, 9 ] * (1 - xi ) * (1 - zeta)
                              + edge_values[j, 10] * (1 + xi ) * (1 - zeta)
                              + edge_values[j, 11] * (1 + xi ) * (1 + zeta)
                              + edge_values[j, 12] * (1 - xi ) * (1 + zeta) ) )

    # Subtracted interior twice, so add back the straight-sided hexahedral mapping
    coordinate[j] += coordinate_straight[j]
  end

  return coordinate
end

# Construct the (x, y, z) node coordinates in the volume of a curved sided hexahedral element
function calc_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, element, nodes,
                                face_curves::AbstractVector{<:CurvedFace})

  for k in eachindex(nodes), j in eachindex(nodes), i in eachindex(nodes)
    node_coordinates[:, i, j, k, element] .= transfinite_hex_map(nodes[i], nodes[j], nodes[k], face_curves)
  end

  return node_coordinates
end

end # @muladd


text:
"""
    MD_TOKENS

Dictionary of tokens for Markdown. Note that for each, there may be several
possibilities to consider in which case the order is important: the first case
that works will be taken.

Dev: F_* are greedy match, see `md_utils.jl`.

Try: https://spec.commonmark.org/dingus
"""
const MD_TOKENS = Dict{Char, Vector{Pair{TokenFinder, Symbol}}}(
    '\n' => [
        F_LINE_RETURN       => :LINE_RETURN,
        forward_match("\n") => :LINE_RETURN,
        ],
    '?' => [
        forward_match("???") => :RAW
    ],
    '~' => [
        forward_match("~~~") => :RAW_HTML
        ],
    '%' => [
        forward_match("%%%") => :RAW_LATEX
        ],
    '(' => [
        forward_match("(") => :BRACKET_OPEN
        ],
    ')' => [
        forward_match(")") => :BRACKET_CLOSE
        ],
    '{' => [
        forward_match("{") => :CU_BRACKET_OPEN
        ],
    '}' => [
        forward_match("}") => :CU_BRACKET_CLOSE,
        ],
    '<' => [
        forward_match("<!--")           => :COMMENT_OPEN,
        forward_match("<", ALPHA_LATIN) => :AUTOLINK_OPEN
        ],
    '>' => [
        forward_match(">") => :AUTOLINK_CLOSE
        ],
    '[' => [
        forward_match("[") => :SQ_BRACKET_OPEN,
        ],
    ']' => [
        forward_match("]") => :SQ_BRACKET_CLOSE,
        ],
    ':' => [
        F_EMOJI => :CAND_EMOJI,
        ],
    '|' => [
        forward_match("|") => :PIPE,
        ],
    '-' => [
        forward_match("-->") => :COMMENT_CLOSE,
        ],
    '+' => [
        forward_match("+++", ['\n']) => :MD_DEF_BLOCK
        ],
    '\\' => [
        # -- special characters (https://www.amp-what.com/unicode/search)
        # commonmark specs: https://spec.commonmark.org/0.30/#backslash-escapes
        # OK \#\$\~\{\}\*\@\\\!\"\%\&\'\+\,\-\.\/\:\;\<\=\>\?\^\_\`\|
        # NO \[\]\(\)
        forward_match("\\\\") => :LINEBREAK,         # --> <br/>
        forward_match("\\", SPACE_CHAR) => :CHAR_92,
        # -- maths
        forward_match("\\[")  => :MATH_DISPL_B_OPEN,      # \[ ...
        forward_match("\\]")  => :MATH_DISPL_B_CLOSE,     #    ... \]
        # -- latex
        forward_match("\\newenvironment", ['{']) => :LX_NEWENVIRONMENT,
        forward_match("\\newcommand",     ['{']) => :LX_NEWCOMMAND,
        forward_match("\\begin",          ['{']) => :LX_BEGIN,
        forward_match("\\end",            ['{']) => :LX_END,
        F_LX_COMMAND                             => :LX_COMMAND,  # \command⎵*
        # -- other special characters
        forward_match("\\*")  => :CHAR_42,
        forward_match("\\_")  => :CHAR_95,
        forward_match("\\`")  => :CHAR_96,
        forward_match("\\@")  => :CHAR_64,
        forward_match("\\#")  => :CHAR_35,
        forward_match("\\{")  => :CHAR_123,
        forward_match("\\}")  => :CHAR_125,
        forward_match("\\\$") => :CHAR_36,
        forward_match("\\~")  => :CHAR_126,
        forward_match("\\!")  => :CHAR_33,
        forward_match("\\\"") => :CHAR_34,
        forward_match("\\%")  => :CHAR_37,
        forward_match("\\&")  => :CHAR_38,
        forward_match("\\'")  => :CHAR_39,
        forward_match("\\+")  => :CHAR_43,
        forward_match("\\,")  => :CHAR_44,
        forward_match("\\-")  => :CHAR_45,
        forward_match("\\.")  => :CHAR_46,
        forward_match("\\/")  => :CHAR_47,
        forward_match("\\:")  => :CHAR_58,
        forward_match("\\;")  => :CHAR_59,
        forward_match("\\<")  => :CHAR_60,
        forward_match("\\=")  => :CHAR_61,
        forward_match("\\>")  => :CHAR_62,
        forward_match("\\?")  => :CHAR_63,
        forward_match("\\^")  => :CHAR_94,
        forward_match("\\|")  => :CHAR_124,
        ],
    '@' => [
        forward_match("@def", [' '])    => :MD_DEF_OPEN,    # @def var = ...
        forward_match("@@", SPACE_CHAR) => :DIV_CLOSE,      # @@⎵*
        F_DIV_OPEN                      => :DIV_OPEN,       # @@dname
        ],
    '#' => [
        forward_match("#",      [' ']) => :H1_OPEN,
        forward_match("##",     [' ']) => :H2_OPEN,
        forward_match("###",    [' ']) => :H3_OPEN,
        forward_match("####",   [' ']) => :H4_OPEN,
        forward_match("#####",  [' ']) => :H5_OPEN,
        forward_match("######", [' ']) => :H6_OPEN,
        ],
    '&' => [
        F_HTML_ENTITY => :CHAR_HTML_ENTITY,
        ],
    '$' => [
        forward_match("\$", ['$'], false) => :MATH_INLINE,   # $⎵*
        forward_match("\$\$")             => :MATH_DISPL_A,  # $$⎵*
        ],
    '_' => [
        forward_match("___", ['_'], false) => :EM_STRONG,
        forward_match("__",  ['_'], false) => :STRONG,
        forward_match("_",   ['_'], false) => :EM,
        ],
    '`' => [
        forward_match("`",   ['`'], false) => :CODE_SINGLE,  # `⎵
        forward_match("``",  ['`'], false) => :CODE_DOUBLE,  # ``⎵*
        forward_match("```", ['`'], false) => :CODE_TRIPLE,  # ```⎵*
        forward_match("`"^4, ['`'], false) => :CODE_QUAD,    # ````⎵*
        forward_match("`"^5, ['`'], false) => :CODE_PENTA,   # `````⎵*
        ],
    '*' => [
        forward_match("***", ['*'], false) => :EM_STRONG,
        forward_match("**",  ['*'], false) => :STRONG,
        forward_match("*",   ['*'], false) => :EM,
        ]
    )  # end dict

"""
MD_MATH_TOKENS

Tokens that should be considered within a math environment.
"""
const MD_MATH_TOKENS = Dict{Char, Vector{Pair{TokenFinder, Symbol}}}(
    '{' => [
        forward_match("{") => :CU_BRACKET_OPEN
        ],
    '}' => [
        forward_match("}") => :CU_BRACKET_CLOSE,
        ],
    '\\' => [
        forward_match("\\begin", ['{']) => :LX_BEGIN,
        forward_match("\\end", ['{'])   => :LX_END,
        F_LX_COMMAND                    => :LX_COMMAND,  # \command⎵*
        ],
    ) # end dict

"""
    END_OF_LINE

All tokens that indicate the end of a line.
"""
const END_OF_LINE = (:LINE_RETURN, :EOS)

"""
    MD_IGNORE

Tokens that may be left over after partition but should be ignored in text blocks.
"""
const MD_IGNORE = (:SOS, :LINE_RETURN, :EOS)

"""
    MD_HEADERS

Tokens for headers.
"""
const MD_HEADERS = (
    :H1_OPEN,
    :H2_OPEN,
    :H3_OPEN,
    :H4_OPEN,
    :H5_OPEN,
    :H6_OPEN,
)


text:
#!/usr/bin/env julia

using Random

# mod must always be positive

function mod(a, n)
    if (t = a % n) < 0
        n + t
    else
        t
    end
end

#=
a^b (mod n) using the method of repeated squares.

The key here is that every integer can be written as a sum of powers of 2 (binary numbers)
and that includes the exponent. By repeated squaring we get a raised to a power of 2. Also
recall that a^b * a^c = a^(b + c), so rather than adding we multiply since we are dealing
with the exponent.
=#

function powerMod(a, d, n)
    v = 1 # Value
    p = a # Powers of a
    while d > 0
        if isodd(d) # 1 bit in the exponent
           v = mod(v * p, n)
        end
        p = mod(p^2, n) # Next power of two
        d >>>= 1
    end
    v
end

#=
Greatest common divisor, Euclidean version.
=#

function gcd(a, b)
    while b ≠ 0
        a, b = b, mod(a, b)
    end
    a
end

#=
Witness loop of the Miller-Rabin probabilistic primality test.
=#

function witness(a, n)
    u, t = n - 1, 0
    while iseven(u) # n = u * 2^t + 1
        t += 1   # Increase exponent
        u >>>= 1 # Decrease the multiplier
    end
    x = powerMod(a, u, n)
    for i in 1:t
        y = powerMod(x, 2, n)
        if y == 1 && x ≠ 1 && x ≠ n - 1
            return true
        end
        x = y
    end
    x ≠ 1
end

#=
Miller-Rabin probabilistic primality test: the chance of being wrong is ≈ 1/4 each pass through
the loop.
=#

function isPrime(n, k)
    if n < 2 || (n ≠ 2 && iseven(n)) # 0, 1, and even except for 2 are not prime.
        return false
    elseif n < 4 # 3 is prime
        return true
    end # We must test all others
    for j in 1:k
        a = rand(2:n - 2) # Choose a random witness
        if witness(a, n)
            return false
        end
    end
    true
end

#=
We need a random prime number in [low, high] and for now a 4^–100 chance of a composite is
good enough.
=#

function randomPrime(low, high)
    guess = 0 # Certainly not prime!
    while !isPrime(guess, 100)
        guess = rand(low:high) # Half will be even, the rest have Pr[prime] ≈ 1/log(N).
    end
    guess
end

#=
A safe prime is the one following a Sophie German prime. If prime(p) and prime(2p + 1) then
2p + 1 is a safe prime.
=#

function safePrime(low, high)
    p = randomPrime(low, high)
    while !isPrime(2 * p + 1,100)
        p = randomPrime(low, high)
    end
    return 2 * p + 1
end

# Rabin prime

function RabinPrime(safe, low, high)
    f = safe ? safePrime : randomPrime
    p = f(low, high)
    while mod(p, 4) != 3
        p = f(low, high)
    end
    return p
end

# Euclidean extended greatest common divisor

function extendedGCD(a, b)
    (r, rP) = (a, b)
    (s, sP) = (1, 0)
    (t, tP) = (0, 1)
    while rP != 0
        q = r ÷ rP
        (r, rP) = (rP, r - q * rP)
        (s, sP) = (sP, s - q * sP)
        (t, tP) = (tP, t - q * tP)
        end
    return (r, (s, t))
end

#=
Add comment
=#

function key(safe, k)
    x = k + 32 # Make room for the tag
    p = RabinPrime(safe, big"2"^(x - 1), big"2"^x - 1)
    q = RabinPrime(safe, big"2"^(x - 1), big"2"^x - 1)
    while p == q
        q = RabinPrime(safe, big"2"^(x - 1), big"2"^x - 1)
    end
    return (p * q, (p, q))
end

#=
We need dicriminate amongst the four possible roots.
=#

using CRC # requires the CRC packages
crc32 = crc(CRC_32)

h = BigInt(crc32("Michael O. Rabin"))

#=
Add comment
=#

encrypt(m, n) = mod((m << 32 + h)^2, n)

#=
Add comment
=#

function decrypt(m, key)
    (p, q) = key
    n = p * q
    (g, (yP, yQ)) = extendedGCD(p, q)
    mP = powerMod(m, (p + 1) ÷ 4, p)
    mQ = powerMod(m, (q + 1) ÷ 4, q)
    x = mod(yP * p * mQ + yQ * q * mP, n)
    y = mod(yP * p * mQ - yQ * q * mP, n)
    msgs = [x, n - x, y, n - y]
    for d in  msgs
        if mod(d, 2^32) == h
            return d ÷ 2^32
        end
    end
    1279869254 # FAIL
end

function encode(s)
    sum::BigInt = 0
    pow::BigInt = 1
    for c in s
        sum += pow * BigInt(c)
        pow *= 256
    end
    BigInt(sum)
end

#=
Transform a BigInt back into a string, subtracting off the 0xAA. We treat it as a base-256
integer and just pull off the digits.
=#

function decode(n)
    s = ""
    while n > 0
        s = s * Char(mod(n, 256))
        n ÷= 256
    end
    s
end

print("How many bits? ")

bits = parse(Int64, readline())

(n, k) = key(false, bits)

println("n = $n")
println("k = $k")

print(">> ")
for m in eachline()
    c = encrypt(encode(m), n); println("En[$m] = $c")
    t = decode(decrypt(c, k)); println("De[$c] = $t")
    print(">> ")
end

max_stars_repo_name: darrelllong/Rabin-in-Julia

max_issues_repo_name: darrelllong/Rabin-in-Julia

max_forks_repo_name: darrelllong/Rabin-in-Julia


text:
export ScikitClassifier

@with_kw mutable struct ScikitClassifier <: BaseEstimator
    pyo::PyObject = PyNULL()
    name::String = "logistic"
    alpha::Float32 = 1.0
    penalty::String = "l2"
    loss::String = "hinge"
    gamma::Float32 = 0
    n_neighbors::Int = 10
    num_layers::Int = 1
    hidden_size::Int = 100
    lr::Float32 = 0.001
end

is_classifier(::ScikitClassifier) = true

function gridparams(m::ScikitClassifier)
    grid = [
        [
            "name" => ["logistic"],
            "alpha" => [0.001, 0.01, 0.1, 1, 5],
            "penalty" => ["l2", "l1"]
        ],
        [
            "name" => ["linearsvc"],
            "alpha" => [0.001, 0.01, 0.1, 1, 5],
            [
                ["penalty" => ["l2"], "loss" => ["hinge", "squared_hinge"]],
                ["penalty" => ["l1"], "loss" => ["squared_hinge"]]
            ]
        ],
        [
            "name" => ["svc"],
            "alpha" => [0.001, 0.01, 0.1, 1, 5, 10, 20, 50],
            "gamma" => [0, 0.01, 0.1, 1, 10, 100]
        ],
        [
            "name" => ["mlp"],
            "num_layers" => [1, 2],
            "hidden_size" => [10, 100],
            "lr" => [0.001, 0.01]
        ]
    ]
    gridparams(grid)
end

function fit!(m::ScikitClassifier, x, y, w = nothing; columns = nothing)
    columns = something(columns, string.(1:size(x, 1)))
    @unpack name, alpha, penalty, loss, gamma = m
    @unpack n_neighbors, num_layers, hidden_size, lr = m
    x, y, w = pymat(x), vec(y), vec(w)
    if name == "logistic"
        @from sklearn.linear_model imports LogisticRegression
        pyo = LogisticRegression(C = 1 / alpha, penalty = penalty, solver = "saga", n_jobs = 20, verbose = true)
    elseif name == "linearsvc"
        @from sklearn.svm imports LinearSVC
        pyo = LinearSVC(dual = loss == "hinge", C = 1 / alpha, penalty = penalty, loss = loss, verbose = true)
    elseif name == "svc"
        use_thunder = true
        try
            @from thundersvm imports SVC
            pyo = SVC(C = 0.5 / alpha, gamma = gamma == 0 ? "auto" : gamma, n_jobs = 20, verbose = true)
        catch e
            use_thunder = false
            @from sklearn.svm imports SVC
            pyo = SVC(C = 0.5 / alpha, gamma = gamma == 0 ? "auto" : gamma, verbose = true)
        end
    elseif name == "mlp"
        @from sklearn.neural_network imports MLPClassifier
        pyo = MLPClassifier(alpha = alpha, learning_rate_init = lr, verbose = true,
                        hidden_layer_sizes = ntuple(i -> hidden_size, num_layers))
    end
    if name == "svc" && use_thunder || name == "mlp"
        pyo.fit(x, y)
    else
        pyo.fit(x, y, sample_weight = w)
    end
    @pack! m = pyo
    BSON.bson("model.bson", model = m)
    visualize(m, columns)
    return m
end

function predict_proba(m::ScikitClassifier, x)
    @unpack name, pyo = m
    if occursin("svc", name)
        pyo.decision_function(pymat(x))
    else
        pyo.predict_proba(pymat(x))
    end
end

predict(m::ScikitClassifier, x) = m.pyo.predict(pymat(x))

modelhash(m::ScikitClassifier) = hash(m.pyo)

function visualize(m::ScikitClassifier, columns)
    @unpack name, pyo = m
    if name == "logistic" || name == "linearsvc"
        write_feaimpt(pyo.coef_, columns)
    end
end


text:
struct _bi01IP
    C::Matrix{Int} # objective functions, k=1..2, j=1..n
    A::Matrix{Int} # matrix of constraints, i=1..m, j=1..n
    b::Vector{Int} # right-hand side, i=1..m
end

mutable struct solution
    val_objectif::Vector{Int64}
    sol::Vector{Int64}
    cout::Vector{Int64}
end


text:
module TestMajorityVote
	using Base.Test
	using kNN

	votes = ["A", "B", "B", "A", "C", "B"]
    @test isequal(kNN.majority_vote(votes), "B")

    ["A", "B", "B", "A", "C", "B", "C", "C"]
    @test isequal(kNN.majority_vote(votes), "B")

    votes = ["A", "B", "B", "A", "C", "B", "C", "C", "C"]
    @test isequal(kNN.majority_vote(votes), "C")
end


text:
module Lu_bench
#
  importall Shared_types
  using Shared_funcs
  importall Arg
  import Lu
#
  import Lu_julia_dense
  import Lu_julia_sparse
  import Lu_elim_I
#
  type DbProblem
    creator::Function
    descr::String
  end
#
  module_db = [Lu_julia_dense, Lu_julia_sparse, Lu_elim_I]
  module_name_map = merge([x.ArgName => x for x in module_db], [string(i) => module_db[i] for i in 1:length(module_db)])
  prob_db = DbProblem[]
  prob_last = Array(Any, 0)
#
  function problem_t1(params::Params)
    return Lu.create_Problem(params, [1 1; 1 0])
  end; push!(prob_db, DbProblem(problem_t1, "t1"))
  function problem_file(params::Params)
    println("[", params["file"], "]")
    A = Shared_funcs.read_matrix_mtx(params["file"])
    return Lu.create_Problem(params, A)
  end; push!(prob_db, DbProblem(problem_file, "file"))
  function problem_random(params::Params)
  return Lu.create_Problem(params, Shared_funcs.random_matrix(params))
  end; push!(prob_db, DbProblem(problem_random, "random"))
  function solve_problem(code_module::Module, params::Params, lu_prob)
    T = Lu.get_t(lu_prob)
    dat = code_module.construct_dat(T)
    code_module.fill_dat(lu_prob, dat)
    sol = Lu.construct_solution(T, params)
    code_module.solve_dat(dat, sol)
    return sol
  end
  function solve_problem_timed(code_module::Module, params::Params, lu_prob)
    @time sol = solve_problem(code_module, params, lu_prob)
    return sol
  end
  function solve_problem(params::Params, lu_prob)
    chosen = []
    if params["module"] == "all"
      chosen = [x.ArgName for x in module_db]
    elseif params["module"] == "choose"
      chosen = [module_db[i].ArgName for i in Shared_funcs.print_choose(["$(x.ArgName): $(x.Descr)"for x in module_db])]
    else
      chosen = split(params["module"], ",")
    end
    return [solve_problem_timed(module_name_map[x], params, lu_prob) for x in chosen]
  end
  function solve(arg_str::String = "/prob:t1")
    def_params = { "prob"=>"t1", "type"=>"Float32", "module"=>"choose" }
    params = deepcopy(def_params)
    arg_get(arg_create(arg_str), params)
    #
    prob_key = params["prob"]
    if (length(prob_key) == 0 || prob_key == "all")
      for i = 1:length(prob_db) solve(i, arg_str); end;
      return
    end
    #
    found_prob = false; dbprob = 0;
    for i = 1:length(prob_db)
      iprob = prob_db[i]
      if (contains(iprob.descr, prob_key))
        found_prob = true; dbprob = iprob;
        break
      end
    end
    #
    if (found_prob)
      println()
      println("Problem: '", dbprob.descr, "'")
      println("Module: '", params["module"], "'")
      println("------------")
      lu_prob = dbprob.creator(params)
      #
      if (length(Lu_bench.prob_last) > 0) pop!(Lu_bench.prob_last) end
      push!(Lu_bench.prob_last, deepcopy(lu_prob))
      #
      println("+++++++",
        " n:", Lu.get_r(lu_prob), ", m:", Lu.get_c(lu_prob),
        ", density:", Shared_funcs.format_percent(Lu.comp_density(lu_prob)),
        ", type:", Lu.get_t(lu_prob),
        ", form:", Lu.get_form(lu_prob),
        " +++++++")
      end
      sols = solve_problem(params, lu_prob)
      println("************")
      println("Distance: ", [Lu.calc_solution_distance(lu_prob, x) for x in sols])
      println("")
      return sols
    end
end

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
import DataStructures: OrderedDict

# This set of tests directly uses the high-level `convert` functions
# And checks the behaviour is as expected.

J.def_GLOBAL_LXDEFS!()
cmd   = st -> J.convert_md(st, collect(values(J.GLOBAL_LXDEFS)))
chtml = t -> J.convert_html(t...)
conv  = st -> st |> cmd |> chtml

# convenience function that squeezes out all whitespaces and line returns out of a string
# and checks if the resulting strings are equal. When expecting a specific string +- some
# spaces, this is very convenient. Use == if want to check exact strings.
isapproxstr(s1::String, s2::String) =
    isequal(map(s->replace(s, r"\s|\n"=>""), (s1, s2))...)

# this is a slightly ridiculous requirement but apparently the `eval` blocks
# don't play well with Travis nor windows while testing, so you just need to forcibly
# specify that LinearAlgebra and Random are used (even though the included block says
# the same thing).
if get(ENV, "CI", "false") == "true" || Sys.iswindows()
    import Pkg; Pkg.add("LinearAlgebra"); Pkg.add("Random");
    using LinearAlgebra, Random;
end

# takes md input and outputs the html (good for integration testing)
function seval(st)
    J.def_GLOBAL_PAGE_VARS!()
    J.def_GLOBAL_LXDEFS!()
    m, v = J.convert_md(st, collect(values(J.GLOBAL_LXDEFS)))
    h = J.convert_html(m, v)
    return h
end

function set_globals()
    J.def_GLOBAL_PAGE_VARS!()
    J.def_GLOBAL_LXDEFS!()
    empty!(J.LOCAL_PAGE_VARS)
    J.def_LOCAL_PAGE_VARS!()
    J.CUR_PATH_WITH_EVAL[] = ""
end

function explore_md_steps(mds)
    J.def_GLOBAL_PAGE_VARS!()
    J.def_GLOBAL_LXDEFS!()
    empty!(J.RSS_DICT)

    J.def_LOCAL_PAGE_VARS!()
    J.def_PAGE_HEADERS!()
    J.def_PAGE_EQREFS!()
    J.def_PAGE_BIBREFS!()
    J.def_PAGE_FNREFS!()
    J.def_PAGE_LINK_DEFS!()

    steps = OrderedDict{Symbol,NamedTuple}()

    # tokenize
    tokens  = J.find_tokens(mds, J.MD_TOKENS, J.MD_1C_TOKENS)
    fn_refs = J.validate_footnotes!(tokens)
    J.validate_headers!(tokens)
    tokens  = J.find_indented_blocks(tokens, mds)
    steps[:tokenization] =  (tokens=tokens,)
    steps[:fn_validation] = (fn_refs=fn_refs,)

    # ocblocks
    blocks, tokens = J.find_all_ocblocks(tokens, J.MD_OCB_ALL)
    J.merge_indented_blocks!(blocks, mds)
    J.filter_indented_blocks!(blocks)
    steps[:ocblocks] = (blocks=blocks, tokens=tokens)

    # filter
    filter!(τ -> τ.name ∉ J.L_RETURNS, tokens)
    steps[:filter] = (tokens=tokens, blocks=blocks)

    J.validate_and_store_link_defs!(blocks)

    # latex commands
    lxdefs, tokens, braces, blocks = J.find_md_lxdefs(tokens, blocks)
    lxcoms, _ = J.find_md_lxcoms(tokens, lxdefs, braces)
    steps[:latex] = (lxdefs=lxdefs, tokens=tokens, braces=braces,
                     blocks=blocks, lxcoms=lxcoms)

    sp_chars = J.find_special_chars(tokens)
    steps[:spchars] = (spchars=sp_chars,)

    blocks2insert = J.merge_blocks(lxcoms, J.deactivate_divs(blocks), sp_chars)
    steps[:blocks2insert] = (blocks2insert=blocks2insert,)

    inter_md, mblocks = J.form_inter_md(mds, blocks2insert, lxdefs)
    steps[:inter_md] = (inter_md=inter_md, mblocks=mblocks)

    inter_html = J.md2html(inter_md; stripp=false)
    steps[:inter_html] = (inter_html=inter_html,)

    lxcontext = J.LxContext(lxcoms, lxdefs, braces)
    hstring   = J.convert_inter_html(inter_html, mblocks, lxcontext)
    steps[:hstring] = (hstring=hstring,)

    return steps
end

function explore_h_steps(hs, allvars=J.PageVars())
    steps = OrderedDict{Symbol,NamedTuple}()

    tokens = J.find_tokens(hs, J.HTML_TOKENS, J.HTML_1C_TOKENS)
    steps[:tokenization] = (tokens=tokens,)

    hblocks, tokens = J.find_all_ocblocks(tokens, J.HTML_OCB)
    filter!(hb -> hb.name != :COMMENT, hblocks)
    steps[:ocblocks] = (hblocks=hblocks, tokens=tokens)

    qblocks = J.qualify_html_hblocks(hblocks)
    steps[:qblocks] = (qblocks=qblocks,)

    fhs = J.process_html_qblocks(hs, allvars, qblocks)
    steps[:fhs] = (fhs=fhs,)

    return steps
end