text:
map_str = ""
width, height = Int32(0), Int32(0)
open("day03/day03.in") do io
    global map_str, width, height
    for line in readlines(io)
        map_str = map_str * line
        width = length(line)
        height += 1
    end
end

function encounters(right::Int64, down::Int64)::Int64
    posx, posy = Int64(1+right), Int64(1+down)
    tree_count = Int64(0)
    while posy <= height
        global map_str
        if map_str[(posy-1)*width+posx] == '#'
            tree_count += 1
        end
        posx = (posx+right-1) % width + 1
        posy = posy+down

    end
    return tree_count
end

c = encounters(1, 1)
c *= encounters(3, 1)
c *= encounters(5, 1)
c *= encounters(7, 1)
c *= encounters(1, 2)
println(c)

max_stars_repo_name: detrin/advent-of-code

max_issues_repo_name: detrin/advent-of-code

max_forks_repo_name: detrin/advent-of-code


text:
function thermalflowlimits(m::JuMP.AbstractModel, system_formulation::Type{S}, devices::Array{B,1}, time_periods::Int64) where {B <: PowerSystems.Branch, S <: PM.AbstractDCPForm}

    fbr = m[:fbr]
    name_index = m[:fbr].axes[1]
    time_index = m[:fbr].axes[2]

    device_index = Dict(value => key for (key, value) in Dict(collect(enumerate([d.name for d in devices]))))

    Flow_max_tf = JuMP.JuMPArray(Array{ConstraintRef}(undef,length(name_index), time_periods), name_index, time_index)
    Flow_max_ft = JuMP.JuMPArray(Array{ConstraintRef}(undef,length(name_index), time_periods), name_index, time_index)

    for t in time_index, (ix, name) in enumerate(name_index)
        if name in keys(device_index)
            if name == devices[device_index[name]].name
                Flow_max_tf[name, t] = @constraint(m, fbr[name, t] <= devices[device_index[name]].rate)
                Flow_max_ft[name, t] = @constraint(m, fbr[name, t] >= -1*devices[device_index[name]].rate)
            else
                @error "Branch name in Array and variable do not match"
            end
        else
            @warn "No flow limit constraint populated for $(name)"
        end
    end

    JuMP.register_object(m, :Flow_max_ToFrom, Flow_max_tf)
    JuMP.register_object(m, :Flow_max_FromTo, Flow_max_ft)

    return m
end

function thermalflowlimits(m::JuMP.AbstractModel, system_formulation::Type{S}, devices::Array{B,1}, time_periods::Int64) where {B <: PowerSystems.Branch, S <: PM.AbstractDCPLLForm}

    fbr_fr = m[:fbr_fr]
    fbr_to = m[:fbr_to]
    name_index = m[:fbr_fr].axes[1]
    time_index = m[:fbr_to].axes[2]

    device_index = Dict(value => key for (key, value) in Dict(collect(enumerate([d.name for d in devices]))))

    Flow_max_tf = JuMP.JuMPArray(Array{ConstraintRef}(undef,length(name_index), time_periods), name_index, time_index)
    Flow_max_ft = JuMP.JuMPArray(Array{ConstraintRef}(undef,length(name_index), time_periods), name_index, time_index)

    for t in time_index, (ix, name) in enumerate(name_index)
        if name in keys(device_index)
            if name == devices[device_index[name]].name
                Flow_max_tf[name, t] = @constraint(m, fbr_fr[name, t] <= devices[device_index[name]].rate)
                Flow_max_ft[name, t] = @constraint(m, fbr_to[name, t] >= -1*devices[device_index[name]].rate)
            else
                @error "Branch name in Array and variable do not match"
            end
        else
            @warn "No flow limit constraint populated for $(name)"
        end
    end

    JuMP.register_object(m, :Flow_max_ToFrom, Flow_max_tf)
    JuMP.register_object(m, :Flow_max_FromTo, Flow_max_ft)

    return m
end

#TODO: Implement Limits in AC. Use Norm from JuMP Implemented norms.

ext: jl


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

name = "libsass"
version = v"3.5.5"

