text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# Modified version of https://github.com/genkuroki/public/blob/main/0016/apricot/julia_translation_of_python_reimpl.ipynb

# %%
#using Seaborn
using ScikitLearn: @sk_import
@sk_import datasets: fetch_covtype
#using Random
#using StatsBase: sample

# %%
digits_data = fetch_covtype()

# %%
X_digits = permutedims(abs.(digits_data["data"]))
summary(X_digits)

# %%
"""`calculate_gains!(X, gains, current_values, idxs, current_concave_values_sum)` mutates `gains` only"""
function calculate_gains!(X, gains, current_values, idxs, current_concave_values_sum)
    Threads.@threads for i in eachindex(idxs)
        @inbounds idx = idxs[i]
        @inbounds gains[i] = sum(j -> sqrt(current_values[j] + X[j, idx]), axes(X, 1))
    end
    gains .-= current_concave_values_sum
end

@doc calculate_gains!

# %%
"""
The revised version of using popat!(idxs, idx) without mask
"""
function fit_popat(X, k; calculate_gains! = calculate_gains!)
    d, n = size(X)

    cost = 0.0

    ranking = Int[]
    total_gains = Float64[]

    current_values = zeros(d)
    current_concave_values_sum = sum(sqrt, current_values)

    idxs = collect(1:n)

    gains = zeros(n)
    while cost < k
        calculate_gains!(X, gains, current_values, idxs, current_concave_values_sum)

        idx = argmax(gains)
        best_idx = idxs[idx]
        curr_cost = 1.0

        cost + curr_cost > k && break

        cost += curr_cost
        # Calculate gains
        gain = gains[idx] * curr_cost

        # Select next
        current_values .+= @view X[:, best_idx]
        current_concave_values_sum = sum(sqrt, current_values)

        push!(ranking, best_idx)
        push!(total_gains, gain)

        popat!(idxs, idx)
    end
    return ranking, total_gains
end

# %%
"""
The revised version of using mask::BitVector with findall(mask)
"""
function fit_bitvector(X, k; calculate_gains! = calculate_gains!)
    d, n = size(X)

    cost = 0.0

    ranking = Int[]
    total_gains = Float64[]

    mask = trues(n) # `false` stands for "masked".
    current_values = zeros(d)
    current_concave_values_sum = sum(sqrt, current_values)

    idxs = collect(1:n)

    gains = zeros(n)
    while cost < k
        calculate_gains!(X, gains, current_values, idxs, current_concave_values_sum)

        idx = argmax(gains)
        best_idx = idxs[idx]
        curr_cost = 1.0

        cost + curr_cost > k && break

        cost += curr_cost
        # Calculate gains
        gain = gains[idx] * curr_cost

        # Select next
        current_values .+= @view X[:, best_idx]
        current_concave_values_sum = sum(sqrt, current_values)

        push!(ranking, best_idx)
        push!(total_gains, gain)

        mask[best_idx] = 0
        idxs = findall(mask)
    end
    return ranking, total_gains
end

# %%
"""
The revised version of using mask::Vector{Float64} with findall(mask .== 0)
"""
function fit_f64vector(X, k; calculate_gains! = calculate_gains!)
    d, n = size(X)

    cost = 0.0

    ranking = Int[]
    total_gains = Float64[]

    mask = zeros(n)
    current_values = zeros(d)
    current_concave_values_sum = sum(sqrt, current_values)

    idxs = collect(1:n)

    gains = zeros(n)
    while cost < k
        calculate_gains!(X, gains, current_values, idxs, current_concave_values_sum)

        idx = argmax(gains)
        best_idx = idxs[idx]
        curr_cost = 1.0

        cost + curr_cost > k && break

        cost += curr_cost
        # Calculate gains
        gain = gains[idx] * curr_cost

        # Select next
        current_values .+= @view X[:, best_idx]
        current_concave_values_sum = sum(sqrt, current_values)

        push!(ranking, best_idx)
        push!(total_gains, gain)

        mask[best_idx] = 1
        idxs = findall(mask .== 0)
    end
    return ranking, total_gains
end

# %%
"""
The original version of using mask::Vector{Float64} with findall(==(0), mask)
"""
function fit_f64vector_org(X, k; calculate_gains! = calculate_gains!)
    d, n = size(X)

    cost = 0.0

    ranking = Int[]
    total_gains = Float64[]

    mask = zeros(n)
    current_values = zeros(d)
    current_concave_values_sum = sum(sqrt, current_values)

    idxs = collect(1:n)

    gains = zeros(n)
    while cost < k
        calculate_gains!(X, gains, current_values, idxs, current_concave_values_sum)

        idx = argmax(gains)
        best_idx = idxs[idx]
        curr_cost = 1.0

        cost + curr_cost > k && break

        cost += curr_cost
        # Calculate gains
        gain = gains[idx] * curr_cost

        # Select next
        current_values .+= @view X[:, best_idx]
        current_concave_values_sum = sum(sqrt, current_values)

        push!(ranking, best_idx)
        push!(total_gains, gain)

        mask[best_idx] = 1
        idxs = findall(==(0), mask)
    end
    return ranking, total_gains
end

# %%
k = 1000

# %%
@time ranking0_pa, gains0_pa = fit_popat(X_digits, k; calculate_gains! = calculate_gains!);

# %%
@time ranking0_pa, gains0_pa = fit_popat(X_digits, k; calculate_gains! = calculate_gains!);

# %%
@time ranking0_pa, gains0_pa = fit_popat(X_digits, k; calculate_gains! = calculate_gains!);

# %%
@time ranking0_bv, gains0_bv = fit_bitvector(X_digits, k; calculate_gains! = calculate_gains!);

# %%
@time ranking0_bv, gains0_bv = fit_bitvector(X_digits, k; calculate_gains! = calculate_gains!);

# %%
@time ranking0_bv, gains0_bv = fit_bitvector(X_digits, k; calculate_gains! = calculate_gains!);

# %%
@time ranking0_f64v, gains0_f64v = fit_f64vector(X_digits, k; calculate_gains! = calculate_gains!);

# %%
@time ranking0_f64v, gains0_f64v = fit_f64vector(X_digits, k; calculate_gains! = calculate_gains!);

# %%
@time ranking0_f64v, gains0_f64v = fit_f64vector(X_digits, k; calculate_gains! = calculate_gains!);

# %%
@time ranking0_org, gains0_org = fit_f64vector_org(X_digits, k; calculate_gains! = calculate_gains!);

# %%
@time ranking0_org, gains0_org = fit_f64vector_org(X_digits, k; calculate_gains! = calculate_gains!);

# %%
@time ranking0_org, gains0_org = fit_f64vector_org(X_digits, k; calculate_gains! = calculate_gains!);

# %%
@show ranking0_pa == ranking0_bv == ranking0_f64v == ranking0_org
@show gains0_pa == gains0_bv == gains0_f64v == gains0_org;

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using FuzzySets, Test

step_size = 0.005
levels = collect(0:step_size:1)

A₁ = FuzzyNumber(levels, number=1)
A₂ = FuzzyNumber(levels, number=2)
B₁ = FuzzyNumber(levels, number=0)
B₂ = FuzzyNumber(levels, number=0)

A⃗ = FuzzyVector([A₁, A₂])
B⃗ = FuzzyVector([B₁])

@test FuzzySets.d_dsw(A⃗, A⃗) == SingletonFuzzyNumber(levels, number=0)
@test FuzzySets.d_dsw(A⃗, B⃗) == false

@test FuzzySets.d_interval(A⃗, A⃗) == SingletonFuzzyNumber(levels, number=0)
@test FuzzySets.d_interval(A⃗, B⃗) == false

a = [Interval(-3), Interval(-2)]
b = [Interval(-3), Interval(0)]
@test FuzzySets.d_dsw(a, b; squared=false) == Interval(((-3 - (-3))^2 + (-2 - 0)^2)^0.5)
@test FuzzySets.d_dsw(a, b; squared=true) == Interval(((-3 - (-3))^2 + (-2 - 0)^2))
@test FuzzySets.d_interval(a, b; squared=false) == Interval(((-3 - (-3))^2 + (-2 - 0)^2)^0.5)
@test FuzzySets.d_interval(a, b; squared=true) == Interval(((-3 - (-3))^2 + (-2 - 0)^2))

a = [Interval(-3, -1), Interval(-2, 0)]
b = [Interval(-3, -1), Interval(0, 2)]
d_interval_sol = (a[1] - b[1])^2 + (a[2] - b[2])^2
@test FuzzySets.d_interval(a, b; squared=true) == d_interval_sol
@test FuzzySets.d_interval(a, b; squared=false) == d_interval_sol^0.5

x = [Interval(-3), Interval(-0)]
c1 = [Interval(-2.1), Interval(0)]
c2 = [Interval(2.1), Interval(0)]
d² = [FuzzySets.d_interval(x, c1; squared=true), FuzzySets.d_interval(x, c2; squared=true)]
@test round(FuzzySets.u_interval(d², 1; m=2), digits=2) == Interval(0.97)
@test round(FuzzySets.u_interval(d², 2; m=2), digits=2) == Interval(0.03)

# Test cases from "Analysis and Efficient Implementation of a Linguistic Fuzzy C-Means" by S Auephanwiriyakul and JM Keller
X = [Interval(1, 4), Interval(2, 3)]
u = [Interval(0, 0.4), Interval(0.6, 1)]
m = 2.0
@test round(FuzzySets.km_iwa(X, u, bound="lower", m=m), digits=2) == 1.69
@test round(FuzzySets.km_iwa(X, u, bound="upper", m=m), digits=2) == 3.31
@test_throws ArgumentError FuzzySets.km_iwa(X, u, bound="something", m=2.0)

# Test cases from "Uncertain Rule-Based Fuzzy Systems Introduction and New Directions" by Jerry M. Mendel
X = [Interval(1), Interval(2), Interval(3), Interval(4), Interval(5)]
u = [Interval(0.6, 0.9), Interval(0.5, 0.7), Interval(0.65, 0.8), Interval(0.2, 0.4), Interval(0.3, 0.75)]
m = 1.0
@test round(FuzzySets.km_iwa(X, u, bound="lower", m=m), digits=3) == 2.382
@test round(FuzzySets.km_iwa(X, u, bound="upper", m=m), digits=3) == 3.069

butterfly_patterns = FuzzyVector[
    FuzzyVector([FuzzyNumber(levels, number=-3), FuzzyNumber(levels, number=-2)]),
    FuzzyVector([FuzzyNumber(levels, number=-3), FuzzyNumber(levels, number=0)]),
    FuzzyVector([FuzzyNumber(levels, number=-3), FuzzyNumber(levels, number=2)]),
    FuzzyVector([FuzzyNumber(levels, number=-2), FuzzyNumber(levels, number=-1)]),
    FuzzyVector([FuzzyNumber(levels, number=-2), FuzzyNumber(levels, number=0)]),
    FuzzyVector([FuzzyNumber(levels, number=-2), FuzzyNumber(levels, number=1)]),
    FuzzyVector([FuzzyNumber(levels, number=-1), FuzzyNumber(levels, number=0)]),
    FuzzyVector([FuzzyNumber(levels, number=0), FuzzyNumber(levels, number=0)]),
    FuzzyVector([FuzzyNumber(levels, number=1), FuzzyNumber(levels, number=0)]),
    FuzzyVector([FuzzyNumber(levels, number=2), FuzzyNumber(levels, number=-1)]),
    FuzzyVector([FuzzyNumber(levels, number=2), FuzzyNumber(levels, number=0)]),
    FuzzyVector([FuzzyNumber(levels, number=2), FuzzyNumber(levels, number=1)]),
    FuzzyVector([FuzzyNumber(levels, number=3), FuzzyNumber(levels, number=-2)]),
    FuzzyVector([FuzzyNumber(levels, number=3), FuzzyNumber(levels, number=0)]),
    FuzzyVector([FuzzyNumber(levels, number=3), FuzzyNumber(levels, number=2)]),
]
cluster_prototypes = FuzzyVector[
    FuzzyVector([FuzzyNumber(levels, number=-2.1), FuzzyNumber(levels, number=0)]),
    FuzzyVector([FuzzyNumber(levels, number=2.1), FuzzyNumber(levels, number=0)]),
]

m = 2.0

grades = [FuzzySets.u_dsw(FuzzySets.cut(butterfly_patterns[1], α), FuzzySets.cut(cluster_prototypes, α), 1, m=m) for α in levels]
A = FuzzyNumber(levels, grades)
draw(A)
grades = [FuzzySets.u_dsw(FuzzySets.cut(butterfly_patterns[1], α), FuzzySets.cut(cluster_prototypes, α), 2, m=m) for α in levels]
B = FuzzyNumber(levels, grades)
draw(B)

