text:
using Compat
using Requests
using JSON
using Base.Test
using Libz

import Requests: get, post, put, delete, options, bytes, text, json, history

# simple calls, no headers, data or query params -------

@test get("http://httpbin.org/get").status == 200
@test post("http://httpbin.org/post").status == 200
@test put("http://httpbin.org/put").status == 200
@test delete("http://httpbin.org/delete").status == 200
@test options("http://httpbin.org/get").status == 200

# check query params -------

data = json(get("http://httpbin.org/get";
                      query = Dict("key1" => "value1",
                                   "key with spaces" => "value with spaces",
                                   "key2" => ["value2.1", "value2.2"])))
@test data["args"]["key1"] == "value1"
@test data["args"]["key with spaces"] == "value with spaces"
@test data["args"]["key2"] == ["value2.1", "value2.2"]

data = json(post("http://httpbin.org/post";
                       query = Dict("key1" => "value1",
                                            "key2" => "value2",
                                            "key with spaces" => "value with spaces")))
@test data["args"]["key1"] == "value1"
@test data["args"]["key2"] == "value2"
@test data["args"]["key with spaces"] == "value with spaces"

data = json(put("http://httpbin.org/put";
                      query = Dict("key1" => "value1",
                                   "key2" => "value2",
                                   "key3" => 3,
                                   "key with spaces" => "value with spaces")))
@test data["args"]["key1"] == "value1"
@test data["args"]["key2"] == "value2"
@test data["args"]["key3"] == "3"
@test data["args"]["key with spaces"] == "value with spaces"

data = json(delete("http://httpbin.org/delete";
                         query = Dict("key1" => "value1",
                                      "key4" => 4.01,
                                      "key with spaces" => "value with spaces")))
@test data["args"]["key1"] == "value1"
@test data["args"]["key4"] == "4.01"
@test data["args"]["key with spaces"] == "value with spaces"

response = options("http://httpbin.org/get";
                          query = Dict("key1" => "value1",
                                       "key2" => "value2",
                                       "key3" => 3,
                                       "key4" => 4.01))
@test length(response.data) == 0

# check data -------

data = json(post("http://httpbin.org/post";
                       json = Dict("key1" => "value1",
                                    "key2" => "value2")))
@test data["json"]["key1"] == "value1"
@test data["json"]["key2"] == "value2"

data = json(put("http://httpbin.org/put";
                      json = Dict("key1" => "value1",
                                  "key2" => "value2",
                                  "key3" => 3)))
@test data["json"]["key1"] == "value1"
@test data["json"]["key2"] == "value2"
@test data["json"]["key3"] == 3

data = json(delete("http://httpbin.org/delete";
                         json = Dict("key1" => "value1",
                                    "key4" => 4.01)))
@test data["json"]["key1"] == "value1"
@test data["json"]["key4"] == 4.01

# form-encoded posts
data = json(post("http://httpbin.org/post", data=Dict("email"=>"a@b.com")))
@test data["form"]["email"] == "a@b.com"

# query + data -------

data = json(post("http://httpbin.org/post";
                       query = (Dict("qkey1" => "value1",
                                     "qkey2" => "value2")),
                       json = (Dict("dkey1" => "data1",
                                    "dkey2" => "data2"))))
@test data["args"]["qkey1"] == "value1"
@test data["args"]["qkey2"] == "value2"
@test data["json"]["dkey1"] == "data1"
@test data["json"]["dkey2"] == "data2"

data = json(put("http://httpbin.org/put";
                      query = (Dict("qkey1" => "value1",
                                    "qkey2" => "value2",
                                    "qkey3" => 3)),
                      json = (Dict("dkey1" => "data1",
                                    "dkey2" => "data2",
                                    "dkey3" => 5))))
@test data["args"]["qkey1"] == "value1"
@test data["args"]["qkey2"] == "value2"
@test data["args"]["qkey3"] == "3"
@test data["json"]["dkey1"] == "data1"
@test data["json"]["dkey2"] == "data2"
@test data["json"]["dkey3"] == 5

data = json(delete("http://httpbin.org/delete";
                         query = (Dict("qkey1" => "value1",
                                       "qkey4" => 4.01)),
                         json = (Dict("dkey1" => "data1",
                                      "dkey2" => 9.01))))
@test data["args"]["qkey1"] == "value1"
@test data["args"]["qkey4"] == "4.01"
@test data["json"]["dkey1"] == "data1"
@test data["json"]["dkey2"] == 9.01

data = json(post(URI("http://httpbin.org/post");
                       data = "√",
                       headers = Dict("Content-Type" => "text/plain")))

@test data["data"] == "√"

# Test custom content type
data = json(post("http://httpbin.org/post",
             data="{\"a\": \"b\"}",
             headers=Dict("Content-Type"=>"application/json")))
@test data["json"]["a"] == "b"

# Test file upload
filename = Base.source_path()

files = [
  FileParam(readstring(filename),"text/julia","file1","runtests.jl"),
  FileParam(open(filename,"r"),"text/julia","file2","runtests.jl",true),
  FileParam(IOBuffer(readstring(filename)),"text/julia","file3","runtests.jl"),
  ]

res = post(URI("http://httpbin.org/post"); files = files)

filecontent = readstring(filename)
data = json(res)
@test data["files"]["file1"] == filecontent
@test data["files"]["file2"] == filecontent
@test data["files"]["file3"] == filecontent

# Test for chunked responses (we expect 100 from split as there are 99 '\n')
@test size(split(text(get("http://httpbin.org/stream/99")), "\n"), 1) == 100

# Test for gzipped responses
@test json(get("http://httpbin.org/gzip"))["gzipped"] == true
@test json(get("http://httpbin.org/deflate"))["deflated"] == true
@test json(get("http://httpbin.org/headers"; compressed=true))["headers"]["Accept-Encoding"] == "gzip, deflate"

# Test for gzipped requests
foo_result = "data:application/octet-stream;base64," * (Vector{UInt8}("foo") |> ZlibDeflateInputStream |> read |> base64encode)
@test json(post("http://httpbin.org/post"; data="foo", gzip_data=true))["data"] == foo_result
@test json(post("http://httpbin.org/post"; data=Vector{UInt8}("foo"), gzip_data=true))["data"] == foo_result

# Test timeout delay
let
    timeout = Dates.Millisecond(500)
    @test_throws Requests.TimeoutException get("http://httpbin.org/delay/3", timeout=timeout)
end

# Test cookies
let
    r = get("http://httpbin.org/cookies/set?a=1&b=2", allow_redirects=false)
    cookies = r.cookies
    @test length(cookies) == 2
    @test cookies["a"].value == "1"
    @test cookies["b"].value == "2"
    @test cookies["a"].attrs["Path"] == "/"
    r = json(get("http://httpbin.org/cookies", cookies=cookies))
    @test r["cookies"]["a"] == "1"
    @test r["cookies"]["b"] == "2"
end

# Test redirects
let
    r = get("http://httpbin.org/absolute-redirect/3")
    @test length(history(r)) == 3
    r = get("http://httpbin.org/relative-redirect/3")
    @test length(history(r)) == 3
    @test_throws Requests.RedirectException get("http://httpbin.org/redirect/3", max_redirects=1)
end

# Test HTTPS
@test statuscode(get("https://httpbin.org")) == 200

#=
# Test output streaming
let
    stream = Requests.post_streaming(
      "http://httpbin.org/post", write_body=false,
      headers=Dict("Transfer-Encoding"=>"chunked"))

    write_chunked(stream, "ab")
    write_chunked(stream, "cde")
    write_chunked(stream, "")

    response = JSON.parse(readstring(stream))
    @test response["data"] == "abcde"
end
=#

# Test input streaming
let
    stream = Requests.get_streaming("http://httpbin.org/stream-bytes/100", query=Dict(:chunk_size=>10))
    N = 0
    while !eof(stream)
        bytes = readavailable(stream)
        N += length(bytes)
    end
    @test N==100
end

# Test that streaming input doesn't produce an early EOF
let
    num_bytes = 0
    total_bytes = 128
    stream = Requests.get_streaming("http://httpbin.org/stream-bytes/$total_bytes", query=Dict(:chunk_size=>1))

    is_eof = false
    num_iterations = 0
    while !is_eof && num_bytes < total_bytes
        num_bytes += length(readavailable(stream))
        is_eof = eof(stream)
        num_iterations += 1
    end

    # This test requires that input could not be read all at one time. We may need to
    # increase the `total_bytes` to ensure multiple read iterations occur.
    @test num_iterations > 1

    @test is_eof == true
    @test nb_available(stream) == 0
    @test num_bytes == total_bytes
end

# Proxy testing. Would be better to use a real proxy instead of using the real site
# as the "proxy", but hard to find a reliable unmetered public proxy for testing.

@test get("http://httpbin.org/get"; proxy=Nullable(URI("http://httpbin.org"))).status == 200

# Test proxy settings via local squid Docker container if REQUESTS_TEST_PROXY is set to 1
# before testing.
# Requires Docker and docker-machine (obtainable via Docker toolbox)
if get(ENV, "REQUESTS_TEST_PROXY", "0") == "1"
    run(`docker-machine create -d virtualbox proxytest`)
    cmds = readstring(`docker-machine env proxytest`)
    for line in split(cmds, '\n')
        m = match(r"export (?<name>.*?)=(?<val>.*)", line)
        m===nothing && continue
        ENV[m[:name]] = m[:val]
    end
    @show ENV
    run(`docker run --name squid -d --restart=always -p 3128:3128 quay.io/sameersbn/squid:3.3.8-2`)
    ip = IPv4(readstring(`docker-machine ip proxytest`))
    proxy_vars = ["http_proxy", "https_proxy"]
    for var in proxy_vars
        ENV[var] = "http://$ip:3128"
    end
    ENV["REQUESTS_TEST_PROXY"] = "0"
    Pkg.test("Requests")
    ENV["REQUESTS_TEST_PROXY"] = "1"
    for var in proxy_vars
        pop!(ENV, var)
    end
    run(`docker stop squid`)
    run(`docker rm -v squid`)
    run(`docker-machine rm proxytest`)
end

# Basic auth
@test get("https://user:pw@httpbin.org/basic-auth/user/pw").status == 200
@test get("https://user:wrongpw@httpbin.org/basic-auth/user/pw").status == 401


text:
using SafeTestsets
using Metatheory
using Test

Metatheory.options.verbose = false
Metatheory.options.printiter = false

@metatheory_init

@timev begin
   @testset "All Tests" begin
      @safetestset "EGraphs Basics" begin include("test_egraphs.jl") end
      @safetestset "EMatch" begin include("test_ematch.jl") end
      @safetestset "EMatch Assertions" begin include("test_ematch_assertions.jl") end
      @safetestset "EGraph Analysis" begin include("test_analysis.jl") end
      @safetestset "EGraph Extraction" begin include("test_extraction.jl") end
      @safetestset "EGraphs Dynamic Rules" begin include("test_dynamic_ematch.jl") end
      @safetestset "Mu Puzzle" begin include("test_mu.jl") end
      @safetestset "Boson" begin include("test_boson.jl") end
      @safetestset "While Interpreter" begin include("test_while_interpreter.jl") end
      @safetestset "Classical Rewriting" begin include("test_reductions.jl") end
      @safetestset "Taylor Series" begin include("test_taylor.jl") end
      @safetestset "While Superinterpreter" begin include("test_while_superinterpreter.jl") end
      @safetestset "EGraphs Inequalities" begin include("test_inequality.jl") end
      @safetestset "PatEquiv" begin include("test_patequiv.jl") end
      @safetestset "Custom Types" begin include("test_custom_types.jl") end

      # @testset "EGraphs Multipattern" begin include("test_multipat.jl") end
      # @testset "PatAllTerm" begin include("test_patallterm.jl") end
      # use cases

      @safetestset "Fibonacci" begin include("fib/test_fibonacci.jl") end
      @safetestset "Calculational Logic" begin include("logic/test_calculational_logic.jl") end
      @safetestset "PROP Logic" begin include("logic/test_logic.jl") end
      @safetestset "CAS Infer" begin include("cas/test_infer.jl") end

      # @testset "CAS" begin include("cas/test_cas.jl") end
      # @safetestset "Categories" begin include("category/test_cat.jl") end

      @safetestset "Knuth Bendix Alternative Hurwitz Groups" begin include("group/test_kb_benchmark.jl") end

      # @testset "Proofs" begin include("proof/test_proof.jl") end
      # TODO n-ary splatvar
   end
