text:
# Functions for checking blanchard_kahn conditions
function blanchard_kahn(m::RiskAdjustedLinearization; verbose::Symbol = :high)
    Γ₅ = m.linearization.Γ₅
    Γ₆ = m.linearization.Γ₆
    A = [Γ₅ Γ₆; Matrix{eltype(Γ₅)}(I, m.Nz, m.Nz) zeros(m.Ny, m.Nz)]
    B = [(-Γ₃(m) - JV(m)) (-Γ₄(m)); Γ₁(m) Γ₂(m)]

    if count(abs.(eigen(A, B).values) .> 1) != m.Nz
        BlanchardKahnError("First-order perturbation around stochastic steady state is not saddle-path stable")
    else
        if verbose in [:low, :high]
            println("Blanchard-Kahn conditions for a unique locally bounded stochastic steady-state perturbation are satisfied")
        end

        return true
    end
end

mutable struct BlanchardKahnError <: Exception
    msg::String
end
BlanchardKahnError() = BlanchardKahnError("First-order perturbation is not saddle-path stable")
Base.showerror(io::IO, ex::BlanchardKahnError) = print(io, ex.msg)


text:
using Test, ReachabilityBenchmarks

@testset "Relative path macro" begin
    file = @relpath "my_data.dat"
    @test file == joinpath(@__DIR__, "my_data.dat")
end

ext: jl

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
using COSMO, Random, Test, Pkg
rng = Random.MersenneTwister(12345)

include("./UnitTests/COSMOTestUtils.jl")
# Define the types to run the unit tests with
UnitTestFloats = [Float32; Float64; BigFloat]

@testset "COSMO Native Testset" begin

  include("./UnitTests/sets.jl")
  include("./UnitTests/nuclear_norm_minimization.jl")
  include("./UnitTests/psd_completion.jl")
  include("./UnitTests/psd_completion_and_merging.jl")
  include("./UnitTests/clique_merging_example.jl")
  include("./UnitTests/reduced_clique_graph.jl")
  include("./UnitTests/chordal_decomposition_triangle.jl")
  include("./UnitTests/InfeasibilityTests/runTests.jl")
  include("./UnitTests/Acceleration/run_acceleration_tests.jl")
  # floating-point precision type agnostic unit tests
  # N.B.: Some infeasibility tests rely on sufficient precision (at least Float64 in some cases)
  # SDPs are currently only supported with Float32/Float64
  include("./UnitTests/simple.jl")
  include("./UnitTests/model_modifications.jl")
  include("./UnitTests/constraints.jl")
  include("./UnitTests/model.jl")
  include("./UnitTests/socp-lasso.jl")
  include("./UnitTests/closestcorr.jl")
  include("./UnitTests/exp_cone.jl")
  include("./UnitTests/pow_cone.jl")
  include("./UnitTests/qp-box.jl")
  include("./UnitTests/algebra.jl")
  include("./UnitTests/splitvector.jl")
  include("./UnitTests/interface.jl")
  include("./UnitTests/interface_python.jl")
  include("./UnitTests/kktsolver.jl")
  include("./UnitTests/print.jl")

    # optional unittests
  if pkg_installed("Pardiso", "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2")
    include("./UnitTests/options_factory.jl")
  end

end
nothing


text:
###############################################################################
# Stochiometry for a given reaction is a vector of pairs mapping species id to
# stochiometric coefficient.
###############################################################################

@fastmath function evalrxrate(speciesvec::AbstractVector{T}, rateconst,
                              stochmat::AbstractVector{Pair{S,V}})::typeof(rateconst) where {T,S,V}
    val = one(T)

    @inbounds for specstoch in stochmat
        specpop = speciesvec[specstoch[1]]
        val    *= specpop
        @inbounds for k = 2:specstoch[2]
            specpop -= one(specpop)
            val     *= specpop
        end
    end

     rateconst * val
end

@inline @fastmath function executerx!(speciesvec::AbstractVector{T},
                                      net_stoch::AbstractVector{Pair{S,V}}) where {T,S,V}
    @inbounds for specstoch in net_stoch
        speciesvec[specstoch[1]] += specstoch[2]
    end
    nothing
end

@inbounds function scalerates!(unscaled_rates, stochmat::Vector{Vector{Pair{S,T}}}) where {S,T}
    for i in eachindex(unscaled_rates)
        coef = one(T)
        for specstoch in stochmat[i]
            coef *= factorial(specstoch[2])
        end
        unscaled_rates[i] /= coef
    end
    nothing
end

###############################################################################
# Stochiometry for a given reaction is a vector of 2xN matrices. 1st row is
# species indices and second row is stoichiometry
###############################################################################

@fastmath function evalrxrate(speciesvec::AbstractVector{T}, rateconst, stochmat::AbstractArray{S,2}) where {T,S}
    val = one(T)

    i = 1
    while i < length(stochmat)
        @inbounds specpop = speciesvec[stochmat[i]]
        val    *= specpop
        i      += 1
        @inbounds for k = 2:stochmat[i]
            specpop -= one(specpop)
            val     *= specpop
        end
        i += 1
    end

    return rateconst * val
end

@inline @fastmath function executerx!(speciesvec::AbstractVector{T}, net_stoch::AbstractArray{S,2}) where {T,S}
    @inbounds for i = 1:2:length(net_stoch)
        speciesvec[net_stoch[i]] += net_stoch[i+1]
    end
    nothing
end

@inbounds function scalerates!(unscaled_rates, stochmat::Vector{Array{T,2}}) where {T}
    for i in eachindex(unscaled_rates)
        coef = one(T)
        stoch = stochmat[i]
        for k in 2:2:length(stoch)
            coef *= factorial(stoch[k])
        end
        unscaled_rates[i] /= coef
    end
    nothing
end

# OLD!!!

# @inbounds @fastmath function evalrxrate(speciesvec::AbstractVector{T}, rateconst, stochmat::AbstractArray{T,2}) where T
#     val = one(T)

#     i = 1
#     while i < (length(stochmat))
#         specpop = speciesvec[stochmat[i]]
#         val    *= specpop
#         i      += 1
#         for k = 2:stochmat[i]
#             specpop -= one(specpop)
#             val     *= specpop
#         end
#         i += 1
#     end

#     return rateconst * val
# end

# @inline @inbounds @fastmath function executerx!(speciesvec::AbstractVector{T}, net_stoch::AbstractArray{T,2}) where T
#     for i = 1:2:length(net_stoch)
#         speciesvec[net_stoch[i]] += net_stoch[i+1]
#     end
#     nothing
# end

# @inbounds @fastmath function evalrxrate(speciesvec::AbstractVector{T}, rateconst,
#                                         stochmat::AbstractVector{Pair{T,T}}) where T
#     val = one(T)

#     for specstoch in stochmat
#         specpop = speciesvec[specstoch[1]]
#         val    *= specpop
#         for k = 2:specstoch[2]
#             specpop -= one(specpop)
#             val     *= specpop
#         end
#     end

#     return rateconst * val
# end

# @inline @inbounds @fastmath function executerx!(speciesvec::AbstractVector{T},
#                                                 net_stoch::AbstractVector{Pair{T,T}}) where T
#     for specstoch in net_stoch
#         speciesvec[specstoch[1]] += specstoch[2]
#     end
#     nothing
# end

max_stars_repo_name: isaacsas/DiffEqJumpExtensions

max_issues_repo_name: isaacsas/DiffEqJumpExtensions

max_forks_repo_name: isaacsas/DiffEqJumpExtensions


text:
# ---
# title: 461. Hamming Distance
# id: problem461
# author: Indigo
# date: 2021-02-04
# difficulty: Easy
# categories: Bit Manipulation
# link: <https://leetcode.com/problems/hamming-distance/description/>
# hidden: true
# ---
#
# The [Hamming distance](https://en.wikipedia.org/wiki/Hamming_distance) between
# two integers is the number of positions at which the corresponding bits are
# different.
#
# Given two integers `x` and `y`, calculate the Hamming distance.
#
# **Note:**
# 0 ≤ `x`, `y` < 231.
#
# **Example:**
#
#
#
#     Input: x = 1, y = 4
#
#     Output: 2
#
#     Explanation:
#     1   (0 0 0 1)
#     4   (0 1 0 0)
#            ↑   ↑
#
#     The above arrows point to positions where the corresponding bits are different.
#
#
#
## @lc code=start
using LeetCode

function hamming_distance(x::Int, y::Int)
    res = 0
    while x != 0 || y != 0
        res += (x & 1) ⊻ (y & 1)
        x >>= 1
        y >>= 1
    end
    res
end
## @lc code=end


text:
using Visualisations
using Test

include( "../src/uses.jl")
include( "../src/examples.jl")
include( "../src/runner_libs.jl" )
include( "../src/display_constants.jl")
include( "../src/static_texts.jl")
include( "../src/table_libs.jl")
include( "../src/text_html_libs.jl")
include( "../src/server_libs.jl")

sys = load_system()
settings = deepcopy(BASE_SETTINGS)

@testset "Examples" begin
    exres = calc_examples( sys, sys, settings )
    println( exres )
    extext = make_examples( exres )
    println( extext )
end

@testset "Parse Params" begin
    req=Dict{Any, Any}(
        :query => "it_basic_rate=20&it_higher_rate=41&it_top_rate=46&it_pa=12750&bi_adult=4800&bi_pensioner=8780&bi_pens_age=66&bi_child=3000&bi_adult_age=17&ubi_mtbens_keep_as_is=false&ubi_mtbens_abolish=false&ubi_mtbens_keep_housing=true&ubi_abolish_sick=true&ubi_abolish_pensions=false&ubi_abolish_esa=true&ubi_abolish_others=false&ubi_as_mt_income=true&ubi_taxable=true",
        :method => "GET",
        :parsed_querystring => Dict{AbstractString, Any}(
            "bi_adult" => 4800,
            "bi_pens_age" => 66,
            "ubi_taxable" => true,
            "ubi_as_mt_income" => true,
            "it_top_rate" => 46,
            "bi_pensioner" => 8780,
            "ubi_abolish_sick" => true,
            "ubi_abolish_pensions" => false,
            "ubi_abolish_others" => false,
            "ubi_mtbens_keep_housing" => true,
            "ubi_mtbens_abolish" => false,
            "it_pa" => 12750,
            "ubi_abolish_esa" => true,
            "it_basic_rate" => 20,
            "bi_adult_age" => 17,
            "it_higher_rate" => 41,
            "ubi_mtbens_keep_as_is" => false,
            "bi_child" => 3000 ))
    msys = web_map_params( req )
    println( msys )
end

@testset "Queuing" begin

        for i in 1:1
            uuid = submit_job( "run$i", sys, settings )
            #=
            @info "iter $i; submitted job $uuid "
            @info "iter $i IN_QUEUE=$IN_QUEUE"
            @info "iter $i OUT_QUEUE=$OUT_QUEUE"
            =#
        end

end


text:
module Knapsack

using ..Diderot
using ..Diderot: Arc, Node, Layer

struct Instance
    values::Vector{Float64}
    weights::Vector{Int}
    capacity::Int

    function Instance(values, weights, capacity)
        @assert length(values) == length(weights)
        new(values, weights, capacity)
    end
end

Base.length(instance::Instance) = length(instance.values)
Diderot.domain_type(instance::Instance) = Bool
Diderot.value_type(instance::Instance) = Float64

Diderot.initial_state(instance::Instance) = instance.capacity

struct DecreasingWeight end

function Diderot.next_variable(instance, diagram, ::DecreasingWeight)
    n = length(instance)
    fixed = Diderot.fixed_variables(diagram)
    # TODO: efficient!
    perm = sortperm(1:n, by=i->instance.weights[i], rev=true)
    for i in perm
        if !(i in fixed)
            return i
        end
    end
    return nothing
end

function Diderot.transitions(instance::Instance, state, variable)
    results = Dict{Arc{Int, Bool, Float64}, Int}()

    # true
    slack = state - instance.weights[variable]
    if slack >= 0
        arc = Arc(state, true, instance.values[variable])
        results[arc] = slack
    end

    # false
    results[Arc(state, false, 0.0)] = state # unchanged

    return results
end

### Relaxation

# For relaxation:
# - Keep states with largest remaining capacity.
# - Merge states by maximum value.
sort_by(tuple) = -tuple.first
merge = max
relax(max_width) = Diderot.RelaxAllInOne(max_width, merge, sort_by)

end # module