m = 2.0
U = [[0.86, 0.14], [0.97, 0.03], [0.86, 0.14], [0.95, 0.05], [1.00, 0.00], [0.95, 0.05], [0.89, 0.11], [0.50, 0.50], [0.11, 0.89], [0.05, 0.95], [0.00, 1.00], [0.05, 0.95], [0.14, 0.86], [0.03, 0.97], [0.14, 0.86]]
b = FuzzySets.cut(cluster_prototypes, 1)
for i = 1:length(U)
    a = FuzzySets.cut(butterfly_patterns[i], 1)
    @test round(FuzzySets.u_dsw(a, b, 1, m=m).left, digits=2) == round(FuzzySets.u_dsw(a, b, 1, m=m).right, digits=2) == U[i][1]
    @test round(FuzzySets.u_dsw(a, b, 2, m=m).left, digits=2) == round(FuzzySets.u_dsw(a, b, 2, m=m).right, digits=2) == U[i][2]
end

m = 1.25
U = [[1.00, 0.00], [1.00, 0.00], [1.00, 0.00], [1.00, 0.00], [1.00, 0.00], [1.00, 0.00], [1.00, 0.00], [0.50, 0.50], [0.00, 1.00], [0.00, 1.00], [0.00, 1.00], [0.00, 1.00], [0.00, 1.00], [0.00, 1.00], [0.00, 1.00]]
b = FuzzySets.cut(cluster_prototypes, 1)
for i = 1:length(U)
    a = FuzzySets.cut(butterfly_patterns[i], 1)
    @test round(FuzzySets.u_dsw(a, b, 1, m=m).left, digits=2) == round(FuzzySets.u_dsw(a, b, 1, m=m).right, digits=2) == U[i][1]
    @test round(FuzzySets.u_dsw(a, b, 2, m=m).left, digits=2) == round(FuzzySets.u_dsw(a, b, 2, m=m).right, digits=2) == U[i][2]
end

m = 6
U = [[0.59, 0.41], [0.67, 0.33], [0.59, 0.41], [0.64, 0.36], [0.82, 0.18], [0.64, 0.36], [0.60, 0.40], [0.50, 0.50], [0.40, 0.60], [0.36, 0.64], [0.18, 0.82], [0.36, 0.64], [0.41, 0.59], [0.33, 0.67], [0.41, 0.59]]
b = FuzzySets.cut(cluster_prototypes, 1)
for i = 1:2
    for k = 1:length(U)
        a = FuzzySets.cut(butterfly_patterns[k], 1)
        print(round(FuzzySets.u_dsw(a, b, i, m=m).left, digits=2), " ")
    end
    println()
end
for i = 1:length(U)
    a = FuzzySets.cut(butterfly_patterns[i], 1)
    @test round(FuzzySets.u_dsw(a, b, 1, m=m).left, digits=2) == round(FuzzySets.u_dsw(a, b, 1, m=m).right, digits=2) == U[i][1]
    @test round(FuzzySets.u_dsw(a, b, 2, m=m).left, digits=2) == round(FuzzySets.u_dsw(a, b, 2, m=m).right, digits=2) == U[i][2]
end

m = 6
h = 1 / (1 - m)
c1 = [-2.1, 0]
c2 = [2.1, 0]
x = [-3, 0]
term1 = sum((x - c1).^2)^h
term2 = sum((x - c2).^2)^h
denom = term1 + term2
u1_crisp = term1 / denom
u2_crisp = term2 / denom

n_points = 30
@test FuzzySets.get_sequential_endpoints(n_points, 1) == ones(n_points)
@test FuzzySets.get_sequential_endpoints(n_points, 2^n_points) == ones(n_points) .+ 1

ext: jl


text:
function ArcArc(arc1::Geometry, arc2::Geometry)
    if arc1.GType != "arc" || arc2.GType != "arc"
        error("Input Type error： arc needed")
    end
    Cp1 = ArcCenter(arc1)
    Cp2 = ArcCenter(arc2)
    r1 = 1/abs(arc1.Sp.cur)
    r2 = 1/abs(arc2.Sp.cur)
    Dis = Distance2D(Cp1, Cp2)
    if Dis > r1 + r2 || Dis < abs(r1 - r2) # separate or contain
        return (intersect=false, type="none")
    elseif Dis <= 1.0e-6 && r1 == r2 # coincident || overlap
        error("arcs are in the same circle, may be handled later.")
    else
        a = (r1^2 - r2^2 + Dis^2)/(2*Dis)
        h = sqrt(r1^2 - a^2)
        xm = Cp1.x + a*(Cp2.x - Cp1.x)/Dis
        ym = Cp1.y + a*(Cp2.y - Cp1.y)/Dis
        Ipx1 = xm + h*(Cp2.y - Cp1.y)/Dis
        Ipy1 = ym - h*(Cp2.x - Cp1.x)/Dis
        Ipx2 = xm - h*(Cp2.y - Cp1.y)/Dis
        Ipy2 = ym + h*(Cp2.x - Cp1.x)/Dis
        Ip1 = Point2D{Float64}(Ipx1, Ipy1)
        Ip2 = Point2D{Float64}(Ipx2, Ipy2)
        Infos = Array{NamedTuple{(:Ip, :Dist1, :Dist2),Tuple{Point2D{Float64}, Float64, Float64}}, 1}()
        for Ipi in [Ip1, Ip2]
            Result1 = PointInSector(Ipi, arc1)
            Result2 = PointInSector(Ipi, arc2)
            if Result1[1] && Result2[1]
                push!(Infos, (Ip=Ipi, Dist1=Result1[2], Dist2=Result2[2]))
            end
        end
        if length(Infos) == 0
            return (intersect=false, type="none")
        else
            return (intersect=true, type="cross", Infos[1]...) # actually both points should be used to calculate TTC. But now I just use one and this will be changed later
        end
    end
end


text:
module EFT
export two_sum, fast_two_sum, two_prod

import ..SIMDops
using ..SIMDops: vabs, vfma, vifelse, vless

include("errorfree.jl")

end


text:
#data type for the investment universe
struct Universe
  df::DataFrame #source dataframe must be properly formatted with a :date column but can include extra observations
  datebegin::Date
  dateend::Date
  returns::Matrix{MFloat64} #should contain the numerical returns data in the DataFrame
  #NOTE for performance reasons, this is M by N!!!

  μ::Vector{Float64} #we can calc this in advance
  N::Int #number of observations in the sample period
  M::Int #number of return series

  missingindices::Vector{Vector{Int}} # a collection of all non-missing indices in returns
end

#creates an investment universe which will be used to generate portfolio returns
function Universe(df::AbstractDataFrame; datebegin::Date = minimum(df[:date]),
  dateend::Date = maximum(df[:date]))::Universe
  subdf::SubDataFrame = view(df,(df[:date] .≥ datebegin) .& (df[:date] .≤ dateend),:)

  local N::Int = size(subdf,1) #get the size in number of observations
  local M::Int = size(subdf, 2) - 1 #exclude date column

  local xcols::Vector{Symbol} = setdiff(names(df), [:date])
  local returns::Matrix{MFloat64} = Matrix{MFloat64}(subdf[xcols])'

  #compute the mean for use later
  local μ::Vector{Float64} = (r::Int->mean(skipmissing(returns[r,:]))).(1:M)

  #map the missing values
  local missingindices::Vector{Vector{Int}} = Vector{Vector{Float64}}()
  for c ∈ 1:N
    push!(missingindices, collect(1:M)[(!ismissing).(returns[:, c])])
  end

  return Universe(df,datebegin, dateend, returns, μ, N, M, missingindices)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using NLopt

# test model from MathProgBase, using undocumented
# function Base.find_package for now (julia#27592):
nlp = joinpath(dirname(Base.find_package(NLopt, "MathProgBase")), "..", "test", "nlp.jl")
include(nlp)

nlptest(NLoptSolver(algorithm=:LD_SLSQP))
# test derivative-free
rosenbrocktest(NLoptSolver(algorithm=:LN_PRAXIS))


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

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Random, StatsBase, DataFrames, Plots
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.12"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
begin
	triangDist = TriangularDist(0,2,1)
	xGrid = 0:0.1:2
		N = 10^6
	inverseSampledData = quantile.(triangDist,rand(N))
end

# ╔═╡ b690890a-2135-11eb-1abc-933c757a93f8
begin
	histogram( inverseSampledData, bins=30, normed=true,
		ylims=(0,1.1), label="Inverse transform data")
	plot!( xGrid, pdf.(triangDist,xGrid), c=:red, lw=4,
		xlabel="x", label="PDF", ylabel = "Density", legend=:topright)
end

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.12"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═b690890a-2135-11eb-1abc-933c757a93f8
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
using Test
using CartesianGeneticProgramming
import Cambrian
import Random

test_filename = string(@__DIR__, "/test.yaml")

@testset "Mutation" begin
    cfg = get_config(test_filename)
    parent = CGPInd(cfg)

    # Uniform mutation
    child = uniform_mutate(cfg, parent)
    @test any(parent.chromosome .!= child.chromosome)
    @test any(parent.genes .!= child.genes)

    # Goldman mutation : ensure structural difference
    child = goldman_mutate(cfg, parent)
    @test any(parent.chromosome .!= child.chromosome)
    @test any(parent.genes .!= child.genes)

    # Profiling mutation: ensure output different for provided inputs
    inputs = rand(cfg.n_in, 1)
    child = profiling_mutate(cfg, parent, inputs)
    @test any(parent.chromosome .!= child.chromosome)
    @test any(parent.genes .!= child.genes)

    out_parent = CartesianGeneticProgramming.process(parent, inputs[:, 1])
    out_child = CartesianGeneticProgramming.process(child, inputs[:, 1])
    @test any(out_parent .!= out_child)
end

# TODO: convert or remove
"""
A simple function module handling integers.
"""
module IntFunctions
    global arity = Dict()
    function fgen(name::Symbol, ar::Int, s1::Union{Symbol, Expr})
        @eval function $name(x::Int64, y::Int64, p::Array{Float64}=Float64[])::Int64
            $s1
        end
        arity[String(name)] = ar
    end
    fgen(:f_add, 2, :(x + y))
    fgen(:f_subtract, 2, :(x - y))
    fgen(:f_mult, 2, :(x * y))
    fgen(:f_div, 1, :(convert(Int64, floor(x / 2))))
    fgen(:f_abs, 2, :(abs(x)))
end

"""
Constructor for custom CGPInd based on configuration.
"""
function CustomCGPInd(cfg::NamedTuple)
    buffer = zeros(Int64, cfg.rows * cfg.columns + cfg.n_in)
    CartesianGeneticProgramming.CGPInd(cfg; buffer=buffer)
end

"""
Constructor for custom CGPInd based on configuration and chromosome.
"""
function CustomCGPInd(cfg::NamedTuple, chromosome::Array{Float64})
    buffer = zeros(Int64, cfg.rows * cfg.columns + cfg.n_in)
    CartesianGeneticProgramming.CGPInd(cfg, chromosome; buffer=buffer)
end

# TODO: randomly fails
#
@testset "Mutation with custom CGPInd" begin
    cfg = get_config(test_filename; functions=["f_add", "f_subtract", "f_mult", "f_div", "f_abs"])
    parent = CustomCGPInd(cfg)

    # Uniform mutation
    child = uniform_mutate(cfg, parent, init_function=CustomCGPInd)
    @test any(parent.chromosome .!= child.chromosome)
    @test any(parent.genes .!= child.genes)

    # Goldman mutation : ensure structural difference
    child = goldman_mutate(cfg, parent, init_function=CustomCGPInd)
    @test any(parent.chromosome .!= child.chromosome)
    @test any(parent.genes .!= child.genes)

    # Profiling mutation: ensure output different for provided inputs
    inputs = rand(Int64, cfg.n_in, 1)
    child = profiling_mutate(cfg, parent, inputs, init_function=CustomCGPInd)
    @test any(parent.chromosome .!= child.chromosome)
    @test any(parent.genes .!= child.genes)

    out_parent = CartesianGeneticProgramming.process(parent, inputs[:, 1])
    out_child = CartesianGeneticProgramming.process(child, inputs[:, 1])
    @test any(out_parent .!= out_child)
end


text:
version = uint8(131)

# See http://erlang.org/doc/apps/erts/erl_ext_dist.html
newfloattag = uint8(70)
bitbintag = uint8(77)
compressedtag = uint8(80)
atomcachereftag = uint8(82)
smallinttag = uint8(97)
inttag = uint8(98)
floattag = uint8(99)
atomtag = uint8(100)
reftag = uint8(101)
porttag = uint8(102)
pidtag = uint8(103)
smalltupletag = uint8(104)
largetupletag = uint8(105)
niltag = uint8(106)
stringtag = uint8(107)
listtag = uint8(108)
bintag = uint8(109)
smallbiginttag = uint8(110)
largebiginttag = uint8(111)
newfuntag = uint8(112)
exporttag = uint8(113)
newreftag = uint8(114)
smallatomtag = uint8(115)
maptag = uint8(116)
funtag = uint8(117)
atomutf8tag = uint8(118)
smallatomutf8tag = uint8(119)