end

# exported consistency test
for m ∈ [Metatheory, Metatheory.Util, Metatheory.EGraphs, Metatheory.EGraphs.Schedulers]
   for i ∈ propertynames(m)
      xxx = getproperty(m, i)
   end
end

ext: jl


text:
const YELLOWKEYS =["Comdty", "Corp", "Curncy", "Equity", "Govt", "Index", "M-Mkt", "Mtge", "Muni", "Pfd"]

function sedolsum(x::AbstractString)::Int
    weights = [1, 3, 1, 7, 3, 9, 1]
    s = 0
    for (w, c) in zip(weights, x)
        s += w * parse(Int, c, base=36)
    end
    return s
end

function issedol(x::AbstractString)::Bool
    return length(x) == 7 && all(isascii, x) && sedolsum(x) % 10 == 0
end

function cusipsum(x::AbstractString)::Int
    s = 0
    for (i, c) in enumerate(x)
        if isdigit(c) || isletter(c)
            v = parse(Int, c, base=36)
        elseif c == '*'
            v = 36
        elseif c == '@'
            v = 37
        elseif c == '#'
            v = 38
        end
        if iseven(i); v *= 2 end
        s += div(v, 10) + rem(v, 10)
    end
    return s
end

function iscusip(x::AbstractString)::Bool
    return length(x) == 9 && all(isascii, x) && cusipsum(x) % 10 == 0
end

function luhntest(x::Integer)::Bool
    (sum(digits(x)[1:2:end]) + sum(map(x->sum(digits(x)), 2 * digits(x)[2:2:end]))) % 10 == 0
end

luhntest(x::AbstractString) = luhntest(parse(Int, x))

function isisin(x::AbstractString)::Bool
    return length(x) == 12 && all(map(c -> isdigit(c) || isletter(c), collect(x[3:end]))) && all(isletter, x[1:2]) && parse.(Int, collect(x), base=36) |> join |> luhntest
end

function isticker(x::AbstractString)::Bool
    xs = split(x)
    return length(xs) > 1 && titlecase(xs[end]) in YELLOWKEYS
end

function isfigi(x::AbstractString)::Bool
    return length(x) == 12 && all(isletter, x[1:3]) && x[3] == 'G'
end

function identifiertype(x::AbstractString)::DataType
    x = strip(x)
    if issedol(x); return Sedol
    elseif iscusip(x); return Cusip
    elseif isfigi(x); return Figi
    elseif isisin(x); return Isin
    else; return Ticker
    end
end

max_stars_repo_name: robfs/figi

max_issues_repo_name: robfs/figi

max_forks_repo_name: robfs/figi


text:
# Color mapping

import .LibCImPlot.GetColormapSize

function SetColormap(colormap, samples::Integer = 1)
    LibCImPlot.SetColormapPlotColormap(colormap, Cint(samples))
end

function SetColormap(colors::ImVec4, num_colors::Integer)
    LibCImPlot.SetColormapVec4Ptr(colors, Cint(num_colors))
end

function GetColormapColor(index::Integer)
    out = Ref{ImVec4}()
    LibCImPlot.GetColormapColor(out, Cint(index))
    return out[]
end

function LerpColormap(t::Integer)
    out = Ref{ImVec4}()
    LibCImPlot.LerpColormap(out, Cint(t))
    return out[]
end

function PushColormap(colormap)
    LibCImPlot.PushColormapPlotColormap(colormap)
end

function PushColormap(colormap::Vector{ImVec4}, size::Integer)
    LibCImPlot.PushColormapVec4Ptr(colormap, size)
end

function PopColormap(count::Integer = 1)
    LibCImPlot.PopColormap(count)
end


text:
import Rotations: lmult, vmat, hmat

"""
    Satellite

A simple satellite model, modeling only the attitude dynamics with no external forces. Has
direct control over torques about the body axes.

# Constructor
    Satellite()
    Satellite(J)

where `J` is a `Diagonal` inertia matrix.
"""
@autodiff struct Satellite <: LieGroupModel
    J::Diagonal{Float64,SVector{3,Float64}}
end

Satellite() = Satellite(Diagonal(@SVector ones(3)))

RobotDynamics.control_dim(::Satellite) = 3
Base.position(::Satellite, x) = @SVector zeros(3)
orientation(::Satellite, x) = UnitQuaternion(x[4], x[5], x[6], x[7])

RobotDynamics.LieState(::Satellite) = RobotDynamics.LieState(UnitQuaternion{Float64}, (3,0))

function dynamics(model::Satellite, x::SVector, u::SVector)
    ω = @SVector [x[1], x[2], x[3]]
    q = normalize(@SVector [x[4], x[5], x[6], x[7]])
    J = model.J

    ωdot = J\(u - ω × (J*ω))
    qdot = 0.5*lmult(q)*hmat()*ω
    return [ωdot; qdot]
end


text:
using Underscores
using Utils

function parse_file(f)
    re::Array{Union{Rotate,Move}} = []
    for line in eachline(f)
        number = parse(UInt16, line[2:end])
        d = if line[1] == 'N'
            Move(north, number)
        elseif line[1] == 'E'
            Move(east, number)
        elseif line[1] == 'S'
            Move(south, number)
        elseif line[1] == 'W'
            Move(west, number)
        elseif line[1] == 'L'
            Rotate(left, number)
        elseif line[1] == 'R'
            Rotate(right, number)
        else
            @assert line[1] == 'F'
            Move(forward, number)
        end
        push!(re, d)
    end
    re
end

function load()
    open(joinpath(@__DIR__, "12_ship-movement.txt"), "r") do f
        parse_file(f)
    end
end

@enum LeftRight::UInt16 left=0 right=180 full=360

function Base.:!(lr::LeftRight)
    @assert lr != full
    lr == left ? right : left
end

@enum NESW::UInt16 north=0 east=270 south=180 west=90 forward=1

function Base.:+(dir::NESW, degrees::UInt16)::NESW
    dir == forward && return forward
    return NESW((Integer(dir) + degrees) % Integer(full))
end

struct Pos
    lat::Int
    long::Int
    heading::NESW
end
Pos() = Pos(0, 0, east)
Pos(heading::NESW) = Pos(0, 0, heading)

struct Waypoint
    lat::Int
    long::Int
end

struct PosWithWaypoint
    lat::Int
    long::Int
    waypoint::Waypoint
end
PosWithWaypoint(w::Waypoint) = PosWithWaypoint(0, 0, w)

abssum(pos) = abs(pos.lat) + abs(pos.long)

struct Move
    direction::NESW
    by::Int
end

function Base.:+(pos::Pos, nesw::Move)
    lat = pos.lat
    long = pos.long
    direction = nesw.direction == forward ? pos.heading : nesw.direction
    if direction == north
        lat += nesw.by
    elseif direction == south
        lat -= nesw.by
    elseif direction == east
        long += nesw.by
    elseif direction == west
        long -= nesw.by
    end
    Pos(lat, long, pos.heading)
end

function Base.:+(pos::PosWithWaypoint, nesw::Move)
    by = nesw.by
    if nesw.direction == forward
        PosWithWaypoint(pos.lat + by * pos.waypoint.lat, pos.long + by * pos.waypoint.long, pos.waypoint)
    elseif nesw.direction == north
        PosWithWaypoint(pos.lat, pos.long, Waypoint(pos.waypoint.lat + by, pos.waypoint.long))
    elseif nesw.direction == south
        PosWithWaypoint(pos.lat, pos.long, Waypoint(pos.waypoint.lat - by, pos.waypoint.long))
    elseif nesw.direction == east
        PosWithWaypoint(pos.lat, pos.long, Waypoint(pos.waypoint.lat, pos.waypoint.long + by))
    elseif nesw.direction == west
        PosWithWaypoint(pos.lat, pos.long, Waypoint(pos.waypoint.lat, pos.waypoint.long - by))
    end
end

struct Rotate
    direction::LeftRight
    by::UInt16
end

function Base.:+(heading::NESW, turn::Rotate)
    if turn.by % 180 == 0
        heading + turn.by
    else
        heading + (turn.by + Integer(turn.direction))
    end
end

function Base.:+(pos::Pos, rot::Rotate)
    Pos(pos.lat, pos.long, pos.heading + rot)
end

function Base.:+(pos::PosWithWaypoint, rot::Rotate)
    PosWithWaypoint(pos.lat, pos.long, pos.waypoint + rot)
end

function Base.:+(w::Waypoint, rot::Rotate)
    if rot.by == 180
        Waypoint(-w.lat, -w.long)
    elseif rot.by == 270
        w + Rotate(!rot.direction, 90)
    elseif rot.direction == left
        @assert rot.by == 90
        Waypoint(w.long, -w.lat)
    elseif rot.direction == right
        @assert rot.by == 90
        Waypoint(-w.long, w.lat)
    end
end

function ex(route, init)
    foldl(+, route, init=init) |> abssum
end

ex1(route) = ex(route, Pos(east))
ex2(route) = ex(route, PosWithWaypoint(Waypoint(1, 10)))

route = load()

println("Manhatten distance: ", ex1(route))
println("Manhatten distance by waypoint: ", ex2(route))

using Test

@testset "Rain Risk" begin
    @test length(route) == 783

    @testset "operators rotate" begin
        for start in instances(NESW)
            if start != forward
                angle = Integer(full)
                while angle >= 90
                    times = Integer(full) / angle
                    la, lo = rand(Int, 2)
                    pos = Pos(la, lo, start)
                    for _ in 1:times
                        pos = pos + Rotate(left, angle)
                    end
                    @test pos == Pos(la, lo, start)
                    for _ in 1:times
                        pos = pos + Rotate(right, angle)
                    end
                    @test pos == Pos(la, lo, start)
                    angle ÷= 2
                end
            end
        end

        @test Waypoint(4, 10) + Rotate(right, 90) == Waypoint(-10, 4)
    end

    @testset "operators move" begin
        @test Pos(1, 2, north) + Move(forward, 10) == Pos(11, 2, north)
        @test Pos(1, 2, east) + Move(forward, 10) == Pos(1, 12, east)
        @test Pos(1, 2, south) + Move(forward, 10) == Pos(-9, 2,south)
        @test Pos(1, 2, west) + Move(forward, 10) == Pos(1, -8, west)
    end

    input = """
    F10
    N3
    F7
    R90
    F11
    """

    example = parse_file(IOBuffer(input))

    @testset "example 1" begin
        @test ex1(example) == 25
    end

    @testset "example 2" begin
        @test ex2(example) == 286
    end

    @testset "results" begin
        @test ex1(route) == 1152
        @test ex2(route) == 58637
    end
end

max_stars_repo_name: CmdQ/AoC


text:
struct ModelParam{p, T, V <: AbstractVector{T}, P <: Tuple} <: parameter{T}
    v::V
    Θ::P
end
function ModelParam(::Type{Val{p}}, x::V, Θ::P) where {p, T, V <: AbstractVector{T}, P <: Tuple}
    ModelParam{p, T, V, P}(x, Θ)
end

log_jacobian(A::ModelParam) = sum(log_jacobian.(A.Θ))
update!(A::ModelParam) = update!.(A.Θ)
Base.size(A::ModelParam{p}) where p = ( p, )
Base.length(A::ModelParam{p}) where p = p
Base.@pure param_length(x::Tuple) = Val{sum(length.(x))}
@generated cum_inds(x::Tuple{P}) where P = tuple(cumsum([ConstrainedParameters.type_length.(P)...])...)
construct(Θ::Tuple, ::Type{Val{p}} = param_length(Θ)) where p = construct(SlidingVec(zeros(p,1),0,p), Θ, Val{p})
sconstruct(Θ::Tuple, ::Type{Val{p}} = param_length(Θ)) where p = construct(SVector{p}(zeros(p)), Θ, Val{p})
function construct(x::V, Θ::P, ::Type{Val{p}} = param_length(Θ)) where {p, T, V <: AbstractVector{T}, P <: Tuple}
    ModelParam(Val{p}, x, construct.(Θ, (x,), cum_inds(Θ)))