# Collection of sources required to build SassBuilder
sources = [
    "https://github.com/sass/libsass.git" =>
    "39e30874b9a5dd6a802c20e8b0470ba44eeba929",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libsass
autoreconf --force --install
./configure --prefix=$prefix --host=$target
make -j${nproc}
make install
"""

platforms = [
    # glibc Linuces
    Linux(:i686),
    Linux(:x86_64),
    Linux(:aarch64),
    Linux(:armv7l),
    Linux(:powerpc64le),

    # musl Linuces
    Linux(:i686, :musl),
    Linux(:x86_64, :musl),
    Linux(:aarch64, :musl),
    Linux(:armv7l, :musl),

    # BSDs
    MacOS(:x86_64),
    FreeBSD(:x86_64),

    # Windows
    Windows(:i686),
    Windows(:x86_64, compiler_abi = CompilerABI(:gcc7))
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libsass", :libsass_so)
]

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

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100
function compute(d::Int, n::Int)::Int
    function mod1(k::Float64)::Float64
        k %= 1
        return k < 0.5 ? k : k - 1
    end

    check(k::Float64)::Bool = (lower && lower_limit ≤ k < upper_limit) || (!lower && lower_limit ≥ k > upper_limit)

    limit = (log(d + 1) - log(d)) / log(10)
    lower_limit, upper_limit = log(d) / log(10) % 1, log(d + 1) / log(10) % 1
    lower = lower_limit < upper_limit
    error = log(2) / log(10)
    x = error % 1
    t, convergent, convergent_error = 0, 1, 1
    while abs(convergent_error) > limit
        x, integer = modf(1 / x)
        convergent, t = convergent * integer + t, convergent
        convergent_error = mod1(convergent * error)
    end
    x, integer = modf(1 / x)
    semi_convergent, semi_convergent_error = t, 1
    while abs(semi_convergent_error) > limit
        semi_convergent += convergent
        semi_convergent_error = mod1(semi_convergent * error)
    end
    differences = (convergent, semi_convergent, convergent + semi_convergent)
    if convergent_error > 0
        convergent_limit = convergent * ceil(lower_limit / convergent_error)
    else
        convergent_limit = semi_convergent * ceil(lower_limit / semi_convergent_error)
    end
    convergents = Vector{Float64}()
    while convergent < convergent_limit
        append!(convergents, convergent)
        convergent, t = integer * convergent + t, convergent
        x, integer = modf(1 / x)
    end
    convergents = convergents[2:end]
    for difference ∈ reverse(convergents)
        while convergent_limit > difference
            x = (convergent_limit - difference) * error % 1
            if check(x)
                convergent_limit -= difference
            else
                break
            end
        end
    end
    while true
        checked = false
        for difference ∈ differences
            if difference > convergent_limit
                continue
            end
            x = (convergent_limit - difference) * error % 1
            if check(x)
                checked = true
                convergent_limit -= difference
                break
            end
        end
        if !checked break end
    end
    while true
        for difference ∈ differences
            x = (convergent_limit + difference) * error % 1
            if check(x)
                n -= 1
                if n == 0
                    return trunc(Int, convergent_limit)
                end
                convergent_limit += difference
                break
            end
        end
    end
end

compute(12, 1)

compute(12, 2)

compute(123, 45)

compute(123, 678910)

@benchmark compute(123, 678910)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using Test

Test.@testset "rngs" begin
    include("test_rng_init.jl")
end


text:
using ADMPS
using ADMPS: model_tensor, tensorfromclassical
using Test

@testset "exampletensor" begin
    β = rand()
    @test model_tensor(Ising(β)) ≈ tensorfromclassical([β -β; -β β])
    @test mag_tensor(Ising(β)) !== nothing
    @test energy_tensor(Ising(β)) !== nothing
end


text:
module Turkie

using AbstractPlotting: Scene, Point2f0
using AbstractPlotting: barplot!, lines!, scatter! # Plotting tools
using AbstractPlotting: Observable, Node, lift, on # Observable tools
using AbstractPlotting: recordframe! # Recording tools
using AbstractPlotting.MakieLayout # Layouting tool
using Colors, ColorSchemes # Colors tools
using KernelDensity # To be able to give a KDE
using OnlineStats # Estimators
using DynamicPPL: VarInfo, Model

export TurkieParams, TurkieCallback

export addIO!, record

include("online_stats_plots.jl")

const std_colors = ColorSchemes.seaborn_colorblind

name(s::Symbol) = string(s)
name(s::OnlineStat) = nameof(typeof(s))

"""
    TurkieCallback(model::DynamicPPL.Model, plots::Series/AbstractVector = )

    ## Keyword arguments
    - `showtrace=true` : Show the trace of the variable
    - `window=0` : Use a window for plotting the trace, 0 will not use a window

"""
TurkieCallBack

struct TurkieCallback
    scene::Scene
    data::Dict{Symbol, MovingWindow}
    axis_dict::Dict
    vars::Dict{Symbol, Any}
    params::Dict{Any, Any}
    iter::Observable{Int64}
end

function TurkieCallback(model::Model, plots::Union{Series, AbstractVector} = [:histkde, Mean(), Variance(), AutoCov(20)]; kwargs...)
    variables = VarInfo(model).metadata
    return TurkieCallback(Dict(Pair.(keys(variables), Ref(plots))),
    Dict(kwargs...))
end

function TurkieCallback(varsdict::Dict; kwargs...)
    return TurkieCallback(varsdict, Dict(kwargs...))
end

function TurkieCallback(vars::Dict, params::Dict)
# Create a scene and a layout
    outer_padding = 5
    scene, layout = layoutscene(outer_padding, resolution = (1200, 700))
    display(scene)

    window = get!(params, :window, 1000)

    n_rows = length(keys(vars))
    n_cols = maximum(length.(values(vars)))
    n_plots = n_rows * n_cols
    iter = Node(0)
    data = Dict{Symbol, MovingWindow}(:iter => MovingWindow(window, Int64))
    obs = Dict{Symbol, Any}()
    axis_dict = Dict()
    for (i, (variable, plots)) in enumerate(vars)
        data[variable] = MovingWindow(window, Float32)
        axis_dict[(variable, :varname)] = layout[i, 1, Left()] = Label(scene, string(variable), textsize = 30)
        axis_dict[(variable, :varname)].padding = (0, 50, 0, 0)
        onlineplot!(scene, layout, axis_dict, plots, iter, data, variable, i)
    end
    on(iter) do i
        if i > 10 # To deal with autolimits a certain number of samples are needed
            for (variable, plots) in vars
                for p in plots
                    autolimits!(axis_dict[(variable, p)])
                end
            end
        end
    end
    MakieLayout.trim!(layout)
    TurkieCallback(scene, data, axis_dict, vars, params, iter)
end

function addIO!(cb::TurkieCallback, io)
    cb.params[:io] = io
end

function (cb::TurkieCallback)(rng, model, sampler, transition, iteration)
    fit!(cb.data[:iter], iteration)
    for (vals, ks) in values(transition.θ)
        for (k, val) in zip(ks, vals)
            if haskey(cb.data, Symbol(k))
                fit!(cb.data[Symbol(k)], Float32(val))
            end
        end
    end
    cb.iter[] += 1
    if haskey(cb.params, :io)
        recordframe!(cb.params[:io])
    end
end

end


text:
struct IndividualState #TODO change to immutable
    health::HealthState
    freedom::FreedomState
    detected::DetectionStatus
    quarantine_level::SafeUInt8
  end

  IndividualState() = IndividualState(
    Healthy,
    Free,
    Undetected,
    0
  )

  show(io::IO, s::IndividualState) = print(io, "(",s.health, ", ", s.freedom, ", ", s.detected, ", ", s.quarantine_level, ")")


text:
@testset "SparseRadialMap Multi-threading test with optimization and evaluation of the resulting map" begin

      X = Matrix([6.831108232125667   6.831108465893829  10.748176564682273  17.220877261555913;
      7.341399771164814   7.341399033871628  11.440770534993137  18.015116463525178;
      6.877822405336048   6.877822863906546  10.793381797612751  17.379545045967195;
      6.971116690210811   6.971117274254750  10.941318835214862  17.571340756941474;
      7.800189373932513   7.800189247976281  11.929084958885008  18.852191934527564;
      6.706165535739543   6.706164938380805  10.596868885136731  17.100665224055952;
      7.162689725595848   7.162690676866085  11.236535720904731  17.691824497167673;
      7.700262639909202   7.700262776963255  11.870709137884843  18.589199307190668;
      6.941836510989174   6.941836045832242  10.822143995520504  17.705463568541791;
      7.284400905037024   7.284402833890172  11.438244192350945  17.713180862215779;
      6.114110264862537   6.114110399157427   9.744642327011279  16.427244631267527;
      6.884587321648968   6.884587349033709  10.880879551471059  17.360443337222691;
      6.964593674795581   6.964592742775720  11.015485574206627  17.379079978320874;
      7.174466645033212   7.174468203234895  11.165089908005559  17.802087293640490;
      7.625148752606352   7.625148374481078  11.784628799983244  18.470370177818385;
      7.878948627898106   7.878948112320454  12.118366730171301  18.770876038639738;
      7.406860218451342   7.406861450644357  11.470902695133232  18.293040279245396;
      7.308212626757207   7.308212927902813  11.381148951581917  17.962547003679237;
      7.519983057000856   7.519981781231936  11.597728250168640  18.407564589808729;
      7.813106349722561   7.813107415959238  11.973288464744240  18.802480702520111]')

      order = [[-1], [1; 1], [-1; 1; 0], [-1; 1; 1; 0]]
      T       = SparseRadialMap(4, order; λ = 0.0)
      Tserial = SparseRadialMap(4, order; λ = 0.0)
      Tthread = SparseRadialMap(4, order; λ = 0.0)

      # Run serial code
      optimize(T, X, nothing; start = 2)

      # Run serial code
      optimize(Tserial, X, nothing; start = 2, P = serial)

      # Run multi-threading code
      @show Threads.nthreads()
      optimize(Tthread, X, nothing; start = 2, P = thread)

      @test norm(T(X[:,1]) -     [6.831108232125667;
                                         0.07966241537360474;
                                        -0.18141931240603526;
                                         -2.0970508045942893])<1e-8

      @test norm(Tserial(X[:,1]) -   [6.831108232125667;
                                             0.07966241537360474;
                                            -0.18141931240603526;
                                             -2.0970508045942893])<1e-8
      @test norm(Tthread(X[:,1]) -   [6.831108232125667;
                                             0.07966241537360474;
                                            -0.18141931240603526;
                                             -2.0970508045942893])<1e-8
end


text:
################################################################################
#  Copyright 2020, Marta Vanin, Tom Van Acker                                  #
################################################################################
# PowerModelsDistributionStateEstimation.jl                                    #
# An extention package of PowerModels(Distribution).jl for Static Power System #
# State Estimation.                                                            #
################################################################################

# using pkgs
using Distributions, GaussianMixtures
using HDF5
using Ipopt
using PowerModels, PowerModelsDistribution
using PowerModelsDistributionStateEstimation
#using SCS #removed while SDP tests are not active
using Test

# pkg const
const _DST = Distributions
const _GMM = GaussianMixtures
const _PMD = PowerModelsDistribution
const _PMDSE = PowerModelsDistributionStateEstimation

#network and feeder from ENWL for tests
ntw, fdr = 4, 2

season     = "summer"
time_step  = 144
elm        = ["load", "pv"]
pfs        = [0.95, 0.90]
rm_transfo = true
rd_lines   = true

# set solvers
ipopt_solver = optimizer_with_attributes(Ipopt.Optimizer,"max_cpu_time"=>300.0,
                                                              "tol"=>1e-9,
                                                              "print_level"=>0)

# scs_solver = optimizer_with_attributes(SCS.Optimizer, "max_iters"=>20000, "eps"=>1e-5,
#                                                             "alpha"=>0.4, "verbose"=>0) #deactivated while SDP tests not active

@testset "PowerModelsDistributionStateEstimation" begin

    include("bad_data.jl")
    include("distributions.jl")
    include("estimation_criteria.jl")
    include("mixed_measurements.jl")
    include("non_exact_forms.jl")
    include("power_flow.jl")
    include("pseudo_measurements.jl")
    include("single_conductor_branches.jl")
    include("utils_and_start_val.jl")
    include("with_errors.jl")

end


text:
module Day20

using AdventOfCode2019

function day20(input::String = readInput(joinpath(@__DIR__, "..", "data", "day20.txt")))
    return solve(input, (true, true))
end

function solve(input, parts)
    solution = Array{Int,1}(undef, sum(parts))
    cm = char_matrix(input)
    imat = int_matrix(cm)
    specialPoints = find_portals(cm)
    portals = create_portal_lookup(specialPoints)
    start = specialPoints["AA"][1]
    finish = specialPoints["ZZ"][1]
    if parts[1]
        solution[1] = flood!(copy(imat), portals, start, finish)
    end
    if parts[2]
        solution[2] = flood_levels(imat, portals, start, finish)
    end
    return solution
end

function char_matrix(input::String)
    m = split.(split(input, "\n"), "")
    nrows = length(m)
    nrows > 1 || throw(ArgumentError("Invalid input"))
    ncols = length(m[1])
    if length(m[end]) <= 1
        nrows -= 1
    end
    for i = 1:nrows
        length(m[i]) == ncols || throw(ArgumentError("Invalid input in row $i"))
    end
    charmat = Array{Char, 2}(undef, nrows, ncols)
    for i = 1:nrows
        for (j, s) in enumerate(m[i])
            charmat[i, j] = m[i][j][1]
        end
    end
    return charmat
end

function int_matrix(charmat::Array{Char,2})
    m, n = size(charmat)
    imat = Array{Int,2}(undef, m, n)
    for i = 1:m
        for j = 1:n
            if charmat[i, j] == '.'  # walkable square
                imat[i, j] = -1
            elseif isuppercase(charmat[i, j])  # letter
                imat[i, j] = -2
            else  # non-walkable square
                imat[i, j] = -3
            end
        end
    end
    return imat
end

function find_portals(charmat::Array{Char,2})
    letterDict = Dict{String,Array{CartesianIndex,1}}()
    letterPos = findall(x -> isuppercase(x), charmat)
    m, n = size(charmat)
    for ci in letterPos
        i, j = ci.I
        if j + 1 <= n && isuppercase(charmat[i, j + 1])
            s = join(charmat[i, j:j+1])
            if !haskey(letterDict, s)
                letterDict[s] = Array{CartesianIndex,1}()
            end
            if j + 2 <= n && charmat[i, j + 2] == '.'
                push!(letterDict[s], CartesianIndex(i, j + 2))
            elseif j - 1 >= 0 && charmat[i, j - 1] == '.'
                push!(letterDict[s], CartesianIndex(i, j - 1))
            end
        end
        if i + 1 <= m && isuppercase(charmat[i + 1, j])
            s = join(charmat[i:i+1, j])
            if !haskey(letterDict, s)
                letterDict[s] = Array{CartesianIndex,1}()
            end
            if i + 2 <= m && charmat[i + 2, j] == '.'
                push!(letterDict[s], CartesianIndex(i + 2, j))
            elseif i - 1 >= 0 && charmat[i - 1, j] == '.'
                push!(letterDict[s], CartesianIndex(i - 1, j))
            end
        end
    end
    return letterDict
end

function create_portal_lookup(letterDict::Dict{String,Array{CartesianIndex,1}})
    lookup = Dict{CartesianIndex,CartesianIndex}()
    for (k, v) in letterDict
        (k == "AA" || k == "ZZ") && continue
        length(v) == 2 || throw(AssertionError("Too many indices for portal $k"))
        lookup[v[1]] = v[2]
        lookup[v[2]] = v[1]
    end
    return lookup
end

function flood!(imat::Array{Int,2}, portals::Dict{CartesianIndex,CartesianIndex}, start::CartesianIndex, finish::CartesianIndex)
    imat[start] = 0
    queue = Array{CartesianIndex,1}()
    push!(queue, start)
    while length(queue) > 0
        current = popfirst!(queue)
        if current == finish
            return imat[current]
        end
        for (k, l) in ((0, 1), (0, -1), (1, 0), (-1, 0))
            if imat[(current.I .+ (k, l))...] == -1
                p = CartesianIndex(current.I .+ (k, l))
                imat[p] = imat[current] + 1
                push!(queue, p)
            end
        end
        if haskey(portals, current)
            p = portals[current]
            if imat[p] == -1
                imat[p] = imat[current] + 1
                push!(queue, p)
            end
        end
    end
end

function is_outer(c::CartesianIndex, m::Int, n::Int)
    (c[1] == 3 || c[1] == m - 2) && return true
    (c[2] == 3 || c[2] == n - 2) && return true
    return false
end

function flood_levels(imat::Array{Int,2}, portals::Dict{CartesianIndex,CartesianIndex}, start::CartesianIndex, finish::CartesianIndex)
    maze = Dict{Int,Array{Int,2}}()
    maze[0] = copy(imat)
    maze[0][start] = 0
    m, n = size(imat)

    queue = Array{Tuple{CartesianIndex,Int},1}()
    level = 0
    push!(queue, (start, level))
    while length(queue) > 0
        current, level = popfirst!(queue)
        if level == 0 && current == finish
            return maze[level][current]
        end
        for (k, l) in ((0, 1), (0, -1), (1, 0), (-1, 0))
            if maze[level][(current.I .+ (k, l))...] == -1
                p = CartesianIndex(current.I .+ (k, l))
                maze[level][p] = maze[level][current] + 1
                push!(queue, (p, level))
            end
        end
        if haskey(portals, current)
            if !is_outer(current, m, n)
                p = portals[current]
                if !haskey(maze, level + 1)
                    maze[level + 1] = copy(imat)
                end
                if maze[level + 1][p] == -1
                    maze[level + 1][p] = maze[level][current] + 1
                    push!(queue, (p, level + 1))
                end
            elseif level > 0  # only inner "portals" allowed for level 0
                p = portals[current]
                if maze[level - 1][p] == -1
                    maze[level - 1][p] = maze[level][current] + 1
                    push!(queue, (p, level - 1))
                end
            end
        end
    end
end

end # module


text:
using DataFrames
using Gadfly
# Read the datasource
include(Pkg.dir("BismarkSummary","src","bismark-report.jl"))
#using BismarkSummary

sampleinfo_path = joinpath(Pkg.dir(), "BismarkSummary","testdata","datasource.tsv")
pipeline_name = nothing
run_number = nothing

# Get filenames of bismark reports

test_bismark_report1 = joinpath(Pkg.dir(), "BismarkSummary","testdata","bismark_report.txt")
test_bismark_report2 = joinpath(Pkg.dir(), "BismarkSummary","testdata","bismark_report1.txt")
test_bismark_report3 = joinpath(Pkg.dir(), "BismarkSummary","testdata","bismark_report2.txt")
bismark_report_filenames = [test_bismark_report1,test_bismark_report2,test_bismark_report3]

bismark_report_filenames
sampleinfo = readtable(sampleinfo_path)
report_dict=parse_bismark_reports( bismark_report_filenames)

methods(parse_bismark_reports)

append_report_info_to_sampleinfo!(sampleinfo,report_dict)
report_dict
sampleinfo
Gadfly.plot(sampleinfo,x="sex",y="Mapping efficiency", Geom.boxplot)


text:
# function Crst_direct(r::Int, s::Int, t::Int)
#     # direct method to calculate Crst
#     # see Products of Laguerre Polynomials
#     # Joseph Gillis and  George Weiss
#     # Mathematics of Computation, Vol. 14, No. 69 (Jan., 1960), pp. 60-63
#     # this function implements equation (7)
#     # this is computationally expensive.
#     # use recurrence relation instead.
#     (r >= 0 && s >= 0 && t >= 0) || throw("negative orders not allowed.")
#     # t needs to be in |r-s| and r+s
#     # otherwise the coef is zero
#     if t < abs(r - s) || t > r + s
#         return 0
#     end

#     p = r + s - t
#     # define nmax and nmin and find them
#     nmax = r
#     if s < nmax
#         nmax =  s
#     end
#     if p < nmax
#         nmax = p
#     end
#     nmin = ceil(Int, p / 2)
#     # use rational numbers for now to keep it exact.
#     C = 0 // 1
#     for n in nmin:nmax
#         C += (2 // 1)^(2 * n) * factorial(r + s - n) / ( factorial(r - n) * factorial(s - n) * factorial(2n - p) * factorial(p - n))
#     end
#     # prefactor
#     C *= (-1 // 2)^p
# end

function Crst(r::Int, s::Int, t::Int)
    # this implements equation (22) and (23) in the same paper.
    if r < 0 || s < 0 || t < 0
        return 0 // 1
    end
    # (r >= 0 && s >= 0 && t >= 0) || throw("negative orders not allowed.")
    # notice that C_{rst} is symmetric with respect to any swapping of indices.
    # take r to be the larer one of r and s
    # the recurrence is given for t+1
    # shift it to calculate for t;
    if r < s
        r, s = s, r # swap.
    end
    # t needs to be in |r-s| and r+s
    # otherwise the coef is zero
    if t < r - s || t > r + s
        return 0 // 1
    end

    sigma = r + s + 1 // 1
    delta = r - s
    if t == delta
        return Rational(binomial(r, s))
    end

    # handle general case.
    f1 = t * (delta^2 - t^2)
    f2 = (2t - 1) * (delta^2 - (t - 1)^2 + t * (sigma - t ))
    f3 = -(t - 1) * (delta^2 - (t - 2)^2 + (4t - 3) * (sigma - t + 1))
    f4 = 2(t - 1) * (t - 2) * (sigma - t + 2)
    return (f2 * Crst(r, s, t - 1) + f3 * Crst(r, s, t - 2) + f4 * Crst(r, s, t - 3) ) / f1
end

function laprodexpand(r::Int, s::Int)
    # expand the product of two Laguerre polynomials (with coefficients being 1) into a linear
    # combination of Laguerre polynomials.
    # Lr(x)Ls(x) into sum_{t}C_{rst}Lt(x)
    orders = abs(r - s):(r + s) |> collect
    coeffs = [Crst(r, s, t) for t in orders]
    return orders, coeffs
end

function Base.:*(l1::LaguerrePolynomial{T}, l2::LaguerrePolynomial{S}) where {T,S}
    # take two linear combinations of Laguerre polynomials.
    # and reexpand into a single Laguerre polynomial series.
    lst = []
    R = promote_type(T, S)
    for (r, c1) in l1
        for (s, c2) in l2
            # construct the product
            orders, coeffs = laprodexpand(r, s)
            # convert to type R
            coeffs = R.(coeffs)
            # add this to the list
            la = LaguerrePolynomial(orders, coeffs)
            la = la * c1 * c2
            push!(lst, la)
        end
    end
    # add elements in the list.
    la = lst[1]
    if length(lst) > 1
        for i = 2:length(lst)
            la += lst[i]
        end
    end
    return la
end


text:
using Flux: convfilter, Zeros, outdims

# Build a VGG block
#  ifilters: number of input filters
#  ofilters: number of output filters
#  batchnorm: add batchnorm
function vgg_block(ifilters, ofilters, depth, batchnorm)
  k = (3,3)
  p = (1,1)
  layers = []
  for l in 1:depth
    if batchnorm
      w = convfilter(k, ifilters=>ofilters)
      b = Zeros()
      push!(layers, Conv(weight=w, bias=b, pad=p))
      push!(layers, BatchNorm(ofilters, relu))
    else
      push!(layers, Conv(k, ifilters=>ofilters, relu, pad=p))
    end
    ifilters = ofilters
  end
  return layers
end

# Build convolutionnal layers
#  config: :A (vgg11) :B (vgg13) :D (vgg16) :E (vgg19)
#  inchannels: number of channels in input image (3 for RGB)
function convolutional_layers(config, batchnorm, inchannels)
  layers = []
  ifilters = inchannels
  for c in config
    push!(layers, vgg_block(ifilters, c..., batchnorm)...)
    push!(layers, MaxPool((2,2)))
    ifilters, _ = c
  end
  return layers
end

# Build classification layers
#  imsize: image size (w, h, c)
#  nclasses: number of classes
#  fcsize: size of fully connected layers (usefull for smaller nclasses than ImageNet)
#  dropout: dropout obviously
function classifier_layers(imsize, nclasses, fcsize, dropout)
  layers = []
  push!(layers, flatten)
  push!(layers, Dense(Int(prod(imsize)), fcsize, relu))
  push!(layers, Dropout(dropout))
  push!(layers, Dense(fcsize, fcsize, relu))
  push!(layers, Dropout(dropout))
  push!(layers, Dense(fcsize, nclasses))
  push!(layers, softmax)
  return layers
end

function vgg(imsize; config, inchannels, batchnorm=false, nclasses, fcsize, dropout)
  conv = convolutional_layers(config, batchnorm, inchannels)
  imsize = foldl((insize, layer) -> outdims(layer, insize), conv; init = imsize)
  class = classifier_layers((imsize..., config[end][1]), nclasses, fcsize, dropout)
  return Chain(conv..., class...)
end

const configs = Dict(:A => [(64,1), (128,1), (256,2), (512,2), (512,2)],
                     :B => [(64,2), (128,2), (256,2), (512,2), (512,2)],
                     :D => [(64,2), (128,2), (256,3), (512,3), (512,3)],
                     :E => [(64,2), (128,2), (256,4), (512,4), (512,4)])

vgg11(imsize; inchannels=3, nclasses=1000, fcsize=4096, dropout=0.5) =
  vgg(imsize, config=configs[:A], inchannels=inchannels, nclasses=nclasses, fcsize=fcsize, dropout=dropout)

vgg11bn(imsize; inchannels=3, nclasses=1000, fcsize=4096, dropout=0.5) =
  vgg(imsize, config=configs[:A], batchnorm=true, inchannels=inchannels, nclasses=nclasses, fcsize=fcsize, dropout=dropout)

vgg13(imsize; inchannels=3, nclasses=1000, fcsize=4096, dropout=0.5) =
  vgg(imsize, config=configs[:B], inchannels=inchannels, nclasses=nclasses, fcsize=fcsize, dropout=dropout)

vgg13bn(imsize; inchannels=3, nclasses=1000, fcsize=4096, dropout=0.5) =
  vgg(imsize, config=configs[:B], batchnorm=true, inchannels=inchannels, nclasses=nclasses, fcsize=fcsize, dropout=dropout)

vgg16(imsize; inchannels=3, nclasses=1000, fcsize=4096, dropout=0.5) =
  vgg(imsize, config=configs[:D], inchannels=inchannels, nclasses=nclasses, fcsize=fcsize, dropout=dropout)

vgg16bn(imsize; inchannels=3, nclasses=1000, fcsize=4096, dropout=0.5) =
  vgg(imsize, config=configs[:D], batchnorm=true, inchannels=inchannels, nclasses=nclasses, fcsize=fcsize, dropout=dropout)

vgg19(imsize; inchannels=3, nclasses=1000, fcsize=4096, dropout=0.5) =
  vgg(imsize, config=configs[:E], inchannels=inchannels, nclasses=nclasses, fcsize=fcsize, dropout=dropout)

vgg19bn(imsize; inchannels=3, nclasses=1000, fcsize=4096, dropout=0.5) =
  vgg(imsize, config=configs[:E], batchnorm=true, inchannels=inchannels, nclasses=nclasses, fcsize=fcsize, dropout=dropout)


text:
################################################################################
# CutPruners
# A package to manage polyhedral convex functions
################################################################################

module CutPruners

using LinearAlgebra, SparseArrays
using JuMP

# Redudancy checking
include("redund.jl")

include("abstract.jl")
include("avg.jl")
include("decay.jl")
include("dematos.jl")
include("exact.jl")

end # module


text:
using DrWatson
@quickactivate "FractalDimension"
include(srcdir("style.jl"))
using PyPlot, DataFrames, CSVFiles

# %% computation times for different lengths

data = wload(datadir("benchmark", "benchmark_number_NofN=10_base=10_lowerN=3_upperN=5.csv")) |> DataFrame
fig, axs = subplots(2,1; sharex = true, figsize = (10, 10))
for (i, name) in enumerate(unique(data.Method))
    color = "C$(i-1)"

    lorenz_data = data[(data.Model .== "lorenz") .& (data.Method .== name), :]
    axs[1].plot(lorenz_data.N, lorenz_data.Time, label = name, color = color)

    hénon_data = data[(data.Model .== "henon") .& (data.Method .== name), :]
    axs[2].plot(hénon_data.N, hénon_data.Time, color = color)
end
axs[2].set_xlabel("\$N\$")
axs[2].set_ylabel("\$t_{Hénon}[\\mathrm{ns}]\$")
axs[1].set_ylabel("\$t_{Lorenz63}[\\mathrm{ns}]\$")

axs[1].set_yscale("log")
axs[2].set_yscale("log")
axs[1].set_xscale("log")
axs[2].set_xscale("log")
for ax in axs; ax.legend(loc = "upper left"); ax.grid(); end
axs[1].tick_params(labelbottom=false)
axs[1].set_title("minimum computation time")
wsave(plotsdir("benchmarks", "comparing_btime_N.png"), fig)

# %% computation times for different system dimensions

data = wload(datadir("benchmark", "benchmark_dim_N=10000_ddim=1_lowerdim=4_upperdim=16.csv")) |> DataFrame
fig, axs = subplots(1,1; sharex = true, figsize = (10, 16))
for (i, method) in enumerate(unique(data.Method))
    color = "C$(i-1)"

    x = data[data.Method .== method, :Dimension]
    y = data[data.Method .== method, :Time]
    axs.plot(x, y, color = color, label = method)
end
axs.set_xlabel("\$\\mathrm{Dimension}\$")
axs.set_ylabel("\$t [\\mathrm{ns}]\$")

axs.set_yscale("log")

axs.legend(loc = "upper left")
axs.grid()

axs.set_title("minimum computation time (N=10000)")
wsave(plotsdir("benchmarks", "comparing_btime_dim.png"), fig)

# %% comparison of computation times for length and dimensionality

path1 = "benchmark_dim_N=50000_ddim=1_lowerdim=4_upperdim=16.csv"
path2 = "benchmark_number_NofN=10_base=10_lowerN=3_upperN=5.csv"
titlename = "\$ \\mathrm{Computation}\\;\\mathrm{Time}\\;\\mathrm{for}\\;\\mathrm{Calculation}\\;\\mathrm{of}\\;\\mathrm{Probabilities} \$"
labels = [
    "\$ C_2^\\mathrm{classic} \$",
    "\$ C_2^\\mathrm{tree} \$",
    "\$ C_2^\\mathrm{box} \$",
    "\$ C_2^\\mathrm{prism} \$",
    "\$ H_1^\\mathrm{Molteno} \$",
    "\$ H_1^\\mathrm{classic} \$",
]

function plot_dim_length(title, labels, path1, path2; save_figure = true)
    data_dim = wload(
        datadir(
                "benchmark",
                path1,
            )) |> DataFrame
    data_var = wload(
            datadir(
                "benchmark",
                path2,
            )) |> DataFrame
    data_length = data_var[(data_var.Model .== "lorenz"), :]

    fig, axs = subplots(2,1; sharex = false)

    c_names = Dict()
    ls_names = Dict()
    agents = []
    for (i, method) in enumerate(unique(data_dim.Method))
        color = "C$(i-1)"
        ls = LINESTYLES[i]
        c_names[method] = color
        ls_names[method] = ls

        x = data_dim[data_dim.Method .== method, :Dimension]
        y = data_dim[data_dim.Method .== method, :Time]
        agent, = axs[1].plot(x, y, color = color, label = method, ls = ls)
        push!(agents, agent)
    end
    for method in unique(data_length.Method)
        if method == "Correlation"
            color = c_names["Classic Correlation"]
            ls = ls_names["Classic Correlation"]
        elseif method == "Generalized Dimension"
            color = c_names["Generalized Entropy"]
            ls = ls_names["Generalized Entropy"]
        else
            color = c_names[method]
            ls = ls_names[method]
        end

        x = data_length[data_length.Method .== method, :N]
        y = data_length[data_length.Method .== method, :Time]
        axs[2].plot(x, y, color = color, label = nothing, ls = ls)
    end

    for ax in axs; ax.set_yscale("log"); ax.grid(true); end

    axs[1].set_xlabel("\$ \\mathrm{Dimension} \$")
    axs[1].set_xticks(4:4:16)
    axs[1].set_xlim((4,16))
    axs[1].set_title("\$ \\mathrm{Lorenz-96} \\quad N=50000 \$")

    axs[1].set_ylabel("\$t [\\mathrm{ns}]\$")
    axs[1].set_yticks(10 .^ (7:12))
    leg = axs[1].legend(
        handles = agents,
        labels = labels,
        bbox_to_anchor=(0., 1.25, 1., .102), loc="lower left",
        ncol=3, mode="expand", borderaxespad=0, handlelength=2,
        title = title,
    )

    axs[2].set_xlabel("\$\\mathrm{N}\$")
    axs[2].set_yticks(10 .^ (5:2:12))
    axs[2].set_xlim((10^3, 10^5))
    axs[2].set_xscale("log")
    axs[2].set_title("\$ \\mathrm{Lorenz-63} \$")
    #axs[2].legend().remove()

    axs[2].set_ylabel("\$t [\\mathrm{ns}]\$")

    axs[1].add_artist(leg)
    fig.tight_layout()
    fig.subplots_adjust(top = 1.0, hspace = 0.5, left = 0.2)

    fig
end

fig = plot_dim_length(title, labels, path1, path2)
fig.tight_layout(;pad = 0.3)
wsave(plotsdir("benchmarks", "comparing_btime_dim_N.png"), fig)

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
using FSM
using CSV
using DataFrames
using Plots

data_force = CSV.File("data/met_CdP_0506.csv") |> DataFrame

data_ref = CSV.File("fortran/output/out_CdP_0506_11111.txt", header=["year", "month", "day", "hour", "alb", "Roff", "snowdepth", "SWE", "Tsurf", "Tsoil"], delim=" ", ignorerepeated=true) |> DataFrame

input = Input{Float64}(
    data_force.year,
    data_force.month,
    data_force.day,
    data_force.hour,
    data_force.SW,
    data_force.LW,
    data_force.Sf,
    data_force.Rf,
    data_force.Ta,
    data_force.RH,
    data_force.Ua,
    data_force.Ps,
    )

ebm = EBM{Float64}(
        am=1,
        cm=1,
        dm=1,
        em=1,
        hm=1,
        zT=1.5,
        zvar=false,
        Tsoil=[282.98, 284.17, 284.70, 284.70]
    )

cn = Constants{Float64}()

snowdepth = similar(input.Ta)
SWE = similar(input.Ta)
Tsurf = similar(input.Ta)

run!(ebm, cn, snowdepth, SWE, Tsurf, input)

plot(data_ref.SWE - SWE)
plot(data_ref.Tsurf - Tsurf .+ 273.15)


text:
if !("." in LOAD_PATH)
    push!(LOAD_PATH, ".")
end
import Coverage
import Plots

# https://github.com/jheinen/GR.jl/issues/278#issuecomment-587090846
ENV["GKSwstype"] = "nul"

all_cases = [
    # "ARM_SGP",
    # "Bomex",
    # "DryBubble",
    # "DYCOMS_RF01",
    # "GABLS",
    # "GATE_III",
    # "life_cycle_Tan2018",
    # "Nieuwstadt",
    "Rico",
    # "Soares",
    # "SP",
    "TRMM_LBA",
    "LES_driven_SCM",
]

filter!(x -> x ≠ "GATE_III", all_cases) # no mse tables for GATE_III
filter!(x -> x ≠ "SP", all_cases) # not currently running SP
allocs = Dict()
for case in all_cases
    ENV["ALLOCATION_CASE_NAME"] = case
    run(`julia --project=test/ --track-allocation=user perf/alloc_per_case.jl`)
    allocs[case] = Coverage.analyze_malloc(".")

    # Clean up files
    all_files = [joinpath(root, f) for (root, dirs, files) in Base.Filesystem.walkdir(".") for f in files]
    all_mem_files = filter(x -> endswith(x, ".mem"), all_files)
    for f in all_mem_files
        rm(f)
    end
end

@info "Post-processing allocations"

function plot_allocs(case_name, allocs_per_case, n_unique_bytes)
    p = Plots.plot()
    @info "Allocations for $case_name"

    filename_only(fn) = first(split(fn, ".jl")) * ".jl"
    function compile_tc(fn, linenumber)
        c1 = endswith(filename_only(fn), "TurbulenceConvection.jl")
        c2 = linenumber == 1
        return c1 && c2
    end

    filter!(x -> x.bytes ≠ 0, allocs_per_case)
    filter!(x -> !compile_tc(x.filename, x.linenumber), allocs_per_case)

    for alloc in allocs_per_case
        println(alloc)
    end
    println("Number of allocating sites: $(length(allocs_per_case))")
    case_bytes = getproperty.(allocs_per_case, :bytes)[end:-1:1]
    case_filename = getproperty.(allocs_per_case, :filename)[end:-1:1]
    case_linenumber = getproperty.(allocs_per_case, :linenumber)[end:-1:1]
    all_bytes = Int[]
    filenames = String[]
    linenumbers = Int[]
    loc_ids = String[]
    for (bytes, filename, linenumber) in zip(case_bytes, case_filename, case_linenumber)
        compile_tc(filename, linenumber) && continue # Skip loading module
        loc_id = "$(filename_only(filename))" * "$linenumber"
        if !(bytes in all_bytes) && !(loc_id in loc_ids)
            push!(all_bytes, bytes)
            push!(filenames, filename)
            push!(linenumbers, linenumber)
            push!(loc_ids, loc_id)
            if length(all_bytes) ≥ n_unique_bytes
                break
            end
        end
    end

    all_bytes = all_bytes ./ 10^3
    max_bytes = maximum(all_bytes)
    @info "$case_name: $all_bytes"
    xtick_name(filename, linenumber) = "$filename, line number: $linenumber"
    markershape = (:square, :hexagon, :circle, :star, :utriangle, :dtriangle)
    for (bytes, filename, linenumber) in zip(all_bytes, filenames, linenumbers)
        Plots.plot!(
            [0],
            [bytes];
            seriestype = :scatter,
            label = xtick_name(filename_only(filename), linenumber),
            markershape = markershape[1],
            markersize = 1 + bytes / max_bytes * 10,
        )
        markershape = (markershape[end], markershape[1:(end - 1)]...)
    end
    p1 = Plots.plot!(ylabel = "Allocations (KB)", title = case_name)
    p2 = Plots.plot(
        1:length(allocs_per_case),
        getproperty.(allocs_per_case, :bytes)[end:-1:1] ./ 1000;
        xlabel = "i-th allocating line (sorted)",
        ylabel = "Allocations (KB)",
        markershape = :circle,
    )
    Plots.plot(p1, p2, layout = Plots.grid(2, 1))
    Plots.savefig(joinpath(folder, "allocations_$case_name.png"))
end

folder = "perf/allocations_output"
mkpath(folder)

@info "Allocated bytes for single tendency per case:"
for case in all_cases
    plot_allocs(case, allocs[case], 10)
end


text:
"""
    Rimu.EmbarrassinglyDistributed
Module that provides an embarrassingly parallel option for breaking up long
time-series with `lomc!()` into chunks performed in parallel using the
`Distributed` package.

### Exports:
* [`d_lomc!()`](@ref) - run [`lomc!()`](@ref) in embarrassingly parallel mode
* [`combine_dfs()`](@ref) - combine the resulting `DataFrame`s to a single one
* [`setup_workers()`](@ref) - set up workers for distributed computing
* [`seedCRNGs_workers!()`](@ref) - seed random number generators for distributed computing
"""
module EmbarrassinglyDistributed

using Random, Parameters, DataFrames
using Rimu, Rimu.ConsistentRNG
using Distributed

export d_lomc!, setup_workers, seedCRNGs_workers!, combine_dfs

"""
    seedCRNGs_workers!([seed])
Seed the random number generators `CRNG` from [`ConsistentRNG`](@ref) on all
available processes in a `Distributed` environment deterministically from `seed`
but such that their pseudo-random number sequences are statistically
independent.

If no `seed` is given, obtain one from system entropy
(with [`Random.RandomDevice()`](@ref)).
"""
function seedCRNGs_workers!(seed=rand(Random.RandomDevice(),UInt))
    @everywhere seedCRNG!($seed + hash(myid()))
end

"""
    d_lomc!(ham, v; eqsteps, kwargs...)
    -> (; dfs = DataFrame[], eqsteps)
Perform linear operator Monte Carlo with [`lomc!()`](@ref) in embarrassingly parallel
mode using `Distributed` computing. Returns all dataframes.

### Keyword arguments in addition to those of [`lomc!()`](@ref):
* `eqsteps` - Number of time steps used for equilibration. Each worker will run an independent simulation with `eqstep + (laststep - step) ÷ nworkers()` time steps.

### Example:
```julia
using Rimu, Rimu.EmbarrassinglyDistributed

setup_workers(4) # set up to run on 4 workers
seedCRNGs_workers!(127) # seed random number generators for deterministic evolution

add = BoseFS((1,1,0,1))
ham = HubbardReal1D(add, u=4.0)
v = DVec(add => 2, capacity = 200)
# run `lomc!()` for 20_100 time steps by
# performing 4 parallel runs of `lomc!()` with 5_100 time steps each and
# stiching the results together into a single dataframe:
df, eqsteps = d_lomc!(ham, v; eqsteps = 100, laststep = 20_100) |> combine_dfs
# or
energies = d_lomc!(ham, v; eqsteps = 100, laststep = 20_100) |> combine_dfs |> autoblock
```

### See also:
* [`setup_workers()`](@ref)
* [`seedCRNGs_workers!()`](@ref)
* [`combine_dfs()`](@ref)
"""
function d_lomc!(ham, v;
    eqsteps,
    laststep = nothing,
    params::Rimu.FciqmcRunStrategy = Rimu.RunTillLastStep(),
    kwargs...
    )

    nw = nworkers()
    nw < 2 && @warn "Not enough workers available for parallel execution" nw

    if !isnothing(laststep)
        params.laststep = laststep
    end
    # unpack the parameters:
    @unpack step, laststep = params
    @assert laststep - step - eqsteps > nw "not enough time steps to run for"
    psteps = laststep - step - eqsteps # steps to be run after equilibration
    stepseach = psteps ÷ nw # to be run on each worker after equilibration

    # @everywhere do_it() = Rimu.lomc!($ham, sizehint!($v,($vc*3)>>1); $kwargs..., laststep = $step+$stepseach, params = $params, threading = false).df
    # start shorter jobs in parallel
    # futures = [@spawnat(p, Main.do_it()) for p in workers()]
    futures = [@spawnat(p, Rimu.lomc!(ham, v; kwargs..., laststep = step+eqsteps+stepseach, params = params, threading = false).df) for p in workers()]
    # futures = [@spawnat(p, Rimu.lomc!(ham, v; kwargs..., laststep = step+stepseach, params = params, threading = false).df) for p in workers()]
    dfs = fetch.(futures) # will now be an array of dataframes
    return (; dfs, eqsteps) # returns NamedTuple
end

"""
    combine_dfs(dfs::AbstractVector, eqsteps)
    combine_dfs((dfs, eqsteps))
Return a dataframe with compounded data from
the parallel runs assuming that after `eqsteps` time steps the runs are
equilibrated and uncorrelated.
"""
combine_dfs(tuple) = combine_dfs(tuple...) # for easy chaining

function combine_dfs(dfs::AbstractVector, eqsteps)
    # add worker id to dfs
    for i in eachindex(dfs)
        df = dfs[i]
        df.workerid = [workers()[i] for n in 1:size(df)[1]]
        df.stepsorg = copy(df.steps) # remember original time step
    end
    # check that the dataframes are all the same size:
    @assert mapreduce(isequal(size(dfs[1])), &, size.(dfs))

    # last time step recorded in the first dataframe
    s = dfs[1].steps[end]

    # create views of all the remaining dataframes starting after the eqsteps
    lelem = size(dfs[1])[1]
    df_adds = [view(dfs[i], (eqsteps+2):lelem, :) for i in 2:length(dfs)]
    # now number the time steps in the views consecutively
    for df in df_adds
        df[:,:steps] .= [s+i for i in 1:length(df.steps)]
        s = df.steps[end]
    end
    dfm = vcat(dfs[1],df_adds...)
    return (; df=dfm, eqsteps)
end # d_lomc!()

"""
    setup_workers([nw])
Set up and prepare `nw` workers for distributed computing. If `nw` is not given,
but multiple threads are available then as many workers will be prepared. The
`Rimu` package code will be made available on all workers.
"""
function setup_workers(nw = nothing)
    # organise the right number of workers
    if isnothing(nw)
        # if `nw` is not specified, use at least as many workers are threads are
        # available
        while Threads.nthreads() > nworkers()
            addprocs(1)
        end
        nw = nworkers()
    else
        while nw > nworkers()
            addprocs(1)
        end
    end
    @assert nw ≥ 1 "not enough workers to run in parallel"
    while nworkers() > nw
        # reduce the number of workers if it is larger than `nw`
        rmprocs(workers()[end])
    end
    # just to check that we have the right number of workers now:
    @assert nw == nworkers()

    @eval @everywhere using Rimu # load code on all walkers

    return nw
end

# function __init__()
#     if myid() == 1
#         nw = setup_workers()
#         @info "`REmbarrassinglyDistributed`: $nw workers set up and ready to go"
#     end
# end

end # module


text:
Base.one(::Type{P}) where {C, T, P <: Polynomial{C, T}} = Polynomial(one(T))
Base.one(p::Polynomial) = one(typeof(p))

Base.zero(::Type{Polynomial{C, T, A}}) where {C, T, A} = Polynomial(A())
Base.zero(t::PolynomialLike) = zero(typeof(t))

combine(t1::Term, t2::Term) = combine(promote(t1, t2)...)
combine(t1::T, t2::T) where {T <: Term} = Term(t1.coefficient + t2.coefficient, t1.monomial)
compare(t1::Term, t2::Term) = monomial(t1) > monomial(t2)

# Graded Lexicographic order
# First compare total degree, then lexicographic order
function isless(m1::Monomial{V}, m2::Monomial{V}) where {V}
    d1 = degree(m1)
    d2 = degree(m2)
    if d1 < d2
        return true
    elseif d1 > d2
        return false
    else
        return exponents(m1) < exponents(m2)
    end
end

jointerms(terms1::AbstractArray{<:Term}, terms2::AbstractArray{<:Term}) = mergesorted(terms1, terms2, compare, combine)

(+)(p1::Polynomial, p2::Polynomial) = Polynomial(jointerms(terms(p1), terms(p2)))
(-)(p1::Polynomial, p2::Polynomial) = Polynomial(jointerms(terms(p1), (-).(terms(p2))))

(==)(::Variable{N}, ::Variable{N}) where {N} = true
(==)(::Variable, ::Variable) = false
(==)(m1::Monomial{V}, m2::Monomial{V}) where {V} = exponents(m1) == exponents(m2)

# Multiplication is handled as a special case so that we can write these
# definitions without resorting to promotion:
MP.multconstant(α, v::Monomial) = Term(α, v)
MP.multconstant(α, v::Variable) = Term(α, Monomial(v))

(*)(v1::V, v2::V) where {V <: Variable} = Monomial{(V(),), 1}((2,))
(*)(v1::Variable, v2::Variable) = (*)(promote(v1, v2)...)

function MP.divides(m1::Monomial{V, N}, m2::Monomial{V, N}) where {V, N}
    reduce((d, exp) -> d && (exp[1] <= exp[2]), zip(m1.exponents, m2.exponents), init=true)
end
MP.divides(m1::Monomial, m2::Monomial) = divides(promote(m1, m2)...)
function MP.mapexponents(op, m1::M, m2::M) where M<:Monomial
    M(map(op, m1.exponents, m2.exponents))
end
MP.mapexponents(op, m1::Monomial, m2::Monomial) = mapexponents(op, promote(m1, m2)...)

# TODO: this could be faster with an in-place summation
function (*)(p1::Polynomial{S}, p2::Polynomial{T}) where {S, T}
    C = Base.promote_op(*, S, T)
    M = promote_type(monomialtype(p1), monomialtype(p2))
    result = Polynomial(termtype(M, C)[])
    for t1 in terms(p1)
        for t2 in terms(p2)
            result += t1 * t2
        end
    end
    result
end

^(v::V, x::Integer) where {V <: Variable} = Monomial{(V(),), 1}((x,))

# dot(v1::AbstractVector{<:TermLike}, v2::AbstractVector) = dot(v1, v2)
# dot(v1::AbstractVector, v2::AbstractVector{<:TermLike}) = dot(v1, v2)
# dot(v1::AbstractVector{<:TermLike}, v2::AbstractVector{<:TermLike}) = dot(v1, v2)

# All of these types are immutable, so there's no need to copy anything to get
# a shallow copy.
Base.copy(x::TermLike) = x

Base.copy(p::Polynomial) = Polynomial(copy(terms(p)))

adjoint(v::Variable) = v
adjoint(m::Monomial) = m
adjoint(t::Term) = Term(adjoint(coefficient(t)), monomial(t))
adjoint(x::Polynomial) = Polynomial(adjoint.(terms(x)))


text:
using LinearAlgebra: tril!, triu!
function batched_tril!(A::CuArray{T, N}, d) where {T, N}
  if N < 2
    error("MethodError: no method matching tril!(::Array{Float64,1})")
  elseif N == 2
    return tril!(x, d)
  else
    s = size(A)
    m, n = s[1], s[2]
    l = m*n
    bs = Int(length(A) // l)
    function batch_tril_kernel!(_A, _d)
      li = (blockIdx().x - 1) * blockDim().x + threadIdx().x
      b = (blockIdx().y - 1) * blockDim().y + threadIdx().y
      @inbounds if 0 < li <= l && b <= bs
        id = Tuple(CartesianIndices(_A)[Base._to_linear_index(_A, mod1(li, m), fld1(li, m), b)])
        i, j = id
        if i < j - _d
          _A[id...] = 0
        end
      end
      return nothing
    end
    max_threads = 256
    thread_x = min(max_threads, l)
    thread_y = min(max_threads ÷ thread_x, bs)
    threads = (thread_x, thread_y)
    blocks = ceil.(Int, (l, bs) ./ threads)
    @cuda threads=threads blocks=blocks batch_tril_kernel!(A, d)
    return A
  end
end

function batched_triu!(A::CuArray{T, N}, d) where {T, N}
  if N < 2
    error("MethodError: no method matching triu!(::Array{Float64,1})")
  elseif N == 2
    return tril!(x, d)
  else
    s = size(A)
    m, n = s[1], s[2]
    l = m*n
    bs = Int(length(A) // l)
    function batch_triu_kernel!(_A, _d)
      li = (blockIdx().x - 1) * blockDim().x + threadIdx().x
      b = (blockIdx().y - 1) * blockDim().y + threadIdx().y
      @inbounds if 0 < li <= l && b <= bs
        id = Tuple(CartesianIndices(_A)[Base._to_linear_index(_A, mod1(li, m), fld1(li, m), b)])
        i, j = id
        if j < i + _d
          _A[id...] = 0
        end
      end
      return nothing
    end
    max_threads = 256
    thread_x = min(max_threads, l)
    thread_y = min(max_threads ÷ thread_x, bs)
    threads = (thread_x, thread_y)
    blocks = ceil.(Int, (l, bs) ./ threads)
    @cuda threads=threads blocks=blocks batch_triu_kernel!(A, d)
    return A
  end
end


text:
const SMALL_SIGNATURES = get(ENV, "PBC_SMALL_SIGNATURES", "n") == "y"
const SMALL_IDENTITIES = get(ENV, "PBC_SMALL_IDENTITIES", "y") == "y"
const NPROCS = get(ENV, "PBC_NPROCS", "") == "auto" ? Sys.CPU_THREADS : parse(Int, get(ENV, "PBC_NPROCS", "1"))
const BATCH_SIZE = Threads.nthreads() * parse(Int, get(ENV, "PBC_BATCH_SCALE_FACTOR", "1024"))

macro EP()
    return SMALL_SIGNATURES ? :(Curve.EP) : :(Curve.EP2)
end
macro EP2()
    return SMALL_SIGNATURES ? :(Curve.EP2) : :(Curve.EP)
end
macro ID()
    return SMALL_IDENTITIES ? :(Int64) : :(Int128)
end

const G1 = Curve.curve_gen(@EP)
const G2 = Curve.curve_gen(@EP2)
const PUBLIC_KEY_SIZE = sizeof(G2.x)
const SIGNATURE_SIZE = sizeof(G1.x)

const PRIME = Curve.fp_prime_get()
const ORDER = Curve.curve_order(@EP)

# Effective number of bits required to store the private key
const PRIVATE_KEY_SIZE_BITS = ceil(Int, log2(BigInt(ORDER)))

# Effective number of bytes required to store the private key
const PRIVATE_KEY_SIZE = ceil(Int, PRIVATE_KEY_SIZE_BITS // 8)

# whether or not we have bits over for encoding LSB(y)
const CAN_STUFF_Y = ceil(Int, log2(BigInt(PRIME))) < (8 * Curve.FP_ST_SIZE)


text:
int_rules_1_1_1_2 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.1.1.2*(a+b*x)^m*(c+d*x)^n =#
    @apply_utils Antiderivative((~a + ~(b') * ~x) ^ ~(m') * (~c + ~(d') * ~x), ~x) => (~d * ~x * (~a + ~b * ~x) ^ (~m + 1)) / (~b * (~m + 2)) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && EqQ(~a * ~d - ~b * ~c * (~m + 2), 0)
    @apply_utils Antiderivative(1 / ((~a + ~(b') * ~x) * (~c + ~(d') * ~x)), ~x) => Antiderivative(1 / (~a * ~c + ~b * ~d * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ(~b * ~c + ~a * ~d, 0)
    @apply_utils Antiderivative(1 / ((~(a') + ~(b') * ~x) * (~(c') + ~(d') * ~x)), ~x) => (~b / (~b * ~c - ~a * ~d)) * Antiderivative(1 / (~a + ~b * ~x), ~x) - (~d / (~b * ~c - ~a * ~d)) * Antiderivative(1 / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ(~b * ~c - ~a * ~d, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~(m') * (~(c') + ~(d') * ~x) ^ ~n, ~x) => ((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n + 1)) / ((~b * ~c - ~a * ~d) * (~m + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ(~m + ~n + 2, 0) && NeQ(~m, -1)))
    @apply_utils Antiderivative((~a + ~(b') * ~x) ^ ~m * (~c + ~(d') * ~x) ^ ~m, ~x) => (~x * (~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~m) / (2 * ~m + 1) + ((2 * ~a * ~c * ~m) / (2 * ~m + 1)) * Antiderivative((~a + ~b * ~x) ^ (~m - 1) * (~c + ~d * ~x) ^ (~m - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && IGtQ(~m + 1 / 2, 0))
    @apply_utils Antiderivative(1 / ((~a + ~(b') * ~x) ^ (3 / 2) * (~c + ~(d') * ~x) ^ (3 / 2)), ~x) => ~x / (~a * ~c * sqrt(~a + ~b * ~x) * sqrt(~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ(~b * ~c + ~a * ~d, 0)
    @apply_utils Antiderivative((~a + ~(b') * ~x) ^ ~m * (~c + ~(d') * ~x) ^ ~m, ~x) => (-(~x) * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~m + 1)) / (2 * ~a * ~c * (~m + 1)) + ((2 * ~m + 3) / (2 * ~a * ~c * (~m + 1))) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~m + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && ILtQ(~m + 3 / 2, 0))
    @apply_utils Antiderivative((~a + ~(b') * ~x) ^ ~(m') * (~c + ~(d') * ~x) ^ ~(m'), ~x) => Antiderivative((~a * ~c + ~b * ~d * (~x) ^ 2) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (IntegerQ(~m) || GtQ(~a, 0) && GtQ(~c, 0)))
    @apply_utils Antiderivative((~a + ~(b') * ~x) ^ ~m * (~c + ~(d') * ~x) ^ ~m, ~x) => (((~a + ~b * ~x) ^ FracPart(~m) * (~c + ~d * ~x) ^ FracPart(~m)) / (~a * ~c + ~b * ~d * (~x) ^ 2) ^ FracPart(~m)) * Antiderivative((~a * ~c + ~b * ~d * (~x) ^ 2) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && Not(IntegerQ(2 * ~m)))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n, ~x) => ((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n) / (~b * (~m + 1)) - ((~d * ~n) / (~b * (~m + 1))) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (ILtQ(~m, -1) && (Not(IntegerQ(~n)) && GtQ(~n, 0))))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n, ~x) => ((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n + 1)) / ((~b * ~c - ~a * ~d) * (~m + 1)) - ((~d * (~m + ~n + 2)) / ((~b * ~c - ~a * ~d) * (~m + 1))) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (ILtQ(~m, -1) && (Not(IntegerQ(~n)) && LtQ(~n, 0))))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~(m') * (~(c') + ~(d') * ~x) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (IGtQ(~m, 0) && (Not(IntegerQ(~n)) || (EqQ(~c, 0) && LeQ(7 * ~m + 4 * ~n + 4, 0) || (LtQ(9 * ~m + 5 * (~n + 1), 0) || GtQ(~m + ~n + 2, 0))))))
    @apply_utils Antiderivative((~a + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (ILtQ(~m, 0) && (IntegerQ(~n) && Not(IGtQ(~n, 0) && LtQ(~m + ~n + 2, 0)))))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n, ~x) => ((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n + 1)) / ((~b * ~c - ~a * ~d) * (~m + 1)) - ((~d * Simplify(~m + ~n + 2)) / ((~b * ~c - ~a * ~d) * (~m + 1))) * Antiderivative((~a + ~b * ~x) ^ Simplify(~m + 1) * (~c + ~d * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (ILtQ(Simplify(~m + ~n + 2), 0) && (NeQ(~m, -1) && (Not(LtQ(~m, -1) && (LtQ(~n, -1) && (EqQ(~a, 0) || NeQ(~c, 0) && (LtQ(~m - ~n, 0) && IntegerQ(~n))))) && (SumSimplerQ(~m, 1) || Not(SumSimplerQ(~n, 1)))))))
    @apply_utils Antiderivative(1 / ((~a + ~(b') * ~x) ^ (9 / 4) * (~c + ~(d') * ~x) ^ (1 / 4)), ~x) => -4 / (5 * ~b * (~a + ~b * ~x) ^ (5 / 4) * (~c + ~d * ~x) ^ (1 / 4)) - (~d / (5 * ~b)) * Antiderivative(1 / ((~a + ~b * ~x) ^ (5 / 4) * (~c + ~d * ~x) ^ (5 / 4)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && NegQ((~a) ^ 2 * (~b) ^ 2))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n, ~x) => ((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n) / (~b * (~m + 1)) - ((~d * ~n) / (~b * (~m + 1))) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (GtQ(~n, 0) && (LtQ(~m, -1) && (Not(IntegerQ(~n) && Not(IntegerQ(~m))) && (Not(ILeQ(~m + ~n + 2, 0) && (FractionQ(~m) || GeQ(2 * ~n + ~m + 1, 0))) && IntLinearQ(~a, ~b, ~c, ~d, ~m, ~n, ~x))))))
    @apply_utils Antiderivative(1 / ((~a + ~(b') * ~x) ^ (5 / 4) * (~c + ~(d') * ~x) ^ (1 / 4)), ~x) => -2 / (~b * (~a + ~b * ~x) ^ (1 / 4) * (~c + ~d * ~x) ^ (1 / 4)) + ~c * Antiderivative(1 / ((~a + ~b * ~x) ^ (5 / 4) * (~c + ~d * ~x) ^ (5 / 4)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && NegQ((~a) ^ 2 * (~b) ^ 2))
    @apply_utils Antiderivative((~a + ~(b') * ~x) ^ ~m * (~c + ~(d') * ~x) ^ ~n, ~x) => ((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n) / (~b * (~m + ~n + 1)) + ((2 * ~c * ~n) / (~m + ~n + 1)) * Antiderivative((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (IGtQ(~m + 1 / 2, 0) && (IGtQ(~n + 1 / 2, 0) && LtQ(~m, ~n))))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n, ~x) => ((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n) / (~b * (~m + ~n + 1)) + ((~n * (~b * ~c - ~a * ~d)) / (~b * (~m + ~n + 1))) * Antiderivative((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (GtQ(~n, 0) && (NeQ(~m + ~n + 1, 0) && (Not(IGtQ(~m, 0) && (Not(IntegerQ(~n)) || GtQ(~m, 0) && LtQ(~m - ~n, 0))) && (Not(ILtQ(~m + ~n + 2, 0)) && IntLinearQ(~a, ~b, ~c, ~d, ~m, ~n, ~x))))))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n, ~x) => ((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n + 1)) / ((~b * ~c - ~a * ~d) * (~m + 1)) - ((~d * (~m + ~n + 2)) / ((~b * ~c - ~a * ~d) * (~m + 1))) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (LtQ(~m, -1) && (Not(LtQ(~n, -1) && (EqQ(~a, 0) || NeQ(~c, 0) && (LtQ(~m - ~n, 0) && IntegerQ(~n)))) && IntLinearQ(~a, ~b, ~c, ~d, ~m, ~n, ~x))))
    @apply_utils Antiderivative(1 / (sqrt(~a + ~(b') * ~x) * sqrt(~c + ~(d') * ~x)), ~x) => acosh((~b * ~x) / ~a) / ~b <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~a + ~c, 0) && (EqQ(~b - ~d, 0) && GtQ(~a, 0)))
    @apply_utils Antiderivative(1 / (sqrt(~a + ~(b') * ~x) * sqrt(~(c') + ~(d') * ~x)), ~x) => Antiderivative(1 / sqrt((~a * ~c - ~b * (~a - ~c) * ~x) - (~b) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~b + ~d, 0) && GtQ(~a + ~c, 0))
    @apply_utils Antiderivative(1 / (sqrt(~(a') + ~(b') * ~x) * sqrt(~(c') + ~(d') * ~x)), ~x) => (2 / sqrt(~b)) * Subst(Antiderivative(1 / sqrt((~b * ~c - ~a * ~d) + ~d * (~x) ^ 2), ~x), ~x, sqrt(~a + ~b * ~x)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (GtQ(~b * ~c - ~a * ~d, 0) && GtQ(~b, 0))
    @apply_utils Antiderivative(1 / ((~(a') + ~(b') * ~x) * (~(c') + ~(d') * ~x) ^ (1 / 3)), ~x) => With([q = Rt((~b * ~c - ~a * ~d) / ~b, 3)], (-(log(RemoveContent(~a + ~b * ~x, ~x))) / (2 * ~b * q) - (3 / (2 * ~b * q)) * Subst(Antiderivative(1 / (q - ~x), ~x), ~x, (~c + ~d * ~x) ^ (1 / 3))) + (3 / (2 * ~b)) * Subst(Antiderivative(1 / (q ^ 2 + q * ~x + (~x) ^ 2), ~x), ~x, (~c + ~d * ~x) ^ (1 / 3))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && PosQ((~b * ~c - ~a * ~d) / ~b)
    @apply_utils Antiderivative(1 / ((~(a') + ~(b') * ~x) * (~(c') + ~(d') * ~x) ^ (1 / 3)), ~x) => With([q = Rt(-((~b * ~c - ~a * ~d)) / ~b, 3)], (log(RemoveContent(~a + ~b * ~x, ~x)) / (2 * ~b * q) - (3 / (2 * ~b * q)) * Subst(Antiderivative(1 / (q + ~x), ~x), ~x, (~c + ~d * ~x) ^ (1 / 3))) + (3 / (2 * ~b)) * Subst(Antiderivative(1 / ((q ^ 2 - q * ~x) + (~x) ^ 2), ~x), ~x, (~c + ~d * ~x) ^ (1 / 3))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NegQ((~b * ~c - ~a * ~d) / ~b)
    @apply_utils Antiderivative(1 / ((~(a') + ~(b') * ~x) * (~(c') + ~(d') * ~x) ^ (2 / 3)), ~x) => With([q = Rt((~b * ~c - ~a * ~d) / ~b, 3)], (-(log(RemoveContent(~a + ~b * ~x, ~x))) / (2 * ~b * q ^ 2) - (3 / (2 * ~b * q ^ 2)) * Subst(Antiderivative(1 / (q - ~x), ~x), ~x, (~c + ~d * ~x) ^ (1 / 3))) - (3 / (2 * ~b * q)) * Subst(Antiderivative(1 / (q ^ 2 + q * ~x + (~x) ^ 2), ~x), ~x, (~c + ~d * ~x) ^ (1 / 3))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && PosQ((~b * ~c - ~a * ~d) / ~b)
    @apply_utils Antiderivative(1 / ((~(a') + ~(b') * ~x) * (~(c') + ~(d') * ~x) ^ (2 / 3)), ~x) => With([q = Rt(-((~b * ~c - ~a * ~d)) / ~b, 3)], -(log(RemoveContent(~a + ~b * ~x, ~x))) / (2 * ~b * q ^ 2) + (3 / (2 * ~b * q ^ 2)) * Subst(Antiderivative(1 / (q + ~x), ~x), ~x, (~c + ~d * ~x) ^ (1 / 3)) + (3 / (2 * ~b * q)) * Subst(Antiderivative(1 / ((q ^ 2 - q * ~x) + (~x) ^ 2), ~x), ~x, (~c + ~d * ~x) ^ (1 / 3))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NegQ((~b * ~c - ~a * ~d) / ~b)
    @apply_utils Antiderivative(1 / ((~(a') + ~(b') * ~x) ^ (1 / 3) * (~(c') + ~(d') * ~x) ^ (2 / 3)), ~x) => With([q = Rt(~d / ~b, 3)], (((-(sqrt(3)) * q) / ~d) * atan((2 * q * (~a + ~b * ~x) ^ (1 / 3)) / (sqrt(3) * (~c + ~d * ~x) ^ (1 / 3)) + 1 / sqrt(3)) - (q / (2 * ~d)) * log(~c + ~d * ~x)) - ((3q) / (2 * ~d)) * log((q * (~a + ~b * ~x) ^ (1 / 3)) / (~c + ~d * ~x) ^ (1 / 3) - 1)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && PosQ(~d / ~b))
    @apply_utils Antiderivative(1 / ((~(a') + ~(b') * ~x) ^ (1 / 3) * (~(c') + ~(d') * ~x) ^ (2 / 3)), ~x) => With([q = Rt(-(~d) / ~b, 3)], ((sqrt(3) * q) / ~d) * atan(1 / sqrt(3) - (2 * q * (~a + ~b * ~x) ^ (1 / 3)) / (sqrt(3) * (~c + ~d * ~x) ^ (1 / 3))) + (q / (2 * ~d)) * log(~c + ~d * ~x) + ((3q) / (2 * ~d)) * log((q * (~a + ~b * ~x) ^ (1 / 3)) / (~c + ~d * ~x) ^ (1 / 3) + 1)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && NegQ(~d / ~b))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~c + ~(d') * ~x) ^ ~m, ~x) => (((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~m) / (~a * ~c + (~b * ~c + ~a * ~d) * ~x + ~b * ~d * (~x) ^ 2) ^ ~m) * Antiderivative((~a * ~c + (~b * ~c + ~a * ~d) * ~x + ~b * ~d * (~x) ^ 2) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (LtQ(-1, ~m, 0) && (LeQ(3, Denominator(~m), 4) && AtomQ(~b * ~c + ~a * ~d))))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~c + ~(d') * ~x) ^ ~m, ~x) => (((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~m) / ((~a + ~b * ~x) * (~c + ~d * ~x)) ^ ~m) * Antiderivative((~a * ~c + (~b * ~c + ~a * ~d) * ~x + ~b * ~d * (~x) ^ 2) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (LtQ(-1, ~m, 0) && LeQ(3, Denominator(~m), 4)))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n, ~x) => With([p = Denominator(~m)], (p / ~b) * Subst(Antiderivative((~x) ^ (p * (~m + 1) - 1) * ((~c - (~a * ~d) / ~b) + (~d * (~x) ^ p) / ~b) ^ ~n, ~x), ~x, (~a + ~b * ~x) ^ (1 / p))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (LtQ(-1, ~m, 0) && (LeQ(-1, ~n, 0) && (LeQ(Denominator(~n), Denominator(~m)) && IntLinearQ(~a, ~b, ~c, ~d, ~m, ~n, ~x)))))
    @apply_utils Antiderivative((~(b') * ~x) ^ ~m * (~c + ~(d') * ~x) ^ ~n, ~x) => (((~c) ^ ~n * (~b * ~x) ^ (~m + 1)) / (~b * (~m + 1))) * HypergeometricFunctions._₂F₁(-(~n), ~m + 1, ~m + 2, (-(~d) * ~x) / ~c) <-- FreeQ([~b, ~c, ~d, ~m, ~n], ~x) && (Not(IntegerQ(~m)) && (IntegerQ(~n) || GtQ(~c, 0) && Not(EqQ(~n, -1 / 2) && (EqQ((~c) ^ 2 - (~d) ^ 2, 0) && GtQ(-(~d) / (~b * ~c), 0)))))
    @apply_utils Antiderivative((~(b') * ~x) ^ ~m * (~c + ~(d') * ~x) ^ ~n, ~x) => ((~c + ~d * ~x) ^ (~n + 1) / (~d * (~n + 1) * (-(~d) / (~b * ~c)) ^ ~m)) * HypergeometricFunctions._₂F₁(-(~m), ~n + 1, ~n + 2, 1 + (~d * ~x) / ~c) <-- FreeQ([~b, ~c, ~d, ~m, ~n], ~x) && (Not(IntegerQ(~n)) && (IntegerQ(~m) || GtQ(-(~d) / (~b * ~c), 0)))
    @apply_utils Antiderivative((~(b') * ~x) ^ ~m * (~c + ~(d') * ~x) ^ ~n, ~x) => (((~c) ^ IntPart(~n) * (~c + ~d * ~x) ^ FracPart(~n)) / (1 + (~d * ~x) / ~c) ^ FracPart(~n)) * Antiderivative((~b * ~x) ^ ~m * (1 + (~d * ~x) / ~c) ^ ~n, ~x) <-- FreeQ([~b, ~c, ~d, ~m, ~n], ~x) && (Not(IntegerQ(~m)) && (Not(IntegerQ(~n)) && (Not(GtQ(~c, 0)) && (Not(GtQ(-(~d) / (~b * ~c), 0)) && (RationalQ(~m) && Not(EqQ(~n, -1 / 2) && EqQ((~c) ^ 2 - (~d) ^ 2, 0)) || Not(RationalQ(~n)))))))
    @apply_utils Antiderivative((~(b') * ~x) ^ ~m * (~c + ~(d') * ~x) ^ ~n, ~x) => ((((-(~b) * ~c) / ~d) ^ IntPart(~m) * (~b * ~x) ^ FracPart(~m)) / ((-(~d) * ~x) / ~c) ^ FracPart(~m)) * Antiderivative(((-(~d) * ~x) / ~c) ^ ~m * (~c + ~d * ~x) ^ ~n, ~x) <-- FreeQ([~b, ~c, ~d, ~m, ~n], ~x) && (Not(IntegerQ(~m)) && (Not(IntegerQ(~n)) && (Not(GtQ(~c, 0)) && Not(GtQ(-(~d) / (~b * ~c), 0)))))
    @apply_utils Antiderivative((~a + ~(b') * ~x) ^ ~m * (~c + ~(d') * ~x) ^ ~n, ~x) => (((~b * ~c - ~a * ~d) ^ ~n * (~a + ~b * ~x) ^ (~m + 1)) / ((~b) ^ (~n + 1) * (~m + 1))) * HypergeometricFunctions._₂F₁(-(~n), ~m + 1, ~m + 2, (-(~d) * (~a + ~b * ~x)) / (~b * ~c - ~a * ~d)) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (Not(IntegerQ(~m)) && IntegerQ(~n)))
    @apply_utils Antiderivative((~a + ~(b') * ~x) ^ ~m * (~c + ~(d') * ~x) ^ ~n, ~x) => ((~a + ~b * ~x) ^ (~m + 1) / (~b * (~m + 1) * (~b / (~b * ~c - ~a * ~d)) ^ ~n)) * HypergeometricFunctions._₂F₁(-(~n), ~m + 1, ~m + 2, (-(~d) * (~a + ~b * ~x)) / (~b * ~c - ~a * ~d)) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (Not(IntegerQ(~m)) && (Not(IntegerQ(~n)) && (GtQ(~b / (~b * ~c - ~a * ~d), 0) && (RationalQ(~m) || Not(RationalQ(~n) && GtQ(-(~d) / (~b * ~c - ~a * ~d), 0)))))))
    @apply_utils Antiderivative((~a + ~(b') * ~x) ^ ~m * (~c + ~(d') * ~x) ^ ~n, ~x) => ((~c + ~d * ~x) ^ FracPart(~n) / ((~b / (~b * ~c - ~a * ~d)) ^ IntPart(~n) * ((~b * (~c + ~d * ~x)) / (~b * ~c - ~a * ~d)) ^ FracPart(~n))) * Antiderivative((~a + ~b * ~x) ^ ~m * Simp((~b * ~c) / (~b * ~c - ~a * ~d) + (~b * ~d * ~x) / (~b * ~c - ~a * ~d), ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (Not(IntegerQ(~m)) && (Not(IntegerQ(~n)) && (RationalQ(~m) || Not(SimplerQ(~n + 1, ~m + 1))))))
    @apply_utils Antiderivative((~(a') + ~(b') * ~u) ^ ~(m') * (~(c') + ~(d') * ~u) ^ ~(n'), ~x) => (1 / Coefficient(~u, ~x, 1)) * Subst(Antiderivative((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n, ~x), ~x, ~u) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (LinearQ(~u, ~x) && NeQ(Coefficient(~u, ~x, 0), 0))
    #= IntLinearQ(a,b,c,d,m,n,x)*returns*True*iff*(a+b*x)^m*(c+d*x)^n*is  integrable*wrt*x*in*terms*of*non-hypergeometric*functions. =#IntLinearQ((~a), (~b), (~c), (~d), (~m), (~n), (~x)) := IGtQ(m, 0) || IGtQ(n, 0) || IntegersQ(3*m, 3*n) || IntegersQ(4*m, 4*n) || IntegersQ(2*m, 6*n) || IntegersQ(6*m, 2*n) || ILtQ(m + n, -1) || IntegerQ(m + n) && RationalQ(m)
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.2 (a+b x)^m (c+d x)^n.jl


text:
using DataFrames: DataFrame, DataFrame!
using CSV: File
using DataStructures
using Parameters
using Convex, SCS
using Plots
using MathOptInterface
const MOI = MathOptInterface

using EllipsisNotation

#=

Notation:

pᵢ      - inflow proportions
CFᵢ     - cumulative inflows
CFₒ     - cumulative outflows
ECᵢ      - expected costs of vehicles coming in
ECₒ      - expected costs of vehicles going out

DTC       - departure time choices
SR      - splitting rates

=#

ROUND_DIGITS = 13

include("types.jl")
include("crange.jl")
include("utils.jl")
include("models.jl")

ugap, bgap = 3, 4 # used in computing sending and receiving functions
u, w = 1/ugap, 1/bgap
Q, N, δ = 24.0, 80, w/u
demand_level, T, Tm = 1800, 132*ugap, 80*ugap
M = 10 * T
α, β, γ, trgt = 1., 0.5, 2., 60*ugap

linkdata = DataFrame!(File("nguyen.csv"; delim=", "))
linkdata.l *= ugap;
net = Network(linkdata);

nclasses, nsinks, nlinks = 1, 2, numlinks(net);
srcs, snks, mrgs, divs = sources(net), sinks(net), merges(net), diverges(net);
trips = ones(length(srcs), length(snks), nclasses) .* demand_level;

include("initchoices.jl");
include("simulation.jl");
include("costcomputation.jl");
include("equilibrium.jl");

DTC, SR = initchoices();

i = 5; DTC, SR = DTC_list[i], SR_list[i];
CFᵢ, CFₒ, states = simulate();
ECᵢ, ECₒ, rtracs = computecosts();
relgap(ECᵢ, DTC)

DTC, SR = updatechoices!()

relgap(ECᵢ, DTC)

DTC_list = [DTC]
SR_list = [SR]

for i in 1:30
    global DTC, SR, CFᵢ, CFₒ, ECᵢ, ECₒ
    #@show i
    CFᵢ, CFₒ, states = simulate()#(DTC, SR)
    ECᵢ, ECₒ, rtracs = computecosts()#(states)
    @show relgap(ECᵢ, DTC)

    #@assert approxpositive(Fᵢ) # all((Fᵢ .>= 0.) .| (Fᵢ .≈ 0.))
    #@assert approxpositive(Fₒ) # all(Fₒ .>= 0.)
    @assert all(ECᵢ .>= 0.)
    @assert all(ECₒ .>= 0.)

    newDTC, newSR = updatechoices!() #(ECᵢ, DTC, SR)
    push!(DTC_list, DTC)
    push!(SR_list, SR)
    DTC, SR = newDTC, newSR
end


text:
using MatrixCompletion

@testset "$(format("Exponential Family: forward_map[poisson]"))" begin
    let
        tc1 = rand(100)
        @test check(:l2diff, forward_map(Poisson(),tc1),exp.(tc1),0)
        @test check(:l2diff, forward_map(:Poisson,tc1),exp.(tc1),0)
        tc2 = rand(100,100)
        @test check(:l2diff, forward_map(Poisson(),tc2),exp.(tc2),0)
        @test check(:l2diff, forward_map(:Poisson,tc2),exp.(tc2),0)
    end
end

@testset "$(format("Exponential Family: forward_map[gamma]"))" begin
    let
        tc1 = rand(100)
        @test check(:l2diff, forward_map(Gamma(),tc1),1 ./ tc1,0)
        @test check(:l2diff, forward_map(:Gamma,tc1),1 ./ tc1,0)
        tc2 = rand(100,100)
        @test check(:l2diff, forward_map(Gamma(),tc2),1 ./ tc2,0)
        @test check(:l2diff, forward_map(:Gamma,tc2),1 ./ tc2,0)
    end
end

@testset "$(format("Exponential Family: forward_map[gaussian]"))" begin
    let
        tc1 = rand(100)
        @test check(:l2diff, forward_map(Gaussian(),tc1),tc1,0)
        @test check(:l2diff, forward_map(:Gaussian,tc1), tc1,0)
        tc2 = rand(100,100)
        @test check(:l2diff, forward_map(Gaussian(),tc2),tc2,0)
        @test check(:l2diff, forward_map(:Gaussian,tc2),tc2,0)
    end
end

@testset "$(format("Exponential Family: forward_map[bernoulli]"))" begin
    let
        logit = (x) -> log.(x./(1 .- x))
        tc1 = rand(100)
        tc1_logit = logit(tc1)
        @test check(:l2diff, forward_map(Bernoulli(),tc1_logit),tc1,0)
        @test check(:l2diff, forward_map(:Bernoulli,tc1_logit), tc1,0)
        tc2 = rand(100,100)
        tc2_logit = logit(tc2)
        @test check(:l2diff, forward_map(Bernoulli(),tc2_logit),tc2,0)
        @test check(:l2diff, forward_map(:Bernoulli,tc2_logit), tc2,0)
    end
end

@testset "$(format("Exponential Family: predict[poisson]"))" begin
    let
        tc1 = rand(2:20,100)
        log_tc1 = log.(tc1)
        @test check(:l2diff, predict(Poisson(),forward_map(Poisson(),log_tc1)),tc1,0.0)
        @test check(:l2diff, predict(:Poisson,forward_map(:Poisson,log_tc1)),tc1,0.0)
        @test predict(:Poisson,[1,1];custom_prediction_function=1) == -1 #
        tc2 = rand(2:20,100,100)
        log_tc2 = log.(tc2)
        @test check(:l2diff, predict(Poisson(),forward_map(Poisson(),log_tc2)),tc2,0.0)
        @test check(:l2diff, predict(:Poisson,forward_map(:Poisson,log_tc2)),tc2,0.0)
        @test predict(:Poisson,[1 1;1 1];custom_prediction_function=1) == -1 #
    end
end

@testset "$(format("Exponential Family: predict[bernoulli]"))" begin
    let
        logit = (x) -> log.(x./(1 .- x))
        tc1 = rand(100)
        tc1_int = Int.(tc1 .> 0.5)
        logit_tc1 = logit(tc1)
        @test check(:l2diff, predict(Bernoulli(), forward_map(Bernoulli(),logit_tc1)), tc1_int, 0.0)
        @test check(:l2diff, predict(:Bernoulli, forward_map(:Bernoulli, logit_tc1)), tc1_int, 0.0)
        @test predict(:Bernoulli,[1,1];custom_prediction_function=1) == -1 #
        tc2 = rand(100,100)
        tc2_int = Int.(tc2 .> 0.5)
        logit_tc2 = logit(tc2)
        @test check(:l2diff, predict(Bernoulli(),forward_map(Bernoulli(),logit_tc2)),tc2_int,0.0)
        @test check(:l2diff, predict(:Bernoulli,forward_map(:Bernoulli,logit_tc2)),tc2_int,0.0)
        @test predict(:Bernoulli,[1 1;1 1];custom_prediction_function=1) == -1 #
    end
end

@testset "$(format("Exponential Family: predict[gaussian]"))" begin
    let
        tc1 = rand(100)
        @test check(:l2diff, predict(Gaussian(),forward_map(Gaussian(),tc1)),tc1,0.0)
        @test check(:l2diff, predict(:Gaussian,forward_map(:Gaussian,tc1)),tc1,0.0)
        @test predict(:Gaussian,[1,1];custom_prediction_function=1) == -1 #
        tc2 = rand(100,100)
        @test check(:l2diff, predict(Gaussian(),forward_map(Gaussian(),tc2)),tc2,0.0)
        @test check(:l2diff, predict(:Gaussian,forward_map(:Gaussian,tc2)),tc2,0.0)
        @test predict(:Poisson,[1 1;1 1];custom_prediction_function=1) == -1 #
    end
end

@testset "$(format("Exponential Family: predict[gaussian]"))" begin
    let
        tc1 = rand(100)
        inv_tc1 = 1 ./ tc1
        @test check(:l2diff, predict(Gamma(),forward_map(Gamma(),1 ./ tc1)),tc1,0.0)
        @test check(:l2diff, predict(:Gamma,forward_map(:Gamma,1 ./ tc1)),tc1,0.0)
        @test predict(:Gamma,[1,1];custom_prediction_function=1) == -1 #
        tc2 = rand(100,100)
        inv_tc2 = 1 ./ tc2
        @test check(:l2diff, predict(Gamma(),forward_map(Gamma(),inv_tc2)),tc2,0.0)
        @test check(:l2diff, predict(:Gamma,forward_map(:Gamma,inv_tc2)),tc2,0.0)
        @test predict(:Poisson,[1 1;1 1];custom_prediction_function=1) == -1 #
    end
end


text:
##############################################################################
##
## Type with stored qr
##
##############################################################################
struct DenseQRAllocatedSolver{Tqr <: StridedMatrix, Tu <: AbstractVector} <: AbstractAllocatedSolver
    qrm::Tqr
    u::Tu
    function DenseQRAllocatedSolver{Tqr, Tu}(qrm, u) where {Tqr <: StridedMatrix, Tu <: AbstractVector}
        length(u) == size(qrm, 1) || throw(DimensionMismatch("u must have length size(J, 1)"))
        new(qrm, u)
    end
end

function DenseQRAllocatedSolver(qrm::Tqr, u::Tu) where {Tqr <: StridedMatrix, Tu <: AbstractVector}
    DenseQRAllocatedSolver{Tqr, Tu}(qrm, u)
end

##############################################################################
##
## solve J'J \ J'y by QR
##
##############################################################################

function AbstractAllocatedSolver(nls::LeastSquaresProblem{Tx, Ty, Tf, TJ, Tg}, optimizer::Dogleg{QR}) where {Tx, Ty, Tf, TJ <: StridedVecOrMat, Tg}
    return DenseQRAllocatedSolver(similar(nls.J), _zeros(nls.y))
end

function LinearAlgebra.ldiv!(x::AbstractVector, J::StridedMatrix,  y::AbstractVector, A::DenseQRAllocatedSolver)
    u, qrm = A.u, A.qrm
    copyto!(qrm, J)
    copyto!(u, y)
    if VERSION >= v"1.7.0"
        ldiv!(qr!(qrm, ColumnNorm()), u)
    else
        ldiv!(qr!(qrm, Val(true)), u)
    end
    for i in 1:length(x)
        x[i] = u[i]
    end
    return x, 1
end

##############################################################################
##
## solve (J'J + diagm(damp)) \ J'y by QR
##
##############################################################################

function AbstractAllocatedSolver(nls::LeastSquaresProblem{Tx, Ty, Tf, TJ, Tg}, optimizer::LevenbergMarquardt{QR}) where {Tx, Ty, Tf, TJ <: StridedVecOrMat, Tg}
    qrm = zeros(eltype(nls.J), length(nls.y) + length(nls.x), length(nls.x))
    u = zeros(eltype(nls.y), length(nls.y) + length(nls.x))
    return DenseQRAllocatedSolver(qrm, u)
end

function LinearAlgebra.ldiv!(x::AbstractVector, J::StridedMatrix, y::AbstractVector,
                damp::AbstractVector, A::DenseQRAllocatedSolver, verbose::Bool = false)
    u, qrm = A.u, A.qrm

    # transform dammp
    length(u) ==  length(y) + length(x) || throw(DimensionMismatch("length(u) should equal length(x) + length(y)"))

    # update qr as |J; diagm(damp)|
    fill!(qrm, zero(eltype(qrm)))
    for j in 1:size(J, 2)
        for i in 1:size(J, 1)
            qrm[i, j] = J[i, j]
        end
    end

    leny = length(y)
    for i in 1:length(damp)
        qrm[leny + i, i] = sqrt(damp[i])
    end

    # update u as |J; 0|
    fill!(u, zero(eltype(u)))
    for i in 1:length(y)
        u[i] = y[i]
    end
    if verbose
        @show mean(u)
        sleep(1)
    end

    # solve
    if VERSION >= v"1.7.0"
        ldiv!(qr!(qrm, ColumnNorm()), u)
    else
        ldiv!(qr!(qrm, Val(true)), u)
    end
    if verbose
        @show mean(u)
        sleep(1)
    end
    for i in 1:length(x)
        x[i] = u[i]
    end
    return x, 1
end

max_stars_repo_name: matthieugomez/LeastSquares

max_issues_repo_name: matthieugomez/LeastSquares

max_forks_repo_name: matthieugomez/LeastSquares


text:
@i function isolve(out!::T, sg::AbstractSpinglass{SquareLattice, T}, reg::ArrayReg{B,TT}, A_STACK, B_STACK) where {B,T,TT<:Tropical{T}}
    @invcheckoff begin
    @routine begin
        Lx ← sg.lattice.Nx
        Ly ← sg.lattice.Ny
        Js ← sg.Js
        hs ← sg.hs
    end
    @safe println("Layer 1/$Lx, stack size: $(A_STACK.top) & $(B_STACK.top)")
    k ← 0
    for j=1:Ly
        apply_Gh!(reg, j, hs[j], A_STACK)
    end
    for j=1:Ly-1
        k += 1
        apply_Gvb!(reg, (@const (j, j+1)), Js[k], A_STACK)
    end
    for i=2:Lx
        @safe println("Layer $i/$Lx, stack size: $(A_STACK.top) & $(B_STACK.top)")
        @routine begin
            for j=1:Ly
                k += 1
                apply_Ghb!(reg, j, Js[k], A_STACK)
            end
        end
        # bookup current state and loss
        incstack!(B_STACK)
        store_state!(B_STACK, reg.state)

        # clean up `NiLang.GLOBAL_STACK`
        ~@routine
        # but wait, `k` should not be uncomputed
        k += Ly

        # restore the state
        swap_state!(B_STACK, reg.state)

        for j=1:Ly
            apply_Gh!(reg, j, hs[(i-1)*Ly+j], A_STACK)
        end
        for j=1:Ly-1
            k += 1
            apply_Gvb!(reg, (@const (j, j+1)), Js[k], A_STACK)
        end
    end
    summed ← one(TT)
    i_sum(summed, reg.state)
    NiLang.SWAP(summed.n, out!)
    k → length(Js)
    summed → one(TT)
    ~@routine
    end
end

@i function swap_state!(B_STACK, state)
    @invcheckoff for j=1:length(state)
        @inbounds NiLang.SWAP(state[j], B_STACK[j])
    end
end

@i function store_state!(B_STACK, state)
    @invcheckoff for j = 1:length(state)
        @inbounds TropicalYao.Reversible.unsafe_store!(B_STACK[j], state[j])
    end
end

@i function isolve_largemem(out!, sg::Spinglass{SquareLattice, T}, reg::ArrayReg{B,TT}, REG_STACK) where {B,T,TT<:Tropical{T}}
    @invcheckoff begin
    @routine begin
        Lx ← sg.lattice.Nx
        Ly ← sg.lattice.Ny
        Js ← sg.Js
        hs ← sg.hs
    end
    @safe println("Layer 1/$Lx")
    k ← 0
    for j=1:Ly
        apply_Gh!(reg, j, hs[j], REG_STACK)
    end
    for j=1:Ly-1
        k += 1
        apply_Gvb!(reg, (@const (j, j+1)), Js[k], REG_STACK)
    end
    for i=2:Lx
        @safe println("Layer $i/$Lx")
        for j=1:Ly
            k += 1
            apply_Ghb!(reg, j, Js[k], REG_STACK)
        end
        for j=1:Ly
            apply_Gh!(reg, j, hs[(i-1)*Ly+j], REG_STACK)
        end
        for j=1:Ly-1
            k += 1
            apply_Gvb!(reg, (@const (j, j+1)), Js[k], REG_STACK)
        end
    end
    summed ← one(TT)
    i_sum(summed, reg.state)
    NiLang.SWAP(summed.n, out!)
    k → length(Js)
    summed → one(TT)
    ~@routine
    end
end

cachesize_A(lt::SquareLattice) = lt.Ny
cachesize_B(lt::SquareLattice) = lt.Nx-1
cachesize_largemem(lt::SquareLattice) = (lt.Nx-1) * lt.Ny


text:
using WordCloud
using HTTP

url = "http://en.wikipedia.org/wiki/Special:random"
try
    resp = HTTP.request("GET", url, redirect=true)
    println(resp.request)
    content = resp.body |> String
    wc = wordcloud(content |> html2text |> processtext) |> generate!
    println("results are saved to fromweb.png")
    paint(wc, "fromweb.png")
    wc
catch e
    println(e)
end
#eval# runexample(:fromweb)
#md# ![](fromweb.png)


text:
using Random
using Dates
using JuMP, BilevelJuMP

using MathOptInterface
MOI = MathOptInterface
MOIU = MOI.Utilities

MAX_TIME = 600#600

include("svr.jl")
include("rand.jl")
include("forecast.jl")
include("toll.jl")

mode = BilevelJuMP.SOS1Mode()

with_att = JuMP.optimizer_with_attributes

FA = BilevelJuMP.FortunyAmatMcCarlMode

using Xpress
using CPLEX
using Gurobi
using SCIP
using Cbc
using GLPK
using Ipopt
using KNITRO
using AmplNLWriter
using Mosek, MosekTools
using Juniper
using QuadraticToBinary

QB = QuadraticToBinary.Optimizer{Float64}
cache(opt) = MOIU.CachingOptimizer(
    MOIU.UniversalFallback(MOIU.Model{Float64}()), opt)

function cpx()
    s=CPLEX.Optimizer()
    MOI.set(s, MOI.RawParameter("CPXPARAM_TimeLimit"),MAX_TIME*1)
    s
end

# const KN_OPT = KNITRO.Optimizer(maxtimecpu = MAX_TIME*1.0)
# function new_knitro()
#     MOI.empty!(KN_OPT)
#     return KN_OPT
# end

SOLVERS = [
    #=
    =#
    #=
        SOS1
    =#
    # (with_att(Gurobi.Optimizer, "TimeLimit" => MAX_TIME*1), BilevelJuMP.SOS1Mode(), "gurobi_sos1"),
    # (with_att(CPLEX.Optimizer, "CPXPARAM_TimeLimit" => MAX_TIME*1), BilevelJuMP.SOS1Mode(), "cplex_sos1"),
    # (with_att(Xpress.Optimizer, "MAXTIME" => -MAX_TIME*1, "logfile" => "output.log"), BilevelJuMP.SOS1Mode(), "xpress_sos1"),
    # (with_att(Cbc.Optimizer, "seconds" => MAX_TIME*1.0), BilevelJuMP.SOS1Mode(), "cbc_sos1"),
    # (with_att(SCIP.Optimizer, "limits/time" => MAX_TIME*1), BilevelJuMP.SOS1Mode(), "scip_sos1"),
    #=
        indicator
    =#
    # (with_att(CPLEX.Optimizer, "CPXPARAM_TimeLimit" => MAX_TIME), BilevelJuMP.IndicatorMode(), "cplex_indc"),
    # # (with_att(Gurobi.Optimizer, "TimeLimit" => MAX_TIME), BilevelJuMP.IndicatorMode(), "gurobi_indc"), # not supporting indicator
    # (with_att(Xpress.Optimizer, "MAXTIME" => -MAX_TIME, "logfile" => "output.log"), BilevelJuMP.IndicatorMode(), "xpress_indc"),
    # (with_att(Cbc.Optimizer, "seconds" => MAX_TIME*1.0), BilevelJuMP.IndicatorMode(), "cbc_indc"),
    # # (with_att(SCIP.Optimizer, "limits/time" => MAX_TIME), BilevelJuMP.IndicatorMode(), "scip_indc"), # weird offset error
    #=
        Fortuny-Amat 10
    =#
    # # (with_att(GLPK.Optimizer, "tm_lim" => MAX_TIME * 1_000), FA(primal_big_M = 10, dual_big_M = 10), "glpk_fa10"),
    # # (with_att(Mosek.Optimizer, "MIO_MAX_TIME" => MAX_TIME * 1.0, "OPTIMIZER_MAX_TIME" => MAX_TIME * 1.0), FA(primal_big_M = 10, dual_big_M = 10), "mosek_fa10"), # no sos1
    # # (with_att(Gurobi.Optimizer, "TimeLimit" => MAX_TIME*1), FA(primal_big_M = 10, dual_big_M = 10), "gurobi_fa10"),
    # (with_att(CPLEX.Optimizer, "CPXPARAM_TimeLimit" => MAX_TIME*1), FA(primal_big_M = 10, dual_big_M = 10), "cplex_fa10"), #TODO
    # (with_att(Xpress.Optimizer, "MAXTIME" => -MAX_TIME*1, "logfile" => "output.log"), FA(primal_big_M = 10, dual_big_M = 10), "xpress_fa10"),
    # (with_att(Cbc.Optimizer, "seconds" => MAX_TIME*1.0), FA(primal_big_M = 10, dual_big_M = 10), "cbc_fa10"),
    # (with_att(SCIP.Optimizer, "limits/time" => MAX_TIME*1), FA(primal_big_M = 10, dual_big_M = 10), "scip_fa10"),
    #=
        Fortuny-Amat 100
    =#
    # (with_att(GLPK.Optimizer, "tm_lim" => MAX_TIME * 1_000), FA(primal_big_M = 100, dual_big_M = 100), "glpk_fa100"),
    # (with_att(Mosek.Optimizer, "MIO_MAX_TIME" => MAX_TIME * 1.0, "OPTIMIZER_MAX_TIME" => MAX_TIME * 1.0), FA(primal_big_M = 100, dual_big_M = 100), "mosek_fa100"), # no sos1
    # (with_att(Gurobi.Optimizer, "TimeLimit" => MAX_TIME*1), FA(primal_big_M = 100, dual_big_M = 100), "gurobi_fa100"),
    # (with_att(CPLEX.Optimizer, "CPXPARAM_TimeLimit" => MAX_TIME*1), FA(primal_big_M = 100, dual_big_M = 100), "cplex_fa100"),
    # (with_att(Xpress.Optimizer, "MAXTIME" => -MAX_TIME*1, "logfile" => "output.log"), FA(primal_big_M = 100, dual_big_M = 100), "xpress_fa100"),
    # (with_att(Cbc.Optimizer, "seconds" => MAX_TIME*1.0), FA(primal_big_M = 100, dual_big_M = 100), "cbc_fa100"),
    # (with_att(SCIP.Optimizer, "limits/time" => MAX_TIME*1), FA(primal_big_M = 100, dual_big_M = 100), "scip_fa100"),
    #=
        Product NLP
    =#
    # (with_att(Ipopt.Optimizer, "max_cpu_time" => MAX_TIME*1.0), BilevelJuMP.ProductMode(1e-7), "ipopt_prod"),
    # (new_knitro(), BilevelJuMP.ProductMode(1e-7), "knitro_prod"),
    #=
        Product BIN 10
    =#
    # (()->QB(GLPK.Optimizer(tm_lim=MAX_TIME*1_000),lb=-10,ub=10), BilevelJuMP.ProductMode(1e-7), "glpk_prod10"),
    # (()->QB(Mosek.Optimizer(MIO_MAX_TIME=MAX_TIME*1.0,OPTIMIZER_MAX_TIME=MAX_TIME*1.0),lb=-10,ub=10), BilevelJuMP.ProductMode(1e-7), "mosek_prod10"),
    # (()->QB(Gurobi.Optimizer(TimeLimit=MAX_TIME*1),lb=-10,ub=10), BilevelJuMP.ProductMode(1e-7), "gurobi_prod10"),
    # (()->QB(cpx(),lb=-10,ub=10), BilevelJuMP.ProductMode(1e-7), "cplex_prod10"), #TODO
    # (()->QB(Xpress.Optimizer(MAXTIME=-MAX_TIME*1),lb=-10,ub=10), BilevelJuMP.ProductMode(1e-7), "xpress_prod10"),
    # (()->QB(cache(Cbc.Optimizer(seconds=MAX_TIME*1.0)),lb=-10,ub=10), BilevelJuMP.ProductMode(1e-7), "cbc_prod10"),
    # (()->QB(SCIP.Optimizer(limits_time=MAX_TIME*1),lb=-10,ub=10), BilevelJuMP.ProductMode(1e-7), "scip_prod10"),
    #=
        Product BIN 100
    =#
    # (()->QB(GLPK.Optimizer(tm_lim=MAX_TIME*1_000),lb=-100,ub=100), BilevelJuMP.ProductMode(1e-7), "glpk_prod100"),
    # (()->QB(Mosek.Optimizer(MIO_MAX_TIME=MAX_TIME*1.0,OPTIMIZER_MAX_TIME=MAX_TIME*1.0),lb=-100,ub=100), BilevelJuMP.ProductMode(1e-7), "mosek_prod100"),
    # (()->QB(Gurobi.Optimizer(TimeLimit=MAX_TIME*1),lb=-100,ub=100), BilevelJuMP.ProductMode(1e-7), "gurobi_prod100"),
    # (()->QB(cpx(),lb=-100,ub=100), BilevelJuMP.ProductMode(1e-7), "cplex_prod100"),
    # (()->QB(Xpress.Optimizer(MAXTIME=-MAX_TIME*1),lb=-100,ub=100), BilevelJuMP.ProductMode(1e-7), "xpress_prod100"),
    # (()->QB(cache(Cbc.Optimizer(seconds=MAX_TIME*1.0)),lb=-100,ub=100), BilevelJuMP.ProductMode(1e-7), "cbc_prod100"),
    # (()->QB(SCIP.Optimizer(limits_time=MAX_TIME*1),lb=-100,ub=100), BilevelJuMP.ProductMode(1e-7), "scip_prod100"),
    #=
        PrimalDual BIN 10
    =#
    # (()->QB(GLPK.Optimizer(tm_lim=MAX_TIME*1_000),lb=-10,ub=10), BilevelJuMP.StrongDualityEqualityMode(), "glpk_sd10"),
    # (()->QB(Mosek.Optimizer(MIO_MAX_TIME=MAX_TIME*1.0,OPTIMIZER_MAX_TIME=MAX_TIME*1.0),lb=-10,ub=10), BilevelJuMP.StrongDualityEqualityMode(), "mosek_sd10"),
    # (()->QB(Gurobi.Optimizer(TimeLimit=MAX_TIME*1),lb=-10,ub=10), BilevelJuMP.StrongDualityEqualityMode(), "gurobi_sd10"),
    # (()->QB(cpx(),lb=-10,ub=10), BilevelJuMP.StrongDualityEqualityMode(), "cplex_sd10"),
    # (()->QB(Xpress.Optimizer(MAXTIME=-MAX_TIME*1),lb=-10,ub=10), BilevelJuMP.StrongDualityEqualityMode(), "xpress_sd10"),
    # (()->QB(cache(Cbc.Optimizer(seconds=MAX_TIME*1.0)),lb=-10,ub=10), BilevelJuMP.StrongDualityEqualityMode(), "cbc_sd10"), #TODO
    (()->QB(SCIP.Optimizer(limits_time=MAX_TIME*1),lb=-10,ub=10), BilevelJuMP.StrongDualityEqualityMode(), "scip_sd10"),
    #=
        PrimalDual BIN 100
    =#
    (()->QB(GLPK.Optimizer(tm_lim=MAX_TIME*1_000),lb=-100,ub=100), BilevelJuMP.StrongDualityEqualityMode(), "glpk_sd100"),
    (()->QB(Mosek.Optimizer(MIO_MAX_TIME=MAX_TIME*1.0,OPTIMIZER_MAX_TIME=MAX_TIME*1.0),lb=-100,ub=100), BilevelJuMP.StrongDualityEqualityMode(), "mosek_sd100"),
    (()->QB(Gurobi.Optimizer(TimeLimit=MAX_TIME*1),lb=-100,ub=100), BilevelJuMP.StrongDualityEqualityMode(), "gurobi_sd100"),
    (()->QB(cpx(),lb=-100,ub=100), BilevelJuMP.StrongDualityEqualityMode(), "cplex_sd100"),
    (()->QB(Xpress.Optimizer(MAXTIME=-MAX_TIME*1),lb=-100,ub=100), BilevelJuMP.StrongDualityEqualityMode(), "xpress_sd100"),
    (()->QB(cache(Cbc.Optimizer(seconds=MAX_TIME*1.0)),lb=-100,ub=100), BilevelJuMP.StrongDualityEqualityMode(), "cbc_sd100"),
    (()->QB(SCIP.Optimizer(limits_time=MAX_TIME*1),lb=-100,ub=100), BilevelJuMP.StrongDualityEqualityMode(), "scip_sd100"),
    #=
        Complemets
    =#
    # #(with_att(KNITRO.Optimizer, "maxtimecpu" => MAX_TIME*1.0), BilevelJuMP.ComplementMode(), "knitro_comp"),
    #=
        Product global
    =#
    # # (() -> AmplNLWriter.Optimizer("bonmin"), BilevelJuMP.ProductMode(1e-5)),
    # # (() -> AmplNLWriter.Optimizer("couenne"), BilevelJuMP.ProductMode(1e-5)),
]

PROBLEMS = [
    :SVR,
    # :RAND,
    :TOLL,
    :FORECAST,
]

SEEDS = [
    1234,
    # 2345,
    # 3456,
    # 4567,
    # 5678,
    # 6789,
    # 7890,
    # 8901,
    # 9012,
    # 0123,
]

SVR = [
    # (features, sample_size)
    (  1,  10),
    (  1,  10),
    (  2,  10),
    (  5,  10),
    (  1, 100),
    (  2, 100), # 600
    (  5, 100),
    ( 10, 100),
    ( 20, 100),
    ( 50, 100),
    # (  1,1000), # hard for prod10
    # (  2,1000),
    # (  5,1000),
    # # # ( 10,1000),
    # # # ( 20,1000),
    # # # ( 50,1000),
    # # # (100,1000),
    # # # (200,1000),
    # # # (500,1000),
    # # () for i in [1,2,5,10,20,50,100,200], j in [10, 100, 1000]
]

RAND = [
    # (rows, cols)
    (   5,   5),
    (  10,   5),
    (   5,  10),
    (  10,  10),
    (  50,  10),
    (  10,  50),
    (  50,  50), # 600
    ( 100,  50),
    (  50, 100),
    ( 100, 100),
    # # ( 500, 100),
    # # ( 100, 500),
    # # ( 500, 500),
    # # (1000, 500),
    # # ( 500,1000),
    # # (1000,1000),
    # # (5000,1000),
    # # (1000,5000),
    # # (5000,5000),
]

TOLL = [
    # nodes
    5,
    10,
    20,
    # 50, # hard for prod10
    # 100,
    # 200, # also massive on memory
    # 500, # 600 - broke gurobi
    # # 1000,
    # # 2000,
    # # 5000,
]

FORECAST = [
    # (products, sample_size)
    (  1,  10),
    (  2,  10),
    (  5,  10),
    (  1, 100),
    # (  2, 100), # hard for prod10
    # (  5, 100), # 600
    # ( 10, 100),
    # ( 20, 100),
    # ( 50, 100),
    # (  1,1000),
    # (  2,1000),
    # (  5,1000),
    # # ( 10,1000),
    # # ( 20,1000),
    # # ( 50,1000),
    # # (100,1000),
    # # (200,1000),
    # # (500,1000),
    # # () for i in [1,2,5,10,20,50,100,200], j in [10, 100, 1000]
]

function separator()
    println()
    println()
    println("============================================================")
    println("============================================================")
    println()
    println()
end

function new_file()
    cd(dirname(@__FILE__))
    FILE = open("bench$(replace("$(now())",":"=>"_")).log", "w")
    println(FILE, "opt_mode, prob, inst, primal_status, termination_status, solve_time, build_time, lower_obj, upper_obj")
    flush(FILE)
    return FILE
end
function newline(FILE, data, opt, prb, inst, seed)
    println(FILE, "$opt, $prb, $inst, $seed, $(data[1]),$(data[2]),$(data[3]),$(data[4]),$(data[5]),$(data[6]),$(data[7])")
    flush(FILE)
end
FILE = new_file()
for seed in SEEDS
    for (optimizer, mode, name) in SOLVERS
        if :SVR in PROBLEMS
            for (features, samples) in SVR
                separator()
                @show features, samples, seed, name
                separator()
                ret = bench_svr(features, samples, optimizer, mode, seed)
                newline(FILE, ret, name, :SVR, (features, samples), seed)
            end
        end
        if :RAND in PROBLEMS
            for (rows, cols) in RAND
                separator()
                @show rows, cols, seed, name
                separator()
                ret = bench_rand(rows, cols, 0.5, optimizer, mode, seed)
                newline(FILE, ret, name, :RAND, (rows, cols), seed)

            end
        end
        if :TOLL in PROBLEMS
            for nodes in TOLL
                separator()
                @show nodes, seed, name
                separator()
                ret = bench_toll(nodes, optimizer, mode, seed)
                newline(FILE, ret, name, :TOLL, (nodes,), seed)
            end
        end
        if :FORECAST in PROBLEMS
            for (products, samples) in FORECAST
                separator()
                @show products, samples, seed, name
                separator()
                ret = bench_forecast(products, samples, optimizer, mode, seed)
                newline(FILE, ret, name, :FORE, (products, samples), seed)
            end
        end
    end
end

close(FILE)

exit(0)


text:
using MultivariateOrthogonalPolynomials, ArrayLayouts, BandedMatrices, BlockArrays, Test
import MultivariateOrthogonalPolynomials: ModalInterlace, ModalInterlaceLayout, ModalTrav

@testset "modalTrav" begin
    a = ModalTrav(randn(2,5))
    b = PseudoBlockArray(a)
    v = Vector(a)

    @test zero(a) isa ModalTrav
    @test zero(a) == zero(b)

    @test exp.(a) isa ModalTrav
    @test 2a isa ModalTrav
    @test a+a isa ModalTrav
    @test a+b isa PseudoBlockArray
    @test a+v isa BlockArray
    @test a .+ exp.(a .+ 1) isa ModalTrav

    @test exp.(a) == exp.(b)
    @test a + a == 2a == a+b

    @test a .+ exp.(a .+ 1) == b .+ exp.(b .+ 1)

    for k = 1:6
        a[k] = k
    end
    @test a == 1:6

    m = ModalTrav(reshape(1:10, 2, 5))
    @test m[Block(3)] == [2,7,9]
    @test m == [1,3,5,2,7,9]

    @test copy(m) isa ModalTrav{Int,Matrix{Int}}
    @test copy(m) == m
end

@testset "ModalInterlace" begin
    ops = [brand(2,3,1,2), brand(1,2,1,1), brand(1,2,1,2)]
    A = ModalInterlace(ops, (3,5), (2,4))
    @test MemoryLayout(A) isa ModalInterlaceLayout
    @test A[[1,4],[1,4,11]] == ops[1]
    @test A[[2],[2,7]] == ops[2]
    @test A[[5],[5,12]] == A[[6],[6,13]] == ops[3]

    b = ModalTrav(1:15)
    @test A*b ≈ Matrix(A) * Vector(b)

    ops = [brand(3,3,1,2), brand(2,2,1,1), brand(2,2,1,2), brand(1,1,1,1), brand(1,1,1,2)]
    B = ModalInterlace(ops, (5,5), (2,4))
    @test B\b ≈ Matrix(B) \ Vector(b)
end