text:
# TODO
# Only for 1dof joints
mutable struct PID{T,N} <: Controller
    P::SVector{N,T}
    I::SVector{N,T}
    D::SVector{N,T}

    eqcids::SVector{N,Int64}
    goals::SVector{N,T}
    integratederrors::SVector{N,T}
    lasterrors::SVector{N,T}

    control!::Function

    function PID(mechanism::Mechanism{T}, eqcid::Int64, goal; P = zero(T), I = zero(T), D = zero(T), controlfunction::Function = control_pid!) where T
        eqc = geteqconstraint(mechanism, eqcid)
        Nb = 6 * length(unique(eqc.childids))
        Nc = ConstrainedDynamics.length(eqc)
        @assert Nb-Nc == 1 "Only 1 DOF joints are supported"

        new{T,1}([P], [I], [D], [eqcid], [goal], [0], [0], controlfunction)
    end

    function PID(mechanism, eqcids::AbstractVector{Int64}, goals::AVec;
            P::AVec = zeros(T,length(eqcids)), I::AVec = zeros(T,length(eqcids)), D::AVec = zeros(T,length(eqcids)),
            controlfunction::Function = control_pid!
        ) where {T, AVec<:AbstractVector{T}}

        N = length(eqcids)
        for eqcid in eqcids
            eqc = geteqconstraint(mechanism, eqcid)
            Nb = 6 * length(unique(eqc.childids))
            Nc = 0
            Nc = ConstrainedDynamics.length(eqc)
            @assert Nb-Nc == 1 "Only 1 DOF joints are supported"
        end
        new{T,N}(P, I, D, eqcids, goals, zeros(T,N), zeros(T,N), controlfunction)
    end
end

# Angle correction only works for revolute joints
function stateError_pid(mechanism, eqc, goal; anglecorrection = true)
    if anglecorrection && (typeof(eqc.constraints[1]) <: ConstrainedDynamics.Translational3 && typeof(eqc.constraints[2]) <: ConstrainedDynamics.Rotational2)
        mincoords = minimalCoordinates(mechanism, eqc)[1]
        tempdiff = goal - mincoords
        if tempdiff > pi
            return tempdiff - 2*pi
        elseif tempdiff < -pi
            return tempdiff + 2*pi
        else
            return tempdiff
        end
    else
        return goal - minimalCoordinates(mechanism, eqc)[1]
    end
end

# @generated function error_pid(mechanism, eqcids::SVector{N,Int64}, goals) where {N}
#     vec = [:(goals[$i] - minimalCoordinates(mechanism, geteqconstraint(mechanism, eqcids[$i]))[1]) for i = 1:N]
#     return :(svcat($(vec...)))
# end

function error_pid(mechanism, eqcids::SVector{N,Int64}, goals; anglecorrection = true) where {N}
    vec = [stateError_pid(mechanism, geteqconstraint(mechanism, eqcids[i]), goals[i]; anglecorrection = anglecorrection) for i = 1:N]
    return svcat(vec...)
end

function control_pid!(mechanism, pid::PID{T,N}, k) where {T,N}
    Δt = mechanism.Δt

    currenterrors = error_pid(mechanism, pid.eqcids, pid.goals, anglecorrection = true)
    k==1 && (pid.lasterrors = currenterrors)

    perrors = currenterrors
    pid.integratederrors += currenterrors * Δt
    differentialerrors = (currenterrors - pid.lasterrors) / Δt

    u = pid.P .* perrors + pid.I .* pid.integratederrors + pid.D .* differentialerrors

    pid.lasterrors = currenterrors

    for i=1:N
        setForce!(mechanism, geteqconstraint(mechanism, pid.eqcids[i]), SVector(u[i]))
    end

    return
end


text:
"""
    Compute ln p(y[i] | y[-i], z[-i], z[i] = k), given `i` and `k`,
    see Bouchard-Côté et al. (2017, eq. 32) and Sudderth (2006, eq. 2.163)
"""
function log_pl(
    s::SuffStats{A, B, C},
    c::ChainState,
    i::Int,
    k::Int
) where {A, B, C}
    @unpack y, x, K̄, D, ν, r, u, S = s
    @unpack γ, z = c
    yi = y[i]
    zi = z[i]
    j  = x[i]^γ[x[i]]

    if zi == k
        r1 = r[j, k];   r0 = r1 - 1
        ν1 = ν[j, k];   ν0 = ν1 - 1
        S1 = S[j, k]
        S0 = S[j, K̄]
        S0.factors .= S1.factors
        u[j, K̄] .= (yi - u[j, k]) * √(r1 / r0)
        lowrankdowndate!(S0, u[j, K̄])
    else
        r0 = r[j, k];   r1 = r0 + 1
        ν0 = ν[j, k];   ν1 = ν0 + 1
        S0 = S[j, k]
        S1 = S[j, K̄]
        S1.factors .= S0.factors
        u[j, K̄] .= (r0 * u[j, k] + yi) / r1
        u[j, K̄] .= (yi - u[j, K̄]) * √(r1 / r0)
        lowrankupdate!(S1, u[j, K̄])
    end

    val =
        0.5ν0 * logdet(S0) -
        0.5ν1 * logdet(S1) +
        logmvgamma(D, ν1 / 2) -
        logmvgamma(D, ν0 / 2) +
        0.5D * log(r0 / r1) -
        0.5D * log(π)
end

"""
    Compute ln p({y[i]: x[i] = j, z[i] = k} | x, z), given `j` and `k`,
    see Bouchard-Côté et al. (2017, eq. 31)
"""
function log_ml(s::SuffStats{A, B, C}, j::Int, k::Int) where {A, B, C}
    @unpack m, D, r, ν, S = s
    @unpack r0, ν0, S0 = m
    r1 = r[j, k]
    ν1 = ν[j, k]
    S1 = S[j, k]
    ll =
        0.5ν0 * logdet(S0) -
        0.5ν1 * logdet(S1) +
        logmvgamma(D, ν1 / 2) -
        logmvgamma(D, ν0 / 2) +
        0.5D * log(r0 / r1) -
        0.5D * log(π) * (r1 - r0)
end


text:
#InspectDR: Secondary annotation (Cairo layer)
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#
const ALIGN_MAP = Dict{Symbol, CAlignment}(
	:tl => ALIGN_TOP | ALIGN_LEFT,
	:tc => ALIGN_TOP | ALIGN_HCENTER,
	:tr => ALIGN_TOP | ALIGN_RIGHT,

	:cl => ALIGN_VCENTER | ALIGN_LEFT,
	:cc => ALIGN_VCENTER | ALIGN_HCENTER,
	:cr => ALIGN_VCENTER | ALIGN_RIGHT,

	:bl => ALIGN_BOTTOM | ALIGN_LEFT,
	:bc => ALIGN_BOTTOM | ALIGN_HCENTER,
	:br => ALIGN_BOTTOM | ALIGN_RIGHT,
)

#==Rendering text annotation
===============================================================================#
function render(ctx::CairoContext, rstrip::RStrip2D, a::TextAnnotation)
	align = get(ALIGN_MAP, a.align, ALIGN_BOTTOM | ALIGN_LEFT)
	angle = deg2rad(a.angle)
	pt = axis2dev(a.pos, rstrip.xf, rstrip.ixf, rstrip.bb)
	render(ctx, a.text, pt, a.font, angle=angle, align=align)
	return
end

#==Rendering markers
===============================================================================#

function render(ctx::CairoContext, rstrip::RStrip2D, mkr::HVMarker)
	graphbb = rstrip.bb #WANTCONST
	if :none == mkr.line.style
		return
	end

	setlinestyle(ctx, LineStyle(mkr.line))
	pt = axis2aloc(mkr.pos, rstrip.ixf)
	pt = apply(rstrip.xf, pt)

	if mkr.vmarker
		drawline(ctx, Point2D(pt.x, graphbb.ymin), Point2D(pt.x, graphbb.ymax))
	end
	if mkr.hmarker
		drawline(ctx, Point2D(graphbb.xmin, pt.y), Point2D(graphbb.xmax, pt.y))
	end

	return
end

#==Rendering polyline for annotation
===============================================================================#
function render(ctx::CairoContext, rstrip::RStrip2D, a::PolylineAnnotation)
	x = a.x; y = a.y
	setlinestyle(ctx, LineStyle(a.line))
	pt = apply(rstrip.xf, Point2D(x[1], y[1]))
	Cairo.move_to(ctx, pt.x, pt.y)
	for i in 2:length(x)
		pt = apply(rstrip.xf, Point2D(x[i], y[i]))
		Cairo.line_to(ctx, pt.x, pt.y)
	end
	if a.closepath
		Cairo.close_path(ctx)
		renderfill(ctx, a.fillcolor)
	end
	Cairo.stroke(ctx)
	return
end

#==Generic rendering algorithm for annotations
===============================================================================#

#Basic render function for PlotAnnotation types:
function render(ctx::CairoContext, rstrip::RStrip2D, a::PlotAnnotation)
	if 0 == a.strip || a.strip == rstrip.istrip
		render(ctx, rstrip, a)
	end
	return
end

function render(ctx::CairoContext, rstrip::RStrip2D, alist::Vector{T}) where T<:PlotAnnotation
	for a in alist
		render(ctx, rstrip, a)
	end
	return
end

#Last line


text:
#######################################################################
# Coverage.jl
# Take Julia test coverage results and bundle them up in JSONs
# https://github.com/JuliaCI/Coverage.jl
#######################################################################

using Coverage, Test, LibGit2

if VERSION < v"1.1"
isnothing(x) = false
isnothing(x::Nothing) = true
end

@testset "Coverage" begin
withenv("DISABLE_AMEND_COVERAGE_FROM_SRC" => nothing) do

@testset "iscovfile" begin
    # test our filename matching. These aren't exported functions but it's probably
    # a good idea to have explicit tests for them, as they're used to match files
    # that get deleted
    @test Coverage.iscovfile("test.jl.cov")
    @test Coverage.iscovfile("test.jl.2934.cov")
    @test Coverage.iscovfile("/home/somebody/test.jl.2934.cov")
    @test !Coverage.iscovfile("test.ji.2934.cov")
    @test !Coverage.iscovfile("test.jl.2934.cove")
    @test !Coverage.iscovfile("test.jicov")
    @test !Coverage.iscovfile("test.c.cov")
    @test Coverage.iscovfile("test.jl.cov", "test.jl")
    @test !Coverage.iscovfile("test.jl.cov", "other.jl")
    @test Coverage.iscovfile("test.jl.8392.cov", "test.jl")
    @test Coverage.iscovfile("/somedir/test.jl.8392.cov", "/somedir/test.jl")
    @test !Coverage.iscovfile("/otherdir/test.jl.cov", "/somedir/test.jl")
end

@testset "isfuncexpr" begin
    @test Coverage.isfuncexpr(:(f() = x))
    @test Coverage.isfuncexpr(:(function() end))
    @test Coverage.isfuncexpr(:(function g() end))
    @test Coverage.isfuncexpr(:(function g() where {T} end))
    @test !Coverage.isfuncexpr("2")
    @test !Coverage.isfuncexpr(:(f = x))
    @test Coverage.isfuncexpr(:(() -> x))
    @test Coverage.isfuncexpr(:(x -> x))
    @test Coverage.isfuncexpr(:(f() where A = x))
    @test Coverage.isfuncexpr(:(f() where A where B = x))
end