end

function log_density(x::AbstractVector, θ::Tuple, data, ::Type{p} = param_length(θ)) where p
    param = construct(x, θ, p)
    log_jacobian(param) + Main.log_density(data, param.Θ...)
end
function log_density(Θ::ModelParam, data)
    log_jacobian(Θ) + Main.log_density(data, Θ.Θ...)
end
function log_density(x::AbstractVector, Θ::ModelParam, data)
    Θ.v .= x
    update!(Θ)
    log_jacobian(Θ) + Main.log_density(data, Θ.Θ...)
end


text:
using BenchmarkTools, DynamicalBilliards, IterTools

const SUITE = BenchmarkGroup(["DynamicalBilliards"])
bt = billiard_mushroom()
bt2 = billiard_sinai(;setting="periodic")

################################################################################
## COLLISION & PROPAGATION
################################################################################

particles = [Particle(0.05, 0.05, -0.1), MagneticParticle(0.05,0.05,-0.1,1.0)]
obstacles = [bt[1], bt[6], bt2[1], bt2[5]] #distinct obstacles for resolvecollision! tests
proptime = 4.2
ptypes = ["straight", "magnetic"]
colf = (collision,
        next_collision,
        bounce!,
        resolvecollision!,
        propagate!,
        lyapunovspectrum
        )
name = (f) -> split(string(f), '.')[end]

for f in colf
    SUITE[name(f)] = BenchmarkGroup(["propagation", "collision"])
    for ptype ∈ ptypes
        SUITE[name(f)][ptype] = BenchmarkGroup(["propagation", "collision", ptype])
    end
end

for (f, p) in zip(["straight", "magnetic"], particles)
    for o in chain(bt, bt2)
        SUITE["collision"][f][o.name] =
            @benchmarkable collision($p, $o)
    end
end

for (f, p) in zip(["straight", "magnetic"], particles)
    for (bname, bil) in zip(["mushroom", "psinai"], (bt, bt2))
        SUITE["next_collision"][f][bname] =
            @benchmarkable next_collision($p, $bil)
    end
end

for (f, p) in zip(["straight", "magnetic"], particles)
    for (bname, bil) in zip(["mushroom", "psinai"], (bt, bt2))
        ploc = copy(p) #location mutation screws up tests for different bts
        SUITE["bounce!"][f][bname] =
            @benchmarkable bounce!($ploc, $bil)
    end
end

let (f, p) = ("straight", particles[1])
    #resolvecollision! is indepent of particle type
    for o in obstacles
        ploc = copy(p)
        SUITE["resolvecollision!"][f][o.name] =
            @benchmarkable resolvecollision!($ploc, $o)
    end
end

for (f, p) in zip(["straight", "magnetic"], particles)
    ploc = copy(p)
    SUITE["propagate!"][f] =
        @benchmarkable propagate!($ploc, $proptime)
end

for (f, p) in zip(["straight", "magnetic"], particles)
    for (bname, bil) in zip(["mushroom", "psinai"], (bt, bt2))
        SUITE["lyapunovspectrum"][f][bname] =
            @benchmarkable lyapunovspectrum($p, $bil, 10000.0)
    end
end

################################################################################
## randominside – rewrite to add more non-collision benchmarks
################################################################################

SUITE["randominside"] = BenchmarkGroup() #I don't know any tags for randominside...
SUITE["randominside"]["straight"] = BenchmarkGroup(["straight"])

let f = "straight"
    for (bname, bil) in zip(["mushroom", "psinai"], (bt, bt2))
        SUITE["randominside"][f][bname] = @benchmarkable randominside($bil)
    end
end


text:
import GEAS:qmsim_base, sim_pt_QTL, breeding_program, bp_summary
using Distributions

"""
    function my_sim_pt_qtl()
---
1. Simulate a base population with `QMSim`.
2. go through a breeding program
3. summarize
"""
function my_sim_pt_qtl()
    dicpar = Dict(
        :hg_size => "1200 [0]\n1200 [2000]",
        :male => 600,
        :female => 600,         # 600 + 600 => 1200
        :nchr => 29,
        :lchr => 100,
        :nsnp => 1800,          # to simulate  1725*29 => 50025 SNP
    )
    bpar = (; dicpar...)
    base = qmsim_base(bpar, dir="dat/test/QMSim_Linux")
    #d = MvNormal(zeros(2), [1 -.1; -.1 1])
    #qtl = sim_pt_QTL(base, 1000, d)
    #
    #nqtl = [length(qtl[1].pos), length(qtl[2].pos)]
    #par = begin
    #    Parameters = Dict(:nSire=> 100,
    #                      :nDam => 200, # male:female = 1:2
    #                      :nSib => 60,
    #                      :nC7e => 30, # number of sib to be challenged
    #                      :nG8n => 10, # number of generations
    #                      :nQTL => nqtl, # usually [1000, 1000]
    #                      :h²   => [.5, .5],
    #                      :p8e  => .5, # percentage of affected in the binary trait
    #                      :e19e => .8, # edit_successsful_rate
    #                      :w4t  => 1 # weight on the binary trait EBV
    #                      )
    #    (; Parameters...)       # named tuple.  contents as above
    #end
    #prd, snp = breeding_program(base, par, qtl)
    #bp_summary(prd, snp, qtl, dir="dat/test/QMSim_Linux")
end


text:
@testset "Gradient Plan" begin
    io = IOBuffer()
    M = ManifoldsBase.DefaultManifold(2)
    x = [4.,2.]
    o = GradientDescentOptions(x, StopAfterIteration(20), ConstantStepsize(1.))
    o.∇ = [1., 0.]
    f = y -> distance(M,y,x).^2
    ∇f = y -> -2*log(M,y,x)
    p = GradientProblem(M,f,∇f)
    @test get_initial_stepsize(p,o) == 1.
    @test get_stepsize(p,o,1) == 1.
    @test get_last_stepsize(p,o,1) == 1.
    # Check Fallbacks of Problen
    @test get_cost(p,o.x) == 0.
    @test get_gradient(p,o.x) == zero_tangent_vector(M,x)
    @test_throws ErrorException getProximalMap(p,1.,o.x,1)
    @test_throws ErrorException get_subgradient(p,o.x)
    # Additional Specific Debugs
    a1 = DebugGradient(false, x -> print(io,x))
    a1(p,o,1)
    @test String(take!(io)) == "∇F(x):[1.0, 0.0]"
    a1a = DebugGradient("s:", x -> print(io,x))
    a1a(p,o,1)
    @test String(take!(io)) == "s:[1.0, 0.0]"
    a2 = DebugGradientNorm(false, x -> print(io,x))
    a2(p,o,1)
    @test String(take!(io)) == "|∇F(x)|:1.0"
    a2a = DebugGradientNorm("s:", x -> print(io,x))
    a2a(p,o,1)
    @test String(take!(io)) == "s:1.0"
    a3 = DebugStepsize(false, x -> print(io,x))
    a3(p,o,1)
    @test String(take!(io)) == "s:1.0"
    a3a = DebugStepsize("S:", x -> print(io,x))
    a3a(p,o,1)
    @test String(take!(io)) == "S:1.0"
    # Additional Specific Records
    b1 = RecordGradient(o.∇)
    b1(p,o,1)
    @test b1.recordedValues == [o.∇]
    b2 = RecordGradientNorm()
    b2(p,o,1)
    @test b2.recordedValues == [1.]
    b3 = RecordStepsize()
    b3(p,o,1)
    b3(p,o,2)
    b3(p,o,3)
    @test b3.recordedValues == [1., 1., 1.]
end


text:
"""
    function justAfunction(a::Int)

Just my personalized `function`.
"""
function justAfunction(a::Int)
    return a
end

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
# Example: Compartmental models in epidemiology
#=
- https://github.com/epirecipes/sir-julia/blob/master/markdown/function_map/function_map.md
- https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#Deterministic_versus_stochastic_epidemic_models
=#
using ModelingToolkit, Test

@inline function rate_to_proportion(r,t)
    1-exp(-r*t)
end;

# Independent and dependent variables and parameters
@parameters t c nsteps δt β γ
D = Difference(t; dt=0.1)
@variables S(t) I(t) R(t) next_S(t) next_I(t) next_R(t)

infection = rate_to_proportion(β*c*I/(S+I+R),δt)*S
recovery = rate_to_proportion(γ,δt)*I

# Equations
eqs = [D(S) ~ S-infection,
       D(I) ~ I+infection-recovery,
       D(R) ~ R+recovery]

# System
sys = DiscreteSystem(eqs,t,[S,I,R],[c,nsteps,δt,β,γ]; controls = [β, γ])

# Problem
u0 = [S => 990.0, I => 10.0, R => 0.0]
p = [β => 0.05, c => 10.0, γ => 0.25, δt => 0.1, nsteps => 400]
tspan = (0.0,ModelingToolkit.value(substitute(nsteps,p))) # value function (from Symbolics) is used to convert a Num to Float64
prob_map = DiscreteProblem(sys,u0,tspan,p)

# Solution
using OrdinaryDiffEq
sol_map = solve(prob_map,FunctionMap());

# Direct Implementation

function sir_map!(u_diff,u,p,t)
    (S,I,R) = u
    (β,c,γ,δt) = p
    N = S+I+R
    infection = rate_to_proportion(β*c*I/N,δt)*S
    recovery = rate_to_proportion(γ,δt)*I
    @inbounds begin
        u_diff[1] = S-infection
        u_diff[2] = I+infection-recovery
        u_diff[3] = R+recovery
    end
    nothing
end;
u0 = [990.0,10.0,0.0];
p = [0.05,10.0,0.25,0.1];
prob_map = DiscreteProblem(sir_map!,u0,tspan,p);
sol_map2 = solve(prob_map,FunctionMap());

@test Array(sol_map) ≈ Array(sol_map2)


text:
module TestBlockLinkedLists

using EventTrackerBase.BlockLinkedLists: BlockLinkedList, alloclast!
using Test

function test_push_collect(xs)
    lst = foldl(push!, xs; init = BlockLinkedList(Vector{eltype(xs)}, 4))
    @test collect(lst) == collect(xs)
end

function test_push_collect()
    @testset "$label" for (label, xs) in [
        "1:3" => 1:3,
        "1:8" => 1:8,
        "1:9" => 1:9,
        "1:10" => 1:10,
        "a-to-z" => sprint(io -> foreach(x -> print(io, x), 'a':'z')),
    ]
        test_push_collect(xs)
    end
end

end  # module