text:
module m2test3
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using LinearAlgebra: norm
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3

    m1 = FinEtoolsDeforNonlinear.MatDeforStVKModule.MatDeforStVK(mr, E, nu)
    m2 = FinEtoolsDeforNonlinear.MatDeforStVKADModule.MatDeforStVKAD(mr, E, nu)
    # @show m
    update! = FinEtoolsDeforNonlinear.MatDeforNonlinearModule.update!
    tangentmoduli! = FinEtoolsDeforNonlinear.MatDeforNonlinearModule.tangentmoduli!

    stress = fill(0.0, 6)
    D = fill(0.0, 6, 6)
    Dtrue = [9.708605576923078e6 4.078846153846155e6 4.0788461538461554e6 0.0 0.0 0.0; 4.078846153846155e6 9.32977913175933e6 3.9984767707539997e6 0.0 0.0 0.0; 4.0788461538461554e6 3.9984767707539997e6 9.329779131759332e6 0.0 0.0 0.0; 0.0 0.0 0.0 2.719230769230769e6 0.0 0.0; 0.0 0.0 0.0 0.0 2.719230769230769e6 0.0; 0.0 0.0 0.0 0.0 0.0 2.665651180502665e6]
    output = FFlt[]
    statev = FFlt[]
    tn = 0.0
    dtn = 0.0
    loc = [0.0 0.0 0.0]
    label = 0
    quantity=:nothing

    Fn1 = [1.01 0 0; 0 1.0 0; 0 0 1.0]
    Fn = [1.0 0 0; 0 1.0 0; 0 0 1.0]
    update!(m1, statev, stress, output, Fn1, Fn, tn, dtn, loc, label, quantity)
    # @show stress
    @test norm(stress-[95648.94230769236, 40184.69154607771, 40184.69154607771, 0.0, 0.0, 0.0]) / E < 1.0e-7
    tangentmoduli!(m1, D, statev, Fn1, Fn, tn, dtn, loc, label)
    @test norm(D-Dtrue) / E < 1.0e-7

    Fn1 = [1.01 0 0; 0 1.0 0; 0 0 1.0]
    Fn = [1.0 0 0; 0 1.0 0; 0 0 1.0]
    update!(m2, statev, stress, output, Fn1, Fn, tn, dtn, loc, label, quantity)
    # @show stress
    @test norm(stress-[95648.94230769236, 40184.69154607771, 40184.69154607771, 0.0, 0.0, 0.0]) / E < 1.0e-7
    tangentmoduli!(m2, D, statev, Fn1, Fn, tn, dtn, loc, label)
    @test norm(D-Dtrue) / E < 1.0e-7

    # Fn1 = [1.0 0 0; 0 1.01 0; 0 0 1.0]
    # Fn = [1.0 0 0; 0 1.0 0; 0 0 1.0]
    # update!(m, statev, stress, output, Fn1, Fn, tn, dtn, loc, label, quantity)
    # # @show stress
    # @test norm(stress-[40184.6915460777, 95648.94230769236, 40184.69154607771, 0.0, 0.0, 0.0] ) / E < 1.0e-7

    # Fn1 = [1.0 0 0; 0 1.0 0.01; 0 0 1.0]
    # Fn = [1.0 0 0; 0 1.0 0; 0 0 1.0]
    # update!(m, statev, stress, output, Fn1, Fn, tn, dtn, loc, label, quantity)
    # # @show stress
    # @test norm(stress-[201.92307692305477, 740.4317307692086, 471.1538461537944, 0.0, 0.0, 26927.78846153846]) / E < 1.0e-7
end
end
using .m2test3
m2test3.test()

module m2test4
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using LinearAlgebra: norm
using BenchmarkTools
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3

    m1 = FinEtoolsDeforNonlinear.MatDeforStVKModule.MatDeforStVK(mr, E, nu)
    m2 = FinEtoolsDeforNonlinear.MatDeforStVKADModule.MatDeforStVKAD(mr, E, nu)
    # @show m
    update! = FinEtoolsDeforNonlinear.MatDeforNonlinearModule.update!
    tangentmoduli! = FinEtoolsDeforNonlinear.MatDeforNonlinearModule.tangentmoduli!

    stress = fill(0.0, 6)
    D = fill(0.0, 6, 6)
    Dtrue = [9.708605576923078e6 4.078846153846155e6 4.0788461538461554e6 0.0 0.0 0.0; 4.078846153846155e6 9.32977913175933e6 3.9984767707539997e6 0.0 0.0 0.0; 4.0788461538461554e6 3.9984767707539997e6 9.329779131759332e6 0.0 0.0 0.0; 0.0 0.0 0.0 2.719230769230769e6 0.0 0.0; 0.0 0.0 0.0 0.0 2.719230769230769e6 0.0; 0.0 0.0 0.0 0.0 0.0 2.665651180502665e6]
    output = FFlt[]
    statev = FFlt[]
    tn = 0.0
    dtn = 0.0
    loc = [0.0 0.0 0.0]
    label = 0
    quantity=:nothing

    Fn1 = [1.01 0 0; 0 1.0 0; 0 0 1.0]
    Fn = [1.0 0 0; 0 1.0 0; 0 0 1.0]
    # @btime $update!($m1, $statev, $stress, $output, $Fn1, $Fn, $tn, $dtn, $loc, $label, $quantity)
    update!(m1, statev, stress, output, Fn1, Fn, tn, dtn, loc, label, quantity)
    # @show stress
    @test norm(stress-[95648.94230769236, 40184.69154607771, 40184.69154607771, 0.0, 0.0, 0.0]) / E < 1.0e-7
    # @btime $tangentmoduli!($m1, $D, $statev, $Fn1, $Fn, $tn, $dtn, $loc, $label)
    tangentmoduli!(m1, D, statev, Fn1, Fn, tn, dtn, loc, label)
    @test norm(D-Dtrue) / E < 1.0e-7

    Fn1 = [1.01 0 0; 0 1.0 0; 0 0 1.0]
    Fn = [1.0 0 0; 0 1.0 0; 0 0 1.0]
    # @btime $update!($m2, $statev, $stress, $output, $Fn1, $Fn, $tn, $dtn, $loc, $label, $quantity)
    update!(m2, statev, stress, output, Fn1, Fn, tn, dtn, loc, label, quantity)
    # @show stress
    @test norm(stress-[95648.94230769236, 40184.69154607771, 40184.69154607771, 0.0, 0.0, 0.0]) / E < 1.0e-7
    # @btime $tangentmoduli!($m2, $D, $statev, $Fn1, $Fn, $tn, $dtn, $loc, $label)
    tangentmoduli!(m2, D, statev, Fn1, Fn, tn, dtn, loc, label)
    @test norm(D-Dtrue) / E < 1.0e-7

    # Fn1 = [1.0 0 0; 0 1.01 0; 0 0 1.0]
    # Fn = [1.0 0 0; 0 1.0 0; 0 0 1.0]
    # update!(m, statev, stress, output, Fn1, Fn, tn, dtn, loc, label, quantity)
    # # @show stress
    # @test norm(stress-[40184.6915460777, 95648.94230769236, 40184.69154607771, 0.0, 0.0, 0.0] ) / E < 1.0e-7

    # Fn1 = [1.0 0 0; 0 1.0 0.01; 0 0 1.0]
    # Fn = [1.0 0 0; 0 1.0 0; 0 0 1.0]
    # update!(m, statev, stress, output, Fn1, Fn, tn, dtn, loc, label, quantity)
    # # @show stress
    # @test norm(stress-[201.92307692305477, 740.4317307692086, 471.1538461537944, 0.0, 0.0, 26927.78846153846]) / E < 1.0e-7
end
end
using .m2test4
m2test4.test()


text:
abstract type AstNode end

struct VariableNode <: AstNode
  identifier::String
end

struct AndNode <: AstNode
  lhs::AstNode
  rhs::AstNode
end

struct OrNode <: AstNode
  lhs::AstNode
  rhs::AstNode
end

struct NotNode <: AstNode
  identifier::String
end


text:
using MAT
using Dates
using GLPK
using JuMP
using Statistics
using CartesianGeneticProgramming
using Cambrian
import Cambrian.mutate
using Plots
using Plots.PlotMeasures
using StatsPlots
using DataFrames
using GLM
include("utils.jl")

function get_wave_t_in_range(param, min, max)
    indices_more = findall(>=(min), param)
    indices_less = findall(<(max), param)

    indices_to_keep = intersect(indices_more, indices_less)
    param = param[indices_to_keep]

    param, indices_to_keep
end

function fit_poly_inf_norm(x,y,n)
	m = Model(GLPK.Optimizer)

    @variable(m, alpha[1:n+1])
	@variable(m, t >= 0.)
	@constraint(m, [i=1:length(y)], sum(alpha[j]*x[i]^(j-1) for j in 1:n+1) - y[i] <= t);
	@constraint(m, [i=1:length(y)], -t <= sum(alpha[j]*x[i]^(j-1) for j in 1:n+1) - y[i]);
	@objective(m, Min, t)

	optimize!(m) # added new 2020-04-08

    return JuMP.value.(alpha)
end

f2(x, args) = (args[2] * x) .+ args[1]

const MATLAB_EPOCH = Dates.DateTime(-1,12,31)

date2num(d::Dates.DateTime) = Dates.value(d-MATLAB_EPOCH)/(1000*60*60*24)
num2date(n::Number) =  MATLAB_EPOCH + Dates.Millisecond(round(Int64, n*1000*60*60*24))

function plot_bathy(xs, profile; show_points=false, x_last_positive=NaN, x_first_negative=NaN, value_last_positive=NaN, value_first_negative=NaN, x_at_zero=NaN, title=NaN)
    p1 = plot(xs, profile, label="profile_t0")
    if show_points
        scatter!([x_last_positive, x_first_negative], [value_last_positive, value_first_negative], label="A, B")
        scatter!([x_at_zero], [0], label="Shoreline")
    end
    display(plot(p1, title=title))
end

function find_x_at_zero(profile, xs, i)
    t0_index_last_positive = findlast(profile.>0)
    if t0_index_last_positive == nothing t0_index_last_positive = 1 end

    t0_index_first_negative = findfirst(profile.<0)
    if t0_index_first_negative == nothing t0_index_first_negative = 1 end

    t0_x_last_positive = xs[t0_index_last_positive]
    t0_x_first_negative = xs[t0_index_first_negative]
    t0_value_last_positive = profile[t0_index_last_positive]
    t0_value_first_negative = profile[t0_index_first_negative]

    t0_coeffs = fit_poly_inf_norm([t0_value_last_positive, t0_value_first_negative], [t0_x_last_positive, t0_x_first_negative], 1)
    t0_x_at_zero = f2([0], t0_coeffs)[1]
    t0_x_at_zero
end

dir_path = get_dataset_path(;dataset_name="grand_popo")
kalman_bathy_vars = matread(string(dir_path, "Kalman_bathy.mat"))
bathy_xb = kalman_bathy_vars["xb"][1, :]# 1×213 Matrix{Float64}
# bathy_tn = kalman_bathy_vars["tn"]      # 1×1214 Matrix{Float64}
bathy_hbK = kalman_bathy_vars["hbK"]    # 213×1214 Matrix{Float64}

min_target_matdate = minimum(filter(!isnan, kalman_bathy_vars["tn"]))
max_target_matdate = maximum(filter(!isnan, kalman_bathy_vars["tn"]))

wave_vars = matread(string(dir_path, "Wave_GPP_ERA.mat"))
wave_t = wave_vars["t"]
wave_Tm = wave_vars["Tm"]
wave_hs = wave_vars["hs"]
wave_dir = wave_vars["dir"]

wave_t, indices = get_wave_t_in_range(wave_t, min_target_matdate, max_target_matdate)
wave_Tm = wave_Tm[indices]
wave_hs = wave_hs[indices]
wave_dir = wave_dir[indices]

bathy_dates = map(num2date, collect(min_target_matdate:max_target_matdate))  # collect() instead of using bathy_tn because bathy_tn includes NaN's, even though bathy_hbK has no NaN
wave_dates = map(num2date, wave_t)

days = []
n_offshore = 0
n_landward = 0
n_none = 0

mean_tp = []
min_tp = []
max_tp = []
std_tp = []

mean_hs = []
min_hs = []
max_hs = []
std_hs = []

mean_dir = []
min_dir = []
max_dir = []
std_dir = []

changes = []
for i in 2:length(bathy_dates)
    _d0 = bathy_dates[i-1]
    _d1 = bathy_dates[i]

    _t0 = date2num(_d0)
    _t1 = date2num(_d1)

    _bathy_t0 = bathy_hbK[:, i-1]
    _bathy_t1 = bathy_hbK[:, i]

    _wave_t = deepcopy(wave_t)
    _wave_Tm = deepcopy(wave_Tm)
    _wave_dir = deepcopy(wave_dir)
    _wave_hs = deepcopy(wave_hs)

    _wave_t, _indices = get_wave_t_in_range(_wave_t, _t0, _t1)
    _wave_Tm = _wave_Tm[_indices]
    _wave_hs = _wave_hs[_indices]
    _wave_dir = _wave_dir[_indices]

    t0_x_at_zero = find_x_at_zero(_bathy_t0, bathy_xb, i)
    t1_x_at_zero = find_x_at_zero(_bathy_t1, bathy_xb, i)
    change = (t1_x_at_zero - t0_x_at_zero)# / 100

    # if change > 1 || change < -1
    #     println(i, "    ", change)
    #     exit()
    # end

    if change > 0
        class = (1, 0, 0)  # (offshore, landward, none)
    elseif change < 0
        class = (0, 1, 0)
    else
        class = (0, 0, 1)
    end

    # _wave_Tm = _wave_Tm ./ 100
    # _wave_hs = _wave_hs ./ 10
    # _wave_dir = _wave_dir ./ 100

    push!(mean_tp, mean(_wave_Tm))
    push!(min_tp, minimum(_wave_Tm))
    push!(max_tp, maximum(_wave_Tm))
    push!(std_tp, std(_wave_Tm))

    push!(mean_hs, mean(_wave_hs))
    push!(min_hs, minimum(_wave_hs))
    push!(max_hs, maximum(_wave_hs))
    push!(std_hs, std(_wave_hs))

    push!(mean_dir, mean(_wave_dir))
    push!(min_dir, minimum(_wave_dir))
    push!(max_dir, maximum(_wave_dir))
    push!(std_dir, std(_wave_dir))

    push!(changes, change)
end

train_df = DataFrame(
    Tp_mean=convert(Vector{Float64}, mean_tp[1:800]),
    Tp_min=convert(Vector{Float64}, min_tp[1:800]),
    Tp_max=convert(Vector{Float64}, max_tp[1:800]),
    Tp_std=convert(Vector{Float64}, std_tp[1:800]),

    Hs_mean=convert(Vector{Float64}, mean_hs[1:800]),
    Hs_min=convert(Vector{Float64}, min_hs[1:800]),
    Hs_max=convert(Vector{Float64}, max_hs[1:800]),
    Hs_std=convert(Vector{Float64}, std_hs[1:800]),

    Dir_mean=convert(Vector{Float64}, mean_dir[1:800]),
    Dir_min=convert(Vector{Float64}, min_dir[1:800]),
    Dir_max=convert(Vector{Float64}, max_dir[1:800]),
    Dir_std=convert(Vector{Float64}, std_dir[1:800]),

    Change=convert(Vector{Float64}, changes[1:800])
)