@testset "Processing coverage" begin
    cd(dirname(@__DIR__)) do
        datadir = joinpath("test", "data")
        # Process a saved set of coverage data...
        r = process_file(joinpath(datadir, "Coverage.jl"))

        # ... and memory data
        malloc_results = analyze_malloc(datadir)
        filename = joinpath(datadir, "testparser.jl.9172.mem")
        @test malloc_results == [Coverage.MallocInfo(96669, filename, 2)]

        lcov = IOBuffer()
        # we only have a single file, but we want to test on the Vector of file results
        LCOV.write(lcov, FileCoverage[r])
        expected = read(joinpath(datadir, "tracefiles", "expected.info"), String)
        if Sys.iswindows()
            expected = replace(expected, "\r\n" => "\n")
            expected = replace(expected, "SF:test/data/Coverage.jl\n" => "SF:test\\data\\Coverage.jl\n")
        end
        @test String(take!(lcov)) == expected

        # LCOV.writefile is a short-hand for writing to a file
        lcov = joinpath(datadir, "lcov_output_temp.info")
        LCOV.writefile(lcov, FileCoverage[r])
        expected = read(joinpath(datadir, "tracefiles", "expected.info"), String)
        if Sys.iswindows()
            expected = replace(expected, "\r\n" => "\n")
            expected = replace(expected, "SF:test/data/Coverage.jl\n" => "SF:test\\data\\Coverage.jl\n")
        end
        @test String(read(lcov)) == expected
        # tear down test file
        rm(lcov)

        # test that reading the LCOV file gives the same data
        lcov = LCOV.readfolder(datadir)
        @test length(lcov) == 1
        r2 = lcov[1]
        r2_filename = r2.filename
        if Sys.iswindows()
            r2_filename = replace(r2_filename, '/' => '\\')
        end
        @test r2_filename == r.filename
        @test r2.source == ""
        @test r2.coverage == r.coverage[1:length(r2.coverage)]
        @test all(isnothing, r.coverage[(length(r2.coverage) + 1):end])
        lcov2 = [FileCoverage(r2.filename, "sourcecode", Coverage.CovCount[nothing, 1, 0, nothing, 3]),
                 FileCoverage("file2.jl", "moresource2", Coverage.CovCount[1, nothing, 0, nothing, 2]),]
        lcov = merge_coverage_counts(lcov, lcov2, lcov)
        @test length(lcov) == 2
        r3 = lcov[1]
        @test r3.filename == r2.filename
        @test r3.source == "sourcecode"
        r3cov = Coverage.CovCount[x === nothing ? nothing : x * 2 for x in r2.coverage]
        r3cov[2] += 1
        r3cov[3] = 0
        r3cov[5] = 3
        @test r3.coverage == r3cov
        r4 = lcov[2]
        @test r4.filename == "file2.jl"
        @test r4.source == "moresource2"
        @test r4.coverage == lcov2[2].coverage

        # Test a file from scratch
        srcname = joinpath("test", "data", "testparser.jl")
        covname = srcname*".cov"
        # clean out any previous coverage files. Don't use clean_folder because we
        # need to preserve the pre-baked coverage file Coverage.jl.cov
        clean_file(srcname)
        cmdstr = "include($(repr(srcname))); using Test; @test f2(2) == 4"
        run(`$(Base.julia_cmd()) --startup-file=no --code-coverage=user -e $cmdstr`)
        r = process_file(srcname, datadir)

        target = Coverage.CovCount[nothing, 2, nothing, 0, nothing, 0, nothing, nothing, nothing, nothing, 0, nothing, nothing, nothing, nothing, nothing, 0, nothing, nothing, 0, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing]
        target_disabled = map(x -> (x !== nothing && x > 0) ? x : nothing, target)
        @test r.coverage == target

        covtarget = (sum(x->x !== nothing && x > 0, target), sum(x->x !== nothing, target))
        @test get_summary(r) == covtarget
        @test get_summary(process_folder(datadir)) == (98, 106)

        r_disabled = withenv("DISABLE_AMEND_COVERAGE_FROM_SRC" => "yes") do
            process_file(srcname, datadir)
        end

        @test r_disabled.coverage == target_disabled
        amend_coverage_from_src!(r_disabled.coverage, r_disabled.filename)
        @test r_disabled.coverage == target

        # Handle an empty coverage vector
        emptycov = FileCoverage("", "", [])
        @test get_summary(emptycov) == (0, 0)

        @test isempty(Coverage.process_cov(joinpath("test", "fakefile"), datadir))

        # test clean_folder
        # set up the test folder
        datadir_temp = joinpath("test", "data_temp")
        cp(datadir, datadir_temp)
        # run clean_folder
        clean_folder(datadir_temp)
        # .cov files should be deleted
        @test !isfile(joinpath(datadir_temp, "Coverage.jl.cov"))
        # other files should remain untouched
        @test isfile(joinpath(datadir_temp, "Coverage.jl"))
        # tear down test data
        rm(datadir_temp; recursive=true)
    end
end