text:
@info "Testing mesh_func..."
@testset "Transforms mesh_func" begin

    @testset "sample_points" begin
        m = load_trimesh(
            joinpath(@__DIR__, "../assets/sphere.obj"),
            joinpath(@__DIR__, "../assets/sphere.obj"),
        )

        samples = sample_points(m, 1000)
        _radius = sqrt.(sum(samples .^ 2; dims = 1))
        @test all(isapprox.(_radius, 1.0, rtol = 1e-2, atol = 1e-5))
        @test gradient(x -> sum(sample_points(x, 1000)), m) isa Tuple
    end

    _mesh = load_trimesh(
        joinpath(@__DIR__, "../assets/teapot.obj"),
        joinpath(@__DIR__, "../assets/sphere.obj"),
    )

    for (inplace, FUNC) in [(true, Flux3D.normalize!), (false, Flux3D.normalize)]
        @testset "$(FUNC)" begin
            m = deepcopy(_mesh)
            m2 = FUNC(m)
            if inplace
                @test m2 === m
            else
                @test m2 !== m
            end
            @test all(
                isapprox.(
                    mean(get_verts_list(m2)[1]; dims = 2),
                    0.0,
                    rtol = 1e-5,
                    atol = 1e-5,
                ),
            )
            @test all(
                isapprox.(
                    std(get_verts_list(m2)[1]; dims = 2),
                    1.0,
                    rtol = 1e-4,
                    atol = 1e-5,
                ),
            )
            @test all(
                isapprox.(
                    mean(get_verts_list(m2)[2]; dims = 2),
                    0.0,
                    rtol = 1e-5,
                    atol = 1e-5,
                ),
            )
            @test all(
                isapprox.(
                    std(get_verts_list(m2)[1]; dims = 2),
                    1.0,
                    rtol = 1e-4,
                    atol = 1e-5,
                ),
            )
        end
    end

    for (inplace, FUNC) in [(true, Flux3D.scale!), (false, Flux3D.scale)]
        @testset "$(FUNC)" begin
            m = deepcopy(_mesh)
            m2 = FUNC(FUNC(m, 2.0), 0.5)
            if inplace
                @test m2 === m
            else
                @test m2 !== m
            end
            @test all(
                isapprox.(
                    get_verts_packed(_mesh),
                    get_verts_packed(m2),
                    rtol = 1e-5,
                    atol = 1e-5,
                ),
            )
        end
    end

    for (inplace, FUNC) in [(true, Flux3D.rotate!), (false, Flux3D.rotate)]
        @testset "$(FUNC)" begin
            m = deepcopy(_mesh)
            rotmat = 2 .* one(rand(Float32, 3, 3))
            rotmat_inv = inv(rotmat)
            rotmat_b = cat(rotmat, rotmat, dims = 3)
            rotmat_inv_b = cat(rotmat_inv, rotmat_inv, dims = 3)
            m2 = FUNC(FUNC(m, rotmat), rotmat_inv)
            m3 = FUNC(FUNC(m, rotmat_b), rotmat_inv_b)
            if inplace
                @test m2 === m
                @test m3 === m
            else
                @test m2 !== m
                @test m3 !== m
            end
            @test all(
                isapprox.(
                    get_verts_packed(_mesh),
                    get_verts_packed(m2),
                    rtol = 1e-5,
                    atol = 1e-5,
                ),
            )
            @test all(
                isapprox.(
                    get_verts_packed(_mesh),
                    get_verts_packed(m3),
                    rtol = 1e-5,
                    atol = 1e-5,
                ),
            )
        end
    end

    for (inplace, FUNC) in [(true, Flux3D.realign!), (false, Flux3D.realign)]
        @testset "$(FUNC)" begin
            m = deepcopy(_mesh)
            src = Flux3D.scale(m, 2.0)
            tgt = Flux3D.scale(m, 1.0)
            src_1 = FUNC(src, tgt, 1)
            if inplace
                @test src === src_1
            else
                @test src !== src_1
            end
            # Transformed PointCloud should be inside the bounding box defined by `tgt` PointCloud
            @test all(
                maximum(get_verts_list(_mesh)[1]; dims = 2) .>=
                get_verts_list(src_1)[1] .>=
                minimum(get_verts_list(_mesh)[1]; dims = 2),
            )
            src_1 = FUNC(src, tgt, 2)
            @test all(
                maximum(get_verts_list(_mesh)[2]; dims = 2) .>=
                get_verts_list(src_1)[1] .>=
                minimum(get_verts_list(_mesh)[2]; dims = 2),
            )
        end
    end

    for (inplace, FUNC) in [(true, Flux3D.translate!), (false, Flux3D.translate)]
        @testset "$(FUNC)" begin
            m = deepcopy(_mesh)
            m2 = FUNC(FUNC(m, 1.0), -1.0)
            if inplace
                @test m2 === m
            else
                @test m2 !== m
            end
            @test all(
                isapprox.(
                    get_verts_packed(_mesh),
                    get_verts_packed(m2),
                    rtol = 1e-5,
                    atol = 1e-5,
                ),
            )
        end
    end

    for (inplace, FUNC) in [(true, Flux3D.offset!), (false, Flux3D.offset)]
        @testset "$(FUNC)" begin
            m = deepcopy(_mesh)
            _offset = ones(size(get_verts_packed(m)))
            m2 = FUNC(FUNC(m, _offset), (_offset .* -1))
            if inplace
                @test m2 === m
            else
                @test m2 !== m
            end
            @test all(
                isapprox.(
                    get_verts_packed(_mesh),
                    get_verts_packed(m2),
                    rtol = 1e-5,
                    atol = 1e-5,
                ),
            )
        end
    end
end


text:
using Avatar

generate()
generate(20, 10, 40, 80, 848743, "avatar2.png")
generate(4, 4, 100, 100, 932683, "avatar3.png")


text:
using LanguageDynamics
using Test

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


text:
module QuDirac

    using Compat
    import Base.Iterators.product
    import LinearAlgebra

    if !(v"0.3-" <= VERSION < v"0.4-")
        @warn("QuDirac v0.1 only officially supports the v0.3 release of Julia. Your version of Julia is $VERSION.")
    end

    ####################
    # String Constants #
    ####################
    const lang = "\u27E8"
    const rang = "\u27E9"
    const otimes = "\u2297"
    const vdots ="\u205E"

    ##################
    # Abstract Types #
    ##################
    abstract type AbstractInner end

    struct UndefinedInner <: AbstractInner end
    struct KroneckerDelta <: AbstractInner end

    abstract type AbstractDirac{P<:AbstractInner,N} end
    abstract type DiracOp{P,N} <: AbstractDirac{P,N} end
    abstract type DiracState{P,N} <: AbstractDirac{P,N} end

    abstract type DiracScalar <: Number end

    #############
    # Functions #
    #############
    # These functions will be in
    # QuBase when it releases
    tensor() = error("Cannot call tensor function without arguments")
    tensor(s...) = reduce(tensor, s)

    ######################
    # Include Statements #
    ######################
    include("labels.jl")
    include("innerexpr.jl")
    include("state.jl")
    include("opsum.jl")
    include("outerproduct.jl")
    include("printfuncs.jl")
    include("dictfuncs.jl")
    include("mapfuncs.jl")
    include("str_macros.jl")

    #################
    # default_inner #
    #################
    # DEFAULT_INNER = KroneckerDelta

    DEFAULT_INNER = KroneckerDelta()
    default_inner() = DEFAULT_INNER
    function default_inner!(ptype::T) where T<:AbstractInner
        global DEFAULT_INNER = ptype
    end

    default_inner!(KroneckerDelta())

    OpSum(dict::Dict) = OpSum(default_inner(), dict)
    Ket(dict::Dict) = Ket(default_inner(), dict)
    ket(label::StateLabel) = Ket(default_inner(), label)
    ket(items...) = Ket(default_inner(), StateLabel(items))

    export AbstractInner,
        UndefinedInner,
        KroneckerDelta,
        default_inner,
        default_inner!,
        AbstractDirac,
        DiracState,
        DiracOp,
        # All functions that conflict
        # with QuBase should be exported
        # below:
        tensor,
        commute,
        anticommute,
        normalize,
        normalize!

end # module QuDirac


text:
function poly_horner(x, a...)
b=zero(x)
for i = length(a):-1:1
b = a[i] + b * x
end
return b
end

max_stars_repo_path: Module 2/Chapter04/poly_naive.jl

max_issues_repo_path: Module 2/Chapter04/poly_naive.jl

max_forks_repo_path: Module 2/Chapter04/poly_naive.jl


text:
using KMeanClusttering;
using Test;
@testset "KMeanClusttering.jl" begin

    for k in 3:10
        dt = KMeanClusttering.generate_test_2d_data(k);
        km = KMeanClusttering.KModel(k,dt);
        printstyled("Runing test k = $k\n", color = :green);
        println("Before K-Mean")
        bd = KMeanClusttering.ave_distance(km);
        println("Ave distance $(bd)");
        println("After K-Mean")
        km = KMeanClusttering.find_clusters(k,km.data);
        km = KMeanClusttering.get_clean(km);
        ad = KMeanClusttering.ave_distance(km);
        println("Ave distance $(ad)");
        println();
        @test bd > ad;
    end

    printstyled("""For a graphic view run the jupyter notebook "Example.ipynb",
    located in the root folder of the package.""", color = :green);
    println();
    println();

end


text:
# Import the progress bar and Dates for the conversion of seconds
import Dates

add_verbose_scope(:ClassGroupProof)

# This is a modified showprogress from Pkg.GitTools

PROGRESS_BAR_PERCENTAGE_GRANULARITY = Ref(0.1)

Base.@kwdef mutable struct MiniProgressBar
    max::Float64 = 1.0
    header::String = ""
    color::Symbol = :white
    width::Int = 40
    current::Float64 = 0.0
    prev::Float64 = 0.0
    has_shown::Bool = false
    time_shown::Float64 = 0.0
end

function showprogress(io::IO, p::MiniProgressBar, info)
  perc = p.current / p.max * 100
  prev_perc = p.prev / p.max * 100
  # Bail early if we are not updating the progress bar,
  # Saves printing to the terminal
  if p.has_shown && !((perc - prev_perc) > PROGRESS_BAR_PERCENTAGE_GRANULARITY[])
    return
  end
  if !isinteractive()
    t = time()
    if p.has_shown && (t - p.time_shown) < NONINTERACTIVE_TIME_GRANULARITY[]
      return
    end
    p.time_shown = t
  end
  p.time_shown = time()
  p.prev = p.current
  p.has_shown = true
  n_filled = ceil(Int, p.width * perc / 100)
  n_left = p.width - n_filled
  print(io, "    ")
  printstyled(io, p.header, color=p.color, bold=true)
  print(io, " [")
  print(io, "="^n_filled, ">")
  print(io, " "^n_left, "]  ", )
  @printf io "%2.1f %%" perc
  print(io, info)
  print(io, "\r")
end

function class_group_proof(clg::ClassGrpCtx, lb::fmpz, ub::fmpz; extra :: fmpz=fmpz(0), prec::Int = 100, do_it=1:ub)
  PB = MiniProgressBar(header = "Class group proof")

  #for all prime ideals P with lb <= norm <= ub, find a relation
  #tying that prime to the factor base
  # if extra is useful, assume that the function was already run for all primes
  # up to norm extra

  if extra == 0
    extra = norm(clg.FB.ideals[1])
  end
  lb = max(lb, norm(clg.FB.ideals[1]))
  #println("expect to need ", Int(floor(li(ub*1.0) - li(lb*1.0))), " primes")
  O = order(clg.FB.ideals[1])
  n = degree(O)
  p = next_prime(root(lb, n))
  np = Int(floor(log(abs(discriminant(O)))/log(2)/2))
  no_primes = 0
  no_ideals = 0
  if do_it.start > 1
    p = fmpz(next_prime(do_it.start))
  end
  r = fmpz()
  _no_of_primes = logarithmic_integral(Float64(ub)) - logarithmic_integral(Float64(lb))
  #gc_enable(false)
  rate = 0.0
  length_eta = 0
  if _no_of_primes < 1000
    interval = 10
  else
    interval = 1000
  end

  while p < do_it.stop
    no_primes += 1

    @v_do :ClassGroupProof if no_primes % interval == 0
      #println("did $no_primes prime numbers so far, now $p, need to reach $ub (~$(no_primes/_no_of_primes))")
      last_time = PB.time_shown
      cur_time = time()
      prev = PB.prev
      PB.current = no_primes/_no_of_primes
      # from PB.current to prev it took cur_time - last_time seconds

      if rate == 0.0
        rate = ((PB.current - PB.prev)/(cur_time - last_time))
      else
        rate = (((PB.current - PB.prev)/(cur_time - last_time)) + rate)/2
      end

      duration = (1 - PB.current)/rate

      duration = round(Int, duration)
      ETA = " (ETA: $(Dates.Time(Dates.Nanosecond(duration * 10^9))))"
      if length(ETA) < length_eta
        ETA = ETA * " "^(length_eta - length(ETA))
      else
        length_eta = length(ETA)
      end

      showprogress(stdout, PB, ETA)
      flush(stdout)
    end
    deg_lim = Int(flog(ub, p))
    if deg_lim == 1
      low_lim = 1
    else
      low_lim = Int(flog(lb, p))
    end
    fac = prime_decomposition(O, Int(p), deg_lim, low_lim)
    for _k in fac
      k = _k[1]
      if norm(k) <= lb
        continue
      end
      no_ideals += 1
      #if no_ideals % 10 == 0
      #  println("done $no_ideals ideals so far...")
      #  x = Base.gc_num()
      #  println("alloc $(x.malloc)   free $(x.freed)  diff: $(x.malloc - x.freed)")
      #  gc_enable(true)
      #  gc()
      #  gc_enable(false)
      #  x = Base.gc_num()
      #  println("alloc $(x.malloc)   free $(x.freed)  diff: $(x.malloc - x.freed)")
      #end
      #println("to be more precise: $k")
      E = class_group_small_real_elements_relation_start(clg, k, limit=10, prec=prec)
      while true
        sucess = false
        a = class_group_small_real_elements_relation_next(E)
        n = norm_div(a, norm(k), np)
        if gcd(numerator(n), p) > extra