test_df = DataFrame(
    Tp_mean=convert(Vector{Float64}, mean_tp[800:end]),
    Tp_min=convert(Vector{Float64}, min_tp[800:end]),
    Tp_max=convert(Vector{Float64}, max_tp[800:end]),
    Tp_std=convert(Vector{Float64}, std_tp[800:end]),

    Hs_mean=convert(Vector{Float64}, mean_hs[800:end]),
    Hs_min=convert(Vector{Float64}, min_hs[800:end]),
    Hs_max=convert(Vector{Float64}, max_hs[800:end]),
    Hs_std=convert(Vector{Float64}, std_hs[800:end]),

    Dir_mean=convert(Vector{Float64}, mean_dir[800:end]),
    Dir_min=convert(Vector{Float64}, min_dir[800:end]),
    Dir_max=convert(Vector{Float64}, max_dir[800:end]),
    Dir_std=convert(Vector{Float64}, std_dir[800:end]),

    Change=convert(Vector{Float64}, changes[800:end])
)

println(describe(train_df))

model = lm(@formula(Change ~ Tp_mean + Tp_min + Tp_max + Tp_std + Hs_mean + Hs_min + Hs_max + Hs_std + Dir_mean + Dir_min + Dir_max + Dir_std), train_df)
println(model)
println("R2 score: ", r2(model))

ypredicted_test = predict(model, test_df)
ypredicted_train = predict(model, train_df)

rmse(target, preds) = mean((target .- preds).^2)

println("Train RMSE: ", rmse(changes[1:800], ypredicted_train))
println("Test RMSE: ", rmse(changes[800:end], ypredicted_test))

corr_train = round(cor(changes[1:800], ypredicted_train), digits=2)
corr_test = round(cor(changes[800:end], ypredicted_test), digits=2)
println("Train corr: ", corr_train)
println("Test corr: ", corr_test)

change_xs = collect(1:length(changes))
train_xs = collect(1:800)
test_xs = collect(800:length(changes))

p1 = plot(change_xs, changes, label="Raw data")
plot!(train_xs, ypredicted_train, label="Train set preds")
plot!(test_xs, ypredicted_test, label="Test set preds")
# display(plot(p1))
# readline()

plot!(size=(1920, 1080))
savefig("/Users/mn/JuliaProjects/cgp-morphcast/output/grand_popo/linear_regression.pdf")

max_stars_repo_name: mahmoud-al-najar/cgp-morphcast

max_issues_repo_name: mahmoud-al-najar/cgp-morphcast

max_forks_repo_name: mahmoud-al-najar/cgp-morphcast


text:
addprocs(5)

include("../../daggen.jl")
using DagScheduler
using BenchmarkTools

isdir(".mempool") && rm(".mempool"; recursive=true)
@everywhere begin
    function resetscheduler()
        Dagger.PLUGIN_CONFIGS[:scheduler] = "Dagger.Sch"
    end
    Dagger.use_shared_array[] = true
    resetscheduler()
end

const L = 6^4
const dag2 = gen_straight_dag(ones(Int, L));

function clean_compute(dag2)
    result = compute(dag2)
    gc() # gc required to prevent memory usage rising exponentially
    @everywhere gc()
    collect(result) # to make it similar to what dagscheduer does today
end

result = clean_compute(dag2);
#@time result = clean_compute(dag2);
@btime clean_compute(dag2);


text:
jpeg_components(::AbstractArray{T}) where T = jpeg_components(T)
jpeg_components(::Type{CT}) where CT<:Colorant = length(CT)

jpeg_color_space(::AbstractArray{T}) where T = jpeg_color_space(T)
jpeg_color_space(::Type{CT}) where CT<:Gray = LibJpeg.JCS_GRAYSCALE
jpeg_color_space(::Type{CT}) where CT<:RGB = LibJpeg.JCS_RGB
jpeg_color_space(::Type{CT}) where CT<:BGR = LibJpeg.JCS_EXT_BGR
jpeg_color_space(::Type{CT}) where CT<:RGBA = LibJpeg.JCS_EXT_RGBA
jpeg_color_space(::Type{CT}) where CT<:BGRA = LibJpeg.JCS_EXT_BGRA
jpeg_color_space(::Type{CT}) where CT<:ABGR = LibJpeg.JCS_EXT_ABGR
jpeg_color_space(::Type{CT}) where CT<:ARGB = LibJpeg.JCS_EXT_ARGB

function jpeg_color_space(v::LibJpeg.J_COLOR_SPACE)::Type{<:Colorant}
    if v === LibJpeg.JCS_GRAYSCALE
        return Gray{N0f8}
    elseif v === LibJpeg.JCS_RGB
        return RGB{N0f8}
    elseif v === LibJpeg.JCS_EXT_RGB
        return RGB{N0f8}
    elseif v === LibJpeg.JCS_EXT_BGR
        return BGR{N0f8}
    elseif v === LibJpeg.JCS_EXT_RGBA
        return RGBA{N0f8}
    elseif v === LibJpeg.JCS_EXT_BGRA
        return BGRA{N0f8}
    elseif v === LibJpeg.JCS_EXT_ABGR
        return ABGR{N0f8}
    elseif v === LibJpeg.JCS_EXT_ARGB
        return ARGB{N0f8}
    else
        throw(ArgumentError("Unsupported colorspace $v"))
    end
end