@testset "codecovio.jl" begin
    # these methods are only used for testing the token generation for local repos
    # and CI whilst not breaking the current API
    construct_uri_string_local(dir=pwd(); kwargs...) = Coverage.Codecov.construct_uri_string(
        ;Coverage.Codecov.add_local_to_kwargs(dir; kwargs...)...)

    construct_uri_string_ci(;kwargs...) = Coverage.Codecov.construct_uri_string(
        ;Coverage.Codecov.add_ci_to_kwargs(;kwargs...)...)

    # empty file coverage for testing
    fcs = FileCoverage[]

    # set up base system ENV vars for testing
    withenv(
        "CODECOV_URL" => nothing,
        "CODECOV_TOKEN" => nothing,
        "TRAVIS" => nothing,
        "TRAVIS_BRANCH" => nothing,
        "TRAVIS_COMMIT" => nothing,
        "TRAVIS_PULL_REQUEST" => nothing,
        "TRAVIS_JOB_ID" => nothing,
        "TRAVIS_REPO_SLUG" => nothing,
        "TRAVIS_JOB_NUMBER" => nothing,
        "APPVEYOR" => nothing,
        "APPVEYOR_PULL_REQUEST_NUMBER" => nothing,
        "APPVEYOR_ACCOUNT_NAME" => nothing,
        "APPVEYOR_PROJECT_SLUG" => nothing,
        "APPVEYOR_BUILD_VERSION" => nothing,
        "APPVEYOR_REPO_BRANCH" => nothing,
        "APPVEYOR_REPO_COMMIT" => nothing,
        "APPVEYOR_REPO_NAME" => nothing,
        "APPVEYOR_JOB_ID" => nothing,
        ) do

        # test local submission process (but only if we are in a git repo)

        _dotgit = joinpath(dirname(@__DIR__), ".git")

        if isdir(_dotgit) || isfile(_dotgit)
            LibGit2.with(LibGit2.GitRepoExt(pwd())) do repo
                # default values
                codecov_url = construct_uri_string_local()
                @test occursin("codecov.io", codecov_url)
                @test occursin("commit", codecov_url)
                @test occursin("branch", codecov_url)
                @test !occursin("service", codecov_url)

                # env var url override
                withenv( "CODECOV_URL" => "https://enterprise-codecov-1.com" ) do

                    codecov_url = construct_uri_string_local()
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("commit", codecov_url)
                    @test occursin("branch", codecov_url)
                    @test !occursin("service", codecov_url)

                    # function argument url override
                    codecov_url = construct_uri_string_local(codecov_url = "https://enterprise-codecov-2.com")
                    @test occursin("enterprise-codecov-2.com", codecov_url)
                    @test occursin("commit", codecov_url)
                    @test occursin("branch", codecov_url)
                    @test !occursin("service", codecov_url)

                    # env var token
                    withenv( "CODECOV_TOKEN" => "token_name_1" ) do

                        codecov_url = construct_uri_string_local()
                        @test occursin("enterprise-codecov-1.com", codecov_url)
                        @test occursin("token=token_name_1", codecov_url)
                        @test !occursin("service", codecov_url)

                        # function argument token url override
                        codecov_url = construct_uri_string_local(token="token_name_2")
                        @test occursin("enterprise-codecov-1.com", codecov_url)
                        @test occursin("token=token_name_2", codecov_url)
                        @test !occursin("service", codecov_url)
                    end
                end
            end
        end

        # test faulty non-CI submission

        @test_throws ErrorException Coverage.Codecov.submit(fcs; dry_run = true)

        # test travis-ci submission process

        # set up travis env
        withenv(
            "TRAVIS" => "true",
            "TRAVIS_BRANCH" => "t_branch",
            "TRAVIS_COMMIT" => "t_commit",
            "TRAVIS_PULL_REQUEST" => "t_pr",
            "TRAVIS_JOB_ID" => "t_job_id",
            "TRAVIS_REPO_SLUG" => "t_slug",
            "TRAVIS_JOB_NUMBER" => "t_job_num",
            ) do

            # default values
            codecov_url = construct_uri_string_ci()
            @test occursin("codecov.io", codecov_url)
            @test occursin("service=travis-org", codecov_url)
            @test occursin("branch=t_branch", codecov_url)
            @test occursin("commit=t_commit", codecov_url)
            @test occursin("pull_request=t_pr", codecov_url)
            @test occursin("job=t_job_id", codecov_url)
            @test occursin("slug=t_slug", codecov_url)
            @test occursin("build=t_job_num", codecov_url)

            # env var url override
            withenv( "CODECOV_URL" => "https://enterprise-codecov-1.com" ) do

                codecov_url = construct_uri_string_ci()
                @test occursin("enterprise-codecov-1.com", codecov_url)
                @test occursin("service=travis-org", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("pull_request=t_pr", codecov_url)
                @test occursin("job=t_job_id", codecov_url)
                @test occursin("slug=t_slug", codecov_url)
                @test occursin("build=t_job_num", codecov_url)

                # function argument url override
                codecov_url = construct_uri_string_ci(;codecov_url = "https://enterprise-codecov-2.com")
                @test occursin("enterprise-codecov-2.com", codecov_url)
                @test occursin("service=travis-org", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("pull_request=t_pr", codecov_url)
                @test occursin("job=t_job_id", codecov_url)
                @test occursin("slug=t_slug", codecov_url)
                @test occursin("build=t_job_num", codecov_url)

                # env var token
                withenv( "CODECOV_TOKEN" => "token_name_1" ) do

                    codecov_url = construct_uri_string_ci()
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("token=token_name_1", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("pull_request=t_pr", codecov_url)
                    @test occursin("job=t_job_id", codecov_url)
                    @test occursin("slug=t_slug", codecov_url)
                    @test occursin("build=t_job_num", codecov_url)

                    # function argument token url override
                    codecov_url = construct_uri_string_ci(token="token_name_2")
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("token=token_name_2", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("pull_request=t_pr", codecov_url)
                    @test occursin("job=t_job_id", codecov_url)
                    @test occursin("slug=t_slug", codecov_url)
                    @test occursin("build=t_job_num", codecov_url)
                end
            end
        end

        # test appveyor submission process

        # set up appveyor env
        withenv(
            "APPVEYOR" => "true",
            "APPVEYOR_PULL_REQUEST_NUMBER" => "t_pr",
            "APPVEYOR_ACCOUNT_NAME" => "t_account",
            "APPVEYOR_PROJECT_SLUG" => "t_slug",
            "APPVEYOR_BUILD_VERSION" => "t_version",
            "APPVEYOR_REPO_BRANCH" => "t_branch",
            "APPVEYOR_REPO_COMMIT" => "t_commit",
            "APPVEYOR_REPO_NAME" => "t_repo",
            "APPVEYOR_JOB_ID" => "t_job_num",
            ) do

                # default values
                codecov_url = construct_uri_string_ci()
                @test occursin("codecov.io", codecov_url)
                @test occursin("service=appveyor", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("pull_request=t_pr", codecov_url)
                @test occursin("job=t_account%2Ft_slug%2Ft_version", codecov_url)
                @test occursin("build=t_job_num", codecov_url)

                # env var url override
                withenv( "CODECOV_URL" => "https://enterprise-codecov-1.com" ) do

                    codecov_url = construct_uri_string_ci()
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("service=appveyor", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("pull_request=t_pr", codecov_url)
                    @test occursin("job=t_account%2Ft_slug%2Ft_version", codecov_url)
                    @test occursin("build=t_job_num", codecov_url)

                    # function argument url override
                    codecov_url = construct_uri_string_ci(codecov_url = "https://enterprise-codecov-2.com")
                    @test occursin("enterprise-codecov-2.com", codecov_url)
                    @test occursin("service=appveyor", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("pull_request=t_pr", codecov_url)
                    @test occursin("job=t_account%2Ft_slug%2Ft_version", codecov_url)
                    @test occursin("build=t_job_num", codecov_url)

                    # env var token
                    withenv( "CODECOV_TOKEN" => "token_name_1" ) do

                        codecov_url = construct_uri_string_ci()
                        @test occursin("enterprise-codecov-1.com", codecov_url)
                        @test occursin("token=token_name_1", codecov_url)
                        @test occursin("branch=t_branch", codecov_url)
                        @test occursin("commit=t_commit", codecov_url)
                        @test occursin("pull_request=t_pr", codecov_url)
                        @test occursin("job=t_account%2Ft_slug%2Ft_version", codecov_url)
                        @test occursin("build=t_job_num", codecov_url)

                        # function argument token url override
                        codecov_url = construct_uri_string_ci(token="token_name_2")
                        @test occursin("enterprise-codecov-1.com", codecov_url)
                        @test occursin("token=token_name_2", codecov_url)
                        @test occursin("branch=t_branch", codecov_url)
                        @test occursin("commit=t_commit", codecov_url)
                        @test occursin("pull_request=t_pr", codecov_url)
                        @test occursin("job=t_account%2Ft_slug%2Ft_version", codecov_url)
                        @test occursin("build=t_job_num", codecov_url)
                    end
                end
            end

        # test circle ci submission process

        # set up circle ci env
        withenv(
            "CIRCLECI" => "true",
            "CIRCLE_PR_NUMBER" => "t_pr",
            "CIRCLE_PROJECT_USERNAME" => "t_proj",
            "CIRCLE_BRANCH" => "t_branch",
            "CIRCLE_SHA1" => "t_commit",
            "CIRCLE_PROJECT_REPONAME" => "t_repo",
            "CIRCLE_BUILD_URL" => "t_url",
            "CIRCLE_BUILD_NUM" => "t_num",
            ) do

            # default values
            codecov_url = construct_uri_string_ci()
            @test occursin("codecov.io", codecov_url)
            @test occursin("service=circleci", codecov_url)
            @test occursin("branch=t_branch", codecov_url)
            @test occursin("commit=t_commit", codecov_url)
            @test occursin("pull_request=t_pr", codecov_url)
            @test occursin("build_url=t_url", codecov_url)
            @test occursin("build=t_num", codecov_url)

            # env var url override
            withenv( "CODECOV_URL" => "https://enterprise-codecov-1.com" ) do

                codecov_url = construct_uri_string_ci()
                @test occursin("enterprise-codecov-1.com", codecov_url)
                @test occursin("service=circleci", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("pull_request=t_pr", codecov_url)
                @test occursin("build_url=t_url", codecov_url)
                @test occursin("build=t_num", codecov_url)

                # function argument url override
                codecov_url = construct_uri_string_ci(codecov_url="https://enterprise-codecov-2.com")
                @test occursin("enterprise-codecov-2.com", codecov_url)
                @test occursin("service=circleci", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("pull_request=t_pr", codecov_url)
                @test occursin("build_url=t_url", codecov_url)
                @test occursin("build=t_num", codecov_url)

                # env var token
                withenv( "CODECOV_TOKEN" => "token_name_1" ) do

                    codecov_url = construct_uri_string_ci()
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("token=token_name_1", codecov_url)
                    @test occursin("service=circleci", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("pull_request=t_pr", codecov_url)
                    @test occursin("build_url=t_url", codecov_url)
                    @test occursin("build=t_num", codecov_url)

                    # function argument token url override
                    codecov_url = construct_uri_string_ci(token="token_name_2")
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("service=circleci", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("pull_request=t_pr", codecov_url)
                    @test occursin("build_url=t_url", codecov_url)
                    @test occursin("build=t_num", codecov_url)
                end
            end
        end

        # test Jenkins ci submission process

        # set up Jenkins ci env
        withenv(
            "JENKINS" => "true",
            "GIT_BRANCH" => "t_branch",
            "GIT_COMMIT" => "t_commit",
            "JOB_NAME" => "t_job",
            "BUILD_ID" => "t_num",
            "BUILD_URL" => "t_url",
            "JENKINS_URL" => "t_jenkins_url",
            ) do

            # default values
            codecov_url = construct_uri_string_ci()
            @test occursin("codecov.io", codecov_url)
            @test occursin("service=jenkins", codecov_url)
            @test occursin("branch=t_branch", codecov_url)
            @test occursin("commit=t_commit", codecov_url)
            @test occursin("build_url=t_url", codecov_url)
            @test occursin("build=t_num", codecov_url)

            # env var url override
            withenv( "CODECOV_URL" => "https://enterprise-codecov-1.com" ) do

                codecov_url = construct_uri_string_ci()
                @test occursin("enterprise-codecov-1.com", codecov_url)
                @test occursin("service=jenkins", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("build_url=t_url", codecov_url)
                @test occursin("build=t_num", codecov_url)

                # function argument url override
                codecov_url = construct_uri_string_ci(codecov_url="https://enterprise-codecov-2.com")
                @test occursin("enterprise-codecov-2.com", codecov_url)
                @test occursin("service=jenkins", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("build_url=t_url", codecov_url)
                @test occursin("build=t_num", codecov_url)

                # env var token
                withenv( "CODECOV_TOKEN" => "token_name_1" ) do

                    codecov_url = construct_uri_string_ci()
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("token=token_name_1", codecov_url)
                    @test occursin("service=jenkins", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("build_url=t_url", codecov_url)
                    @test occursin("build=t_num", codecov_url)

                    # function argument token url override
                    codecov_url = construct_uri_string_ci(token="token_name_2")
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("service=jenkins", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("build_url=t_url", codecov_url)
                    @test occursin("build=t_num", codecov_url)
                end
            end
        end

        # test codecov token masking
        withenv(
            "APPVEYOR" => "true",
            "APPVEYOR_PULL_REQUEST_NUMBER" => "t_pr",
            "APPVEYOR_ACCOUNT_NAME" => "t_account",
            "APPVEYOR_PROJECT_SLUG" => "t_slug",
            "APPVEYOR_BUILD_VERSION" => "t_version",
            "APPVEYOR_REPO_BRANCH" => "t_branch",
            "APPVEYOR_REPO_COMMIT" => "t_commit",
            "APPVEYOR_REPO_NAME" => "t_repo",
            "APPVEYOR_JOB_ID" => "t_job_num",
            "CODECOV_URL" => "https://enterprise-codecov-1.com",
            "CODECOV_TOKEN" => "token_name_1"
            ) do
                codecov_url = construct_uri_string_ci()
                masked = Coverage.Codecov.mask_token(codecov_url)
                @test !occursin("token_name_1", masked)
                @test occursin("token=<HIDDEN>", masked)
        end

        # in the case above, the token is at the end. Let's test explicitly,
        # that this also works if the token occurs earlier in the url
        url = "https://enterprise-codecov-1.com/upload/v2?token=token_name_1&build=t_job_num"
        masked = Coverage.Codecov.mask_token(url)
        @test masked == "https://enterprise-codecov-1.com/upload/v2?token=<HIDDEN>&build=t_job_num"
    end
end

@testset "coveralls" begin
    # NOTE: this only returns actual content if this package is devved.
    # Hence the test is basically on this function returning something
    # (rather than erroring)
    git = Coverage.Coveralls.query_git_info()
    @test git["remotes"][1]["name"] == "origin"
    @test haskey(git["remotes"][1], "url")

    # for testing submit_***()
    fcs = FileCoverage[]

    # an error should be raised if there is no coveralls token set
    withenv("COVERALLS_TOKEN" => nothing,
            "REPO_TOKEN" => nothing,  # this is deprecrated, use COVERALLS_TOKEN
            "APPVEYOR" => "true",  # use APPVEYOR as an example to make the test reach the repo token check
            "APPVEYOR_JOB_ID" => "my_job_id") do
            @test_throws ErrorException Coverage.Coveralls.prepare_request(fcs, false)
    end

    withenv("COVERALLS_TOKEN" => "token_name_1",
            "APPVEYOR" => nothing,
            "APPVEYOR_JOB_ID" => nothing,
            "TRAVIS" => nothing,
            "TRAVIS_JOB_ID" => nothing,
            "JENKINS" => nothing,
            "BUILD_ID" => nothing,
            "CI_PULL_REQUEST" => nothing,
            "GIT_BRANCH" => nothing) do

        # test error if not local and no CI platform can be detected from ENV
        @test_throws ErrorException Coverage.Coveralls.prepare_request(fcs, false)

        # test local submission, when we are local
        _dotgit = joinpath(dirname(@__DIR__), ".git")
        if isdir(_dotgit) || isfile(_dotgit)
                request = Coverage.Coveralls.prepare_request(fcs, true)
                @test request["repo_token"] == "token_name_1"
                @test isempty(request["source_files"])
                @test haskey(request, "git")
                @test request["git"]["remotes"][1]["name"] == "origin"
        end

        # test APPVEYOR
        withenv("APPVEYOR" => "true",
                "APPVEYOR_JOB_ID" => "my_job_id") do
                request = Coverage.Coveralls.prepare_request(fcs, false)
                @test request["repo_token"] == "token_name_1"
                @test request["service_job_id"] == "my_job_id"
                @test request["service_name"] == "appveyor"
        end

        # test Travis
        withenv("TRAVIS" => "true",
                "TRAVIS_JOB_ID" => "my_job_id") do
                request = Coverage.Coveralls.prepare_request(fcs, false)
                @test request["repo_token"] == "token_name_1"
                @test request["service_job_id"] == "my_job_id"
                @test request["service_name"] == "travis-ci"
        end

        # test Jenkins
        withenv("JENKINS" => "true",
                "BUILD_ID" => "my_job_id",
                "CI_PULL_REQUEST" => true) do
                my_git_info = Dict("remote_name" => "my_origin")
                request = Coverage.Coveralls.prepare_request(fcs, false)
                @test request["repo_token"] == "token_name_1"
                @test request["service_job_id"] == "my_job_id"
                @test request["service_name"] == "jenkins-ci"

                withenv("CI_PULL_REQUEST" => "false",
                        "GIT_BRANCH" => "my_remote/my_branch") do
                        request = Coverage.Coveralls.prepare_request(fcs, false)
                        @test haskey(request, "git")
                        @test request["git"]["branch"] == "my_branch"
                end
        end

        # test git_info (only works with Jenkins & local at the moment)
        withenv("JENKINS" => "true",
                "BUILD_ID" => "my_job_id",
                "CI_PULL_REQUEST" => true) do
                # we can pass in our own function, that returns a dict
                my_git_info() = Dict("test" => "test")
                request = Coverage.Coveralls.prepare_request(fcs, false, my_git_info)
                @test haskey(request, "git")

                # or directly a dict
                request = Coverage.Coveralls.prepare_request(fcs, false, Dict("test" => "test"))
                @test haskey(request, "git")
            end
    end
end

end # of withenv("DISABLE_AMEND_COVERAGE_FROM_SRC" => nothing)

end # of @testset "Coverage"


text:
"""
    SerialExecution

Functor object for using serial execution in a progressive-hedging algorithm. Create by supplying a [`Serial`](@ref) object through `execution` in the `ProgressiveHedgingSolver` factory function and then pass to a `StochasticPrograms.jl` model.

"""
struct SerialExecution{T <: AbstractFloat,
                       A <: AbstractVector,
                       S <: LQSolver,
                       PT <: PenaltyTerm} <: AbstractExecution
    subproblems::Vector{SubProblem{T,A,S,PT}}

    function SerialExecution(::Type{T}, ::Type{A}, ::Type{S}, ::Type{PT}) where {T <: AbstractFloat, A <: AbstractVector, S <: LQSolver, PT <: PenaltyTerm}
        return new{T,A,S,PT}(Vector{SubProblem{T,A,S,PT}}())
    end
end

function initialize_subproblems!(ph::AbstractProgressiveHedgingSolver, subsolver::QPSolver, penaltyterm::PenaltyTerm, execution::SerialExecution)
    for i = 1:ph.nscenarios
        push!(execution.subproblems,SubProblem(WS(ph.stochasticprogram, scenario(ph.stochasticprogram,i); solver = subsolver),
                                               i,
                                               probability(ph.stochasticprogram,i),
                                               decision_length(ph.stochasticprogram),
                                               subsolver,
                                               copy(penaltyterm)))
    end
    update_iterate!(ph)
    return ph
end

function resolve_subproblems!(ph::AbstractProgressiveHedgingSolver, execution::SerialExecution{T,A}) where {T <: AbstractFloat, A <: AbstractVector}
    Qs = A(undef, length(execution.subproblems))
    # Update subproblems
    reformulate_subproblems!(execution.subproblems, ph.ξ, penalty(ph))
    # Solve sub problems
    for (i, subproblem) ∈ enumerate(execution.subproblems)
        Qs[i] = subproblem(ph.ξ)
    end
    # Return current objective value
    return sum(Qs)
end

function update_iterate!(ph::AbstractProgressiveHedgingSolver, execution::SerialExecution)
    # Update the estimate
    ξ_prev = copy(ph.ξ)
    ph.ξ[:] = sum([subproblem.π*subproblem.x for subproblem in execution.subproblems])
    # Update δ₁
    ph.data.δ₁ = norm(ph.ξ-ξ_prev, 2)^2
    return nothing
end

function update_subproblems!(ph::AbstractProgressiveHedgingSolver, execution::SerialExecution)
    # Update dual prices
    update_subproblems!(execution.subproblems, ph.ξ, penalty(ph))
    return nothing
end

function update_dual_gap!(ph::AbstractProgressiveHedgingSolver, execution::SerialExecution)
    # Update δ₂
    ph.data.δ₂ = sum([s.π*norm(s.x-ph.ξ,2)^2 for s in execution.subproblems])
    return nothing
end

function calculate_objective_value(ph::AbstractProgressiveHedgingSolver, execution::SerialExecution)
    return sum([get_objective_value(s) for s in execution.subproblems])
end

function fill_first_stage!(ph::AbstractProgressiveHedgingSolver, stochasticprogram::StochasticProgram, nrows::Integer, ncols::Integer, execution::SerialExecution)
    subproblem = execution.subproblems[1]
    μ = try
        MPB.getreducedcosts(subproblem.solver.lqmodel)[1:ncols]
    catch
        fill(NaN, ncols)
    end
    StochasticPrograms.set_first_stage_redcosts!(stochasticprogram, μ)
    λ = try
        MPB.getconstrduals(subproblem.solver.lqmodel)[1:nrows]
    catch
        fill(NaN, nrows)
    end
    StochasticPrograms.set_first_stage_duals!(stochasticprogram, λ)
    return nothing
end

function fill_submodels!(ph::AbstractProgressiveHedgingSolver, scenarioproblems::StochasticPrograms.ScenarioProblems, nrows::Integer, ncols::Integer, execution::SerialExecution)
    for (i, submodel) in enumerate(scenarioproblems.problems)
        fill_submodel!(submodel, execution.subproblems[i], nrows, ncols)
    end
    return nothing
end

function fill_submodels!(ph::AbstractProgressiveHedgingSolver, scenarioproblems::StochasticPrograms.DScenarioProblems, nrows::Integer, ncols::Integer, execution::SerialExecution)
    j = 0
    @sync begin
        for w in workers()
            n = remotecall_fetch((sp)->length(fetch(sp).problems), w, scenarioproblems[w-1])
            for i in 1:n
                k = i+j
                @async remotecall_fetch((sp,i,x,μ,λ,nrows,ncols) -> fill_submodel!(fetch(sp).problems[i],x,μ,λ,nrows,ncols),
                                        w,
                                        scenarioproblems[w-1],
                                        i,
                                        get_solution(execution.subproblems[k])...,
                                        nrows,
                                        ncols)
            end
            j += n
        end
    end
    return nothing
end
# API
# ------------------------------------------------------------
function (execution::Serial)(::Type{T}, ::Type{A}, ::Type{S}, ::Type{PT}) where {T <: AbstractFloat, A <: AbstractVector, S <: LQSolver, PT <: PenaltyTerm}
    return SerialExecution(T,A,S,PT)
end

function str(::Serial)
    return "serial execution"
end


text:
# Gleb: the code taken from the ode elimination project, tests are there

using MacroTools
using Oscar

include("power_series_utils.jl")

# P is the type of polynomials in the rhs of the ODE system
struct ODE{P}
    poly_ring::MPolyRing
    x_vars::Array{P, 1}
    u_vars::Array{P, 1}
    parameters::Array{P, 1}
    equations::Dict{P, <: Union{P, Generic.Frac{P}}}

    function ODE{P}(eqs::Dict{P, <: Union{P, Generic.Frac{P}}}, inputs::Array{P, 1}) where {P <: MPolyElem{<: FieldElem}}
        #Initialize ODE
        #equations is a dictionary x_i => f_i(x, u, params)

        num, den = unpack_fraction(collect(values(eqs))[1])
        poly_ring = parent(num)
        x_vars = collect(keys(eqs))
        u_vars = inputs
        parameters = filter(v -> (!(v in x_vars) && !(v in u_vars)), gens(poly_ring))
        new(poly_ring, x_vars, u_vars, parameters, eqs)
    end

    function ODE{P}(eqs::Dict{P, <: Union{P, Generic.Frac{P}}}) where {P <: MPolyElem{<: FieldElem}}
        return ODE{P}(eqs, Array{P, 1}())
    end
end

#------------------------------------------------------------------------------

function power_series_solution(
        ode::ODE{P},
        param_values::Dict{P, T},
        initial_conditions::Dict{P, T},
        input_values::Dict{P, Array{T, 1}},
        prec::Int
    ) where {T <: FieldElem, P <: MPolyElem{T}}
    """
    Input:
        - ode, an ode to solve
        - param_values, initial_conditions - parameter values and initial conditions to plug in
          both are dictionaries variable => value
        - input_values - power series for the inpiuts presented as a dictionary
          variable => list of coefficients
        - prec, the precision
    Output: computes a power series solution with precision prec presented as a dictionary
            variable => corresponding coordiante of the solution
    """
    new_varnames = map(var_to_str, vcat(ode.x_vars, ode.u_vars))
    append!(new_varnames, map(v -> var_to_str(v) * "_dot", ode.x_vars))

    new_ring, new_vars = PolynomialRing(base_ring(ode.poly_ring), new_varnames)
    equations = Array{P, 1}()
    evaluation = Dict(k => new_ring(v) for (k, v) in param_values)
    for v in vcat(ode.x_vars, ode.u_vars)
        evaluation[v] = switch_ring(v, new_ring)
    end
    for (v, eq) in ode.equations
        num, den = map(p -> eval_at_dict(p, evaluation), unpack_fraction(eq))
        push!(equations, den * str_to_var(var_to_str(v) * "_dot", new_ring) - num)
    end
    new_inputs = Dict(switch_ring(k, new_ring) => v for (k, v) in input_values)
    new_ic = Dict(switch_ring(k, new_ring) => v for (k, v) in initial_conditions)
    result = ps_ode_solution(equations, new_ic, new_inputs, prec)
    return Dict(v => result[switch_ring(v, new_ring)] for v in vcat(ode.x_vars, ode.u_vars))
end

#------------------------------------------------------------------------------

function power_series_solution(
        ode::ODE{P},
        param_values::Dict{P, Int},
        initial_conditions::Dict{P, Int},
        input_values::Dict{P, Array{Int, 1}},
        prec::Int
    ) where P <: MPolyElem{<: FieldElem}
    bring = base_ring(ode.poly_ring)
    return power_series_solution(
        ode,
        Dict(p => bring(v) for (p, v) in param_values),
        Dict(x => bring(v) for (x, v) in initial_conditions),
        Dict(u => map(v -> bring(v), vv) for (u, vv) in input_values),
        prec
    )
end
#------------------------------------------------------------------------------

function _reduce_poly_mod_p(poly::MPolyElem{Nemo.fmpq}, p::Int)
    """
    Reduces a polynomial over Q modulo p
    """
    den = denominator(poly)
    num = change_base_ring(ZZ, den * poly)
    if GF(p)(den) == 0
        throw(Base.ArgumentError("Prime $p divides the denominator of $poly"))
    end
    return change_base_ring(GF(p), num) * (1 // GF(p)(den))
end

#--------------------------------------

function reduce_ode_mod_p(ode::ODE{<: MPolyElem{Nemo.fmpq}}, p::Int)
    """
    Input: ode is an ODE over QQ, p is a prime number
    Output: the reduction mod p, throws an exception if p divides one of the denominators
    """
    new_ring, new_vars = PolynomialRing(GF(p), map(var_to_str, gens(ode.poly_ring)))
    new_type = typeof(new_vars[1])
    new_inputs = map(u -> switch_ring(u, new_ring), ode.u_vars)
    new_eqs = Dict{new_type, Union{new_type, Generic.Frac{new_type}}}()
    for (v, f) in ode.equations
        new_v = switch_ring(v, new_ring)
        if applicable(numerator, f)
            # if f is a rational function
            num, den = map(poly -> _reduce_poly_mod_p(poly, p), [numerator(f), denominator(f)])
            if den == 0
                throw(Base.ArgumentError("Prime $p divides the denominator of $poly"))
            end
            new_eqs[new_v] = num // den
        else
            new_eqs[new_v] = _reduce_poly_mod_p(f, p)
        end
    end
    return ODE{new_type}(new_eqs, new_inputs)
end

#------------------------------------------------------------------------------

function print_for_SIAN(ode::ODE{P}, outputs::Array{P, 1}) where P <: MPolyElem{<: FieldElem}
    """
    Prints the ODE in the format accepted by SIAN (https://github.com/pogudingleb/SIAN)
    """
    vars_str = Dict(x => var_to_str(x) * "(t)" for x in vcat(ode.x_vars, ode.u_vars))
    merge!(vars_str, Dict(p => var_to_str(p) for p in ode.parameters))
    R_print, vars_print = PolynomialRing(base_ring(ode.poly_ring), [vars_str[v] for v in gens(ode.poly_ring)])
    result = ""

    function _lhs_to_str(lhs)
        num, den = unpack_fraction(lhs)
        result = string(evaluate(num, vars_print))
        if den != 1
             result = "($result) / ($(evaluate(den, vars_print)))"
        end
        return result
    end

    for (x, f) in ode.equations
        result = result * "diff(" * var_to_str(x) * "(t), t) = $(_lhs_to_str(f)), \n"
    end
    for (y_ind, g) in enumerate(outputs)
        result = result * "y_var_$y_ind(t) = $(_lhs_to_str(g)), \n"
    end
    return result
end

#------------------------------------------------------------------------------

function macrohelper_extract_vars(equations::Array{Expr, 1}, summary::Bool=true)
    funcs, x_vars, all_symb = Set(), Set(), Set()
    aux_symb = Set([:(+), :(-), :(=), :(*), :(^), :t, :(/), :(//)])
    for eq in equations
        MacroTools.postwalk(
            x -> begin
                if @capture(x, f_'(t))
                    push!(x_vars, f)
                    push!(all_symb, f)
                elseif @capture(x, f_(t))
                    push!(funcs, f)
                elseif (x isa Symbol) && !(x in aux_symb)
                    push!(all_symb, x)
                end
                return x
            end,
            eq
        )
    end
    u_vars = setdiff(funcs, x_vars)
    params = setdiff(all_symb, union(x_vars, u_vars))
    all_symb = collect(all_symb)
    if summary
        print("Summary of the model:\n")
        print("State variables: ", join(map(string, collect(x_vars)), ", "), "\n")
        print("Parameter: ", join(map(string, collect(params)), ", "), "\n")
        print("Inputs: ", join(map(string, collect(u_vars)), ", "), "\n")
    end
    return collect(u_vars), collect(all_symb)
end

#------------------------------------------------------------------------------

function macrohelper_clean(ex::Expr)
    ex = MacroTools.postwalk(x -> @capture(x, f_'(t)) ? f : x, ex)
    ex = MacroTools.postwalk(x -> @capture(x, f_(t)) ? f : x, ex)
    ex = MacroTools.postwalk(x -> x == :(/) ? :(//) : x, ex)
    return ex
end

#------------------------------------------------------------------------------

macro ODEmodel(ex::Expr...)
    """
    Macros for creating an ODE from a list of equations
    Also injects all variables into the global scope
    """
    extra_params = []
    equations = [ex...]
    if equations[end].head == :vect
        extra_params = equations[end].args
        equations = equations[1:end - 1]
    end

    u_vars, all_symb = macrohelper_extract_vars(equations)
    append!(all_symb, extra_params)

    # creating the polynomial ring
    vars_list = :([$(all_symb...)])
    R = gensym()
    vars_aux = gensym()
    exp_ring = :(($R, $vars_aux) = PolynomialRing(QQ, map(string, $all_symb)))
    assignments = [:($(all_symb[i]) = $vars_aux[$i]) for i in 1:length(all_symb)]

    # preparing equations
    equations = map(macrohelper_clean, equations)
    dict_var = gensym()
    dict_create_expr = :($dict_var = Dict{fmpq_mpoly, Union{fmpq_mpoly, Generic.Frac{fmpq_mpoly}}}())
    eqs_expr = []
    for eq in equations
        if eq.head != :(=)
            throw("Problem with parsing at $eq")
        end
        lhs, rhs = eq.args[1:2]
        loc_all_symb = macrohelper_extract_vars([rhs], false)[2]
        if isempty(loc_all_symb)
            push!(eqs_expr, :($dict_var[$lhs] = $R($rhs)))
        else
            push!(eqs_expr, :($dict_var[$lhs] = ($rhs)))
        end
    end

    # creating the ode object
    ode_expr = :(ODE{fmpq_mpoly}($dict_var, Array{fmpq_mpoly}([$(u_vars...)])))

    result = Expr(
        :block,
        exp_ring, assignments...,
        dict_create_expr, eqs_expr...,
        ode_expr
    )
    return esc(result)
end

#------------------------------------------------------------------------------

function generate_replica(ode::ODE{P}, outputs::Array{P, 1}, r::Int) where P <: MPolyElem
    """
    Returns (ode_r, outputs_r), and r-fold replica of the original pair (ode, outputs).
    States, outputs, and inputs are replicated, parameters are not
    """
    new_varnames = map(string, ode.parameters)
    for v in vcat(ode.x_vars, ode.u_vars)
        append!(new_varnames, [var_to_str(v) * "_r$i" for i in 1:r])
    end
    new_ring, new_vars = PolynomialRing(base_ring(ode.poly_ring), new_varnames)
    new_eqs = Dict{P, Union{P, Generic.Frac{P}}}()
    new_outputs = Array{P, 1}()
    new_us = Array{P, 1}()
    for i in 1:r
        eval = merge(
            Dict(p => switch_ring(p, new_ring) for p in ode.parameters),
            Dict(v => str_to_var(var_to_str(v) * "_r$i", new_ring) for v in vcat(ode.x_vars, ode.u_vars))
        )
        eval_vec = [eval[v] for v in gens(ode.poly_ring)]
        append!(new_outputs, map(p -> evaluate(p, eval_vec), outputs))
        new_eqs = merge(
            new_eqs,
            Dict{P, Union{P, Generic.Frac{P}}}(evaluate(x, eval_vec) => evaluate(f, eval_vec) for (x, f) in ode.equations)
        )
        append!(new_us, [str_to_var(var_to_str(u) * "_r$i", new_ring) for u in ode.u_vars])
    end
    return (ODE{P}(new_eqs, new_us), new_outputs)
end

#------------------------------------------------------------------------------

max_stars_repo_name: pogudingleb/ExperimentsBound

max_issues_repo_name: pogudingleb/ExperimentsBound

max_forks_repo_name: pogudingleb/ExperimentsBound


text:
module QuadraticToBinary

using MathOptInterface
using DataStructures

include("moi.jl")

end # module


text:
using FirthRegression, Base.Test, DataFrames

## Firth regression with admit data
df = readtable(Pkg.dir("GLM","data","admit.csv.gz"))
df[:rank] = pool(df[:rank])
gm = fit(FirthGLM, admit ~ gre + gpa + rank, df, Binomial(), convTol=1e-14)
@test_approx_eq_eps coef(gm) [-3.9076341278382402677, 0.0022173644223770295, 0.7875352644980999628, -0.6657797116943170446, -1.3177582553490327921, -1.5111669280800947845] 1e-8
@test_approx_eq penalized_deviance(gm.model) 431.06681398743257

## With sparse matrices
srand(1)
X = sprand(1000, 10, 0.01)
β = randn(10)
y = Bool[rand() < x for x in logistic(X * β)]

gmdense = fit(FirthGLM, full(X), y, Binomial())
gmsparse = fit(FirthGLM, X, y, Binomial())
@test_approx_eq deviance(gmsparse) deviance(gmdense)
@test_approx_eq coef(gmsparse) coef(gmdense)

## ANOVA
anv = anova(fit(FirthGLM, admit ~ gre + gpa + rank, df, Binomial(), convTol=1e-14).model,
            fit(FirthGLM, admit ~ gre + gpa, df, Binomial(), convTol=1e-14).model, convTol=1e-14)
@test_approx_eq anv.chisq 21.373030365595071
@test anv.df == 3
@test_approx_eq anv.p 8.8071001596912168e-05

## ANOVA w/o prefit models
anv = anova(fit(FirthGLM, admit ~ gre + gpa + rank, df, Binomial(), convTol=1e-14, dofit=false).model,
            fit(FirthGLM, admit ~ gre + gpa, df, Binomial(), convTol=1e-14, dofit=false).model, convTol=1e-14)
@test_approx_eq anv.chisq 21.373030365595071
@test anv.df == 3
@test_approx_eq anv.p 8.8071001596912168e-05

## ANOVA with sparse matrix
anvdense = anova(fit(FirthGLM, full(X[:, 1:5]), y, Binomial(), convTol=1e-14),
                 gmdense, convTol=1e-14)
anvsparse = anova(fit(FirthGLM, X[:, 1:5], y, Binomial(), convTol=1e-14),
                  gmsparse, convTol=1e-14)
@test_approx_eq anvdense.chisq anvsparse.chisq
@test anvdense.df == anvsparse.df
@test_approx_eq anvdense.p anvsparse.p


text:
using Test
import Model_Household: Household, Person, People_Dict, default_bu_allocation
import FRS_Household_Getter
import Example_Household_Getter
using Definitions
import Dates: Date
import IncomeTaxCalculations: old_enough_for_mca, apply_allowance, calc_income_tax
import STBParameters: IncomeTaxSys, get_default_it_system

const RUK_PERSON = 100000001001
const SCOT_HEAD = 100000001002
const SCOT_SPOUSE = 100000001003

function get_tax(; scotland = false ) :: IncomeTaxSys
    it = get_default_it_system( year=2019, scotland=scotland, weekly=false )
    it.non_savings_rates ./= 100.0
    it.savings_rates ./= 100.0
    it.dividend_rates ./= 100.0
    it.personal_allowance_withdrawal_rate /= 100.0
    it.mca_credit_rate /= 100.0
    it.mca_withdrawal_rate /= 100.0
    it
end

@testset "Melville 2019 ch2 examples 1; basic calc Scotland vs RUK" begin
    # BASIC IT Calcaulation on
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    @time names = Example_Household_Getter.initialise()
    income = [11_730,14_493,30_000,33_150.0,58_600,231_400]
    ntests = size(income)[1]
    for i in 1:ntests-1
        income[i] += itsys_scot.personal_allowance # weird way this is expessed in Melville
    end
    taxes_ruk = [2_346.0,2898.60,6_000,6_630.0,15_940.00,89_130.0]
    taxes_scotland = [2_325.51,2_898.60,6_155.07, 7260.57,17_695.07,92_613.07]
    @test size( income ) == size( taxes_ruk) == size( taxes_scotland )
    ruk = Example_Household_Getter.get_household( "mel_c2" )
    scottish = Example_Household_Getter.get_household( "mel_c2_scot" )

    bus = default_bu_allocation( scottish )
    nbus = size(bus)[1]
    println( bus )
    # @test nbus == 1 == size( bus[1])[1]
    pers = bus[1][1]
    intermediate = Dict()
    for i in 1:ntests
        pers.income[wages] = income[i]
        due = calc_income_tax( pers, itsys_scot, intermediate ).total_tax
        println( "Scotland $i : calculated $due expected $(taxes_scotland[i])")
        @test due ≈ taxes_scotland[i]
        due = calc_income_tax( pers, itsys_ruk, intermediate ).total_tax
        println( "rUK $i : calculated $due expected $(taxes_ruk[i])")
        @test due ≈ taxes_ruk[i]
        println( ruk.people[RUK_PERSON].income )
    end
end # example 1

@testset "Melville 2019 ch2 example 2; personal savings allowance" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    names = Example_Household_Getter.initialise()
    ruk = Example_Household_Getter.get_household( "mel_c2" )
    scottish = Example_Household_Getter.get_household( "mel_c2_scot" )
    income = [20_000,37_501,64_000,375_000.0]
    psa = [1_000.0,500.0,500.0,0.0]
    @test size( income ) == size( psa ) # same Scotland and RUK
    pers = ruk.people[RUK_PERSON] # doesn't matter S/RUK
    intermediate = Dict()
    for i in size(income)[1]
        pers.income[wages] = income[i]
        println( "case $i income = $(income[i])")
        due = calc_income_tax( pers, itsys_ruk, intermediate ).total_tax
        @test intermediate["personal_savings_allowance"] == psa[i]
    end
end # example 2

@testset "ch2 example 3; savings calc" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    names = Example_Household_Getter.initialise()
    ruk = Example_Household_Getter.get_household( "mel_c2" )
    pers = ruk.people[RUK_PERSON]
    pers.income[self_employment_income] = 40_000.00
    pers.income[bank_interest] = 1_250.00
    intermediate = Dict()
    tax_due_scotland = 5680.07
    due = calc_income_tax( pers, itsys_scot, intermediate ).total_tax
    println( intermediate )
    @test due ≈ tax_due_scotland
    due = calc_income_tax( pers, itsys_ruk, intermediate ).total_tax
    println( intermediate )
    tax_due_ruk = 5_550.00
    @test due ≈ tax_due_ruk
end # example 3

@testset "ch2 example 4; savings calc" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    intermediate = Dict()
    names = Example_Household_Getter.initialise()
    ruk = Example_Household_Getter.get_household( "mel_c2" )
    pers = ruk.people[RUK_PERSON]
    pers.income[property] = 16_700.00
    pers.income[bank_interest] = 1_100.00
    tax_due_ruk = 840.00
    tax_due_scotland = 819.51
    due = calc_income_tax( pers, itsys_scot, intermediate ).total_tax
    println( intermediate )
    @test due ≈ tax_due_scotland
    due = calc_income_tax( pers, itsys_ruk, intermediate ).total_tax
    println( intermediate )
    @test due ≈ tax_due_ruk
end # example 4

@testset "ch2 example 5; savings calc" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    intermediate = Dict()
    names = Example_Household_Getter.initialise()
    ruk = Example_Household_Getter.get_household( "mel_c2" )
    pers = ruk.people[RUK_PERSON]
    pers.income[self_employment_income] = 58_850.00
    pers.income[bank_interest] = 980.00
    tax_due_ruk = 11_232.00
    tax_due_scotland = 12_864.57
    due = calc_income_tax( pers, itsys_scot, intermediate ).total_tax
    println( intermediate )
    @test due ≈ tax_due_scotland
    due = calc_income_tax( pers, itsys_ruk, intermediate ).total_tax
    println( intermediate )
    @test due ≈ tax_due_ruk
end # example 5

@testset "ch2 example 6; savings calc" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    intermediate = Dict()
    names = Example_Household_Getter.initialise()
    ruk = Example_Household_Getter.get_household( "mel_c2" )
    pers = ruk.people[RUK_PERSON]
    pers.income[self_employment_income] = 240_235.00
    pers.income[bank_interest] = 1_600.00

    tax_due_ruk = 93_825.75
    tax_due_scotland = 97_397.17
    due = calc_income_tax( pers, itsys_scot, intermediate ).total_tax
    println( intermediate )
    @test due ≈ tax_due_scotland
    due = calc_income_tax( pers, itsys_ruk, intermediate ).total_tax
    println( intermediate )
    @test due ≈ tax_due_ruk
end # example 6

@testset "ch2 example 7; savings calc" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    intermediate = Dict()
    names = Example_Household_Getter.initialise()
    ruk = Example_Household_Getter.get_household( "mel_c2" )
    pers = ruk.people[RUK_PERSON]
    pers.income[self_employment_income] = 10_000.00
    pers.income[bank_interest] = 3_380.00
    pers.income[other_investment_income] = 36_680.00/0.8 # gross up at basic
    tax_due_ruk = 10_092.00 # inc already deducted at source
    tax_due_scotland = 10_092.00
    due = calc_income_tax( pers, itsys_scot, intermediate ).total_tax
    println( intermediate )
    @test due ≈ tax_due_scotland
    due = calc_income_tax( pers, itsys_ruk, intermediate ).total_tax
    println( intermediate )
    @test due ≈ tax_due_ruk
end # example 7

#
# stocks_shares
#

@testset "ch2 example 8; simple stocks_shares" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    intermediate = Dict()
    names = Example_Household_Getter.initialise()
    ruk = Example_Household_Getter.get_household( "mel_c2" )
    pers = ruk.people[RUK_PERSON]
    pers.income[property] = 28_590.00
    pers.income[bank_interest] = 1_050.00
    pers.income[stocks_shares] = 204_100.0 # gross up at basic
    tax_due_ruk = 74_834.94 # inc already deducted at source
    tax_due_scotland = 74_834.94+140.97
    due = calc_income_tax( pers, itsys_scot, intermediate ).total_tax
    println( intermediate )
    @test due ≈ tax_due_scotland
    due = calc_income_tax( pers, itsys_ruk, intermediate ).total_tax
    println( intermediate )
    @test due ≈ tax_due_ruk
end # example 8

@testset "ch2 example 9; simple stocks_shares" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    intermediate = Dict()
    names = Example_Household_Getter.initialise()
    ruk = Example_Household_Getter.get_household( "mel_c2" )
    pers = ruk.people[RUK_PERSON]
    pers.income[private_pensions] = 17_750.00
    pers.income[bank_interest] = 195.00
    pers.income[stocks_shares] = 1_600.0 # gross up at basic
    tax_due_ruk = 1_050.00 # inc already deducted at source
    tax_due_scotland = 1_050.00-20.49
    due = calc_income_tax( pers, itsys_scot, intermediate ).total_tax
    println( intermediate )
    @test due ≈ tax_due_scotland
    due = calc_income_tax( pers, itsys_ruk, intermediate ).total_tax
    println( intermediate )
    @test due ≈ tax_due_ruk
end # example 9

@testset "ch3 personal allowances ex 1 - hr allowance withdrawal" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    intermediate = Dict()
    names = Example_Household_Getter.initialise()
    ruk = Example_Household_Getter.get_household( "mel_c2" )
    pers = ruk.people[RUK_PERSON]
    pers.income[self_employment_income] = 110_520.00
    tax_due_ruk = 33_812.00
    due = calc_income_tax( pers, itsys_ruk, intermediate ).total_tax
    println( intermediate )
    @test due ≈ tax_due_ruk
    pers.income[self_employment_income] += 100.0
    due = calc_income_tax( pers, itsys_ruk, intermediate ).total_tax
    println( intermediate )
    tax_due_ruk = 33_812.00+60.0
    @test due ≈ tax_due_ruk

    # tax_due_scotland = 33_812.00+61.5 ## FIXME actually, check this by hand

end # example1 ch3

@testset "ch3 personal allowances ex 2 - marriage allowance" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    intermediate = Dict()
    names = Example_Household_Getter.initialise()
    names = Example_Household_Getter.initialise()
    scot = Example_Household_Getter.get_household( "mel_c2_scot" ) # scots are a married couple
    head = scot.people[SCOT_HEAD]
    spouse = scot.people[SCOT_SPOUSE]
    head.income[self_employment_income] = 11_290.0
    head_tax_due = 0.0
    spouse.income[self_employment_income] = 20_000.0
    spouse_tax_due_ruk = 1_258.0

    result = calc_income_tax( head,spouse, itsys_ruk, intermediate )

    println( result )
    println( intermediate )
    @test result.spouse.total_tax ≈ spouse_tax_due_ruk
end # example 2 ch3

@testset "ch3 blind person" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    intermediate = Dict()
    names = Example_Household_Getter.initialise()
    ruk = Example_Household_Getter.get_household( "mel_c2" )
    pers = ruk.people[RUK_PERSON]
    pers.registered_blind = true
    result = calc_income_tax( pers, nothing, itsys_ruk, intermediate )
    @test result.head.allowance == itsys_scot.personal_allowance + itsys_scot.blind_persons_allowance
    # test that tax is 2450xmr
end

@testset "ch3 tax reducers" begin
    # check that mca is always 10% of amount
    # check marriage transfer is always basic rate tax credit
    # checl MCA only available if 1 spouse born before 6th April  1935
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    intermediate = Dict()
    names = Example_Household_Getter.initialise()
    scot = Example_Household_Getter.get_household( "mel_c2_scot" ) # scots are a married couple
    head = scot.people[SCOT_HEAD]
    spouse = scot.people[SCOT_SPOUSE]
    head_ages = [75,90,90,70] # after 1935
    spouse_ages = [90,70,70,90]
    head_incomes = [19_100.0, 29_710.0,41_080.0,0.0]
    spouse_incomes = [12_450.0,0,13_950.0,49_300.0]
    for i in 1:4
        head.income[private_pensions] = head_incomes[i]
        spouse.income[private_pensions] = spouse_incomes[i]
        # head.income[private_pension] = head_incomes[i]
        head.age = head_ages[i]
        spouse.age = spouse_ages[i]
        result_ruk = calc_income_tax( head, spouse, itsys_ruk, intermediate )
        result_scot = calc_income_tax( head, spouse, itsys_scot, intermediate )
        if i == 1
            @test result_ruk.head.mca ≈ 891.50 ≈ result_scot.head.mca
            @test result_ruk.spouse.mca ≈ 0 ≈ result_scot.spouse.mca
        elseif i == 2
            @test result_ruk.head.mca ≈ 886.00 ≈ result_scot.head.mca
            @test result_ruk.spouse.mca ≈ 0 ≈ result_scot.spouse.mca
        elseif i == 3
            @test result_ruk.head.mca ≈ 345.00 ≈ result_scot.head.mca
            @test result_ruk.spouse.mca ≈ 0 ≈ result_scot.spouse.mca
        elseif i == 4
            @test result_ruk.spouse.mca ≈ 345.00 ≈ result_scot.spouse.mca
            @test result_ruk.head.mca ≈ 0 ≈ result_scot.head.mca
        end

    end
end

@testset "Crude MCA Age Check" begin
    # cut-off for jan 2010 should be age 85
    d = Date( 2020, 1, 28 )
    @test old_enough_for_mca( 85, d )
    @test ! old_enough_for_mca( 84, d )
    @test old_enough_for_mca( 86, d )
end

@testset "Apply Allowance" begin
    allowance = 10_000
    allowance,t1 = apply_allowance( allowance, 5_000 )
    @test t1 == 0
    @test allowance == 5_000
    allowance,t2 = apply_allowance( allowance, 4_000 )
    @test t2 == 0
    @test allowance == 1_000
    allowance,t3 = apply_allowance( allowance, 4_000 )
    @test t3 == 3_000
    @test allowance == 0
end


text:
using RecSys

import RecSys: train, recommend, rmse

if isless(Base.VERSION, v"0.5.0-")
    using SparseVectors
end

type MovieRec
    movie_names::FileSpec
    als::ALSWR
    movie_mat::Nullable{SparseVector{AbstractString,Int64}}

    function MovieRec(trainingset::FileSpec, movie_names::FileSpec)
        new(movie_names, ALSWR(trainingset, ParShmem()), nothing)
    end

    function MovieRec(trainingset::FileSpec, movie_names::FileSpec, thread::Bool)
	new(movie_names, ALSWR(trainingset, ParThread()), nothing)
    end

    function MovieRec(user_item_ratings::FileSpec, item_user_ratings::FileSpec, movie_names::FileSpec)
        new(movie_names, ALSWR(user_item_ratings, item_user_ratings, ParBlob()), nothing)
    end
end

function movie_names(rec::MovieRec)
    if isnull(rec.movie_mat)
        A = read_input(rec.movie_names)
        movie_ids = convert(Array{Int}, A[:,1])
        movie_names = convert(Array{AbstractString}, A[:,2])
        movie_genres = convert(Array{AbstractString}, A[:,3])
        movies = AbstractString[n*" - "*g for (n,g) in zip(movie_names, movie_genres)]
        M = SparseVector(maximum(movie_ids), movie_ids, movies)
        rec.movie_mat = Nullable(M)
    end

    get(rec.movie_mat)
end

train(movierec::MovieRec, args...) = train(movierec.als, args...)
rmse(movierec::MovieRec, args...; kwargs...) = rmse(movierec.als, args...; kwargs...)
recommend(movierec::MovieRec, args...; kwargs...) = recommend(movierec.als, args...; kwargs...)

function print_list(mat::SparseVector, idxs::Vector{Int}, header::AbstractString)
    if isless(Base.VERSION, v"0.5.0-")
        if !isempty(idxs)
            println(header)
            for idx in idxs
                println("[$idx] $(mat[idx])")
            end
        end
    else
        isempty(idxs) || println("$header\n$(mat[idxs])")
    end
end

function print_recommendations(rec::MovieRec, recommended::Vector{Int}, watched::Vector{Int}, nexcl::Int)
    mnames = movie_names(rec)

    print_list(mnames, watched, "Already watched:")
    (nexcl == 0) || println("Excluded $(nexcl) movies already watched")
    print_list(mnames, recommended, "Recommended:")
    nothing
end

function test_thread(dataset_path)
    ratings_file = DlmFile(joinpath(dataset_path, "ratings.csv"); dlm=',', header=true)
    movies_file = DlmFile(joinpath(dataset_path, "movies.csv"); dlm=',', header=true)
    rec = MovieRec(ratings_file, movies_file, true)

    @time train(rec, 10, 10)
end

function test(dataset_path)
    ratings_file = DlmFile(joinpath(dataset_path, "ratings.csv"); dlm=',', header=true)
    movies_file = DlmFile(joinpath(dataset_path, "movies.csv"); dlm=',', header=true)
    rec = MovieRec(ratings_file, movies_file)
    @time train(rec, 10, 10)

    err = rmse(rec)
    println("rmse of the model: $err")

    println("recommending existing user:")
    print_recommendations(rec, recommend(rec, 100)...)

    println("recommending anonymous user:")
    u_idmap = RecSys.user_idmap(rec.als.inp)
    i_idmap = RecSys.item_idmap(rec.als.inp)
    # take user 100
    actual_user = isempty(u_idmap) ? 100 : findfirst(u_idmap, 100)
    rated_anon, ratings_anon = RecSys.items_and_ratings(rec.als.inp, actual_user)
    actual_movie_ids = isempty(i_idmap) ? rated_anon : i_idmap[rated_anon]
    nmovies = isempty(i_idmap) ? RecSys.nitems(rec.als.inp) : maximum(i_idmap)
    sp_ratings_anon = SparseVector(nmovies, actual_movie_ids, ratings_anon)
    print_recommendations(rec, recommend(rec, sp_ratings_anon)...)

    println("saving model to model.sav")
    clear(rec.als)
    localize!(rec.als)
    save(rec, "model.sav")
    nothing
end

# prepare chunks for movielens dataset by running `split_movielens` from `playground/split_input.jl`
function test_chunks(dataset_path, model_path)
    mem = Base.Sys.free_memory()
    mem_model = mem_inputs = round(Int, mem/3)
    user_item_ratings = SparseBlobs(joinpath(dataset_path, "splits", "R_itemwise"); maxcache=mem_model)
    item_user_ratings = SparseBlobs(joinpath(dataset_path, "splits", "RT_userwise"); maxcache=mem_model)
    movies_file = DlmFile(joinpath(dataset_path, "movies.csv"); dlm=',', header=true)
    rec = MovieRec(user_item_ratings, item_user_ratings, movies_file)
    @time train(rec, 10, 10, model_path, mem_inputs)

    err = rmse(rec)
    println("rmse of the model: $err")

    println("recommending existing user:")
    print_recommendations(rec, recommend(rec, 100)...)
    nothing
end

max_stars_repo_path: Chapter 10/movielens.jl

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_path: Chapter 10/movielens.jl

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_path: Chapter 10/movielens.jl

max_forks_repo_name: PacktPublishing/Julia-for-data-science


text:
max_xyz_inv(width, xmask=0, ymask=0, zmask=0) = 1f0/max(width[1]*xmask, width[2]*ymask , width[3]*zmask)

function grid_translation(scale, model_scale, bb, model, i=1, j=1, k=1)
    translationmatrix(Vec3f0(i-1, j-1, k-1).*scale)*scalematrix(model_scale*scale)*translationmatrix(-minimum(bb))*model
end

function visualize{T <: Composable, N}(grid::Array{T, N}, s::Style, data::Dict)
    @gen_defaults! data begin
        scale    = Vec3f0(1) #axis of 3D dimension, can be signed
    end
    for ind=1:length(grid)
        robj     = grid[ind].children[]
        bb_s     = boundingbox(robj)
        w        = const_lift(widths, bb_s)
        model_scale = const_lift(max_xyz_inv, w, Vec{N, Int}(1)...)
        robj[:model] = const_lift(grid_translation, scale, model_scale, bb_s, robj[:model], ind2sub(size(grid), ind)...) # update transformation matrix
    end
    Context(grid...)
end

function list_translation(lastposition, gap, direction, bb)
    directionmask     = unit(Vec3f0, abs(direction))
    alignmask         = abs(1-directionmask)
    move2align        = (alignmask.*lastposition)-minimum(bb) #zeros direction
    move2nextposition = sign(direction)*(directionmask.*widths(bb))*0.5f0
    nextpos           = lastposition + move2nextposition + (directionmask.*gap)
    translationmatrix(lastposition+move2align),nextpos
end

function visualize{T <: Composable}(list::Vector{T}, s::Style, data::Dict)
    @gen_defaults! data begin
        direction    = 2 #axis of 3D dimension, can be signed
        gap          = 0.1f0*unit(Vec3f0, abs(direction))
        lastposition = Vec3f0(0)
    end
    for elem in list
        transl_nextpos = const_lift(list_translation, lastposition, gap, direction, boundingbox(elem))
        transformation(elem, const_lift(first, transl_nextpos))
        lastposition = const_lift(last, transl_nextpos)
    end
    Context(list...)
end


text:
struct Manifold
    G   # Metric
    H   # Co-metric
    U   # Orthogonal basis      # (d, m, n)
    Σ   # H = U Σ U'
    D   # Diffusion mapping
end

dim(mfd::Manifold) = size(mfd.U, 1)
getm(mfd::Manifold) = size(mfd.U, 2)
getn(mfd::Manifold) = size(mfd.U, 3)

function Base.iterate(mfd::Manifold, state = 1)
    if state > 5
        return nothing
    else
        getfield(mfd, state), state + 1
    end
end

"""
    rmetric(DM::DiffMap, d::Integer, lazy=true)

Estimate the Riemannian metric from data. Set lazy=true if only U is needed.
"""
function rmetric(DM::DiffMap; d::Integer, lazy=true)
    m, n = size(DM)

    # ==== BRUTE FORCE
    # H̃ = zeros(m, m, n)
    # for i in 1:n
    #     for k in 1:m
    #         for ℓ in 1:m
    #             for j in 1:n
    #                 H̃[ℓ, k, i] += DM.L[j, i] * (DM.Y[ℓ, j] - DM.Y[ℓ, i]) * (DM.Y[k, j] - DM.Y[k, i])
    #             end
    #         end
    #     end
    # end

    # ==== Somewhat vectorized but could use improvement
    H̃ = Array{Float64}(undef, m, m, n)

    ydiffcache = Dict{Tuple{Int64, Int64}, Array{Float64}}()
    function ydiff(ℓ, i)
        if !haskey(ydiffcache, (ℓ, i))
            ydiffcache[ℓ, i] = DM.Y[ℓ, :] .- DM.Y[ℓ, i]
        end

        return ydiffcache[ℓ, i]
    end

    for i in 1:n
        for k in 1:m
            for ℓ in k:m
                H̃[ℓ, k, i] = ydiff(ℓ, i)' * (DM.L[:, i] .* ydiff(k, i))
            end
        end
    end
    for k in 1:m
        for ℓ in 1:k-1
            H̃[ℓ, k, :] = H̃[k, ℓ, :]
        end
    end

    U = Array{Float64}(undef, d, m, n)
    Σ = Array{Float64}(undef, d, n)
    H = Array{Float64}(undef, m, m, n)
    G = Array{Float64}(undef, m, m, n)

    for i in 1:n
        # Reduced rank SVD
        tmpU, tmpΣ, _ = svd(H̃[:,:,i])
        U[:,:,i] = permutedims(tmpU[:, 1:d])
        Σ[:,i] = tmpΣ[1:d]

        if !lazy
            H[:,:,i] = U[:,:,i]' * diagm(Σ[:,i]) * U[:,:,i]
            G[:,:,i] = U[:,:,i]' * diagm(Σ[:,i].^-1) * U[:,:,i]
        end
    end

    return Manifold(G, H, U, Σ, DM)
end