#          println("a: $a, $(norm(a)), $(norm(k)), $n")
#          println("contains too many conjugates, bad")
          continue
        end
        f, r = issmooth!(clg.FB.fb_int, numerator(n))
        if f
          M = SMat{Int}()
          fl = _factor!(clg.FB, a, false, n)[1]
          if fl
            break
          else
#            println("not smooth, ideal")
          end
        else
#          println("not smooth, int")
        end
      end
    end
    p = next_prime(p)
  end
  @v_do :ClassGroupProof begin
    PB.current = 1.0
    showprogress(stdout, PB, " (ETA: 00:00:00)" * " "^(max(length_eta - 15, 0)))
    println()
  end
  #println("success: used $no_primes numbers and $no_ideals ideals")
  #gc_enable(true)
end


text:
# # Taylor Polynomial Models

#src We generate a documentation page from this source file.
#src Let's tell the user:
#md # This file is automatically generated from source code.
#md # For usage examples refer to [Summary & Quick Examples](@ref taylor_summary).

# We provide vector valued polynomial Taylor models of degree 1 or 2.
# They implement the `SurrogateModel` interface.
#
# We allow the user to either provide gradient and hessian callback handles
# or to request finite difference approximations.
# For using callbacks, we have `TaylorConfigCallbacks`. \
# There are two ways to use finite differences. The old (not recommended way) is to
# use `TaylorConfigFiniteDiff`. This uses `FiniteDiff.jl` and could potentially
# require more evaluations. \
# To make use of the new 2-phase construction procedure, use `TaylorConfig` and
# set the fields `gradients` and `hessians` to an `RFD.FiniteDiffStamp`.
# If they use the same stamp (default: `RFD.CFDStamp(1,3) :: CFDStamp{3,Float64}`),
# it should be the most efficient, because we get the gradients for free from computing the hessians.

include("RecursiveFiniteDifferences.jl")

using .RecursiveFiniteDifferences
const RFD = RecursiveFiniteDifferences

# The actual model is defined only by the gradient vectors at `x₀` and the Hessians (if applicable).
@with_kw struct TaylorModel{
    XT <: AbstractVector{<:Real}, FXT <: AbstractVector{<:Real},
    G <: AbstractVector{<:AbstractVector{<:Real}},
    HT <: Union{Nothing,AbstractVector{<:AbstractMatrix{<:Real}}},
    } <: SurrogateModel

    ## expansion point and value
    x0 :: XT
    fx0 :: FXT

    ## gradient(s) at x0
    g :: G
    H :: HT = nothing
end

fully_linear( :: TaylorModel ) = true

# Note, that the derivative approximations are actually constructed for the function(s)
# ```math
#     f_ℓ ∘ s^{-1}
# ```
# if some internal transformation ``s`` has happened before.
# If the problem is unbounded then ``s = \operatorname{id} = s^{-1}``.

# ## Model Construction

# Because of all the possibilities offered to the user, we actually have several
# (sub-)implementiations of `SurrogateConfig` for Taylor Models.
abstract type TaylorCFG <: SurrogateConfig end

# We make sure, that all subtypes have a field `max_evals`:
max_evals( cfg :: TaylorCFG ) = cfg.max_evals

# ### Recursive Finite Difference Models

# Let's start by defining the recommended way of using Taylor approximations.
# The derivative information is approximated using a dynamic programming approach
# and we take care to avoid unnecessary objective evaluations.

@doc """
    TaylorConfig(; degree, gradients :: RFD.CFDStamp, hessians :: RFD.CFDStamp, max_evals)

Configuration for a polynomial Taylor model using finite difference approximations of the derivatives.
By default we have `degree = 2` and `gradients == hessians == RFD.CFDStamp(1,2)`, that is,
a first order central difference scheme of accuracy order 3 is recursed to compute the Hessians
and the gradients.
In this case, the finite difference scheme is the same for both Hessians and gradients and we profit
from caching intermediate results.
"""
@with_kw struct TaylorConfig{
        S1 <: RFD.FiniteDiffStamp,
        S2 <: Union{Nothing,RFD.FiniteDiffStamp}
    } <: TaylorCFG

    degree :: Int64 = 2

    gradients :: S1 = RFD.CFDStamp(1,2)
    hessians :: S2 = gradients

    max_evals :: Int64 = typemax(Int64)

    @assert 1 <= degree <= 2 "Can only construct linear and quadratic polynomial Taylor models."
end

combinable( :: TaylorConfig ) = true

# The new meta type only stores database indices of sites used for a finite diff approximation
# in the actual construction call and is filled in the `prepare_XXX` methods:
@with_kw struct TaylorIndexMeta{W1, W2} <: SurrogateMeta
    database_indices :: Vector{Int} = Int[]
    grad_setter_indices :: Vector{Int} = Int[]
    hess_setter_indices :: Vector{Int} = Int[]
    hess_wrapper :: W1 = nothing
    grad_wrapper :: W2 = nothing
end

# The end user won't be interested in the wrappers, so we put `nothing` in there:
get_saveable_type( :: TaylorIndexMeta ) = TaylorIndexMeta{Nothing, Nothing}
get_saveable( meta :: TaylorIndexMeta ) = TaylorIndexMeta(;
    grad_setter_indices = meta.grad_setter_indices,
    hess_setter_indices = meta.hess_setter_indices
)

# The new construction process it is a bit complicated.
# We set up a recursive finite diff tree and
# need this little helper:
"Return `unique_elems, indices = unique_with_indices(arr)` such that
`unique_elems[indices] == arr` (and `unique_elems == unique(arr)`)."
function unique_with_indices( x :: AbstractVector{T} ) where T
	unique_elems = T[]
	indices = Int[]
	for elem in x
		i = findfirst( e -> all( isequal.(e,elem) ), unique_elems )
		if isnothing(i)
			push!(unique_elems, elem)
			push!(indices, length(unique_elems) )
		else
			push!(indices, i)
		end
	end
	return unique_elems, indices
end

# Now, if the polynomial degree equals 2 we construct a tree for the Hessian calculation.
# In any case, we need a tree for the gradients/jacobian.
# If the `RFD.FiniteDiffStamp` for the gradients is the same as for the Hessians, we can re-use the
# Hessian tree for this purpose. Else, we need to construct a new one.

function _get_RFD_trees( x, fx, grad_stamp, hess_stamp = nothing, deg = 2)
    if deg >= 2
        @assert !isnothing(hess_stamp)
        ## construct tree for hessian first
        hess_wrapper = RFD.DiffWrapper(; x0 = x, fx0 = fx, stamp = hess_stamp, order = 2 )
    else
        hess_wrapper = nothing
    end

    if !isnothing(hess_wrapper) && grad_stamp == hess_stamp
        grad_wrapper = hess_wrapper
    else
        grad_wrapper = RFD.DiffWrapper(; x0 = x, fx0 = fx, stamp = grad_stamp, order = 1 )
    end

    return grad_wrapper, hess_wrapper
end

function prepare_init_model(cfg :: TaylorConfig, objf :: AbstractObjective,
    mop :: AbstractMOP, iter_data ::AbstractIterData, db :: AbstractDB, algo_cfg :: AbstractConfig; kwargs...)

    return prepare_update_model( nothing, objf, TaylorIndexMeta(), mop, iter_data, db, algo_cfg; kwargs... )
end

# The actual database preparations are delegated to the `prepare_update_model` function.
function prepare_update_model( mod :: Union{Nothing, TaylorModel}, objf, meta :: TaylorIndexMeta, mop,
    iter_data, db, algo_cfg; kwargs... )

    x = get_x( iter_data )
    fx = get_fx( iter_data )
    x_index = get_x_index( iter_data )

    cfg = model_cfg( objf )

    grad_wrapper, hess_wrapper = _get_RFD_trees( x, fx, cfg.gradients, cfg.hessians, cfg.degree )

    XT = typeof(x)
    FXT = typeof(fx)

    lb, ub = full_bounds_internal( mop )

    if cfg.degree >= 2
        RFD.substitute_leaves!(hess_wrapper)
        ## We project into the scaled variable boundaries to avoid violations:
        hess_sites = [ _project_into_box(s,lb,ub) for s in RFD.collect_leave_sites( hess_wrapper ) ]
    else
        hess_sites = XT[]
    end

    ## collect leave sites for gradients
    if grad_wrapper == hess_wrapper
        grad_sites = hess_sites
    else
        RFD.substitute_leaves!( grad_wrapper )
        grad_sites = [ _project_into_box(s, lb,ub) for s in RFD.collect_leave_sites( grad_wrapper ) ]
    end

    combined_sites = [ [x,]; hess_sites; grad_sites ]

    unique_new, unique_indices = unique_with_indices(combined_sites)
    ## now: `combined_sites == unique_new[unique_indices]`

    num_hess_sites = length(hess_sites)
    hess_setter_indices = unique_indices[ 2 : num_hess_sites + 1]
    grad_setter_indices = unique_indices[ num_hess_sites + 2 : end ]
    ## now: `hess_sites == unique_new[ hess_setter_indices ]` and
    ## `grad_sites == unique_new[ grad_setter_indices ]`

    db_indices = [ [x_index,]; [ new_result!(db, ξ, FXT()) for ξ in unique_new[ 2:end ] ] ]
    ## now: `unique_new == get_site.(db, db_indices)`

    ## we return a new meta object in each iteration, so that the node cache is reset in between.
    return TaylorIndexMeta(;
        database_indices = db_indices,
        grad_setter_indices,
        hess_setter_indices,
        grad_wrapper,
        hess_wrapper
    )
end

# If the meta data is set correctly, we only have to set the value vectors for the
# RFD trees and then ask for the right matrices:
function _init_model( cfg :: TaylorConfig, objf :: AbstractObjective, mop :: AbstractMOP,
    iter_data :: AbstractIterData, db :: AbstractDB, algo_config :: AbstractConfig, meta :: TaylorIndexMeta; kwargs... )
    return update_model( nothing, objf, meta, mop, iter_data, db, algo_config; kwargs...)
end

# Note, that we only perform updates if the iterate has changed, `x != mod.x0`, because
# we don't change the differencing parameters.
function update_model( mod :: Union{Nothing, TaylorModel}, objf :: AbstractObjective, meta :: TaylorIndexMeta,
    mop :: AbstractMOP, iter_data :: AbstractIterData, db :: AbstractDB, algo_config :: AbstractConfig; kwargs...)

    x = get_x(iter_data)
    if isnothing(mod) || (x != mod.x0)
        all_leave_vals = get_value.( db, meta.database_indices )

        if !isnothing( meta.hess_wrapper )
            hess_leave_vals = all_leave_vals[ meta.hess_setter_indices ]
            RFD.set_leave_values!( meta.hess_wrapper, hess_leave_vals )
            H = [ RFD.hessian( meta.hess_wrapper; output_index = ℓ ) for ℓ = 1 : num_outputs(objf) ]
        else
            H = nothing
        end

        ## calculate gradients
        if meta.hess_wrapper != meta.grad_wrapper
            grad_leave_vals = all_leave_vals[ meta.grad_setter_indices ]
            RFD.set_leave_values!( meta.grad_wrapper, grad_leave_vals )
        end

        ## if hessians have been calculated before and `grad_wrapper == hess_wrapper` we profit from caching
        J = RFD.jacobian( meta.grad_wrapper )
        g = copy.( eachrow( J ) )

        return TaylorModel(;
            x0 = x,
            fx0 = get_fx( iter_data ),
            g, H
        ), meta
    else
        return mod,meta
    end