text:
module SCOPF
    using Base: String
    using LinearAlgebra:length
    using SparseArrays, LinearAlgebra
    using ..AmplTxt
    using Printf

    mutable struct Bus
        id::Int
        name::String
    end

    mutable struct Branch
        from::Int
        to::Int
        r::Float64
        x::Float64
        name::String
    end
    function get_b(branch::Branch)
        # return branch.x / sqrt(branch.r * branch.r + branch.x * branch.x)
        return 1/branch.x;
    end
    mutable struct Network
        bus_to_i::Dict{Int,Int}
        branch_to_i::Dict{Int,Int}

        branches::Dict{Int,Branch}
        buses::Dict{Int,Bus}
        CC0::Dict{Int,Int}
    end
    function Network()
        return Network(Dict{Int,Int}(), Dict{Int,Int}(), Dict{Int,Branch}(), Dict{Int,Bus}(), Dict{Int,Int}())
    end

    function Network(amplTxt)
        network = Network();
        add_bus!(network, amplTxt);
        add_branches!(network, amplTxt);
        add_generator!(network, amplTxt);
        return network
    end

    function get_or_add(i::Int, d::Dict{Int,Int})
        if ! haskey(d, i)
            next_id = length(d) + 1;
            d[i] = next_id;
        end
        return d[i];
    end
    function add_branches!(network::Network, amplTxt)
        branches = amplTxt["branches"];
        for branch in branches.data
            # println(branch)
            num  = parse(Int, branch[2]);
            num_or = parse(Int, branch[3]);
            num_ex = parse(Int, branch[4]);

            on_CC0 = haskey(network.CC0, num_or) && haskey(network.CC0, num_ex);

            if on_CC0 && num_or != -1  && num_ex != -1
                id = get_or_add(num, network.branch_to_i);
                ior = get_or_add(num_or, network.bus_to_i);
                iex = get_or_add(num_ex, network.bus_to_i);
                r = parse(Float64, branch[8]);
                x = parse(Float64, branch[9]);
                push!(network.branches, id => Branch(ior, iex, r, x, branch[26]));
            end
        end
    end
    function add_bus!(network::Network, amplTxt)
        buses = amplTxt["buses"];
        for bus in buses.data
            num = parse(Int, bus[2]);
            numCC = parse(Int, bus[4]);
            if numCC == 0
                id = get_or_add(num, network.bus_to_i);
                push!(network.buses, id => Bus(num, bus[11]));
                network.CC0[num] = 1;
            end
        end
    end

    function add_generator!(network::Network, amplTxt)
        generators = amplTxt["generators"];
        for genData in generators.data
            gen = parse(Int, genData[2]);
            bus = parse(Int, genData[3]);
            conbus = parse(Int, genData[4]);
            minP = parse(Float64, genData[6]);
            maxP = parse(Float64, genData[7]);
        end
    end
    function get_B(network::Network, DIAG_EPS::Float64)
        B_dict = Dict{Tuple{Int,Int},Float64}();
        for kvp in collect(network.branches)
            b = get_b(kvp[2]);

            ior = kvp[2].from;
            iex = kvp[2].to;

            # ior->iex
            # + dans ior
            # - dans iex

            # diagonal part
            # sur ior
            # +theta_ior
            # -theta_iex
            B_dict[ior, ior] = get(B_dict, (ior, ior), DIAG_EPS) + b;
            B_dict[ior, iex] = get(B_dict, (ior, iex), 0) - b;

            # sur iex
            # -theta_ior
            # +theta_iex
            B_dict[iex, ior] = get(B_dict, (iex, ior), 0) - b;
            B_dict[iex, iex] = get(B_dict, (iex, iex), DIAG_EPS) + b;

            # B_dict[ior, ior] = get(B_dict, (ior, ior), DIAG_EPS) + b;
            # B_dict[iex, iex] = get(B_dict, (iex, iex), DIAG_EPS) - b;

            # B_dict[ior, iex] = get(B_dict, (ior, iex), 0) + b;
            # B_dict[iex, ior] = get(B_dict, (iex, ior), 0) - b;
        end
        Brow = Int[];
        Bcol = Int[];
        Bval = Float64[];
        for kvp in collect(B_dict)
            push!(Brow, kvp[1][1]);
            push!(Bcol, kvp[1][2]);
            push!(Bval, kvp[2]);
        end
        n = length(network.bus_to_i);
        println("n = ", n)
        B = sparse(Brow, Bcol, Bval, n, n);
        Bdense = Matrix(B);
        return Bdense
    end

    function get_B_inv(Bdense::Matrix, ref_bus::Int64)
        println("ref_bus is ", ref_bus)
        n = size(Bdense)[1]
        Bdense2 = Bdense
        Bdense2[ref_bus, 1:n] .= 0
        Bdense2[1:n, ref_bus] .= 0
        Bdense2[ref_bus, ref_bus] = 1
        binv = inv(Bdense2);
        binv[ref_bus,ref_bus] = 0;
        return binv
    end

    function get_PTDF(network::Network, binv::Matrix, ref_bus::Int)
        n = length(network.bus_to_i);
        m = length(network.branches);
        println("n ", n)
        println("m ", m)
        PTDF = zeros(Float64, m, n);
        for kvp in collect(network.branches)
            b = get_b(kvp[2]);

            branchid = kvp[1];
            ior = kvp[2].from;
            iex = kvp[2].to;

            for i in 1:n
                PTDF[branchid, i] += b * binv[ior, i];
                PTDF[branchid, i] -= b * binv[iex, i];
                # if ior != ref_bus
                #     PTDF[branchid, i] += b * binv[ior, i];
                # end
                # if iex != ref_bus
                #     PTDF[branchid, i] -= b * binv[iex, i];
                # end
            end
        end
        return PTDF
    end

    function distribute_slack(PTDF::Matrix, ref_bus::Int)
        m, n = size(PTDF);
        μ = 1/n;
        v = PTDF * ones(n)*μ;
        result = -v*(ones(n)');
        result +=PTDF ;
        return result;
    end

    function write_PTDF(file_path::String, network::Network, ref_bus::Int, PTDF::Matrix, PTDF_TRIMMER::Float64)
        n = length(network.bus_to_i);
        m = length(network.branches);
        open(file_path, "w") do file
            ref_name =  @sprintf("\"%s\"", network.buses[ref_bus].name)
            write(file, @sprintf("#%20s %20s\n", "REF_BUS", ref_name))
            for branch_id in 1:m
                for bus_id in 1:n
                    branch_name =  @sprintf("\"%s\"", network.branches[branch_id].name)
                    bus_name =  @sprintf("\"%s\"", network.buses[bus_id].name)
                    if abs(PTDF[branch_id,bus_id])>PTDF_TRIMMER
                        write(file, @sprintf("%20s %20s %20.6E\n", branch_name, bus_name,PTDF[branch_id,bus_id]))
                    end
                end
            end
        end
    end
    export get_b;
    export add_bus!;
    export add_branches!;
    export get_B;
    export get_B_inv;
    export get_PTDF;
    export write_PTDF;
    export distribute_slack;
end

max_stars_repo_name: azized/pscopf

max_issues_repo_name: azized/pscopf

max_forks_repo_name: azized/pscopf


text:
using ITensors,
      ITensorsGPU,
      Test

@testset "CuMPO Basics" begin
  N = 6
  sites = [Index(2,"Site") for n=1:N]
  @test length(cuMPO()) == 0
  O = cuMPO(sites)
  @test length(O) == N

  str = split(sprint(show, O), '\n')
  @test str[1] == "MPO"
  @test length(str) == length(O) + 2

  O[1] = cuITensor(sites[1], prime(sites[1]))
  @test hasind(O[1],sites[1])
  @test hasind(O[1],prime(sites[1]))
  P = copy(O)
  @test hasind(P[1],sites[1])
  @test hasind(P[1],prime(sites[1]))

  K = randomCuMPO(sites)
  K_ = cuMPO(ITensors.data(K))
  @test all(ITensors.data(K) .== ITensors.data(K_))

  L = randomMPO(sites)
  K = cuMPO(L)
  @test all(ITensors.data(collect(K)) .== ITensors.data(collect(L)))
  @testset "orthogonalize" begin
    phi = randomCuMPS(sites)
    K = randomCuMPO(sites)
    orthogonalize!(phi, 1)
    orthogonalize!(K, 1)
    orig_inner = inner(phi, K, phi)
    orthogonalize!(phi, div(N, 2))
    orthogonalize!(K, div(N, 2))
    @test inner(phi, K, phi) ≈ orig_inner
  end

  @testset "inner <y|A|x>" begin
    phi = randomCuMPS(sites)
    K = randomCuMPO(sites)
    @test maxlinkdim(K) == 1
    psi = randomCuMPS(sites)
    phidag = dag(phi)
    prime!(phidag)
    phiKpsi = phidag[1]*K[1]*psi[1]
    for j = 2:N
      phiKpsi *= phidag[j]*K[j]*psi[j]
    end
    @test phiKpsi[] ≈ inner(phi,K,psi)

    badsites = [Index(2,"Site") for n=1:N+1]
    badpsi = randomCuMPS(badsites)
    @test_throws DimensionMismatch inner(phi,K,badpsi)

    # make bigger random MPO...
    for link_dim in 2:5
        mpo_tensors  = ITensor[cuITensor() for ii in 1:N]
        mps_tensors  = ITensor[cuITensor() for ii in 1:N]
        mps_tensors2 = ITensor[cuITensor() for ii in 1:N]
        mpo_link_inds = [Index(link_dim, "r$ii,Link") for ii in 1:N-1]
        mps_link_inds = [Index(link_dim, "r$ii,Link") for ii in 1:N-1]
        mpo_tensors[1] = randomCuITensor(mpo_link_inds[1], sites[1], sites[1]')
        mps_tensors[1] = randomCuITensor(mps_link_inds[1], sites[1])
        mps_tensors2[1] = randomCuITensor(mps_link_inds[1], sites[1])
        for ii in 2:N-1
            mpo_tensors[ii] = randomCuITensor(mpo_link_inds[ii], mpo_link_inds[ii-1], sites[ii], sites[ii]')
            mps_tensors[ii] = randomCuITensor(mps_link_inds[ii], mps_link_inds[ii-1], sites[ii])
            mps_tensors2[ii] = randomCuITensor(mps_link_inds[ii], mps_link_inds[ii-1], sites[ii])
        end
        mpo_tensors[N] = randomCuITensor(mpo_link_inds[N-1], sites[N], sites[N]')
        mps_tensors[N] = randomCuITensor(mps_link_inds[N-1], sites[N])
        mps_tensors2[N] = randomCuITensor(mps_link_inds[N-1], sites[N])
        K   = MPO(mpo_tensors, 0, N+1)
        psi = MPS(mps_tensors, 0, N+1)
        phi = MPS(mps_tensors2, 0, N+1)
        orthogonalize!(psi, 1; maxdim=link_dim)
        orthogonalize!(K,   1; maxdim=link_dim)
        orthogonalize!(phi, 1; normalize=true, maxdim=link_dim)
        phidag = dag(phi)
        prime!(phidag)
        phiKpsi = phidag[1]*K[1]*psi[1]
        for j = 2:N
          phiKpsi *= phidag[j]*K[j]*psi[j]
        end
        @test scalar(phiKpsi) ≈ inner(phi,K,psi)
    end
  end

  @testset "contract" begin
    phi = randomCuMPS(sites)
    K   = randomCuMPO(sites)
    @test maxlinkdim(K) == 1
    psi = randomCuMPS(sites)
    psi_out = contract(K, psi,maxdim=1)
    @test inner(phi,psi_out) ≈ inner(phi,K,psi)
    @test_throws ArgumentError contract(K, psi, method="fakemethod")

    badsites = [Index(2,"Site") for n=1:N+1]
    badpsi = randomCuMPS(badsites)
    @test_throws DimensionMismatch contract(K,badpsi)

    # make bigger random MPO...
    for link_dim in 2:5
        mpo_tensors  = ITensor[ITensor() for ii in 1:N]
        mps_tensors  = ITensor[ITensor() for ii in 1:N]
        mps_tensors2 = ITensor[ITensor() for ii in 1:N]
        mpo_link_inds = [Index(link_dim, "r$ii,Link") for ii in 1:N-1]
        mps_link_inds = [Index(link_dim, "r$ii,Link") for ii in 1:N-1]
        mpo_tensors[1] = randomCuITensor(mpo_link_inds[1], sites[1], sites[1]')
        mps_tensors[1] = randomCuITensor(mps_link_inds[1], sites[1])
        mps_tensors2[1] = randomCuITensor(mps_link_inds[1], sites[1])
        for ii in 2:N-1
            mpo_tensors[ii] = randomCuITensor(mpo_link_inds[ii], mpo_link_inds[ii-1], sites[ii], sites[ii]')
            mps_tensors[ii] = randomCuITensor(mps_link_inds[ii], mps_link_inds[ii-1], sites[ii])
            mps_tensors2[ii] = randomCuITensor(mps_link_inds[ii], mps_link_inds[ii-1], sites[ii])
        end
        mpo_tensors[N] = randomCuITensor(mpo_link_inds[N-1], sites[N], sites[N]')
        mps_tensors[N] = randomCuITensor(mps_link_inds[N-1], sites[N])
        mps_tensors2[N] = randomCuITensor(mps_link_inds[N-1], sites[N])
        K   = MPO(mpo_tensors, 0, N+1)
        psi = MPS(mps_tensors, 0, N+1)
        phi = MPS(mps_tensors2, 0, N+1)
        orthogonalize!(psi, 1; maxdim=link_dim)
        orthogonalize!(K,   1; maxdim=link_dim)
        orthogonalize!(phi, 1; normalize=true, maxdim=link_dim)
        psi_out = contract(deepcopy(K), deepcopy(psi); maxdim=10*link_dim, cutoff=0.0)
        @test inner(phi, psi_out) ≈ inner(phi, K, psi)
    end
  end
  @testset "add" begin
    shsites = siteinds("S=1/2",N)
    K = randomCuMPO(shsites)
    L = randomCuMPO(shsites)
    M = add(K, L)
    @test length(M) == N
    psi = randomCuMPS(shsites)
    k_psi = contract(K, psi, maxdim=1)
    l_psi = contract(L, psi, maxdim=1)
    @test inner(psi, add(k_psi, l_psi)) ≈ inner(psi, M, psi) atol=5e-3
  end
  @testset "contract(::CuMPO, ::CuMPO)" begin
    psi = randomCuMPS(sites)
    K = randomCuMPO(sites)
    L = randomCuMPO(sites)
    @test maxlinkdim(K) == 1
    @test maxlinkdim(L) == 1
    KL = contract(prime(K), L, maxdim=1)
    psi_kl_out = contract(prime(K), contract(L, psi, maxdim=1), maxdim=1)
    @test inner(psi,KL,psi) ≈ inner(psi, psi_kl_out) atol=5e-3

    # where both K and L have differently labelled sites
    othersitesk = [Index(2,"Site,aaa") for n=1:N]
    othersitesl = [Index(2,"Site,bbb") for n=1:N]
    K = randomCuMPO(sites)
    L = randomCuMPO(sites)
    for ii in 1:N
     replaceind!(K[ii], sites[ii]', othersitesk[ii])
     replaceind!(L[ii], sites[ii]', othersitesl[ii])
    end
    KL = contract(K, L, maxdim=1)
    psik = randomCuMPS(othersitesk)
    psil = randomCuMPS(othersitesl)
    psi_kl_out = contract(K, contract(L, psil, maxdim=1), maxdim=1)
    @test inner(psik,KL,psil) ≈ inner(psik, psi_kl_out) atol=5e-3

    badsites = [Index(2,"Site") for n=1:N+1]
    badL = randomCuMPO(badsites)
    @test_throws DimensionMismatch contract(K,badL)
  end
end


text:
include("init.jl")
include("lib/MatlabDecomp.jl")
mld = MatlabDecomp
include("lib/ProcessLibrariesIO.jl")

using Serialization

import ProgressMeter
import Dates

getr(rs, platform, repeats=10) = begin
    @info "Deserialising adjacency matrices"
    adj_mats = ProcessLibrariesIO.get_adj_mats(platform)
    @info "Loading adjacency matrices into matlab"
    mld.matlab_load_sptensor(adj_mats)
    adj_mats = nothing

    ProgressMeter.@showprogress for r in rs
        @info "Computing $repeats decompositions with r = $r"
        Ds = [mld.matlab_nncp_loaded_spt(r) for i in 1:repeats]
        dir = "../data/processed/pacman/$platform/decomps/$r/"
        run(`mkdir -p $dir`)
        serialize("$dir/$(Dates.now()).jls", Ds)
        @info "Component $r is done $(Dates.now())"
    end
    @info "All $rs × $repeats for $platform done!"
end

repeats = 10

if length(ARGS) > 3
    rstart = parse(Int, ARGS[1])
    rstep = parse(Int, ARGS[2])
    rend = parse(Int, ARGS[3])
    rs = range(rstart, rend, step=rstep)
    if length(ARGS) > 4
        repeats = parse(Int, ARGS[4])
        platform = ARGS[5]
    else
        platform = ARGS[4]
    end
else
    rstart = parse(Int, ARGS[1])
    rend = parse(Int, ARGS[2])
    rs = range(rstart, stop=rend)
    platform = ARGS[3]
end

getr(rs, platform, repeats)


text:
# Compare the times of simulation
using NbodyGradient, DelimitedFiles, BenchmarkTools

elements = readdlm("elements_rebound.txt", ',')
elements[:,end] .= 0.0

# Remove the old times file, if exists
try
  rm("nbg_times.txt")
catch
end

# Set up integrator
t0 = 0.0
h = 0.05
# Rebound integration is for 5000 code units; inner planet has
# a period of 2π, while in our case 1-day, so differs by 5000/2π
tmax = 2500/π

intr = Integrator(h, 0.0, tmax)

# For integrating with only Kepler+drift for adjacent planets & planets+star
# (more distant planets are handled with fast kicks):
function integrate_adj(s::State{T},h::T,tmax::T;grad::Bool=false) where T<:Real
  # Only include Kepler+drift solver for adjacent planets:
  s.pair .= true
  s.pair[1,:] .= false
  s.pair[:,1] .= false
  for i=2:s.n-1
    s.pair[i,i+1] = false
    s.pair[i+1,i] = false
  end
  if grad; d = NbodyGradient.Derivatives(T,s.n); end
  N = abs(round(Int64,tmax/h))
  for _ in 1:N
    if grad
      ahl21!(s,d,h)
    else
      ahl21!(s,h)
    end
  end
  s.t[1] = h*N
  if grad; return d; else return; end
end

for n = 2:11
  ic = ElementsIC(t0, n, elements[1:n,:])
  s = State(ic)
  tt = TransitTiming(intr.tmax, ic)

  # Time without derivatives
  t = @belapsed intr(s, grad=false) setup=(s = State($ic))

  # Time with derivatives
  grad_t = @belapsed intr(s) setup=(s = State($ic))

  # Time with non-adjacent fast-kicks:
  t_adj = @belapsed integrate_adj(s,h,tmax,grad=false) setup=(s=State($ic))
  grad_adj = @belapsed integrate_adj(s,h,tmax,grad=true) setup=(s=State($ic))

  # Now while calculating transit times and derivatives
  grad_tt = @belapsed intr(s, tt) setup=(s = State($ic); tt = TransitTiming($intr.tmax, $ic))

  println("Nplanet: ", n-1, " No gradient: ", t, " gradient: ", grad_t, " t adj: ",t_adj," grad t_adj: ",grad_adj," grad w tts: ", grad_tt," ratio: ", grad_t/t)
  open("nbg_times.txt", "a") do io
    writedlm(io, [n-1 t grad_t t_adj grad_adj grad_tt grad_t/t], ',')
  end
end

max_forks_repo_name: langfzac/nbg-papers


text:
using Test
using Distributed, .Threads

if nprocs() == 1 && nthreads() == 1
    # Lazy stuff should also test distributed
    println("Adding a proc")
    addprocs(1; exeflags="--project")
end

@testset "basics" begin include("basics.jl"); end
@testset "taxi" begin include("taxi.jl"); end

if nprocs() > 1
    rmprocs(workers())
end


text:
import Base.identity
import Base.*
import Base.inv
import Base.show
import Base.one

"""
Artin's braid group of order `n`.
"""
struct BraidGroup
    n::Int
end

"""
An element of Artin's braid group.
"""
struct BraidElement
    group::BraidGroup
    crossings::Array{Tuple{Int,Int},1}
end

"""
Alias for the identity on `g`.
"""
function identity(g::BraidGroup)
    one(g)
end

"""
The identity on `g`.
"""
function one(g::BraidGroup)
    BraidElement(g, [])
end

"""
A generator `sigma_i^sgn` of the group `g`.
"""
function generator(g::BraidGroup, i::Int, sgn::Bool)
    sgnn = sgn ? 1 : -1
    BraidElement(g, [(i % g.n, sgnn)])
end

"""
Get braid from integer encoding of generators.
"""
function braidfromintrep(g::BraidGroup, rep::Array{Int,1})
    crossings = [(abs(x), x > 0 ? 1 : -1) for x in rep]
    BraidElement(g, crossings)
end

"""
Free reduction on the braid word crossings.
"""
function simplifybraidcrossings(crossings::Array{Tuple{Int,Int},1})
    xings0 = crossings
    xings1 = []
    nchanged = 1
    # Simplify exprssion to extent possible
    while nchanged > 0
        lchanged = 0
        if length(xings0) > 1
            for i = 1:length(xings0)
                if (i < length(xings0)) && (xings0[i][1] == xings0[i+1][1]) && (xings0[i][2] == -xings0[i+1][2])
                    # In theory these should always cancel but let's check
                    newpow = xings0[i][2] + xings0[i+1][2]
                    if newpow != 0
                        push!(xings1, (xings0[i][1], newpow))
                    end
                    lchanged += 1
                elseif (i > 1) && ((xings0[i][1] != (xings0[i-1][1])) || (xings0[i][2] != -xings0[i-1][2]))
                    push!(xings1, xings0[i])
                elseif (i == 1) && xings0[i][1] != xings0[i+1][1]
                    push!(xings1, xings0[i])
                elseif (i == 1) && xings0[i][1] == xings0[i+1][1]
                    push!(xings1, xings0[i])
                end
            end
            nchanged = lchanged
            xings0 = xings1
            xings1 = []
        else
            nchanged = 0
        end
    end
    xings = filter(x -> x[2] != 0, xings0) # Remove identities and ^0s
    xings
end

"""
Free reduction on the braid word.
"""
function simplifybraid(b::BraidElement)
    BraidElement(b.group, simplifybraidcrossings(b.crossings))
end

"""
Simplify a braid word.
"""
function simplify(b::BraidElement)
    prettybraid(BraidElement(b.group, simplifybraidcrossings(b.crossings)))
end

"""
Multiplication of braid words.
"""
function *(b1::BraidElement, b2::BraidElement)
    xings0 = cat(1, b1.crossings, b2.crossings)
    xings = simplifybraidcrossings(xings0)
    prettybraid(BraidElement(b1.group, xings))
end

"""
Finds crossing of lowest index in `b`.
"""
function lowestindex(b::BraidElement)
    if length(b.crossings) == 0
        0
    else
        minimum([x[1] for x in b.crossings])
    end
end

"""
Determines whether `b` is reduced in the sense of Dehornoy.
"""
function isreduced(b::BraidElement)
    minidx = lowestindex(b)
    genpows = [q[2] for q in find(z -> z[1] == minidx, b.crossings)]
    (minimum(genpows) < 0) == (maximum(genpows) < 0)
end

"""
Get the inverse of `b`.
"""
function inv(b::BraidElement)
    newxings = [(q[1], -1 * q[2]) for q in b.crossings[length(b.crossings):-1:1]]
    BraidElement(b.group, newxings)
end

"""
Convert `b` to a string.
"""
function Base.show(io::IO, b0::BraidElement)
    b = prettybraid(b0)
    getexp(i) = (i == 1 ? "" : (i == -1 ? "′" : string("^", i)))
    if length(b.crossings) > 0
        sx = [string("σ_", string(i[1]), getexp(i[2])) for i in b.crossings]
        s = join(sx, "⋅")
        #s = s[length("String")+1:end] # Hackety hack
        print(io, s)
    else
        print(io, "ϵ")
    end
end

"""
Convert representation from pure generators to powers of generators.
"""
function prettybraid(b::BraidElement)
    if length(b.crossings) == 0
        return b
    end
    c = []
    curexp = [b.crossings[1][2]]
    curbase = b.crossings[1][1]
    for i in 2:length(b.crossings)
        x = b.crossings[i]
        if x[1] == curbase
            push!(curexp, x[2])
        else
            push!(c, (curbase, sum(curexp)))
            curbase = x[1]
            curexp = [x[2]]
        end
    end
    push!(c, (curbase, sum(curexp)))
    simplifybraid(BraidElement(b.group, c))
end

"""
Convert representation from powers of generators to pure generators.
"""
function unprettybraid(b::BraidElement)
    c = []
    for x in b.crossings
        sgn = x[2] < 0 ? -1 : 1
        for i = 1:abs(x[2])
            push!(c, (x[1], sgn))
        end
    end
    simplifybraid(BraidElement(b.group, c))
end


text:
module PackedVectorsOfVectors

export
    PackedVectorOfVectors,
    pack,
    allocate_packed,
    packed_zeros,
    packed_ones,
    packed_fill,
    packed_indices

"""
    PackedVectorOfVectors{P,V,E} <: AbstractVector{E}

Vector of vectors, stored as a single long vector `v::V` and a vector of
pointers `p::P` indicating the subvector boundaries.

`PackedVectorOfVectors` is roughly equivalent to `SparseMatrixCSC` without the
`rowval` vector.
"""
struct PackedVectorOfVectors{P,V,E} <: AbstractVector{E}
    p::P
    v::V
end

PackedVectorOfVectors(p,v) = PackedVectorOfVectors{typeof(p),typeof(v),typeof(view(v,1:0))}(p,v)

Base.size(pv::PackedVectorOfVectors) = (length(pv.p)-1,)
Base.getindex(pv::PackedVectorOfVectors, i::Int) = view(pv.v, pv.p[i]:pv.p[i+1]-1)
Base.showarg(io::IO, pv::PackedVectorOfVectors, toplevel) = print(io, "pack(::Vector{Vector{", eltype(eltype(pv)), "}})")
Iterators.flatten(pv::PackedVectorOfVectors) = pv.v

"""
    cumsum1(v) -> p

Compute the cumulative sum of the iterable `v`, starting from `p[1] = 1`.

# Examples
```jldoctest
julia> cumsum1([2,1,3])
4-element Vector{Int64}:
 1
 3
 4
 7
```
"""
cumsum1(v) = cumsum1!(Vector{Int}(undef,length(v)+1), v)

"""
    cumsum1!(p,v) -> p

Same as `cumsum1(v)`, but using the preallocated vector `p`.
"""
function cumsum1!(p,v)
    @assert length(p) == length(v)+1  "Mismatched lengths of p and v in cumsum1!(p,v)"
    s = 1
    p[1] = s
    @inbounds for (k,vk) in enumerate(v)
        s += vk
        p[k+1] = s
    end
    return p
end

"""
    pack([T,] vv) -> PackedVectorOfVectors

Convert the vector of vectors `vv` to its packed representation. If `T` is
provided, then the elements of the nested vector are converted to this type.

# Example
```jldoctest
julia> pack([[1,2],[3]])
2-element pack(::Vector{Vector{Int64}}):
 [1, 2]
 [3]

julia> pack(Float64, [[1,2],[3]])
2-element pack(::Vector{Vector{Float64}}):
 [1.0, 2.0]
 [3.0]
```
"""
pack(vv) = pack(eltype(eltype(vv)),vv)
function pack(T::Type, vv)
    p = cumsum1(length(v) for v in vv)
    v = Vector{T}(undef, p[end]-1)
    pv = PackedVectorOfVectors(p,v)
    for (k,vk) in enumerate(vv)
        pv[k] .= vk
    end
    return pv
end

"""
    allocate_packed(T, init, n) -> PackedVectorOfVectors

Allocate a packed vector of vectors with nested element type `T` such that the
`k`th nested vector has length `nth(n,k)`. See the documentation of
`Vector{T}(init,n)` regarding the meaning of `init`.

# Examples
```jldoctest
julia> allocate_packed(Ref{Int}, undef, [1,2,3])
3-element pack(::Vector{Vector{Ref{Int64}}}):
 [#undef]
 [#undef, #undef]
 [#undef, #undef, #undef]
```
"""
function allocate_packed(T::Type, init, n)
    p = cumsum1(n)
    v = Vector{T}(init, p[end]-1)
    return PackedVectorOfVectors(p,v)
end

for f in (:zeros, :ones)
    pf = Symbol(:packed_,f)
    @eval begin
        """
            $($pf)([T=Float64,] n) -> pv::PackedVectorOfVectors

        Create a packed vector of vectors such that `pv[k] == $($f)(T, nth(n,k))`.
        """
        function $pf(T::Type, n)
            p = cumsum1(n)
            v = $f(T, p[end]-1)
            return PackedVectorOfVectors(p,v)
        end
        $pf(n) = $pf(Float64,n)
    end
end

"""
    packed_fill(x,n) -> pv::PackedVectorOfVectors

Create a packed vector of vectors such that `pv[k] == fill(x, nth(n,k))`.
"""
function packed_fill(x, n)
    p = cumsum1(n)
    v = fill(x, p[end]-1)
    return PackedVectorOfVectors(p,v)
end

"""
    packed_indices(n) -> pv::PackedVectorOfVectors

Packed vector of vectors such that the `k`th nested vector has length `nth(n,k)`
and `pv[k][i]` is the index in the flattened vector of the `i`th element in the
`k`th nested vector.

# Examples
```jldoctest
julia> packed_indices([2,3])
2-element pack(::Vector{Vector{Int64}}):
 1:2
 3:5
```
"""
function packed_indices(n)
    p = cumsum1(n)
    v = 1:p[end]-1
    return PackedVectorOfVectors(p,v)
end

"""
    transpose(pv::PackedVectorOfVectors[, n = maximum(maximum.(pv))) -> PackedVectorOfVectors

Interpret `pv` as the sparsity pattern of a sparse matrix and compute the
sparsity pattern of its transpose.

The nested eltype of `pv` must be `Int`.

# Examples
```jldoctest
julia> transpose(pack([[1,2],[2]]))
2-element pack(::Vector{Vector{Int64}}):
 [1]
 [1, 2]
```
"""
function Base.transpose(
    pv::PackedVectorOfVectors,
    n = maximum(pv.v; init=0)
)
    @assert eltype(eltype(pv)) == Int  "transpose(pv::PackedVectorOfVectors) is defined only if the nested eltype is Int"
    p,v = pv.p,pv.v

    # Compute row counts
    q = zeros(Int,n+1)
    for vk in v
        q[vk+1] += 1
    end

    # Set up q so q[k+1] points to beginning of column k in w
    q[1] = 1
    s = 1
    @inbounds for k = 2:n+1
        t = q[k]
        q[k] = s
        s += t
    end

    # Copy v into w and simultaneously finalise q
    w = Vector{Int}(undef, length(v))
    @inbounds for i = 1:length(pv)
        for j = pv[i]
            w[q[j+1]] = i
            q[j+1] += 1
        end
    end

    return PackedVectorOfVectors(q,w)
end

end # module


text:
export Frobenius
"""
[`Frobenius` matrix](http://en.wikipedia.org/wiki/Frobenius_matrix)

Frobenius matrices or Gaussian elimination matrices of the form

```
[ 1 0 ...     0 ]
[ 0 1 ...     0 ]
[ .........     ]
[ ... 1 ...     ]
[ ... c1 1 ...  ]
[ ... c2 0 1 ...]
[ ............. ]
[ ... ck ...   1]
```

i.e. an identity matrix with nonzero subdiagonal elements along a single column.

```
 In this implementation, the subdiagonal of the nonzero column is stored as a
 dense vector, so that the size can be inferred automatically as j+k where j is
 the index of the column and k is the number of subdiagonal elements.
```

```julia
julia> using SpecialMatrices

julia> F=Frobenius(3, [1.0,2.0,3.0]) #Specify subdiagonals of column 3
6x6 Frobenius{Float64}:
 1.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  0.0  0.0
 0.0  0.0  2.0  0.0  1.0  0.0
 0.0  0.0  3.0  0.0  0.0  1.0

julia> inv(F) #Special form of inverse
6x6 Frobenius{Float64}:
 1.0  0.0   0.0  0.0  0.0  0.0
 0.0  1.0   0.0  0.0  0.0  0.0
 0.0  0.0   1.0  0.0  0.0  0.0
 0.0  0.0  -1.0  1.0  0.0  0.0
 0.0  0.0  -2.0  0.0  1.0  0.0
 0.0  0.0  -3.0  0.0  0.0  1.0

julia> F*F #Special form preserved if the same column has the subdiagonals
6x6 Frobenius{Float64}:
 1.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0
 0.0  0.0  2.0  1.0  0.0  0.0
 0.0  0.0  4.0  0.0  1.0  0.0
 0.0  0.0  6.0  0.0  0.0  1.0

julia> F*Frobenius(2, [5.0,4.0,3.0,2.0]) #Promotes to Matrix
6x6 Array{Float64,2}:
 1.0   0.0  0.0  0.0  0.0  0.0
 0.0   1.0  0.0  0.0  0.0  0.0
 0.0   5.0  1.0  0.0  0.0  0.0
 0.0   9.0  1.0  1.0  0.0  0.0
 0.0  13.0  2.0  0.0  1.0  0.0
 0.0  17.0  3.0  0.0  0.0  1.0

julia> F*[10.0,20,30,40,50,60.0]
6-element Array{Float64,1}:
  10.0
  20.0
  30.0
  70.0
 110.0
 150.0
```
"""
struct Frobenius{T} <: AbstractArray{T, 2}
    colidx :: Int
    c :: Vector{T}
end

#Basic property computations
size(F::Frobenius, r::Int) = (r==1 || r==2) ? F.colidx + length(F.c) :
    throw(ArgumentError("Frobenius matrix is of rank 2"))

function size(F::Frobenius)
    n = F.colidx + length(F.c)
    n, n
end

#XXX Inefficient but works
getindex(F::Frobenius, i, j) = getindex(Matrix(F), i, j)
isassigned(F::Frobenius, i::Integer, j::Integer) = isassigned(Matrix(F), i, j)

function Matrix(F::Frobenius{T}) where T
    M = Matrix{T}(I, size(F))
    M[F.colidx+1:end, F.colidx] = F.c
    M
end

#Linear algebra stuff
function mul!(F::Frobenius{T}, b::Vector{T}) where T
    (n = size(F, 2)) == length(b) || throw(DimensionMismatch("$n $(length(b))"))
    for i=F.colidx+1:n
        b[i] += F.c[i-F.colidx] * b[F.colidx]
    end
    b
end
*(F::Frobenius{T}, b::Vector{T}) where T = mul!(F, copy(b))

function *(F::Frobenius{T}, G::Frobenius{T}) where T
    (n = size(F, 2)) == size(G, 2) || throw(DimensionMismatch(""))
    if F.colidx == G.colidx #Answer is still expressible as a Frobenius
        return Frobenius(F.colidx, F.c + G.c)
    else
        M = Matrix(F)
        M[G.colidx+1:end, G.colidx] = F.colidx < G.colidx ? G.c :
            Frobenius(F.colidx-G.colidx, F.c) * G.c
        return M
    end
end

inv(F::Frobenius) = Frobenius(F.colidx, -F.c)


text:
## Utils
module MapInLayer

using Flux

export MapIn

include("MapIn.jl")

end

max_stars_repo_name: williamjsdavis/ThermoChroNN

max_issues_repo_name: williamjsdavis/ThermoChroNN

max_forks_repo_name: williamjsdavis/ThermoChroNN


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"TOML Table"
struct Table
    values::Dict{String,Any}
    defined::Bool
end

Table(defined::Bool) = Table(Dict{String,Any}(), defined)
function Base.show(io::IO, tbl::Table, level::Int=1)
    Base.print(io, "T($(tbl.defined)){\n")
    for (k,v) in tbl.values
        Base.print(io, "\t"^level, k, " => ")
        if isa(v, Table)
            Base.show(io, v, level+1)
            Base.print(io, "\n")
        elseif isa(v, Vector{Table})
            Base.print(io, "[\n")
            for i in 1:length(v)
                Base.print(io, "\t"^level, "  ")
                Base.show(io, v[i], level+1)
                i != length(v) && Base.print(io, ",\n")
            end
            Base.print("\t"^level, "]\n")
        else
            Base.print(io, " $v\n")
        end
    end
    Base.print(io, "\t"^(level-1), "}\n")
end
Base.getindex(tbl::Table, key::AbstractString) = tbl.values[key]
Base.haskey(tbl::Table, key::AbstractString) = haskey(tbl.values ,key)

"Parser error exception"
struct ParserError <: Exception
    lo::Int
    hi::Int
    msg::String
end

"TOML Parser"
mutable struct Parser{IO_T <: IO}
    input::IO_T
    errors::Vector{ParserError}
    charbuffer::Base.GenericIOBuffer{Array{UInt8,1}}
    currentchar::Char

    Parser(input::IO_T) where {IO_T <: IO}  = new{IO_T}(input, ParserError[], IOBuffer(), ' ')
end
Parser(input::String) = Parser(IOBuffer(input))
Base.error(p::Parser, l, h, msg) = push!(p.errors, ParserError(l, h, msg))
Base.eof(p::Parser) = eof(p.input)
Base.position(p::Parser) = Int(position(p.input))+1
Base.write(p::Parser, x) = write(p.charbuffer, x)
Base.read(p::Parser) = p.currentchar = read(p.input, Char)

NONE() = nothing
NONE(::Type{T}) where {T} = nothing
SOME(v::T) where {T} = v
NONE(::Type{String}, p::Parser) = (take!(p.charbuffer); nothing)
SOME(::Type{String}, p::Parser) = String(take!(p.charbuffer))
isnull(x) = x === nothing
get(x) = (@assert !isnull(x); x)

"Rewind parser input on `n` characters."
function rewind(p::Parser, n=1)
    pos = position(p.input)
    pos == 0 && return 0
    skip(p.input, -n)
    return Int(position(p.input))
end

"Converts an offset to a line and a column in the original source."
function linecol(p::Parser, offset::Int)
    pos = position(p)
    seekstart(p.input)
    line = 0
    cur = 1
    for (i,l) in enumerate(eachline(p.input, chomp=false))
        if cur + length(l) > offset
            return (i, offset - cur + 1)
        end
        cur += length(l)
        line+=1
    end
    seek(p.input, pos)
    return (line, 0)
end

"Determine next position in parser input"
function nextpos(p::Parser)
    pos = position(p)
    return pos + 1
end

"Peeks ahead `n` characters"
function peek(p::Parser) #, i::Int=0
    eof(p) && return nothing
    res = Base.peek(p.input)
    res == -1 && return nothing
    return Char(res)
end

"Returns `true` and consumes the next character if it matches `ch`, otherwise do nothing and return `false`"
function consume(p::Parser, ch::AbstractChar)
    eof(p) && return false
    c = peek(p)
    if get(c) == ch
        read(p)
        return true
    else
        return false
    end
end

function expect(p::Parser, ch::AbstractChar)
    consume(p, ch) && return true
    lo = position(p)
    if eof(p)
        error(p, lo, lo, "expected `$ch`, but found EOF")
    else
        v = peek(p)
        mark(p.input)
        c = read(p)
        error(p, lo, lo+1, "expected `$ch`, but found `$c`")
        reset(p.input)
    end
    return false
end

"Consumes whitespace ('\t' and ' ') until another character (or EOF) is reached. Returns `true` if any whitespace was consumed"
function whitespace(p::Parser)
    ret = false
    while !eof(p)
        c = read(p)
        if c == '\t' || c == ' '
            ret = true
        else
            rewind(p)
            break
        end
    end
    return ret
end

"Consumes the rest of the line after a comment character"
function comment(p::Parser)
    !consume(p, '#') && return false
    while !eof(p)
        ch = read(p)
        ch == '\n' && break
    end
    return true
end

"Consumes a newline if one is next"
function newline(p::Parser)
    if !eof(p)
        n = 1
        ch = read(p.input, UInt8)
        ch == 0x0a && return true
        if ch == 0x0d
            if !eof(p)
                nch = read(p.input, UInt8)
                nch == 0x0a  && return true
                n+=1
            end
        end
        rewind(p, n)
    end
    return false
end

"Consume ignored symbols"
function ignore(p::Parser)
    while true
        whitespace(p)
        !newline(p) && !comment(p) && break
    end
end

"Parse a single key name starting at next position"
function keyname(p)
    s = nextpos(p)
    key = if consume(p, '"')
        basicstring(p, s, false)
    elseif consume(p, '\'')
        literalstring(p, s, false)
    else
        while !eof(p)
            ch = read(p)
            if  'a' <= ch <= 'z' ||
                'A' <= ch <= 'Z' ||
                isdigit(ch)      ||
                ch == '_'        ||
                ch == '-'
                write(p, ch)
            else
                rewind(p)
                break
            end
        end
        SOME(String, p)
    end

    if !isnull(key) && isempty(get(key))
        error(p, s, s, "expected a key but found an empty string")
        key = NONE(String, p)
    end
    return key
end

"Parse a path into a vector of paths"
function lookup(p::Parser)
    ks = String[]
    eof(p) && return SOME(ks)
    while true
        whitespace(p)
        s = keyname(p)
        if !isnull(s)
            push!(ks, get(s))
        else
            s = integer(p, 0)
            if !isnull(s)
                push!(ks, get(s))
            else
                return NONE()
            end
        end
        whitespace(p)
        !expect(p, '.') && return SOME(ks)
    end
end

"Parses a key-value separator"
function separator(p::Parser)
    whitespace(p)
    !expect(p, '=') && return false
    whitespace(p)
    return true
end

"Insert key-value pair to table, if duplicate key found reports error"
function insertpair(p::Parser, tbl::Table, k, v, st)
    if haskey(tbl, k)
        error(p, st, st+length(k), "duplicate key `$k`")
    else
        tbl.values[k] = v
    end
end

"Parses integer with leading zeros and sign"
function integer(p::Parser, st::Int, allow_leading_zeros::Bool=false, allow_sign::Bool=false)
    if allow_sign
        if consume(p, '-')
            write(p, '-')
        elseif consume(p, '+')
            write(p, '+')
        end
    end
    nch = peek(p)
    if isnull(nch)
        pos = nextpos(p)
        error(p, pos, pos, "expected start of a numeric literal")
        return NONE(String, p)
    else
        ch = get(nch)
        if isdigit(ch)
            c = read(p)
            if c == '0' && !allow_leading_zeros
                write(p, '0')
                nch = peek(p)
                if !isnull(nch)
                    ch = get(nch)
                    if isdigit(ch)
                        error(p, st, position(p), "leading zeroes are not allowed")
                        return NONE(String, p)
                    end
                end
            elseif isdigit(c)
                write(p, c)
            end
        else
            # non-digit
            error(p, st, position(p), "expected a digit, found `$ch`")
            return NONE(String, p)
        end
    end

    underscore = false
    while !eof(p)
        ch = read(p)
        if isdigit(ch)
            write(p, ch)
            underscore = false
        elseif ch == '_' && !underscore
            underscore = true
        else
            rewind(p)
            break
        end
    end
    if underscore
        pos = nextpos(p)
        error(p, pos, pos, "numeral cannot end with an underscore")
        NONE(String, p)
    else
        p.charbuffer.ptr == 1 ? NONE(String, p) : SOME(String, p)
    end
end

"Parses boolean"
function boolean(p::Parser, st::Int)
    ch = '\x00'

    i = 0
    word = "true"
    while !eof(p) && i<length(word)
        ch = read(p)
        i+=1
        ch != word[i] && break
    end
    if i == length(word)
        return SOME(true)
    else
        rewind(p, i)
    end
    ti = i
    tch = ch

    i = 0
    word = "false"
    while !eof(p) && i<length(word)
        ch = read(p)
        i+=1
        ch != word[i] && break
    end
    if i == length(word)
        return SOME(false)
    else
        rewind(p, i)
        if i < ti
            i = ti
            ch = tch
        end
        error(p, st, st+i-1, "unexpected character: `$ch`")
    end

    return NONE(Bool)
end

"Parse number or datetime"
function numdatetime(p::Parser, st::Int)
    isfloat = false

    nprefix = integer(p, st, false, true)
    isnull(nprefix) && return NONE()
    prefix = get(nprefix)

    decimal = if consume(p, '.')
        isfloat = true
        ndecimal = integer(p, st, true, false)
        isnull(ndecimal) && return NONE()
        SOME(get(ndecimal))
    else
        NONE(String, p)
    end

    exponent = if consume(p,'e') || consume(p,'E')
        isfloat = true;
        nexponent = integer(p, st, false, true)
        isnull(nexponent) && return NONE()
        SOME(get(nexponent))
    else
        NONE(String, p)
    end

    pend = nextpos(p)
    nch = peek(p)
    ret = if isnull(decimal) &&
             isnull(exponent) &&
             nch != '+' &&
             nch != '+' &&
             st + 4 == pend &&
             consume(p, '-')
        datetime(p, prefix, st)
    else
        input = if isnull(decimal) && isnull(exponent) #!isfloat
            prefix
        elseif !isnull(decimal) && isnull(exponent)
            "$(prefix)."*get(decimal,"")
        elseif isnull(decimal) && !isnull(exponent)
            "$(prefix)E"*get(exponent,"")
        elseif !isnull(decimal) && !isnull(exponent)
            "$(prefix)."*get(decimal,"")*"E"*get(exponent,"")
        end
        input = lstrip(input, '+')
        try
            SOME(Base.parse(isfloat ? Float64 : Int, input))
        catch
            NONE()
        end
    end
    isnull(ret) && error(p, st, pend, "invalid numeric literal")
    return ret
end

"Parses a datetime value"
function datetime(p::Parser, syear::String, st::Int)

    valid = true

    function parsetwodigits(p, valid)
        eof(p) && return 0, false
        ch1 = read(p)
        valid = valid && isdigit(ch1)
        eof(p) && return 0, false
        ch2 = read(p)
        valid = valid && isdigit(ch2)
        (valid ? Base.parse(Int, String([ch1, ch2])) : 0), valid
    end

    year = Base.parse(Int, syear)
    month,  valid = parsetwodigits(p, valid)
    valid = valid && consume(p, '-')
    day,    valid = parsetwodigits(p, valid)
    valid = valid && consume(p, 'T')
    hour,   valid = parsetwodigits(p, valid)
    valid = valid && consume(p, ':')
    minute, valid = parsetwodigits(p, valid)
    valid = valid && consume(p, ':')
    second, valid = parsetwodigits(p, valid)

    # fractional seconds
    msec = 0
    if consume(p, '.')
        fsec = Char[]
        ch = peek(p)
        valid = valid && !isnull(ch) && isdigit(get(ch))
        while true
            ch = peek(p)
            if !isnull(ch) && isdigit(get(ch))
                push!(fsec, read(p))
            else
                break
            end
        end
        if length(fsec)>0
            msec = Base.parse(Int, String(fsec))
        end
    end

    # time zone
    tzhour=0
    tzminute=0
    tzplus = true
    tzminus = false
    tzsign = true
    if valid && !consume(p, 'Z')
        tzplus = consume(p, '+')
        if !tzplus
            tzminus = consume(p, '-')
        end
        valid = valid && (tzplus || tzminus)
        tzhour,   valid = parsetwodigits(p, valid)
        valid =   valid && consume(p, ':')
        tzminute, valid = parsetwodigits(p, valid)

        tzsign = tzplus
    end

    if valid
        dt = DateTime(year, month, day,
                hour + (tzsign ? tzhour : -tzhour),
                minute + (tzsign ? tzminute : -tzminute),
                second, msec)
        return SOME(dt)
    else
        error(p, st, position(p), "malformed date literal")
        return NONE()
    end
end

"Parses a single or multi-line string"
function basicstring(p::Parser, st::Int)
    !expect(p, '"') && return NONE(String, p)

    multiline = false

    if consume(p, '"')
        if consume(p, '"')
            multiline = true
            newline(p)
        else
            return SOME("")
        end
    end

    basicstring(p, st, multiline)
end

"Finish parsing a basic string after the opening quote has been seen"
function basicstring(p::Parser, st::Int, multiline::Bool)
    while true
        while multiline && newline(p)
            write(p, '\n')
        end
        if eof(p)
            pos = position(p)
            error(p, st, pos, "unterminated string literal")
            return NONE(String, p)
        else
            ch = read(p)
            if ch == '"'
                if multiline
                    if !consume(p, '"')
                        write(p, '"')
                        continue
                    end
                    if !consume(p, '"')
                        write(p, '"')
                        write(p, '"')
                        continue
                    end
                end
                return SOME(String, p)
            elseif ch == '\\'
                pos = position(p)
                ec = escape(p, pos, multiline)
                !isnull(ec) && write(p, get(ec))
            elseif ch < '\x1f'
                pos = position(p)
                error(p, st, pos, "control character `$ch` must be escaped")
            else
                write(p, ch)
            end
        end
    end
end

"Reads character(s) after `\\` and transforms them into proper character"
function escape(p::Parser, st::Int, multiline::Bool)
    if multiline && newline(p)
        while whitespace(p) || newline(p) end
        return NONE()
    end
    pos = position(p)
    if eof(p)
        error(p, st, pos, "unterminated escape sequence")
        return NONE()
    else
        ch = read(p)
        if ch == 'b'
            SOME('\b')
        elseif ch == 't'
            SOME('\t')
        elseif ch == 'n'
            SOME('\n')
        elseif ch == 'f'
            SOME('\f')
        elseif ch == 'r'
            SOME('\r')
        elseif ch == '"'
            SOME('\"')
        elseif ch == '\\'
            SOME('\\')
        elseif ch == 'u' || ch == 'U'
            len = ch == 'u' ? 4 : 8
            ucstr = ch == 'u' ? "\\u" : "\\U"
            snum = String(read(p.input, len))
            try
                if length(snum) < len
                    error(p, st, st+len, "expected $len hex digits after a `$ch` escape")
                    NONE()
                end
                if !all(isxdigit, snum)
                    error(p, st, st+len, "unknown string escape: `$snum`")
                    NONE()
                end
                c = unescape_string(ucstr * snum)[1]
                SOME(c)
            catch
                error(p, st, st+len, "codepoint `$snum` is not a valid unicode codepoint")
                rewind(p, len)
                NONE()
            end
        else
            error(p, st, position(p), "unknown string escape: `$("\\x"*hex(ch,2))`")
            NONE()
        end
    end
end

"Parses a single or multi-line literal string"
function literalstring(p::Parser, st::Int)
    !expect(p, '\'') && return NONE(String, p)

    multiline = false

    if consume(p, '\'')
        if consume(p, '\'')
            multiline = true
            newline(p)
        else
            return SOME("")
        end
    end

    literalstring(p, st, multiline)
end

function literalstring(p::Parser, st::Int, multiline::Bool)
    while true
        if !multiline && newline(p)
            npos = nextpos(p)
            error(p, st, npos, "literal strings cannot contain newlines")
            return NONE(String, p)
        end
        if eof(p)
            error(p, st, position(p), "unterminated string literal")
            return NONE(String, p)
        else
            ch = read(p.input, UInt8)
            if ch == 0x27
                if multiline
                    if !consume(p, '\'')
                        write(p, 0x27)
                        continue
                    end
                    if !consume(p, '\'')
                        write(p, 0x27)
                        write(p, 0x27)
                        continue
                    end
                end
                return SOME(String, p)
            else
                write(p, ch)
            end
        end
    end
end

function array(p::Parser, st::Int)
    !expect(p, '[') && return NONE()
    ret = Any[]
    rettype = Any
    expected = Any
    while true

        # Break out early if we see the closing bracket
        ignore(p)
        consume(p, ']') && return SOME(ret)

        # Attempt to parse a value, triggering an error if it's the wrong type.
        pstart = nextpos(p)
        npvalue = value(p)
        isnull(npvalue) && return NONE()
        pvalue = get(npvalue)

        pend = nextpos(p)
        valtype = isa(pvalue, Array) ? Array : typeof(pvalue)
        expected = rettype === Any ? valtype : expected
        if valtype != expected
            error(p, pstart, pend, "expected type `$expected`, found type `$valtype`")
        else
            rettype = expected
            push!(ret, pvalue)
        end

        # Look for a comma. If we don't find one we're done
        ignore(p)
        !consume(p, ',') && break
    end
    ignore(p)
    !expect(p, ']') && return NONE()
    return SOME(convert(Vector{rettype}, ret))
end

function inlinetable(p::Parser, st::Int)
    !expect(p, '{') && return NONE()
    whitespace(p)

    ret = Table(true)
    consume(p, '}') && return SOME(ret)

    while true
        npos = nextpos(p)
        k = keyname(p)
        isnull(k) && return NONE()
        !separator(p) && return NONE()
        v = value(p)
        isnull(v) && return NONE()
        insertpair(p, ret, get(k), get(v), npos)

        whitespace(p)
        consume(p, '}') && break
        !expect(p, ',') && return NONE()
        whitespace(p)
    end

    return SOME(ret)
end

"Parses a value"
function value(p::Parser)
    whitespace(p)
    c = peek(p)
    isnull(c) && return NONE()
    ch = get(c)
    pos = position(p)+1
    if ch == '"'
        return basicstring(p, pos)
    elseif ch == '\''
        return literalstring(p, pos)
    elseif ch == 't' || ch == 'f'
        return boolean(p, pos)
    elseif ch == '['
        return array(p, pos)
    elseif ch == '{'
        return inlinetable(p, pos)
    elseif ch == '-' || ch == '+' || isdigit(ch)
        return numdatetime(p, pos)
    else
        error(p, pos, pos+1, "expected a value")
        return NONE()
    end
end

"Parses the key-values and fills the given dictionary. Returns true in case of success and false in case of error."
function keyvalues(p::Parser, tbl::Table)
    while true
        whitespace(p)
        newline(p) && continue
        comment(p) && continue
        nc = peek(p)
        isnull(nc) && break
        get(nc) == '[' && break

        # get key
        klo = nextpos(p)
        k = keyname(p)
        isnull(k) && return false

        # skip kv separator
        !separator(p) && return false

        # get value
        v = value(p)
        isnull(v) && return false
        vend = nextpos(p)

        # insert kv into result table
        insertpair(p, tbl, get(k), get(v), klo)

        whitespace(p)
        if !comment(p) && !newline(p)
            pc = peek(p)
            isnull(pc) && return true
            error(p, klo, vend, "expected a newline after a key")
            return false
        end
    end
    return true
end

"Construct nested structures from keys"
function nested(p, into, ks, kstart)
    cnode = into
    kend = 0
    for k in ks[1:end-1]
        kend += length(k)+1
        if !haskey(cnode, k)
            cnode.values[k] = Table(false)
            cnode = cnode[k]
        else
            tmp = cnode[k]
            if isa(tmp, Table)
                cnode = tmp
            elseif isa(tmp, Array)
                if length(tmp)>0 && typeof(tmp[end]) === Table
                    cnode = tmp[end]
                else
                    error(p, kstart, kstart+kend, "array `$k` does not contain tables")
                    return NONE(), kend
                end
            else
                error(p, kstart, kstart+kend, "key `$k` was not previously a table")
                return NONE(), kend
            end
        end
    end
    return SOME(cnode), kend
end

function addtable(p::Parser, into::Table, ks::Vector{String}, tbl::Table, kstart)

    cnode, kend = nested(p, into, ks, kstart)
    isnull(cnode) && return
    cur = get(cnode)

    # fill last level with values
    tkey = ks[end]
    if haskey(cur, tkey)
        ctbl = cur[tkey]
        if isa(ctbl, Table)
            ctbl.defined && error(p, kstart, kstart+kend+length(tkey), "redefinition of table `$tkey`")

            for (k,v) in tbl.values
                insertpair(p, ctbl, k, v, kstart)
            end
        else
            error(p, kstart, kstart+kend+length(tkey), "key `$tkey` was not previously a table")
        end
    else
        cur.values[tkey] = tbl
    end

end

function addarray(p::Parser, into::Table, ks::Vector{String}, val, kstart)

    cnode, kend = nested(p, into, ks, kstart)
    isnull(cnode) && return
    cur = get(cnode)

    akey = ks[end]
    if haskey(cur, akey)
        vec = cur[akey]
        if isa(vec, Array)
            if isa(val, eltype(vec))
                push!(vec, val)
            else
                error(p, kstart, kstart+kend+length(akey), "expected type `$(typeof(val))`, found type `$(eltype(vec))`")
            end
        else
            error(p, kstart, kstart+kend+length(akey), "key `$akey` was previously not an array")
        end
    else
        cur.values[akey] = Any[val]
    end

end

"""Executes the parser, parsing the string contained within.

This function will return the `Table` instance if parsing is successful, or it will return `nothing` if any parse error or invalid TOML error occurs.

If an error occurs, the `errors` field of this parser can be consulted to determine the cause of the parse failure.
"""
function parse(p::Parser)
    ret = Table(false)
    while !eof(p)
        whitespace(p)
        newline(p) && continue
        comment(p) && continue

        # parse section
        if consume(p, '[')
            arr = consume(p, '[')
            npos = nextpos(p)

            # parse the name of the section
            ks = String[]
            while true
                whitespace(p)
                s = keyname(p)
                !isnull(s) && push!(ks, get(s))
                whitespace(p)
                if consume(p, ']')
                    arr && !expect(p, ']') && return NONE(Table)
                    break
                end
                !expect(p, '.') && return NONE()
            end
            isempty(ks) && return NONE(Table)

            # build the section
            section = Table(true)
            !keyvalues(p, section) && return NONE(Table)
            if arr
                addarray(p, ret, ks, section, npos)
            else
                addtable(p, ret, ks, section, npos)
            end
        else
            !keyvalues(p, ret) && return NONE(Table)
        end
    end

    length(p.errors) != 0 && return NONE(Table)

    return SOME(ret)
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
# # Inverse Lighting Tutorial
#
# In this tutorial we shall explore the inverse lighting problem.
# Here, we shall try to reconstruct a target image by optimizing
# the parameters of the light source (using gradients).
# From https://github.com/avik-pal/RayTracer.jl/ by Avik Pal

using RayTracer, Images, Zygote, Flux, Statistics, Random

# ## Configuring the Scene
#
# Reduce the screen_size if the optimization is taking a bit long
screen_size = (w = 180, h = 180)

scene = load_obj("./tree.obj")

# Let us set up the [`Camera`](@ref). For a more detailed understanding of
# the rendering process look into [Introduction to rendering using RayTracer.jl](@ref).
cam = Camera(
    lookfrom = Vec3(0.0f0, 6.0f0, -10.0f0),
    lookat   = Vec3(0.0f0, 2.0f0,  0.0f0),
    vup      = Vec3(0.0f0, 1.0f0,  0.0f0),
    vfov     = 45.0f0,
    focus    = 0.5f0,
    width    = screen_size.w,
    height   = screen_size.h
)

origin, direction = get_primary_rays(cam)

# We should define a few convenience functions. Since we are going to calculate
# the gradients only wrt to `light` we have it as an argument to the function. Having
# `scene` as an additional parameters simply allows us to test our method for other
# meshes without having to run `Zygote.refresh()` repeatedly.
function render(light, scene)
    packed_image = raytrace(origin, direction, scene, light, origin, 2)
    array_image = reshape(hcat(packed_image.x, packed_image.y, packed_image.z),
                          (screen_size.w, screen_size.h, 3, 1))
    return array_image
end

showimg(img) = colorview(RGB, permutedims(zeroonenorm(img)[:,:,:,1], (3,2,1)))

# ## [Ground Truth Image](@id inv_light)
#
# For this tutorial we shall use the [`PointLight`](@ref) source.
# We define the ground truth lighting source and the rendered image. We
# will later assume that we have no information about this lighting
# condition and try to reconstruct the image.
target_light = PointLight(
    color     = Vec3(1.0f0, 1.0f0, 1.0f0),
    intensity = 20000.0f0,
    position  = Vec3(1.0f0, 10.0f0, -50.0f0)
)

target_img = render(target_light, scene)

# The presence of [`zeroonenorm`](@ref) is very important here. It rescales the
# values in the image to 0 to 1. If we don't perform this step `Images` will
# clamp the values while generating the image in RGB format.
showimg(target_img)

# ## Initial Guess of Lighting Parameters
#
# We shall make some arbitrary guess of the lighting parameters (intensity and
# position) and try to get back the image in [Ground Truth Image](@ref inv_light)
light_guess = PointLight(
    color     = Vec3(1.0f0, 1.0f0, 1.0f0),
    intensity = 1.0f0,
    position  = Vec3(-1.0f0, -10.0f0, -50.0f0)
)

showimg(render(light_guess, scene))

# We shall store the images in `results_inv_lighting` directory
mkpath("results_inv_lighting")

save("./results_inv_lighting/inv_light_original.png",
     showimg(render(target_light, scene)))
save("./results_inv_lighting/inv_light_initial.png",
     showimg(render(light_guess, scene)))

# loss function - sum of squares of differences between render and target
function loss_fn(lighting)
    loss = sum((render(lighting, scene) .- target_img) .^ 2)
    @show loss
    return loss
end

# ## Optimization Loop
#
# We will use the ADAM optimizer from Flux. (Try experimenting with other
# optimizers as well). We can also use frameworks like Optim.jl for optimization.
# We will show how to do it in a future tutorial
opt = ADAM(1.5)

@time for i in 1:296
    g = gradient(loss_fn, light_guess)
    update!(opt, light_guess.intensity, g[1].intensity)
    update!(opt, light_guess.position, g[1].position)
    if i % 5 == 1
        save("./results_inv_lighting/iteration_$i.png",
             showimg(render(light_guess, scene)))
    end
end


text:
"""
    logged_solve(prob, args...; kwargs...)

Create a `Logged` ODE solution whose logged variables can be accessed through the `.log`
property.

See also: [`get_log`](@ref), [`solve`](@ref)
"""
function logged_solve(prob, args...; callback=nothing, kwargs...)
    p = deepcopy(prob.p)
    sol = solve(deepcopy(prob), args...; callback=callback, kwargs...)
    log = _get_log(isinplace(prob), prob.f, sol.u, p, sol.t, callback)
    return Logged(sol, log)
end

"""
    get_log(sol; callback=nothing)
    get_log(sol, t; callback=nothing)

Get the variables logged by `@log` from an `ODESolution`. If a callback or callback set
was used to change parameters during the simulation, it must be passed in through the keyword
`callback` to obtain correct results. When doing this, be sure to reset any parameters that
changed over the course of the simulation to their starting values.

Alternatively, replace your `solve` call with `logged_solve` to handle this all automatically.

See also: [`logged_solve`](@ref)
"""
get_log(sol; callback=nothing) = _get_log(isinplace(sol.prob), sol.prob.f, sol.u, sol.prob.p, sol.t, callback)
get_log(sol, t; callback=nothing) = _get_log(isinplace(sol.prob), sol.prob.f, sol(t).u, sol.prob.p, t, callback)
get_log(sol::Logged) = getfield(sol, :log)
get_log(sol::Logged, t) = error("Sorry, `get_log(sol::Logged, t)` does not currently work. Use `get_log(logged_solve(prob, ...; saveat=t))` or `get_log(solve(prob, ...), t)`")

function _get_log(iip, f, u, p, t, callback)
    u = iip ? (zero.(u), u) : (u,)

    reset!()
    activate!()

    _run_funcs(f, u, deepcopy(p), t, callback)

    deactivate!()
    log = deepcopy(GLOBAL_LOG)
    reset!()

    return _postprocess(log, length(t))
end

function _run_funcs(f, u, p, t, callbacks::CallbackSet)
    @assert isempty(callbacks.continuous_callbacks) "Sorry, SimulationLogs can't handle continuous callbacks yet"
    _run_funcs(f, u, p, t, callbacks.discrete_callbacks...)
end
function _run_funcs(f, u, p, t, ::Nothing)
    for i in eachindex(t)
        f(getindex.(u, i)..., p, t[i])
    end
end
function _run_funcs(f, du_u, p, t, callbacks::DiscreteCallback...)
    u = _get_u(du_u...)
    integrator = (; u=u, p=p, t=t)
    f(getindex.(_set_u(du_u..., integrator.u), 1)..., p, t[1])
    hit_last = true
    for i in 2:length(t)
        for callback in callbacks
            hit = callback.condition(integrator.u, t[i], integrator)
            # TODO: See if this works
            if hit==callback.save_positions[1] && hit_last==callback.save_positions[1]
                callback.affect!(integrator)
            end
            hit_last = hit
        end
        f(getindex.(_set_u(du_u..., integrator.u), i)..., p, t[i])
    end
end

_get_u(du, u) = u
_get_u(u) = u

_set_u(du, u_old, u_new) = (du, u_new)
_set_u(u_old, u_new) = (u_new,)

function _postprocess(log, n)
    for (key, val) in value_dict(log)
        if length(val) != n
            if length(val) % n == 0
                setproperty!(log, key, collect(permutedims(reshape(val, :, n))))
            else
                @warn """
                Signal $key was logged $(length(val)) times during $n timesteps. SimulationLogs
                is not currently set up to handle a mixed number of `@log` calls per step.
                """
                delete!(value_dict(log), key)
            end
        end
    end

    return log
end

max_stars_repo_name: jonniedie/SimulationLogs

max_issues_repo_name: jonniedie/SimulationLogs

max_forks_repo_name: jonniedie/SimulationLogs


text:
module PolrtestTest

using Test, OrdinalMultinomialModels, Random

Random.seed!(123)

n, p, J, q = 1000, 5, 7, 5
Xtrue = randn(n, p + q)
X = Xtrue[:, 1:p] # null model matrix
Z = Xtrue[:, p+1:p+q] # covariates to be tested
βtrue = [ones(p); 0.1ones(q)]
θ = collect(1.0:J-1)
α = [θ[1]; [log(θ[i] - θ[i-1]) for i in 2:J-1]]
link = LogitLink() # LogitLink(), ProbitLink(), CauchitLink(), CloglogLink()

Y = rpolr(Xtrue, βtrue, θ, link)

# Gradient based: LD_LBFGS, :LD_MMA, :LD_SLSQP, :LD_CCSAQ, :LD_TNEWTON_PRECOND_RESTART, :LD_TNEWTON_PRECOND, :LD_TNEWTON_RESTART, :LD_VAR2, :LD_VAR1
# Gradient free: :LN_COBYLA
# solver = NLoptSolver(algorithm=:LD_LBFGS)
solver = Ipopt.Optimizer() # more stable but take a lot more iterations
@time dd = polr(X, Y, link, solver;)
@show dd.θ, dd.β
@show stderror(dd)

# testing
@show polrtest(dd, Z; test = :score)
@show polrtest(dd, Z; test = :LRT)

end


text:
module RayTracer
    import Base
    using LinearAlgebra
    using UUIDs

    include("tuples.jl")
    include("canvas.jl")
    include("transformations.jl")
    include("rays.jl")

    include("light.jl")
    include("material.jl")

    include("sphere.jl")

end


text:
#=
 * sais
 * Copyright (c) 2008-2010 Yuta Mori All Rights Reserved.
 *
 * Permission is hereby granted, free of charge, to any person
 * obtaining a copy of this software and associated documentation
 * files (the "Software"), to deal in the Software without
 * restriction, including without limitation the rights to use,
 * copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following
 * conditions:
 *
 * The above copyright notice and this permission notice shall be
 * included in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
 * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
 * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
 * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
 * OTHER DEALINGS IN THE SOFTWARE.
 =#

#= Suffixsorting =#
type IntArray
    a::Array{Int,1}
    pos::Int
end
import Base: getindex, setindex!
getindex(a::IntArray,key) = a.a[a.pos + key]
setindex!(a::IntArray,value,key) = a.a[a.pos + key] = value

# "banana" = [5 3 1 0 4 2]
# "banana" = [6, 4, 2, 1, 5, 3]

#TODO
 #refactor code to simplify
 #build user interface for string operations

function getcounts(T,C,n,k)
    for i = 1:k
        C[i] = 0
    end
    for i = 1:n
        C[Int(T[i])+1] += 1
    end
end

function getbuckets(C,B,k,isend)
    sum = 0
    if isend != false
        for i = 1:k
            sum += C[i]
            B[i] = sum
        end
    else
        for i = 1:k
            sum += C[i]
            B[i] = sum - C[i]
        end
    end
end

function sais(T, SA, fs, n, k, isbwt)
    pidx = 0
    flags = 0
    if k <= 256
        C = IntArray(zeros(Int,k),0)
        if k <= fs
            B = IntArray(SA,n + fs - k)
            flags = 1
        else
            B = IntArray(zeros(Int,k),0)
            flags = 3
        end
    elseif k <= fs
        C = IntArray(SA,n + fs - k)
        if k <= fs - k
            B = IntArray(SA,n + fs - 2k)
            flags = 0
        elseif k <= 1024
            B = IntArray(zeros(Int,k),0)
            flags = 2
        else
            B = C
            flags = 8
        end
    else
        C = B = IntArray(zeros(Int,k),0)
        flags = 4 | 8
    end
    # stage 1
    getcounts(T,C,n,k)
    getbuckets(C,B,k,true)
    for i = 1:n
        SA[i] = 0
    end
    b = -1
    i = j = n
    m = 0
    c0 = c1 = T[n]
    i -= 1
    while 1 <= i && ((c0 = T[i]) >= c1)
        c1 = c0
        i -= 1
    end
    while 1 <= i
        c1 = c0
        i -= 1
        while 1 <= i && ((c0 = T[i]) <= c1)
            c1 = c0
            i -= 1
        end
        if 1 <= i
            0 <= b && (SA[b+1] = j)
            b = (B[Int(c1)+1] -= 1)
            j = i-1
            m += 1
            c1 = c0
            i -= 1
            while 1 <= i && ((c0 = T[i]) >= c1)
                c1 = c0
                i -= 1
            end
        end
    end
    if 1 < m
        LMSsort(T,SA,C,B,n,k)
        name = LMSpostproc(T,SA,n,m)
    elseif m == 1
        SA[b+1] = j+1
        name = 1
    else
        name = 0
    end
    # stage 2
    if name < m
        newfs = n + fs - 2m
        if flags & (1 | 4 | 8) == 0
            if (k + name) <= newfs
                newfs -= k
            else
                flags |= 8
            end
        end
        j = 2m + newfs
        for i = (m + (n >> 1)):-1:(m+1)
            if SA[i] != 0
                SA[j] = SA[i] - 1
                j -= 1
            end
        end
        RA = IntArray(SA, m + newfs)
        sais(RA,SA,newfs,m,name,false)

        i = n
        j = 2m
        c0 = c1 = T[n]
        while 1 <= (i -= 1) && ((c0 = T[i]) >= c1)
            c1 = c0
        end
        while 1 <= i
            c1 = c0
            while 1 <= (i -= 1) && ((c0 = T[i]) <= c1)
                c1 = c0
            end
            if 1 <= i
                SA[j] = i
                j -= 1
                c1 = c0
                while 1 <= (i -= 1) && ((c0 = T[i]) >= c1)
                    c1 = c0
                end
            end
        end
        for i = 1:m
            SA[i] = SA[m+SA[i]+1]
        end
        if flags & 4 != 0
            C = B = IntArray(zeros(Int,k),0)
        end
        if flags & 2 != 0
            B = IntArray(zeros(Int,k),0)
        end
    end
    # stage 3
    flags & 8 != 0 && getcounts(T,C,n,k)
    if 1 < m
        getbuckets(C,B,k,true)
        i = m-1
        j = n
        p = SA[m]
        c1 = T[p+1]
        while true
            c0 = c1
            q = B[Int(c0)+1]
            while q < j
                j -= 1
                SA[j+1] = 0
            end
            while true
                j -= 1
                SA[j+1] = p
                i -= 1
                i < 0 && break
                p = SA[i+1]
                c1 = T[p+1]
                c1 != c0 && break
            end
            i < 0 && break
        end
        while 0 < j
            j -= 1
            SA[j+1] = 0
        end
    end
    if isbwt == false
        induceSA(T,SA,C,B,n,k)
    else
        pidx = computeBWT(T,SA,C,B,n,k)
    end
    return SA
end

function LMSsort(T, SA, C, B, n, k)
    C == B && getcounts(T,C,n,k)
    getbuckets(C,B,k,false)
    j = n - 1
    c1 = T[j+1]
    b = B[Int(c1)+1]
    j -= 1
    SA[b+1] = T[j+1] < c1 ? ~j : j
    b += 1
    for i = 1:n
        if 0 < (j = SA[i])
            if (c0 = T[j+1]) != c1
                B[Int(c1)+1] = b
                c1 = c0
                b = B[Int(c1)+1]
            end
            j -= 1
            SA[b+1] = T[j+1] < c1 ? ~j : j
            b += 1
            SA[i] = 0
        elseif j < 0
            SA[i] = ~j
        end
    end
    C == B && getcounts(T,C,n,k)
    getbuckets(C,B,k,true)
    c1 = 0
    b = B[c1+1]
    for i = n:-1:1
        if 0 < (j = SA[i])
            c0 = T[j+1]
            if c0 != c1
                B[c1+1] = b
                c1 = c0
                b = B[Int(c1)+1]
            end
            j -= 1
            b -= 1
            SA[b+1] = T[j+1] > c1 ? ~(j+1) : j
            SA[i] = 0
        end
    end
end

function LMSpostproc(T,SA,n,m)
    i = 1
    while (p = SA[i]) < 0
        SA[i] = ~p
        i += 1
    end
    if i-1 < m
        j = i
        i += 1
        while true
            if (p = SA[i]) < 0
                SA[j] = ~p
                j += 1
                SA[i] = 0
                j-1 == m && break
            end
            i += 1
        end
    end

    i = j = n
    c0 = c1 = T[n]
    while 1 <= (i -= 1) && ((c0 = T[i]) >= c1)
        c1 = c0
    end
    while 1 <= i
        c1 = c0
        while 1 <= (i -= 1) && ((c0 = T[i]) <= c1)
            c1 = c0
        end
        if 1 <= i
            SA[m + (i >> 1) + 1] = j - i
            j = i + 1
            c1 = c0
            while 1 <= (i -= 1) && ((c0 = T[i]) >= c1)
                c1 = c0
            end
        end
    end
    name = 0
    q = n
    qlen = 0
    for i = 1:m
        p = SA[i]
        plen = SA[m + (p >> 1) + 1]
        diff = true
        if plen == qlen && (q + plen < n)
            j = 0
            while j < plen && T[p+j+1] == T[q+j+1]
                j += 1
            end
            j == plen && (diff = false)
        end
        if diff != false
            name += 1
            q = p
            qlen = plen
        end
        SA[m + (p >> 1) + 1] = name
    end
    return name
end

function induceSA(T,SA,C,B,n,k)
    C == B && getcounts(T,C,n,k)
    getbuckets(C,B,k,false)
    j = n - 1
    c1 = T[j+1]
    b = B[Int(c1)+1]
    SA[b+1] = 0 < j && T[j] < c1 ? ~j : j
    b += 1
    for i = 1:n
        j = SA[i]
        SA[i] = ~j
        if 0 < j
            j -= 1
            if (c0 = T[j+1]) != c1
                B[Int(c1)+1] = b
                c1 = c0
                b = B[Int(c1)+1]
            end
            SA[b+1] = 0 < j && T[j] < c1 ? ~j : j
            b += 1
        end
    end
    C == B && getcounts(T,C,n,k)
    getbuckets(C,B,k,true)
    c1 = 0
    b = B[c1+1]
    for i = n:-1:1
        if 0 < (j = SA[i])
            j -= 1
            if (c0 = T[j+1]) != c1
                B[Int(c1)+1] = b
                c1 = c0
                b = B[Int(c1)+1]
            end
            b -= 1
            SA[b+1] = j == 0 || T[j] > c1 ? ~j : j
        else
            SA[i] = ~j
        end
    end
end

function computeBWT(T,SA,C,B,n,k)
    pidx = -1
    C == B && getcounts(T,C,n,k)
    getbuckets(C,B,k,false)
    j = n-1
    c1 = T[j+1]
    b = B[Int(c1)+1]
    SA[b+1] = 0 < j && T[j] < c1 ? ~j : j
    b += 1
    for i = 1:n
        if 0 < (j = SA[i])
            j -= 1
            c0 = T[j+1]
            SA[i] = ~c0
            if c0 != c1
                B[Int(c1)+1] = b
                c1 = c0
                b = B[c1+1]
            end
            SA[b+1] = 0 < j && T[j] < c1 ? ~j : j
            b += 1
        elseif j != 0
            SA[i] = ~j
        end
    end
    C == B && getcounts(T,C,n,k)
    getbuckets(C,B,k,true)
    c1 = 0
    b = B[Int(c1)+1]
    for i = n:-1:1
        if 0 < (j = SA[i])
            j -= 1
            c0 = T[j+1]
            SA[i] = c0
            if c0 != c1
                B[Int(c1)+1] = b
                c1 = c0
                b = B[Int(c1)+1]
            end
            b -= 1
            SA[b+1] = 0 < j && T[j] > c1 ? ~(T[j]) : j
        elseif j != 0
            SA[i] = ~j
        else
            pidx = i-1
        end
    end
    return pidx
end