end

# ### Callback Models with Derivatives, AD or Adaptive Finite Differencing

# The old way of defining Taylor Models was to provide an objective callback function
# and either give callbacks for the derivatives too or ask for automatic differencing.
# This is very similar to the `ExactModel`s, with the notable difference that the
# gradient and Hessian information is only used to construct models
# ``m_ℓ = f_0 + \mathbf g^T \mathbf h + \mathbf h^T \mathbf H \mathbf h`` **once** per iteration
# and then use these ``m_ℓ`` for all subsequent model evaluations/differentiation.

"""
    TaylorCallbackConfig(;degree=1,gradients,hessians=nothing,max_evals=typemax(Int64))

Configuration for a linear or quadratic Taylor model where there are callbacks provided for the
gradients and -- if applicable -- the Hessians.
The `gradients` keyword point to an array of callbacks where each callback evaluates
the gradient of one of the outputs.
"""
@with_kw struct TaylorCallbackConfig{
        G <:Union{Nothing,AbstractVector{<:Function}},
        J <:Union{Nothing,Function},
        H <:Union{Nothing,AbstractVector{<:Function}},
    } <: TaylorCFG

    degree :: Int64 = 1
    gradients :: G
    jacobian :: J = nothing
    hessians :: H = nothing

    max_evals :: Int64 = typemax(Int64)

    @assert 1 <= degree <= 2 "Can only construct linear and quadratic polynomial Taylor models."
    @assert !(isnothing(gradients) && isnothing(jacobian)) "Provide either `gradients` or `jacobian`."
    @assert isa( gradients, AbstractVector ) && !isempty( gradients ) || !isnothing(jacobian) "Provide either `gradients` or `jacobian`."
    @assert !(isnothing(gradients)) && ( isnothing(hessians) || length(gradients) == length(hessians) ) "Provide same number of gradients and hessians."
end

"""
    TaylorApproximateConfig(;degree=1,mode=:fdm,max_evals=typemax(Int64))

Configure a linear or quadratic Taylor model where the gradients and Hessians are constructed
either by finite differencing (`mode = :fdm`) or automatic differencing (`mode = :autodiff`).
"""
@with_kw struct TaylorApproximateConfig <: TaylorCFG
    degree :: Int64 = 1

    mode :: Symbol = :fdm

    max_evals :: Int64 = typemax(Int64)

    @assert 1 <= degree <= 2 "Can only construct linear and quadratic polynomial Taylor models."
    @assert mode in [:fdm, :autodiff] "Use `mode = :fdm` or `mode = :autodiff`."
end

# For these models, it is not advisable to combine objectives:
combinable( :: Union{TaylorCallbackConfig, TaylorApproximateConfig}) = false

# In both cases we transfer the finalized callbacks to the same meta structs.
# In fact, `GW` and `HW` are `DiffFn`s as defined in `src/diff_wrappers.jl` (or `nothing`):

struct TaylorMetaCallbacks{GW, HW} <: SurrogateMeta
    gw :: GW
    hw :: HW
end

# Again, we have a `tfn::TransformerFn` that represents the (un)scaling. \
# If callbacks are provided, then we use the `GradWrapper` and the `HessWrapper`.

function init_meta( cfg :: TaylorCallbackConfig, objf, tfn )
    gw = GradWrapper( tfn, cfg.gradients, cfg.jacobian )
    hw = cfg.degree == 2 ? isa( cfg.hessians, AbstractVector{<:Function} ) ?
        HessWrapper(tfn, cfg.hessians ) : nothing : nothing;
    return TaylorMetaCallbacks( gw, hw )
end

# If no callbacks are provided, we inspect the `mode` and use the corresponding wrappers:

function init_meta( cfg :: TaylorApproximateConfig, objf, tfn )
    if cfg.mode == :fdm
        gw = FiniteDiffWrapper(objf, tfn, nothing)
    else
        gw = AutoDiffWrapper( objf, tfn, nothing )
    end
    hw = cfg.degree == 2 ? HessFromGrads(gw) : nothing
    return TaylorMetaCallbacks( gw, hw )
end

# The initialization for the legacy config types is straightforward as they don't use
# the new 2-phase process:
function prepare_init_model(cfg :: Union{TaylorCallbackConfig, TaylorApproximateConfig}, objf :: AbstractObjective,
    mop :: AbstractMOP, ::AbstractIterData, ::AbstractDB, :: AbstractConfig; kwargs...)
    tfn = TransformerFn(mop)
    return init_meta( cfg, objf, tfn )
end

# The model construction happens in the `update_model` method and makes use of the `get_gradient` and `get_hessian`
# methods for the wrappers stored in `meta`:

function _init_model(cfg :: Union{TaylorCallbackConfig, TaylorApproximateConfig}, objf :: AbstractObjective,
    mop :: AbstractMOP, iter_data ::AbstractIterData, db ::AbstractDB, algo_config :: AbstractConfig,
    meta :: TaylorMetaCallbacks; kwargs...)
    return update_model(nothing, objf, meta, mop, iter_data, db, algo_config; kwargs...)
end

function update_model( mod :: Union{Nothing,TaylorModel}, objf :: AbstractObjective, meta :: TaylorMetaCallbacks,
    mop :: AbstractMOP, iter_data :: AbstractIterData, db :: AbstractDB, algo_config :: AbstractConfig; kwargs...)

    x = get_x(iter_data)
    if isnothing(mod) || (x != mod.x0)
        fx = get_fx( iter_data )

        num_out = num_outputs( objf )
        g = [ get_gradient(meta.gw , x , ℓ ) for ℓ = 1 : num_out ]

        if !isnothing(meta.hw)
            H = [ get_hessian(meta.hw, x, ℓ) for ℓ = 1 : num_out ]
        else
            H = nothing
        end

        return TaylorModel(; x0 = x, fx0 = fx, g, H ), meta
    else
        return mod, meta
    end

end

# ## Model Evaluation

# The evaluation of a Taylor model of form
#
# ```math
#     m_ℓ(\mathbf x) = f_ℓ(\mathbf x_0) +
#     \mathbf g^T ( \mathbf x - \mathbf x_0 ) + ( \mathbf x - \mathbf x_0 )^T \mathbf H_ℓ ( \mathbf x - \mathbf x_0)
# ```
# is straightforward:
"Evaluate (internal) output `ℓ` of TaylorModel `tm`, provided a difference vector `h = x - x0`."
function _eval_models( tm :: TaylorModel, h :: Vec, ℓ :: Int )
    ret_val = tm.fx0[ℓ] + tm.g[ℓ]'h
    if !isnothing(tm.H)
        ret_val += .5 * h'tm.H[ℓ]*h
    end
    return ret_val
end

"Evaluate (internal) output `ℓ` of `tm` at scaled site `x̂`."
function eval_models( tm :: TaylorModel, x̂ :: Vec, ℓ :: Int )
    h = x̂ .- tm.x0
    return _eval_models( tm, h, ℓ)
 end

# For the vector valued model, we iterate over all (internal) outputs:
function eval_models( tm :: TaylorModel, x̂ :: Vec )
    h = x̂ .- tm.x0
    return [ _eval_models(tm, h, ℓ) for ℓ = eachindex(tm.g)]
end

# The gradient of ``m_ℓ`` can easily be determined:
function get_gradient( tm :: TaylorModel, x̂ :: Vec, ℓ :: Int)
    if isnothing(tm.H)
        return tm.g[ℓ]
    else
        h = x̂ .- tm.x0
        return tm.g[ℓ] .+ .5 * ( tm.H[ℓ]' + tm.H[ℓ] ) * h
    end
end

# And for the Jacobian, we again iterate:
function get_jacobian( tm :: TaylorModel, x̂ :: Vec )
    grad_list = [ get_gradient(tm, x̂, ℓ) for ℓ=eachindex( tm.g ) ]
    return transpose( hcat( grad_list... ) )
end

# ## [Summary & Quick Examples](@id taylor_summary)

# 1. The recommended way to use Finite Difference Taylor models is to define them
#    with TaylorConfig, i.e.,
#    ```julia
#    add_objective!(mop, f, TaylorConfig())
#    ```
# 2. To use `FiniteDiff.jl` instead, do
#    ```julia
#    add_objective!(mop, f, TaylorApproximateConfig(; mode = :fdm))
#    ```
# 3. Have callbacks for the gradients and the Hessians? Great!
#    ```julia
#    add_objective!(mop, f, TaylorCallbackConfig(; degree = 1, gradients = [g1,g2]))
#    ```
# 4. No callbacks, but you want the correct matrices anyways? `ForwardDiff` to the rescue:
#    ```julia
#    add_objective!(mop, f, TaylorApproximateConfig(; degree = 2, mode = :autodiff)
#    ```

# ### Complete usage example
# ```julia
# using Morbit
# Morbit.print_all_logs()
# mop = MixedMOP(3)
#
# add_objective!( mop, x -> sum( ( x .- 1 ).^2 ), Morbit.TaylorApproximateConfig(;degree=2,mode=:fdm) )
# add_objective!( mop, x -> sum( ( x .+ 1 ).^2 ), Morbit.TaylorApproximateConfig(;degree=2,mode=:autodiff) )
#
# x_fin, f_fin, _ = optimize( mop, [-π, ℯ, 0])
# ```


text:
export ARGStream

# TODO: support Base.IO interface

"""
CLI Arguement Stream. Acts like a stream and a Array (support indexing),
it takes in a list of `String` and return an `String` each time
you move its `ptr`.
"""
mutable struct ARGStream
    head::Int
    data::Vector{String}
end

ARGStream(data::Vector{String}) = ARGStream(1, data)

# Stream Interface
import Base: eof
eof(stream::ARGStream) = stream.head > endof(stream.data)
position(stream::ARGStream) = stream.head
function seek(stream::ARGStream, pos::Int)
    # @assert 1 <= pos <= length(stream.data)
    stream.head = pos
    return stream
end
close(stream::ARGStream) = 0

"""
Copyed from Base.Docs for ARGStream. It keeps the `ptr` position
if function `f` returns `nothing` or `false`.
"""
function withstream(f::Function, stream::ARGStream)
    pos = position(stream)
    result = f()
    (result === nothing || result == false) && seek(stream, pos)
    return result
end

"""
flush the rest arguemnts
"""
function flush(stream::ARGStream)
    head = stream.head
    stream.head = endof(stream.data)
    return stream.data[head:end]
end

# Array Interface
import Base: eltype, length, size, getindex, eachindex, endof
eltype(stream::ARGStream) = String
length(stream::ARGStream) = length(stream.data) - stream.head
size(stream::ARGStream) = (length(stream), )
getindex(stream::ARGStream, index) = getindex(stream.data, index)
eachindex(stream::ARGStream) = stream.head:endof(stream.data)
endof(stream::ARGStream) = endof(stream.data)

# Read
export readarg

"""
    readarg(stream::ARGStream)

read an argument
"""
function readarg(from::ARGStream)
    if eof(from)
        return nothing
    end

    head = from.head
    from.head += 1
    return from.data[head]
end

"""
    readarg(stream::ARGStream, size::Int)

read `size` arguments
"""
function readarg(from::ARGStream, size::Int)
    head = from.head
    tail = from.head + size - 1
    tail = tail > length(ARGS) ? length(ARGS) : tail
    from.head = tail + 1
    return from.data[head:tail]
end


text:
# For convenience
const Some{T,N} = Tuple{T, Vararg{T,N}}
const TupleN{T,N} = NTuple{N,T}
const Anys{N} = NTuple{N,Any}

struct True end
struct False end
const TypedBool = Union{True, False}

@pure Base.:(!)(::False) = True()
@pure Base.:(!)(::True) = False()

@pure invert(t::TupleN{TypedBool}) = tuple_map(!, t)

@pure tuple_getindex(xs::Anys{N}, I::NTuple{N,TypedBool}) where {N} = _tuple_getindex(xs, I)
@inline function _tuple_getindex(xs::Some{Any,N}, I::Some{TypedBool,N}) where {N}
    rest  = _tuple_getindex(tail(xs), tail(I))
    return first(I) === True() ? (first(xs), rest...) : rest
end
_tuple_getindex(::Tuple{}, ::Tuple{}) = ()

@inline function tuple_setindex(t::Anys{N}, v::Tuple, I::NTuple{N,TypedBool}) where {N}
    _tuple_setindex(t, v, I, tuple_getindex(I, I))
end
@noinline function _tuple_setindex(t::Tuple, v::Tuple, I::Tuple, ::Anys{M}) where {M}
    throw_setindex_mismatch(M, length(v))
end
@pure function _tuple_setindex(t::Tuple, v::Anys{M}, I::Tuple, ::Anys{M}) where {M}
    __tuple_setindex(t, v, I)
end
@inline function __tuple_setindex(t::Some{Any,N}, v::Some{Any}, I::Some{TypedBool,N}) where {N}
    if first(I) === True()
        (first(v), __tuple_setindex(tail(t), tail(v), tail(I))...)
    else
        (first(t), __tuple_setindex(tail(t), v, tail(I))...)
    end
end
__tuple_setindex(t::Tuple{}, v::Tuple{}, I::Tuple{}) = ()
__tuple_setindex(t::Some{Any,N}, v::Tuple{}, I::Some{False,N}) where {N} = t

@inline tuple_map(f, t::Anys{N}) where {N} = ntuple(i -> (@_inline_meta; f(t[i])), Val(N))

@inline function tuple_split(t::NTuple{L,Any}, ::Val{M}) where {L,M}
    ntuple(i -> t[i], Val(M)), ntuple(i -> t[M+i], Val(L-M))
end


text:
using Symbolics

@variables q[1:4]

T(q) = [
    cos(q) -sin(q)
    sin(q) cos(q)
]

r = [
    1
    2
]

# T_all(q) = T(q[1]) * T(q[2])

# x = T_all(q)
T_(q) = T(q[1])*r
x = T_(q)
J = Symbolics.jacobian(x, [q[1], q[2], q[3], q[4]])
println(J)


text:
################################################################################
###################### DEFINE MASTER PARAMETERS FOR THIS FILE ##################
################################################################################

using Mimi
using Distributions
using CSVFiles
using DataFrames
using CSV
using Random
using StatsBase
using Statistics

include("main_model_growth.jl")

# define the model regions in the order that Mimi returns stuff
myregions = ["EU", "USA", "Other OECD","Former USSR","China","Southeast Asia","Africa","Latin America"]
myyears = [2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200, 2250, 2300]

# define the output directory
dir_output = joinpath(@__DIR__, "../../output/")

# define number of Monte Carlo runs
if !@isdefined samplesize
    samplesize = 50000
end

# define the seed
masterseed = 22081994

# define the pulse size
scc_pulse_size = 75000.

################################################################################
################### MC MODEL RUNS    ###########################################
################################################################################

### PART I: results for fixed persistence values

# create a DataFrame where results will be stored
df_sccMC_singleGE = DataFrame(damagePAGE09=false, permafr=false, seaice=false, ge=-999., scen="-999", pulse_size=-999.,
                                 civvalue=-999.,
                                 ge_adapt = -999.,
                                  mean=-999., median=-999., min=-999., max=-999., perc25=-999.,
                                  perc75=-999., sd=-999., varcoeff=-999., perc05=-999., perc95=-999.,
                                  perc10=-999., perc90=-999.,
                                  share_zeroSCC = -999.)

# loop through all relevant parameter combinations and save out results and SCC summary stats
for jj_page09damages in [false] # whether or not to use the market damage function from PAGE09 instead of the new PAGE-ICE function using Burke et al. (2015)
    for jj_permafr in [true] # whether or not to include the permafrost carbon feedback, see Yumashev et al. (2019)
        for jj_seaice in [true] # whether or not to include the surface albedo feedback, see Yumashev et al. (2019)
            for jj_scen in ["RCP4.5 & SSP2"]
                for jj_civvalue in [1., 10.0^20] # whether or not to hold the civilization value damage cap at its default value

                    # skip undesired combinations of parameters
                    if jj_page09damages != false && (jj_permafr != true  || jj_seaice != true || jj_civvalue != 1.)
                        continue
                    elseif jj_permafr != true && (jj_seaice != true || jj_civvalue != 1.)
                        continue
                    elseif jj_permafr != jj_seaice
                        continue
                    end

                    # print the parameter combination to track progress
                    print(string("page09_", jj_page09damages, " Permafr_", jj_permafr, " Scen_", jj_scen, " Seaice_", jj_seaice
                                     ))

                    # loop through the different values for persistence and its decay rate (= adaptation)
                    for jj_ge in [0.0:0.1:1.0;]
                        for jj_geadapt in [0.0:0.005:0.04;]

                            # skip redundant iterations with non-zero adaptation and zero persistence
                            if jj_ge == 0. && jj_geadapt > 0.
                                continue
                            elseif jj_geadapt != 0. && jj_civvalue != 1.
                                continue
                            end

                            # print out the growth effects to track progress
                            print(string("Rho = ", jj_ge, " ge_adapt = ", jj_geadapt))

                            # define the output directory for the Monte Carlo files
                            dir_MCoutput = string(dir_output, "montecarlo-singleGE/scen", jj_scen, "_permafr", jj_permafr, "_seaice", jj_seaice,
                                                            "_page09", jj_page09damages, "/", "ge", jj_ge,
                                                            "_geadapt", jj_geadapt,
                                                            "_civ", jj_civvalue, "/")

                            # fix the seed and calculate the SCC using a triangular distribution collapsing to a single value
                            Random.seed!(masterseed)
                            global scc_mcs_object = get_scc_mcs(samplesize, 2020, dir_MCoutput,
                                                                    scenario=jj_scen,
                                                                    pulse_size=scc_pulse_size,
                                                                    use_permafrost=jj_permafr,
                                                                    use_seaice=jj_seaice,
                                                                    use_page09damages=jj_page09damages,
                                                                    ge_minimum=jj_ge,
                                                                    ge_maximum=jj_ge + 10^(-30),
                                                                    ge_mode=jj_ge,
                                                                    civvalue_multiplier=jj_civvalue,
                                                                    geadrate = jj_geadapt)

                            # write out the full distribution
                            writedlm(string(dir_MCoutput, "SCC_MC.csv"), scc_mcs_object, ",")

                            # write results into the data frame
                            push!(df_sccMC_singleGE, [jj_page09damages, jj_permafr, jj_seaice, jj_ge,
                                                          jj_scen, scc_pulse_size, jj_civvalue,
                                                          jj_geadapt,
                                                          mean(scc_mcs_object[:, 1]),
                                                          median(scc_mcs_object[:, 1]),
                                                          minimum(scc_mcs_object[:, 1]),
                                                          maximum(scc_mcs_object[:, 1]),
                                                          StatsBase.percentile(scc_mcs_object[:, 1], 25),
                                                          StatsBase.percentile(scc_mcs_object[:, 1], 75),
                                                          Statistics.std(scc_mcs_object[:, 1]),
                                                          Statistics.std(scc_mcs_object[:, 1]) / mean(scc_mcs_object[:, 1]),
                                                          StatsBase.percentile(scc_mcs_object[:, 1], 5),
                                                          StatsBase.percentile(scc_mcs_object[:, 1], 95),
                                                          StatsBase.percentile(scc_mcs_object[:, 1], 10),
                                                          StatsBase.percentile(scc_mcs_object[:, 1], 90),
                                                          count(scc_mcs_object .== 0.) / samplesize])

                          # clean out the MCS objects
                            scc_mcs_object = nothing

                        end
                    end
                end
            end
        end
    end
end

# remove the first placeholder row
df_sccMC_singleGE = df_sccMC_singleGE[df_sccMC_singleGE[!, :scen] .!= "-999", :]

# export the results
CSV.write(string(dir_output, "MimiPageGrowthEffectsResults_SCC_fixedGE.csv"), df_sccMC_singleGE)

### PART II: results for the EMPIRICAL persistence distribution

# create a DataFrame where results will be stored
df_sccMC = DataFrame(permafr=false, seaice=false, ge_string="-999", scen="-999",
                                  convergence=-999.,  bound=-999., eqwshare=-999.,
                                  civvalue=-999., pulse=-999., ge_adapt = -999.,
                                  emfeed = -999., switch = -999.,
                                  mean=-999., median=-999., min=-999., max=-999., perc25=-999.,
                                  perc75=-999., sd=-999., varcoeff=-999.,
                                  perc05=-999., perc95=-999., perc10=-999., perc90=-999.,
                                  share_zeroSCC = -999., count_NaN = -999)

# loop through all relevant parameter/scenario combinations and save out results and SCC summary stats
for jj_scen in ["RCP4.5 & SSP2", "RCP2.6 & SSP1", "RCP8.5 & SSP5", "RCP1.9 & SSP1"]
    for jj_gestring in ["EMPIRICAL"] # which persistence distribution to use (use "EMPIRICAL+" for a version where negative persistence draws are trimmed to zero)
        for jj_permafr in [true]
            for jj_seaice in [true]
                for jj_civvalue in [1., 10.0^20]
                    for jj_cbabs in [740.65] # the level of the subsistence boundary for per capita consumption
                        for jj_eqwshare in [0.99, 0.95, 0.999] # the share of per capita consumption until which damages are equity-weighted
                            for jj_convergence in [1., 0.] # whether or not to use the convergence system for the subsistence and equity-weighting thresholds
                                for jj_pulse in [scc_pulse_size] # the size of the CO2 pulse used for SCC calculations
                                    for jj_geadapt in [0.0:0.01:0.05;] # the decay rate of persistence over time
                                        for jj_emfeed in [1., 0.] # whether or not persistence-related GDP reductions cause proportional reductions in emissions
                                            for jj_regionswitch in [0., 1.] # whether or not persistence is only limited to Global South regions

                                        # overwrite jj_geadapt with nothing if its zero (which will triger the deterministic default value of zero)
                                        if jj_geadapt == 0.
                                            jj_geadapt = nothing
                                        end

                                        # skip undesired or infeasible combinations of parameters/scenarios
                                        if jj_scen != "RCP4.5 & SSP2" && (jj_gestring != "EMPIRICAL" || jj_permafr != true ||
                                                                            jj_seaice != true || jj_civvalue != 1. || jj_cbabs != 740.65 ||
                                                                            jj_eqwshare != 0.99 || jj_convergence != 1. || jj_pulse != scc_pulse_size ||
                                                                            jj_geadapt != nothing || jj_emfeed != 1. || jj_regionswitch != 0.)
                                            continue
                                        elseif jj_permafr != jj_seaice
                                            continue
                                        elseif jj_gestring != "EMPIRICAL" && (jj_cbabs != 740.65 || jj_eqwshare != 0.99 || jj_convergence != 1. ||
                                                                                jj_pulse != scc_pulse_size || jj_geadapt != nothing  || jj_emfeed != 1. || jj_regionswitch != 0.)
                                            continue
                                        elseif jj_civvalue != 1. && (jj_cbabs != 740.65 || jj_eqwshare != 0.99 || jj_convergence != 1. ||
                                                                    jj_pulse != scc_pulse_size || jj_geadapt != nothing || jj_emfeed != 1. || jj_regionswitch != 0.)
                                            continue
                                        elseif jj_cbabs != 740.65 && (jj_civvalue != 1. || jj_gestring != "EMPIRICAL" || jj_eqwshare != 0.99 ||
                                                                    jj_convergence != 1. || jj_pulse != scc_pulse_size || jj_geadapt != nothing || jj_emfeed != 1. || jj_regionswitch != 0.)
                                            continue
                                        elseif jj_eqwshare != 0.99 && (jj_civvalue != 1. || jj_gestring != "EMPIRICAL" || jj_convergence != 1. ||
                                                                        jj_pulse != scc_pulse_size || jj_geadapt != nothing || jj_emfeed != 1. || jj_regionswitch != 0.)
                                            continue
                                        elseif jj_convergence != 1. && (jj_civvalue != 1. || jj_gestring != "EMPIRICAL" || jj_pulse != scc_pulse_size ||
                                                                        jj_geadapt != nothing  || jj_emfeed != 1. || jj_regionswitch != 0.)
                                            continue
                                        elseif jj_pulse != scc_pulse_size && (jj_geadapt != nothing  || jj_emfeed != 1. || jj_regionswitch != 0.)
                                            continue
                                        elseif jj_geadapt != nothing  && (jj_emfeed != 1. || jj_regionswitch != 0.)
                                            continue
                                        elseif jj_emfeed != 1. && jj_regionswitch != 0.
                                            continue
                                        end

                                        # set the ge_string parameters
                                        if jj_gestring == "MILD"
                                            ge_string_min = 0.
                                            ge_string_mode = 0.
                                            ge_string_max = 0.5280681 # value based on converging ratio of marginal impacts from Burke et al regression for 1Lag (temperature and precipitation)
                                            ge_use_empirical = 0.
                                        elseif jj_gestring == "MEDIUM"
                                            ge_string_min = 0.
                                            ge_string_mode = 0.5280681
                                            ge_string_max = 1.
                                            ge_use_empirical = 0.
                                        elseif jj_gestring == "EMPIRICAL" # note: for EMPIRICAL and EMPIRICAL+, the triangular will be overwritten by draws from the empirical distribution
                                            ge_string_min = 0.
                                            ge_string_mode = 0.
                                            ge_string_max = 0.0001
                                            ge_use_empirical = 1.
                                        elseif jj_gestring == "EMPIRICAL+"
                                            ge_string_min = 0.
                                            ge_string_mode = 0.
                                            ge_string_max = 0.0001
                                            ge_use_empirical = 2.
                                        end

                                        # print out the parameters to track progress
                                        print(string("Scen_", jj_scen, " Ge-mode_", jj_gestring,
                                                     " Permafr_", jj_permafr, " Seaice_", jj_seaice, " Civvalue_", jj_civvalue,
                                                     " cbabs_", jj_cbabs, " eqwshare_", jj_eqwshare,
                                                     " convergence", jj_convergence, "pulse_", jj_pulse,
                                                     " geadapt", jj_geadapt, " emfeed", jj_emfeed,
                                                     " switch", jj_regionswitch))

                                        # define the output directory for the Monte Carlo files
                                        dir_MCoutput = string(dir_output, "mc_diGE/ge", jj_gestring,
                                                                        "_scen", jj_scen,
                                                                        "_pf", jj_permafr,
                                                                        "_se", jj_seaice,
                                                                        "_co", jj_convergence, "_bd", jj_cbabs, "_eq", jj_eqwshare,
                                                                        "_ci", jj_civvalue,
                                                                        "_p", jj_pulse,
                                                                        "_gead", jj_geadapt,
                                                                        "_emfe", jj_emfeed,
                                                                        "_sw", jj_regionswitch,
                                                                         "/")

                                        # calculate the MC SCC
                                        Random.seed!(masterseed)
                                        global scc_mcs_object = get_scc_mcs(samplesize, 2020, dir_MCoutput,
                                                                            scenario=jj_scen,
                                                                            pulse_size=jj_pulse,
                                                                            use_permafrost=jj_permafr,
                                                                            use_seaice=jj_seaice,
                                                                            use_page09damages=false,
                                                                            ge_minimum=ge_string_min,
                                                                            ge_maximum=ge_string_max,
                                                                            ge_mode=ge_string_mode,
                                                                            ge_use_empirical=ge_use_empirical,
                                                                            civvalue_multiplier=jj_civvalue,
                                                                            use_convergence=jj_convergence,
                                                                            cbabs=jj_cbabs,
                                                                            eqwbound=jj_eqwshare,
                                                                            geadrate = jj_geadapt,
                                                                            emfeedback = jj_emfeed,
                                                                            ge_use_switch = jj_regionswitch
                                                                            )

                                        # write out the full distribution
                                        writedlm(string(dir_MCoutput, "SCC_MC.csv"),
                                                                scc_mcs_object, ",")

                                        # create a NaN-free version of the SCC distribution
                                        scc_mcs_object = filter(!isnan, scc_mcs_object)

                                        # push the results into the data frame
                                        push!(df_sccMC, [jj_permafr, jj_seaice, jj_gestring,  jj_scen,
                                                        jj_convergence,  jj_cbabs, jj_eqwshare,
                                                        jj_civvalue, jj_pulse,
                                                        ifelse(jj_geadapt == nothing, 0., jj_geadapt),
                                                        jj_emfeed, jj_regionswitch,
                                                        mean(scc_mcs_object[:, 1]),
                                                        median(scc_mcs_object[:, 1]),
                                                        minimum(scc_mcs_object[:, 1]),
                                                        maximum(scc_mcs_object[:, 1]),
                                                        StatsBase.percentile(scc_mcs_object[:, 1], 25),
                                                        StatsBase.percentile(scc_mcs_object[:, 1], 75),
                                                        Statistics.std(scc_mcs_object[:, 1]),
                                                        Statistics.std(scc_mcs_object[:, 1]) / mean(scc_mcs_object[:, 1]),
                                                        StatsBase.percentile(scc_mcs_object[:, 1], 5),
                                                        StatsBase.percentile(scc_mcs_object[:, 1], 95),
                                                        StatsBase.percentile(scc_mcs_object[:, 1], 10),
                                                        StatsBase.percentile(scc_mcs_object[:, 1], 90),
                                                        count(scc_mcs_object .== 0.) / length(scc_mcs_object),
                                                        samplesize - length(scc_mcs_object)])

                                        # clean out the scc_mcs_object
                                        global scc_mcs_object = nothing
                                        ge_string_min = nothing
                                        ge_string_max = nothing
                                        ge_string_mode = nothing

                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

# remove the first placeholder row
df_sccMC = df_sccMC[df_sccMC[!, :scen] .!= "-999", :]

# export the results
CSV.write(string(dir_output, "MimiPageGrowthEffectsResultsMonteCarlo.csv"), df_sccMC)

### PART III: run the default scenario/parameter settings with an empirical persistence distribution derived from low-pass filtering climate vars

# rerun the main settings for an empirical distribution that is derived from low-pass filtering climate variables
include("main_model_growth_lowpass.jl")

# set the output directory for MC results
dir_MCoutput = string(dir_output, "lowpassfiltered_persistence/")

# calculate the stochastic mean SCC using the default parameter & scenario settings
Random.seed!(masterseed)
global scc_mcs_object = get_scc_mcs(samplesize, 2020, dir_MCoutput,
                                    scenario="RCP4.5 & SSP2",
                                    pulse_size=scc_pulse_size,
                                    use_permafrost=true,
                                    use_seaice=true,
                                    use_page09damages=false,
                                    ge_minimum=0.,
                                    ge_maximum=0.0001,
                                    ge_mode=0.,
                                    ge_use_empirical=1.,
                                    civvalue_multiplier=1.,
                                    use_convergence= 1.,
                                    cbabs=740.65,
                                    eqwbound=0.99,
                                    geadrate = 0.,
                                    emfeedback = 1.,
                                    ge_use_switch = 0.
                                    )

# # write out the full distribution
writedlm(string(dir_MCoutput, "SCC_MC.csv"),
                        scc_mcs_object, ",")

# clean out
global scc_mcs_object = nothing


text:
uniform(u::Real) = 1//2 * float64(abs(u) <= 1.0)
uniform_unnormalized(u::Real) = 1//2 * float64(abs(u) <= 1.0)
Base.@vectorize_1arg Real uniform
Base.@vectorize_1arg Real uniform_unnormalized

triangular(u::Real) = (1 - abs(u)) * float64(abs(u) <= 1.0)
triangular_unnormalized(u::Real) = (1 - abs(u)) * float64(abs(u) <= 1.0)
Base.@vectorize_1arg Real triangular
Base.@vectorize_1arg Real triangular_unnormalized

epanechnikov(u::Real) = 3//4 * (1 - u^2) * float64(abs(u) <= 1.0)
epanechnikov_unnormalized(u::Real) = (1 - u^2) * float64(abs(u) <= 1.0)
Base.@vectorize_1arg Real epanechnikov
Base.@vectorize_1arg Real epanechnikov_unnormalized

biweight(u::Real) = 15//16 * (1 - u^2)^2 * float64(abs(u) <= 1.0)
biweight_unnormalized(u::Real) = (1 - u^2)^2 * float64(abs(u) <= 1.0)
Base.@vectorize_1arg Real biweight
Base.@vectorize_1arg Real biweight_unnormalized

triweight(u::Real) = 35//32 * (1 - u^2)^3 * float64(abs(u) <= 1.0)
triweight_unnormalized(u::Real) = (1 - u^2)^3 * float64(abs(u) <= 1.0)
Base.@vectorize_1arg Real triweight
Base.@vectorize_1arg Real triweight_unnormalized

tricube(u::Real) = 70//81 * (1 - abs(u)^3)^3 * float64(abs(u) <= 1.0)
tricube_unnormalized(u::Real) = (1 - abs(u)^3)^3 * float64(abs(u) <= 1.0)
Base.@vectorize_1arg Real tricube
Base.@vectorize_1arg Real tricube_unnormalized

gaussian(u::Real) = (1 / sqrt(2 * pi)) *  exp(-1//2 * u^2)
gaussian_unnormalized(u::Real) = exp(-1//2 * u^2)
Base.@vectorize_1arg Real gaussian
Base.@vectorize_1arg Real gaussian_unnormalized

cosine(u::Real) = (pi / 4) * cos((pi / 2) * u) * float64(abs(u) <= 1.0)
cosine_unnormalized(u::Real) = cos((pi / 2) * u) * float64(abs(u) <= 1.0)
Base.@vectorize_1arg Real cosine
Base.@vectorize_1arg Real cosine_unnormalized

logistic(u::Real) = 1 / (exp(u) + 2 + exp(-u))
logistic_unnormalized(u::Real) = 1 / (exp(u) + 2 + exp(-u))
Base.@vectorize_1arg Real logistic
Base.@vectorize_1arg Real logistic_unnormalized


text:
# bogaczki shampine
struct BogaczkiShampine23 <: AbstractAdaptiveExplicitSolver
end

function order(::BogaczkiShampine23)::Int
    return 3
end

function n_intermediates(::BogaczkiShampine23)
    return 4
end

function n_temporaries(::BogaczkiShampine23)
    return 2
end

function stepper!(ode_solver::BogaczkiShampine23, ode_fcn, stepper_options, stepper_temp::StepperTemporaries)

    # renaming for better readability
    t          = stepper_temp.t_old;
    x          = stepper_temp.x_old;
    dt         = stepper_temp.dt;
    n_prob_dim = n_dim(stepper_temp)

    # views into the temporaries
    k1 = view(stepper_temp.der_x, 1:n_prob_dim, 1)
    k2 = view(stepper_temp.der_x, 1:n_prob_dim, 2)
    k3 = view(stepper_temp.der_x, 1:n_prob_dim, 3)
    k4 = view(stepper_temp.der_x, 1:n_prob_dim, 4)

    x2 = view(stepper_temp.x_tmp, 1:n_prob_dim, 1)
    x3 = view(stepper_temp.x_tmp, 1:n_prob_dim, 2)

    # intermediate steps
    t2 = t + 1/2 * dt
    @. x2 = x + 1/2 * dt * k1
    ode_fcn(k2, t2, x2)

    t3 = t + 3/4 * dt
    @. x3 = x + 3/4 * dt * k2
    ode_fcn(k3, t3, x3)

    # new step
    stepper_temp.t_new = t + dt
    @. stepper_temp.x_new = x + dt * (2/9 * k1 + 1/3 * k2 + 4/9 * k3)
    ode_fcn(k4, stepper_temp.t_new, stepper_temp.x_new)

    # error estimation
    @. stepper_temp.x_est = x + dt * (7/24 * k1 + 1/4 * k2 + 1/3 * k3 + 1/8 * k4)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
temporary_variables(cache::Tsit5ConstantCache) = (;)

function temporary_variables(cache::Tsit5Cache)
    @unpack k1, k2, k3, k4, k5, k6, tmp = cache
    return (;k1, k2, k3, k4, k5, k6, tmp)
end