text:
struct BjerksundStenslandApproximationEngine{B <: AbstractBlackScholesProcess} <: PricingEngine
  process::B
end

function get_phi(S::Float64, gamma::Float64, H::Float64, I::Float64, rT::Float64, bT::Float64, variance::Float64)
  lambda = (-rT + gamma * bT + 0.5 * gamma * (gamma - 1.0) * variance)
  d = -(log(S / H) + (bT + (gamma - 0.5) * variance)) / sqrt(variance)
  kappa = 2.0 * bT / variance + (2.0 * gamma - 1.0)
  w = Normal()
  return exp(lambda) * (cdf(w, d) - ^((I / S), kappa) * cdf(w, d - 2.0 * log(I / S) / sqrt(variance)))
end

function american_call_approximation(S::Float64, X::Float64, rfD::Float64, dD::Float64, variance::Float64)
  bT = log(dD / rfD)
  rT = log(1.0 / rfD)

  beta = (0.5 - bT / variance) + sqrt(^((bT / variance - 0.5), 2.0) + 2.0 * rT / variance)
  BInfinity = beta / (beta - 1.0) * X
  B0 = max(X, rT / (rT - bT) * X)
  ht = -(bT + 2.0 * sqrt(variance)) * B0 / (BInfinity - B0)

  # investigate what happens to I for db -> 0.0
  I = B0 + (BInfinity - B0) * (-expm1(ht))
  I >= X || error("Bjerksund - Stensland approximation not applicable to this set of params")

  if S >= I
    return S - X
  else
    # investigate what happens to alpha for dD -> 0.0
    return (I - X) * ^(S / I, beta) * (1.0 - get_phi(S, beta, I, I, rT, bT, variance)) +
        S * get_phi(S, 1.0, I, I, rT, bT, variance) -
        S * get_phi(S, 1.0, X, I, rT, bT, variance) -
        X * get_phi(S, 0.0, I, I, rT, bT, variance) +
        X * get_phi(S, 0.0, X, I, rT, bT, variance)
  end
end

function _calculate!(pe::BjerksundStenslandApproximationEngine, opt::AmericanOption)
  ex = opt.exercise
  payoff = opt.payoff

  variance = black_variance(pe.process.blackVolatility, ex.dates[end], payoff.strike)

  dividendDiscount = discount(pe.process.dividendYield, ex.dates[end])
  riskFreeDiscount = discount(pe.process.riskFreeRate, ex.dates[end])

  spot = state_variable(pe.process).value
  spot > 0.0 || error("negative or null underlying given")
  strike = payoff.strike

  if isa(payoff.optionType, Put)
    # use put-call symmetry
    spot, strike = strike, spot
    riskFreeDiscount, dividendDiscount = dividendDiscount, riskFreeDiscount

    payoff = PlainVanillaPayoff(Call(), strike)
  end

  if dividendDiscount >= 1.0
    # early exercise is never optimal - use Black formula
    forwardPrice = spot *  dividendDiscount / riskFreeDiscount
    black = BlackCalculator(payoff, forwardPrice, sqrt(variance), riskFreeDiscount)

    opt.results.value = value(black)
    opt.results.delta = delta(black, spot)
    opt.results.deltaForward = delta_forward(black)
    opt.results.elasticity = elasticity(black, spot)
    opt.results.gamma = gamma(black, spot)

    rfdc = pe.process.riskFreeRate.dc
    divdc = pe.process.dividendYield.dc
    voldc = pe.process.blackVolatility.dc

    t = year_fraction(rfdc, reference_date(pe.process.riskFreeRate), opt.exercise.dates[end])
    opt.results.rho = rho(black, t)

    t = year_fraction(divdc, reference_date(pe.process.dividendYield), opt.exercise.dates[end])
    opt.results.dividendRho = dividend_rho(black, t)

    t = year_fraction(voldc, reference_date(pe.process.blackVolatility), opt.exercise.dates[end])
    opt.results.vega = vega(black, t)

    opt.results.theta = theta(black, spot, t)
    opt.results.thetaPerDay = theta_per_day(black, spot, t)

    opt.results.strikeSensitivity = strike_sensitivity(black)
    opt.results.itmCashProbability = itm_cash_probability(black)
  else
    # early exericse can be optimal - use approximation
    opt.results.value = american_call_approximation(spot, strike, riskFreeDiscount, dividendDiscount, variance)
  end

  return opt
end


text:
# SimLynx.jl Process-Based Discrete-Event Simulation Example

#=
This is a simple process-based discrete-event simulation of an N teller, single
queue bank.

There are two (2) processes:
    generator(n::Integer)
      Generates n arrivals into the system with exponentially distributed
      inter-arrival time with a mean of 4.0.
    customer(i::Integer)
      The process representing the ith customer in the system. Each customer
      acquires a teller and works a uniformly distributed time between 2.0 and
      10.0. It then releases the teller and exits.

There is a single resource, tellers, that represents the tellers in the bank.
The number if tellers is set by the N_TELLERS global variable.

Once the simulation runs, statistics are printed for the tellers allocation and
queue length as well as a plot of the queue length over time.
=#

using SimLynx
SimLynx.greet()

using Distributions: Exponential, Uniform
using Random

const N_TELLERS = 2
const MEAN_INTERARRIVAL_TIME = 4.0
const MIN_SERVICE_TIME = 2.0
const MAX_SERVICE_TIME = 10.0

tellers = nothing

"Process to generate n customers arriving into the system."
@process generator(n::Integer) begin
    dist = Exponential(MEAN_INTERARRIVAL_TIME)
    for i = 1:n
        @schedule now customer(i)
        if i < n
            wait(rand(dist))
        end
    end
end

"The ith customer into the system."
@process customer(i::Integer) begin
    dist = Uniform(MIN_SERVICE_TIME, MAX_SERVICE_TIME)
    @with_resource tellers begin
        work(rand(dist))
    end
end

"Run the simulation for n customers."
function run_simulation(n::Integer)
    println("SimLynx.jl Process-Based Discrete-Event Simulation Example")
    println("$N_TELLERS Teller, Single Queue Bank Model")
    println("  Inter-arrival time = Exponential($MEAN_INTERARRIVAL_TIME)")
    println("  Service time = Uniform($MIN_SERVICE_TIME, $MAX_SERVICE_TIME)")
    println("  Number of customers = $n")
    @simulation begin
        global tellers = Resource(N_TELLERS, "tellers")
        current_trace!(true)
        @schedule at 0.0 generator(n)
        start_simulation()
        println("--- Simulation results after $(current_time()) minutes ---")
        print_stats(tellers.allocated, title="Teller Allocation Statistics")
        print_stats(tellers.queue_length, title="Teller Queue Length Statistics")
        plot_history(tellers.queue_length, title="Teller Queue Length History")
        print_stats(tellers.wait, title="Teller Queue Wait Time Statistics")
        plot_history(tellers.wait, title="Teller Queue Wait Time History")
    end
end

run_simulation(10)

ext: jl


text:
push!(LOAD_PATH,"../src/")
using Documenter, BFloat

makedocs(
    sitename = "BFloat.jl",
    modules  = [BFloat],
    doctest  = false
)

deploydocs(
    repo = "github.com/adknudson/BFloat.jl.git",
    versions = ["stable" => "v^", "v#.#", "dev" => "master"]
)


text:
using LightGraphs, GraphIO
using SimpleWeightedGraphs
using EzXML

function parse_graphtool_GraphML(doc)

    g = SimpleWeightedGraph(SimpleGraph(0))

    weight_name = "e_weight"
    weight_key_id = ""
    eid_name = "eid"
    eid_key_id = ""
    nodeid2index = Dict{String, Int64}()
    edgeid2index = Dict{String, Int64}()
    eid2index = Dict{String, Int64}()

    # TODO can be improved using:
    # parse.nodeids="canonical" parse.edgeids="canonical" parse.order="nodesfirst"
    for i in eachelement(root(doc))

        if (i.name == "key")
            key_id = i["id"]
            key_name = i["attr.name"]
            if (key_name == weight_name)
                weight_key_id = key_id
            elseif (key_name == eid_name)
                eid_key_id = key_id
            end
        end

        if (i.name == "graph")
            #node
            for elem in eachelement(i)

                # First pass: add nodes
                if (elem.name == "node")
                    n_nodes = length(nodeid2index)
                    node_id = elem["id"]
                    nodeid2index[node_id] = n_nodes + 1
                    add_vertex!(g)
                end
            end

            #edge
            for elem in eachelement(i)

                # Second pass: add edges
                if (elem.name == "edge")
                    src_id = nodeid2index[elem["source"]]
                    trg_id = nodeid2index[elem["target"]]
                    edge_id = elem["id"]

                    weight = 0.
                    eid = 0
                    has_weight = false
                    for data in eachelement(elem)
                        if (data.name == "data")
                            edata_key = data["key"]
                            if (edata_key == weight_key_id)
                                has_weight = true
                                weight = parse(Float64, data.content)
                            elseif (edata_key == eid_key_id)
                                eid = data.content
                            end
                        end
                    end

                    if (has_weight)
                        n_edges = length(edgeid2index)
                        edgeid2index[edge_id] = n_edges + 1
                        eid2index[eid] = n_edges + 1
                        add_edge!(g, src_id, trg_id, weight)
                    end
                end
            end
        end
    end

    return Dict(:graph=>g, :node_id=>nodeid2index, :edge_id=>edgeid2index, :eid=> eid2index)
end

max_stars_repo_name: kymok/local-betweenness

max_issues_repo_name: kymok/local-betweenness

max_forks_repo_name: kymok/local-betweenness


text:
module Bisquit

    # export all local vars from this module
    export Transaction, Blockchain, Block
    export new_transaction, new_block, create_genesis, hash, valid_proof, proof_of_work, last_block

    # include some utils here...
    # seems to be pretty instable
    # include("./include/utils.jl")

    # include some utils for chaining and formats
    using SHA # sha256bit for hashing
    using JSON # json for formatted output
    using Dates # Dates for tamestamp

    struct Transaction
        sender::String
        recipient::String
        amount::Int
    end

    struct Block
        index::Int
        timestamp::Dates.DateTime
        transactions::Array{Transaction}
        proof::Int
        previous_hash::String
    end

    mutable struct  Blockchain
        chain::Array{Block}
        current_transactions::Array{Transaction}
    end

    function new_transaction( blockchain::Blockchain, sender::String, recipient::String, amount::Int, )::int
        # add some values in here wip todo
        transaction = Transaction( sender, recipient, amount )
        push!( blockchain.current_transactions, transaction )
        return blockchain.chain[ end ].index + 1
    end

    function new_block( blockchain::Blockchain, proof::Int, previous_hash::String)::Block
        block = Block( blockchain.chain[ end ].index + 1, Dates.now( ), blockchain.current_transactions, proof, previous_hash )
        blockchain.current_transactions = [ ]
        push!( blockchain.chain, block )
        return block
    end

    function create_genesis( )::Block
        # " 0 " is set to 0 to vaildate if the current node hasn't already mined blocks
        genesis = Block(1, Dates.now( ), [], 100, "0" )
        blockchain = Blockchain( [ genesis ], [ ] )
        return blockchain
    end

    function hash(block::Block)::String

        # set list of block_string to all strings in bisquit and generate a sha256hash
        block_string = string( block.index, block.timestamp, block.transactions, block.proof, block.previous_hash )
        return bytes2hex( sha256( block_string ) )
    end

    function valid_proof(last_proof::Int, proof::Int)::Bool
        guess = "$last_proof$proof"
        guess_hash = bytes2hex( sha256( guess ) )

        # set the nonce value to 4layered 0000
        return guess_hash[ end - 3 : end ] == "0000"
    end

    function proof_of_work( last_proof::Int )::Int
        proof = 0

        # if proof is neg int
        if last_proof < 0
            print("Method error occured!")
        end

        # while proof == 0 add 1 recursively
        while valid_proof( last_proof, proof ) == 0
            proof += 1
        end

        return proof
    end

    function last_block( blockchain::Blockchain)::Block
        return blockchain.chain[ end ]
    end
end

max_stars_repo_name: timo-cmd/Bisquit


text:
function prune(foodweb::Matrix{Int64})

    S = size(foodweb, 1)

    # Identify basal species
    nprey = sum(foodweb, dims = 1)
    basal = findall(nprey[1, :] .== 0)

    # Identify basal species that have only specialized predators
    prune = zeros(0)
    for i in basal

        # Find all predators of basal species
        preds = findall(foodweb[i, :] .== 1)

        # If those predators each only have one prey, flag prey for removal
        if sum(nprey[preds]) == length(preds)
            append!(prune, i)
        end
    end

    # Prune any prey identified for pruning
    if length(prune) > 0
        retained = findall(x -> !(x in prune), 1:S)
        foodweb = foodweb[retained, retained]
    end

    return foodweb
end

max_stars_repo_name: clint-leach/niche-model

max_issues_repo_name: clint-leach/niche-model

max_forks_repo_name: clint-leach/niche-model


text:
function format_col(col, val)
    return col == 1 ? "\t$val" : val
end

function write_pfn(io::IO, network::PowerFlowNetwork)
    text = """function mpc = $(network.name)\n
%% MATPOWER Case Format : Version 2
mpc.version = '2';

%%-----  Power Flow Data  -----%%
%% system MVA base
mpc.baseMVA = $(network.baseMVA);

%% bus data
%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
mpc.bus = [
\t"""
    write(io, text)
    CSV.write(io, network.bus; delim='\t', append=true, newline=";\n\t")
    text = """];

%% generator data
%	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin	Pc1	Pc2	Qc1min	Qc1max	Qc2min	Qc2max	ramp_agc	ramp_10	ramp_30	ramp_q	apf
mpc.gen = [
\t"""
    write(io, text)
    CSV.write(io, network.gen; delim='\t', append=true, newline=";\n\t")
    text = """];

%% branch data
%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax
mpc.branch = [
\t"""
    write(io, text)
    CSV.write(io, network.branch; delim='\t', append=true, newline=";\n\t")
    text = """];

%%-----  OPF Data  -----%%
%% generator cost data
%	1	startup	shutdown	n	x1	y1	...	xn	yn
%	2	startup	shutdown	n	c(n-1)	...	c0
mpc.gencost = [
\t"""
    write(io, text)
    CSV.write(io, network.gencost; delim='\t', append=true, newline=";\n\t")
    text = "];"
    write(io, text)
end

function write_pfn(file::AbstractString, network::PowerFlowNetwork)
    io = open(file, "w")
    write_pfn(io, network)
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
"""
FunctionBlocks
"""

export PulseGeneratorBlock

mutable struct PulseGeneratorBlock <: AbstractFunctionBlock
    inport::Vector{InPort}
    outport::Vector{OutPort}
    amplitude::AbstractSignalLine
    period::AbstractSignalLine
    pulsewidth::AbstractSignalLine
    phasedelay::AbstractSignalLine

    function PulseGeneratorBlock(::Type{Tv} = Float64; amplitude=1, period=10, pulsewidth=5, phasedelay=0) where Tv
        @createblock new(Vector{InPort}(), Vector{OutPort}(),
            signal(amplitude, Tv),
            signal(period, Tv),
            signal(pulsewidth, Tv),
            signal(phasedelay, Tv)
            ) 0 1
    end
end

"""
IO
"""

function Base.show(io::IO, blk::PulseGeneratorBlock)
    println(io, "Pulse($(blk.amp), $(blk.cycle), $(blk.width), $(blk.delay))")
end

"""
to expr
"""

function _toexpr(blk::PulseGeneratorBlock)
    amplitude = _toexpr(blk.amplitude)
    period = _toexpr(blk.period)
    pulsewidth = _toexpr(blk.pulsewidth)
    phasedelay = _toexpr(blk.phasedelay)
    quote
        if t < $phasedelay
            0.0
        else
            u = ((t - $phasedelay) % $period) / $period * 100
            if u < $pulsewidth
                $amplitude
            else
                0
            end
        end
    end
end


text:
module InteractiveFixedEffectModels

##############################################################################
##
## Dependencies
##
##############################################################################

import Distances: chebyshev
import Base: length, copy!, axpy!, broadcast!, scale!, dot, similar, Ac_mul_B!, A_mul_B!, sumabs2!, map!, sumabs2, maxabs, fill!, norm, maxabs, size, length, eltype, rank, convert, view, clamp!, dot, vecdot
using Base.Cartesian
import StatsBase: coef, nobs, coeftable, vcov, predict, residuals, var, RegressionModel, model_response, stderr, confint, fit, CoefTable,  df_residual
import DataArrays: RefArray, PooledDataVector, DataVector, PooledDataArray, DataArray
import DataFrames: DataFrame, AbstractDataFrame, ModelMatrix, ModelFrame, Terms, coefnames, Formula, complete_cases, names!, pool
import LeastSquaresOptim
using Reexport
@reexport using FixedEffectModels
import FixedEffectModels: title, top, Ones

##############################################################################
##
## Exported methods and types
##
##############################################################################

export InteractiveFixedEffectModel,
InteractiveFixedEffectResult

##############################################################################
##
## Load files
##
##############################################################################
include("types.jl")

include("utils/models.jl")

include("methods/gauss_seidel.jl")
include("methods/ls.jl")

include("fit.jl")

end


text:
@testset "Transformer" begin
    t = Transformer(10, 3, 15, 20)
    td = TransformerDecoder(10, 3, 15, 20)
    x = randn(Float32, 10, 7, 3)
    y = randn(Float32, 10, 6, 3)

    m1 = reshape(Float32[1. 1. 1. 1. 0. 0. 0.; 1. 1. 1. 1. 1. 1. 0.; 1. 1. 1. 1. 1. 1. 1.]', 1, 7, 3)
    m2 = reshape(Float32[1. 1. 1. 0. 0. 0.; 1. 1. 1. 1. 1. 0.; 1. 1. 1. 1. 1. 1.]', 1, 6, 3)
    m = getmask(m1, m2)

    @test size(t(x)) == (10, 7, 3)
    @test size(t(x[:, :, 2])) == (10, 7)

    @test size(td(y, x)) == (10, 6, 3)
    @test size(td(y[:,:,2], x[:,:,2])) == (10, 6)

    @test_nowarn td(y, x, m)

    let STDOUT = stdout
        (outRead, outWrite) = redirect_stdout()
        print(outWrite, t)
        close(outWrite)

        output_string = String(readavailable(outRead))
        close(outRead)

        redirect_stdout(STDOUT)

        @test output_string == "Transformer(head=3, head_size=15, pwffn_size=20, size=10, dropout=0.1)"
    end

    let STDOUT = stdout
        (outRead, outWrite) = redirect_stdout()
        print(outWrite, t.mh)
        close(outWrite)

        output_string = String(readavailable(outRead))
        close(outRead)

        redirect_stdout(STDOUT)

        @test output_string == "MultiheadAttention(head=3, head_size=15, 10=>10, dropout=0.1)"
    end

    let STDOUT = stdout
        (outRead, outWrite) = redirect_stdout()
        print(outWrite, td)
        close(outWrite)

        output_string = String(readavailable(outRead))
        close(outRead)

        redirect_stdout(STDOUT)

        @test output_string == "TransformerDecoder(head=3, head_size=15, pwffn_size=20, size=10, dropout=0.1)"
    end

    Flux.testmode!(t, true)
    Flux.testmode!(td, true)

    let STDOUT = stdout
        (outRead, outWrite) = redirect_stdout()
        print(outWrite, t)
        close(outWrite)

        output_string = String(readavailable(outRead))
        close(outRead)

        redirect_stdout(STDOUT)

        @test output_string == "Transformer(head=3, head_size=15, pwffn_size=20, size=10)"
    end

    let STDOUT = stdout
        (outRead, outWrite) = redirect_stdout()
        print(outWrite, t.mh)
        close(outWrite)

        output_string = String(readavailable(outRead))
        close(outRead)

        redirect_stdout(STDOUT)

        @test output_string == "MultiheadAttention(head=3, head_size=15, 10=>10)"
    end

    let STDOUT = stdout
        (outRead, outWrite) = redirect_stdout()
        print(outWrite, td)
        close(outWrite)

        output_string = String(readavailable(outRead))
        close(outRead)

        redirect_stdout(STDOUT)

        @test output_string == "TransformerDecoder(head=3, head_size=15, pwffn_size=20, size=10)"
    end
end


text:
module Pretrain

using Markdown

using Flux
using BSON
using DataDeps
using ZipFile

import ..Transformers
using ..Basic
using ..Datasets: download_gdrive

export @pretrain_str, load_pretrain, pretrains

isbson(s) = endswith(s, ".bson")
iszip(s) = endswith(s, ".zip")
isnpbson(s) = endswith(s, ".npbson")
istfbson(s) = endswith(s, ".tfbson")

zipname(z::ZipFile.Reader) = z.files[1].name
zipfile(z::ZipFile.Reader, name) = (idx = findfirst(zf->isequal(name)(zf.name), z.files)) !== nothing ? z.files[idx] : nothing
findfile(z::ZipFile.Reader, name) = zipfile(z, joinpath(zipname(z), name))

include("./gpt_pretrain.jl")
include("./bert_pretrain.jl")

function __init__()
    gpt_init()
    bert_init()
end

"""
    pretrain"model-description:item"

convenient macro for loading data from pretrain. Use DataDeps to download automatically, if a model is not downlaod. the string should be in `pretrain"<model>-<model-name>:<item>"` format.

see also `Pretrain.pretrains()`.
"""
macro pretrain_str(name)
    :(load_pretrain($(esc(name))))
end

"""
    load_pretrain(name; kw...)

same as `@pretrain_str`, but can pass keyword argument if needed.
"""
function load_pretrain(name; kw...)
    lowered = lowercase(name)
    if startswith(lowered, "bert")
        col = findlast(isequal(':'), name)
        if col === nothing
            model_name = name[6:end]
            item = :all
        else
            model_name = name[6:col-1]
            item = Symbol(name[col+1:end])
        end
        !haskey(_get_bert_config, model_name) && error("unknown bert pretrain name")
        model_path = @datadep_str("BERT-$model_name/$model_name.tfbson")
        Transformers.load_bert_pretrain(model_path, item)
    elseif startswith(lowered, "gpt") && lowered[4] != '2'
        col = findlast(isequal(':'), name)
        if col === nothing
            model_name = name[5:end]
            item = :all
        else
            model_name = name[5:col-1]
            item = Symbol(name[col+1:end])
        end
        !haskey(_get_gpt_config, model_name) && error("unknown gpt pretrain name")
        model_path = @datadep_str("GPT-$model_name/$model_name.npbson")
        Transformers.load_gpt_pretrain(model_path, item; kw...)
    else
        error("unknown pretrain")
    end
end

"""
    pretrains(model::String = "")

Show all available model.
"""
function pretrains(model::String = "")
    rows = [Any["model", "model name", "support items"]]
    if isempty(model)
        for (model, configs) ∈ zip(("GPT", "Bert"), (_get_gpt_config, _get_bert_config))
            for config ∈ values(configs)
                name = config[:name]
                items = join(config[:items], ", ")
                push!(rows, Any[model, name, items])
            end
        end
    else
        lowered_model = lowercase(model)
        for (model, configs) ∈ zip(("GPT", "Bert"), (_get_gpt_config, _get_bert_config))
            if lowered_model == lowercase(model)
                for config ∈ values(configs)
                    name = config[:name]
                    items = join(config[:items], ", ")
                    push!(rows, Any[model, name, items])
                end
            end
        end
    end
    Markdown.MD(Any[Markdown.Table(rows, [:l, :l, :l])])
end

end


text:
using Setfield

export Vegas

abstract type VegasDamping end
struct NoDamping <: VegasDamping end

"""
Damping for quantile estimation as proposed in original Vegas Paper by Lepage.
"""
@qstruct LepageDamping{T}(;
    alpha::T=1) <: VegasDamping

@settable @qstruct Vegas{RNG,REG}(
        neval::Int=10^5;
        niter::Int=10,
        nbins::Int=max(2, neval ÷ (niter*100)), # number of bins along each dimension
        ndrop::Int=niter÷2,
        rng::RNG=GLOBAL_PAR_RNG,
        regularization::REG=LepageDamping(),
    ) do
        @argcheck neval > 2
        @argcheck niter >= 1
end <: MCAlgorithm

@qstruct VegasGrid{N,T}(
    boundaries::NTuple{N,Vector{T}}
   ) do
    for i in 1:N
        @argcheck length(boundaries[i]) >= 2
        @argcheck issorted(boundaries[i])
    end
end

function Base.:(==)(grid1::VegasGrid, grid2::VegasGrid)
    grid1.boundaries == grid2.boundaries
end

function domain(iq::VegasGrid)
    lo = map(first, iq.boundaries)
    hi = map(last , iq.boundaries)
    Domain(SVector(lo), SVector(hi))
end

function volume(iq::VegasGrid)
    volume(domain(iq))
end

function uniform(rng, iq::VegasGrid)
    uniform(rng, domain(iq))
end

function Base.size(s::VegasGrid)
    map(s.boundaries) do b
        length(b) - 1
    end
end
Base.length(s::VegasGrid) = prod(size(s))
Base.eachindex(s::VegasGrid) = CartesianIndices(s)
Base.axes(s::VegasGrid) = map(Base.OneTo, size(s))
function Base.CartesianIndices(s::VegasGrid)
    CartesianIndices(axes(s))
end
function Base.LinearIndices(s::VegasGrid)
    # TODO these will overflow quickly as dimension of problem
    # becomes large
    LinearIndices(axes(s))
end

function rand_index(rng, s::VegasGrid)
    li = rand(rng, LinearIndices(s))
    CartesianIndices(s)[li]
end

function Base.getindex(s::VegasGrid, index::CartesianIndex)
    lower = map(s.boundaries, Tuple(index)) do walls, i
        walls[i]
    end
    upper = map(s.boundaries, Tuple(index)) do walls, i
        walls[i+1]
    end
    Domain(SVector(lower), SVector(upper))
end

function draw_x_index_cell(rng, s::VegasGrid)
    i = rand_index(rng,s)
    cell = s[i]
    x = uniform(rng, cell)
    (position=x, index=i, cell=cell)
end

function draw(rng, iq::VegasGrid)
    s = draw_x_index_cell(rng, iq)
    wt = volume(s.cell) * length(iq)
    (position=s.position, weight = wt)
end

struct VegasHist{N,D,NT,T}
    grid    ::VegasGrid{N,D}
    normsums::NTuple{N,Vector{NT   }}
    counts  ::NTuple{N,Vector{Int64}}
    sum     ::Base.RefValue{T}
    sum2    ::Base.RefValue{T}
end

struct VegasVR{NFX,N,D,FX}
    grid::VegasGrid{N,D}
    values::Vector{FX}
    vars::Vector{FX}
    neval::Int64
end

function ftype(vr::VegasVR{NFX,N,D,FX}) where {NFX,N,D,FX}
    FX
end

function ntype(vr::VegasVR{NFX,N,D,FX}) where {NFX,N,D,FX}
    NFX
end

@qstruct CDF{P,V}(
    positions::Vector{P},
    values::Vector{V},
    isfinite::Bool,
    ) do
    @argcheck first(values) == 0.
    @argcheck issorted(values)
    @argcheck issorted(positions)
    @argcheck length(positions) == length(values)
    @argcheck length(positions) >= 2
end

function linterpol(x, (x_min, x_max), (y_min, y_max))
    @argcheck x_min <= x <= x_max
    if x == x_min
        y_min
    elseif x == x_max
        y_max
    else
        w1 = x_max - x
        w2 = x - x_min
        w = w1 + w2
        @assert w1 + w2 ≈ w
        y_min * w1/w + y_max * w2/w
    end
end

function quantiles(cdf::CDF, nwalls)
    cdf.isfinite || return cdf.positions
    ret = [first(cdf.positions)]
    qvals = range(first(cdf.values), stop=last(cdf.values), length=nwalls)
    i = 1
    for qval in qvals[2:end]
        while cdf.values[i] < qval
            i += 1
        end
        xs = cdf.values[i-1], cdf.values[i]
        ys = cdf.positions[i-1], cdf.positions[i]
        q = linterpol(qval, xs, ys)
        push!(ret, q)
    end
    @assert length(ret) == nwalls
    ret
end

function midpoint(dom::Domain)
    middle.(dom.lower, dom.upper)
end

function midpoint(grid::VegasGrid)
    midpoint(domain(grid))
end

function init_vegas_hist(f, vr::VegasVR{NFX}) where {NFX}
    T = ftype(vr)
    bincounts = map(vr.grid.boundaries) do walls
        length(walls) - 1
    end
    normsums = map(bincounts) do nbins
        fill(zero(NFX), nbins)
    end
    counts = map(bincounts) do nbins
        fill(zero(Int64), nbins)
    end
    sum  = Ref(zero(T))
    sum2 = Ref(zero(T))
    VegasHist(vr.grid, normsums, counts, sum, sum2)
end

function create_vegas_hist(f, vr, alg)
    create_vegas_hist(f, alg.rng, vr, alg)
end

function create_vegas_hist(f, rng::AbstractRNG, vr, alg)
    h = init_vegas_hist(f, vr)
    fill_vegas_hist!(f,h,rng; neval=ceil(Int,alg.neval/alg.niter))
end

function create_vegas_hist(f, p::ParallelRNG, vr, alg)
    rngs = p.rngs
    nthreads = length(rngs)
    h1 = init_vegas_hist(f, vr)
    hists = [deepcopy(h1) for _ in 1:nthreads]
    neval_i = ceil(Int, alg.neval / alg.niter / nthreads)
    @threads for i in 1:nthreads
        fill_vegas_hist!(f, hists[i], rngs[i], neval=neval_i)
    end
    fuseall(hists)
end

function fuse(h1::VegasHist, h2::VegasHist)
    @assert h1.grid == h2.grid
    VegasHist(
        h1.grid,
        h1.normsums .+ h2.normsums,
        h1.counts .+ h2.counts,
        Ref(h1.sum[] + h2.sum[]),
        Ref(h1.sum2[] + h2.sum2[]),
    )
end

function fuseall(hists::AbstractVector{<:VegasHist})::VegasHist
    reduce(fuse, hists)
end

function fill_vegas_hist!(f, h::VegasHist, rng::AbstractRNG; neval)
    iq = h.grid
    sum = h.sum[]
    sum2 = h.sum2[]
    for _ in 1:neval
        s = draw_x_index_cell(rng, iq)
        cell = s.cell
        fx = f(s.position)
        y  = fx * volume(cell) * length(iq)
        sum  += y
        sum2 += y.^2
        for axis in eachindex(iq.boundaries)
            wt = (cell.upper[axis] - cell.lower[axis]) * size(iq)[axis]
            ibin = s.index[axis]
            h.counts[axis][ibin]   += 1
            h.normsums[axis][ibin] += norm(fx)*wt
        end
    end
    h.sum[] = sum
    h.sum2[] = sum2
    h
end

function axisinds(grid::VegasGrid{N}) where {N}
    ntuple(identity, Val{N}())
end

function estimate_integral(h::VegasHist)
    count = sum(h.counts[1])::Int64
    mean, var = mean_var(h.sum[], h.sum2[], count)
    (value = mean, var=var, neval=count)
end

function update_vegasvr(vr::VegasVR, h::VegasHist, alg::Vegas)
    axes = axisinds(h.grid)
    boundaries = map(axes) do axis
        cdf = estimate_cdf(h, axis, alg)
        nq = length(cdf.positions)
        quantiles(cdf, nq)
    end
    grid = VegasGrid(boundaries)
    val, var, neval = estimate_integral(h)
    values = push!(copy(vr.values), val)
    vars   = push!(copy(vr.vars  ), var)
    VR = typeof(vr)
    VR(grid, values, vars, vr.neval + neval)
end

function estimate_pdf(h::VegasHist, axis::Int, ::NoDamping)
    if any(iszero, h.counts[axis])
        # TODO: we could guarantee, that there are evals in each bin
        msg = """Empty bin in vegas histogram.
        gridsize = $(size(h.grid))
        axis = $axis
        counts[$axis] = $(h.counts[axis])
        neval = $(sum(h.counts[1]))
        Consider increasing `neval` or decreasing `nbins`.
        """
        throw(ArgumentError(msg))
    end
    pdf = map(/, h.normsums[axis], h.counts[axis])
    s = sum(pdf)
    if iszero(s)
        isfinite = false
    else
        pdf ./= s
        isfinite = true
    end
    pdf, isfinite
end

function estimate_pdf(h::VegasHist, axis::Int, r::LepageDamping)
    pdf, isfinite = estimate_pdf(h, axis, NoDamping())
    isfinite || return (pdf, false)
    updf = map(pdf) do p
        ret = ((p - 1) / log(p))^r.alpha
        if iszero(p)
            zero(ret)
        elseif isone(p)
            one(ret)
        else
            ret
        end
    end
    pdf = normalize!(updf, 1)
    @assert !any(isnan.(pdf))
    (pdf, true)
end

function estimate_cdf(h, axis::Int, alg::Vegas)
    pdf, isfinite = estimate_pdf(h, axis, alg.regularization)
    cdf_vals = cumsum(pdf)
    z = zero(eltype(cdf_vals))
    pushfirst!(cdf_vals, z)
    walls = h.grid.boundaries[axis]
    CDF(walls, cdf_vals, isfinite)
end

function tuneonce(f, vr::VegasVR, alg::Vegas=Vegas())
    h = create_vegas_hist(f, vr, alg)
    update_vegasvr(vr, h, alg)
end

function default_grid_size(dom::Domain{N}, nbins) where {N}
    ntuple(_ -> nbins, Val(N))
end

function equidistant_grid(dom::Domain, size)

    bdries = map(Tuple(dom.lower), Tuple(dom.upper), size) do lo, hi, nbins
        nwalls = nbins + 1
        r = range(lo, stop=hi, length=nwalls)
        collect(r)
    end
    VegasGrid(bdries)
end

function initvr(f, dom::Domain{N,D}, alg::Vegas) where {N,D}
    size = default_grid_size(dom, alg.nbins)
    iq = equidistant_grid(dom, size)
    s = draw_x_index_cell(alg.rng, iq)
    x = midpoint(dom)
    y = f(x).^2 ./ 2
    n = norm(y)
    FX = typeof(y)
    NFX = typeof(n)
    neval = 0
    VegasVR{NFX,N,D,FX}(iq, FX[],FX[],neval)
end

function tune(f, vr::VegasVR, alg::Vegas)
    for i in 1:alg.niter
        vr = tuneonce(f, vr, alg)
    end
    vr
end

function tune(f, dom::Domain, alg)
    vr = initvr(f, dom, alg)
    tune(f, vr, alg)
end

# function div_pos(x, y)
#     ret = x / y
#     if iszero(y) & (!iszero(x))
#         typeof(ret)(Inf)
#     else
#         ret
#     end
# end
function _combine_scalar(val1, var1, val2, var2)
    p1 = 1/var1
    p2 = 1/var2
    wt1  = p1 / (p1 + p2)
    wt2 = p2 / (p1 + p2)
    val = wt1*val1 + wt2*val2
    var = (wt1^2)*var1 + (wt2^2)*var2
    T = typeof((val, var))
    ret = if iszero(var1) & iszero(var2)
        @check val1 ≈ val2
        (val1, var1)
    elseif iszero(var1)
        (val1, var1)
    elseif iszero(var2)
        (val2, var2)
    else
        (val, var)
    end
    T(ret)
end

function _combine_vec(val1, var1, val2, var2)
    pairs = _combine_scalar.(val1, var1, val2, var2)
    first.(pairs), last.(pairs)
end

function combine(nt1, nt2) # named tuple val var
    val1, var1 = nt1
    val2, var2 = nt2

    value, var = if nt1.value isa AbstractVector
        _combine_vec(val1,var1,val2,var2)
    else
        _combine_scalar(val1,var1,val2,var2)
    end
    (value=value, var=var)
end

function finish_integral(vr::VegasVR, alg::Vegas)
    index = (alg.ndrop+1):length(vr.values)
    values = vr.values[index]
    vars = vr.vars[index]
    pairs = map(values, vars) do val, var
        (value=val, var=var)
    end
    value, var = reduce(combine, pairs)
    (value=value, std=sqrt.(var), neval=vr.neval)
end

function integral_kernel(f, dom::Domain, alg::Vegas)
    vr = tune(f, dom, alg)
    finish_integral(vr, alg)
end

function integral_kernel(f, dom::VegasGrid, alg::MCVanilla)
    mc_kernel(f, alg.rng, dom, alg.neval)
end

function canonicalize(f, dom::VegasGrid, alg::MCVanilla)
    f, dom, alg
end


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

"""
    InputStream(
        stream::IO,
        location::SourceLocation,
        saved_char::String,
        saved_location::SourceLocation,
        tabulations::Int64,
        saved_token::Union{Token, Nothing},
    )

A high-level wrapper around a stream, used to parse scene files
This class implements a wrapper around a stream, with the following
additional capabilities:
- It tracks the line number and column number;
- It permits to "un-read" characters and tokens.

## Arguments

- `stream::IO` : stream to read from

- `location::SourceLocation` : location of the last char read

- `saved_char::String` : the last char read

- `saved_location::SourceLocation` : location where `saved_char` is in the file

- `tabulations::Int64`: number of space a tab command gives

- `saved_token::Union{Token, Nothing}` : the last token found

## Constructors

-    InputStream(
          stream::IO,
          file_name::String = "",
          tabulations::Int64 = 8
          ) = new(
               stream,
               SourceLocation(file_name, 1, 1),
               "",
               SourceLocation(file_name, 1, 1),
               tabulations,
               nothing
               )

-    InputStream(
          s::IO,
          l::SourceLocation,
          sc::String,
          sl::SourceLocation,
          t::Int64,
          st::Union{Token, Nothing}
          ) = new(s,l,sc,sl,t,st)

See also: [`SourceLocation`](@ref), [`Token`](@ref)
"""
mutable struct InputStream
     stream::IO
     location::SourceLocation
     saved_char::String
     saved_location::SourceLocation
     tabulations::Int64
     saved_token::Union{Token, Nothing}

     InputStream(
          stream::IO,
          file_name::String = "",
          tabulations::Int64 = 8
          ) = new(
               stream,
               SourceLocation(file_name, 1, 1),
               "",
               SourceLocation(file_name, 1, 1),
               tabulations,
               nothing
               )

     InputStream(
          s::IO,
          l::SourceLocation,
          sc::String,
          sl::SourceLocation,
          t::Int64,
          st::Union{Token, Nothing}
          ) = new(s,l,sc,sl,t,st)
end

##########################################################################################92

"""
    update_pos(inputstream::InputStream, ch::String)

Update `location` after having read `ch` from the stream.

See also: [`SourceLocation`](@ref)
"""
function update_pos(inputstream::InputStream, ch::String)
     if ch == ""
          nothing
     elseif ch == "\n"
          inputstream.location.line_num += 1
          inputstream.location.col_num = 1
     elseif ch == "\t"
          inputstream.location.col_num += inputstream.tabulations
     else
          inputstream.location.col_num += 1
     end
end

"""
     read_char(inputstream::InputStream) :: String

Read a new character from the stream.
Calls internally [`update_pos`](@ref).

See also: [`InputStream`](@ref), [`unread_char`](@ref)
"""
function read_char(inputstream::InputStream)
     if inputstream.saved_char ≠ ""
          ch = inputstream.saved_char
          inputstream.saved_char = ""
     elseif eof(inputstream.stream)
          ch = ""
     else
          ch = String([read(inputstream.stream, UInt8)])
     end

     inputstream.saved_location = copy(inputstream.location)  # shallow copy ?
     update_pos(inputstream, ch)

     return ch
end

"""
     unread_char(inputstream::InputStream, ch::String)

Push a character back to the stream.

See also: [`InputStream`](@ref), [`read_char`](@ref)
"""
function unread_char(inputstream::InputStream, ch::String)
     @assert inputstream.saved_char == ""
     inputstream.saved_char = ch
     inputstream.location = copy(inputstream.saved_location) # shallow copy ?
end

"""
     skip_whitespaces_and_comments(inputstream::InputStream)

Keep reading characters until a non-whitespace/non-comment character is found.
Calls internally [`read_char`](@ref) and [`unread_char`](@ref), and it's used
inside the main function [`read_token`](@ref).

See also: [`InputStream`](@ref)
"""
function skip_whitespaces_and_comments(inputstream::InputStream)
     ch = read_char(inputstream)
     while ( (ch in WHITESPACE) || (ch == "#"))

          if ch == "#"
               ch = read_char(inputstream)
               !(ch == "") || (return nothing)

               if ch =="="
                    while true
                         if read_char(inputstream) == "=" && read_char(inputstream)=="#"
                              break
                         end
                    end
               end

               # It's a comment! Keep reading until the end of the line
               #(include the case "", the end-of-file)
               while read_char(inputstream) ∉ ["\r", "\n", ""]
                    nothing
               end
          end

          ch = read_char(inputstream)
          !(ch == "") || (return nothing)
     end

     # Put the non-whitespace character back
     unread_char(inputstream, ch)
     nothing
end

"""
     parse_string_token(
          inputstream::InputStream,
          token_location::SourceLocation
          ) ::Token(::SourceLocation, ::StringToken)

Parse a string from the given input `inputstream` and return
that string inside a `Token(::SourceLocation, ::StringToken)` with the given
`token_location`, throwing `GrammarError` in case of exception.
Works calling [`read_char`](@ref), and it's used
inside the main function [`read_token`](@ref).

See also: [`InputStream`](@ref), [`SourceLocation`](@ref)
[`Token`](@ref), [`StringToken`](@ref), [`GrammarError`](@ref)
"""
function parse_string_token(inputstream::InputStream, token_location::SourceLocation)
     token = ""
     while true
          ch = read_char(inputstream)

          if ch == "\""
               break
          end
          if ch == ""
               throw(GrammarError(token_location, "unterminated string"))
          end

          token *= ch
     end

    return Token(token_location, StringToken(token))
end

"""
     parse_float_token(
          inputstream::InputStream,
          first_char::String,
          token_location::SourceLocation
          ) :: Token{SourceLocation, LiteralNumberToken}

Parse a float from the given input `inputstream` and return
that float inside a `Token(::SourceLocation, ::LiteralNumberToken)` with the given
`token_location`, throwing `GrammarError` in case of exception.
Works calling [`read_char`](@ref) and [`unread_char`](@ref), and it's used
inside the main function [`read_token`](@ref).

See also: [`InputStream`](@ref), [`SourceLocation`](@ref)
[`Token`](@ref), [`LiteralNumberToken`](@ref), [`GrammarError`](@ref)
"""
function parse_float_token(inputstream::InputStream, first_char::String, token_location::SourceLocation)
     token = first_char

     while true
          ch = read_char(inputstream)

          if !( isdigit(ch) || (ch == ".") || (ch ∈ ["e", "E"]) )
               unread_char(inputstream, ch)
               break
          end

          token *= ch
     end

     try
          value = parse(Float64, token)
          return Token(token_location, LiteralNumberToken(value))
     catch ValueError
          throw(
               GrammarError(
                    token_location,
                    """ "$(token)" is an invalid floating-point number"""
               )
          )
     end

end

"""
     parse_keyword_or_identifier_token(
          inputstream::InputStream,
          first_char::String,
          token_location::SourceLocation
          ) ::Union{
                    Token(::SourceLocation, ::KeywordToken),
                    Token(::SourceLocation, ::IdentifierToken)
                    }

Parse a keyword or an identifier from the given input `inputstream` and return
that keyword/identifier inside respectively a `Token(::SourceLocation, ::KeyworkdToken)`
or a `Token(::SourceLocation, ::IdentifierToken)` with the given `token_location`,
throwing `GrammarError` in case of exception.
Works calling [`read_char`](@ref) and [`unread_char`](@ref), and it's used
inside the main function [`read_token`](@ref).

See also: [`InputStream`](@ref), [`SourceLocation`](@ref)
[`Token`](@ref), [`LiteralNumberToken`](@ref), [`GrammarError`](@ref)
"""
function parse_keyword_or_identifier_token(
               inputstream::InputStream,
               first_char::String,
               token_location::SourceLocation
          )

     token = first_char

     while true
          ch = read_char(inputstream)

          if !isalnum(ch)
               unread_char(inputstream, ch)
               break
          end

          token *= ch
     end

     try
          # If it is a keyword, it must be listed in the KEYWORDS dictionary
          return Token(token_location, KeywordToken(KEYWORDS[token]))
     catch KeyError
          # If we got KeyError, it is not a keyword and thus it must be an identifier
          return Token(token_location, IdentifierToken(token))
     end
end

"""
     read_token(inputstream::InputStream) :: Token

Read one of the 6 basic tokens from the stream, raising `GrammarError` if a
lexical error is found.
Calls internally the following functions:
- [`skip_whitespaces_and_comments`](@ref)
- [`read_char`](@ref)
- [`isdecimal`](@ref)
- [`isalpha`](@ref)
- [`copy(::SourceLocation)`](@ref)
- [`parse_string_token`](@ref) for [`StringToken`](@ref)
- [`parse_float_token`](@ref) for [`LiteralNumberToken`](@ref)
- [`parse_keyword_or_identifier_token`](@ref) for [`KeywordToken`](@ref)
  and [`IdentifierToken`](@ref)

See also: [`InputStream`](@ref), [`Token`](@ref), [`SymbolToken`](@ref),
[`StopToken`](@ref), [`GrammarError`](@ref),
"""
function read_token(inputstream::InputStream)

     if !isnothing(inputstream.saved_token)
          result = inputstream.saved_token
          inputstream.saved_token = nothing
          return result
     end

     skip_whitespaces_and_comments(inputstream)

     # At this point we're sure that ch does *not* contain a whitespace character
     ch = read_char(inputstream)
     if ch == ""
          # No more characters in the file, so return a StopToken
          return Token(inputstream.location, StopToken())
     end

     # At this point we must check what kind of token begins with the "ch" character
     # (which has been put back in the stream with self.unread_char). First, we save
     # the position in the stream.
     token_location = copy(inputstream.location)  # shallow copy ?

     if ch ∈ SYMBOLS
          # One-character symbol, like '(' or ','
          return Token(token_location, SymbolToken(ch))
     elseif ch == "\""
          # A literal string (used for file names)
          return parse_string_token(inputstream, token_location)
     elseif ( isdecimal(ch) || (ch ∈ ["."]) )
          # A floating-point number
          return parse_float_token(inputstream, ch, token_location)
     elseif isalpha(ch)
          # Since it begins with an alphabetic character, it must either be
          # a keyword or a identifier
          return parse_keyword_or_identifier_token(inputstream, ch, token_location)
     else
          # We got some weird character, like '@` or `&`
          throw(GrammarError(inputstream.location, "Invalid character $(ch)"))
     end
end

"""
    unread_token(inputstream::InputStream, token::Token)

Make as if `token` were never read from `inputstream`.

See also: [`InputStream`](@ref), [`Token`](@ref), [`read_token`](@ref)
"""
function unread_token(inputstream::InputStream, token::Token)
    @assert isnothing(inputstream.saved_token) "$(inputstream.saved_token) ≠ nothing "
    inputstream.saved_token = copy(token)
end


text:
# generators.jl
#
# Core package functionality that implements algorithmic counterfactual.
module Generators

using ..Models
using ..Losses
using Flux
using LinearAlgebra

export AbstractGenerator, AbstractGradientBasedGenerator, GenericGenerator, GreedyGenerator, CounterfactualState,
    generate_perturbations, conditions_satisified, mutability_constraints

include("base.jl")

end


text:
export SeisAddNoise,
SeisLinearEvents,
SeisParabEvents,
SeisHypEvents

include("SeisAddNoise.jl")
include("SeisLinearEvents.jl")
include("SeisParabEvents.jl")
include("SeisHypEvents.jl")
include("Wavelets/Wavelets.jl")
include("Windows/Windows.jl")


text:
# indicator of the free cone

export IndFree

"""
**Indicator of the free cone**

    IndFree()

Returns the indicator function of the whole space, or "free cone", *i.e.*,
a function which is identically zero.
"""
struct IndFree <: ProximableFunction end

is_separable(f::IndFree) = true
is_convex(f::IndFree) = true
is_affine(f::IndFree) = true
is_cone(f::IndFree) = true
is_smooth(f::IndFree) = true
is_quadratic(f::IndFree) = true

const Zero = IndFree

function (f::IndFree)(x::AbstractArray{T}) where {R, T <: RealOrComplex{R}}
    return R(0)
end

function prox!(y::AbstractArray{T}, f::IndFree, x::AbstractArray{T}, args...) where {R, T <: RealOrComplex{R}}
    y .= x
    return R(0)
end

function gradient!(y::AbstractArray{T}, f::IndFree, x::AbstractArray{T}) where {R, T <: RealOrComplex{R}}
    y .= T(0)
    return R(0)
end

fun_name(f::IndFree) = "indicator of the free cone"
fun_dom(f::IndFree) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::IndFree) = "x ↦ 0"
fun_params(f::IndFree) = "none"

function prox_naive(f::IndFree, x::AbstractArray{R}, gamma=R(1)) where {R, T <: RealOrComplex{R}}
    return x, R(0)
end


text:
lines = readlines("input.txt")
heightmap = parse.(Int, reduce(hcat, collect.(lines)))
idx = CartesianIndices(heightmap)

o = CartesianIndex(1,1)
V = CartesianIndex(-1, 0), CartesianIndex(1, 0), CartesianIndex(0, -1), CartesianIndex(0, 1)
padded = fill(10, size(heightmap) .+ 2)
padded[idx .+ o] = heightmap[idx]

minimum = reduce(.&, [padded[idx .+ o .+ v] .> heightmap for v = V])

function DFS(yx_i, cond, s)
    cond_j = yx_j -> cond(yx_j) && (heightmap[yx_j] >= heightmap[yx_i]) && !(yx_j ∈ s)
    return reduce(vcat, [yx_i, [DFS(yx_i + v, cond, [yx_i, s...]) for v ∈ V if cond_j(yx_i + v)]...])
end

cond = yx -> padded[yx + o] < 9
basins = [Set(DFS(start, cond, [])) for start = findall(==(1), minimum)]
println(prod(sort(length.(basins), rev=true)[1:3]))


text:
abstract type Objective end

mutable struct TrackingObjective{Q,U,C,B} <: Objective
    q::Vector{Q}
    u::Vector{U}
    γ::Vector{C}
    b::Vector{B}
end

function TrackingObjective(model, env, H::Int;
    q = [Diagonal(zeros(SizedVector{model.nq})) for t = 1:H],
    u = [Diagonal(zeros(SizedVector{model.nu})) for t = 1:H],
    γ = [Diagonal(zeros(SizedVector{model.nc})) for t = 1:H],
    b = [Diagonal(zeros(SizedVector{model.nc * friction_dim(env)})) for t = 1:H])
    return TrackingObjective(q, u, γ, b)
end

mutable struct TrackingVelocityObjective{Q,V,U,C,B,VT} <: Objective
    q::Vector{Q}
    v::Vector{V}
    u::Vector{U}
    γ::Vector{C}
    b::Vector{B}
    v_target::Vector{VT}
    q_target::Vector{VT}
end

function TrackingVelocityObjective(model, env, H::Int;
    q = [Diagonal(zeros(SizedVector{model.nq})) for t = 1:H],
    v = [Diagonal(zeros(SizedVector{model.nq})) for t = 1:H],
    u = [Diagonal(zeros(SizedVector{model.nu})) for t = 1:H],
    γ = [Diagonal(zeros(SizedVector{model.nc})) for t = 1:H],
    b = [Diagonal(zeros(SizedVector{model.nc * friction_dim(env)})) for t = 1:H],
    v_target = [zeros(SizedVector{model.nq}) for t = 1:H])

    v_target = [SizedVector{model.nq}(v) for v in v_target]
    if v_target != [zeros(SizedVector{model.nq}) for t = 1:H]
        q_target = [zeros(SizedVector{model.nq})]
        for t = 1:H-1
            qt = q_target[end] + v_target[t]
            push!(q_target, qt)
        end
    else
        q_target = [zeros(SizedVector{model.nq}) for t = 1:H]
    end
    return TrackingVelocityObjective(q, v, u, γ, b, v_target, q_target)
end


text:
const ONNX2JULIA_TYPES = Dict(
    TensorProto_DataType.INT64 => Int64,
    TensorProto_DataType.INT32 => Int32,
    TensorProto_DataType.INT8 => Int8,
    TensorProto_DataType.DOUBLE => Float64,
    TensorProto_DataType.FLOAT => Float32,
    TensorProto_DataType.FLOAT16 => Float16,
)

const ONNX2JULIA_DATA_FIELDS = Dict(
    TensorProto_DataType.INT64 => :int64_data,
    TensorProto_DataType.INT32 => :int32_data,
    # TensorProto_DataType.INT8 => no special field
    TensorProto_DataType.DOUBLE => :double_data,
    TensorProto_DataType.FLOAT => :float_data,
    # TensorProto_DataType.FLOAT16 => no special field
)

"""
    array(p::TensorProto)

Return `p` as an `Array` of the correct type. Second argument can be used to change type of the returned array
"""
function array(p::TensorProto, wrap=Array)
    T = ONNX2JULIA_TYPES[p.data_type]
    fld = get(ONNX2JULIA_DATA_FIELDS, p.data_type, :raw_data)
    bytes = getproperty(p, fld)
    data = !isempty(bytes) ? reinterpret(T, bytes) : reinterpret(T, p.raw_data)
    # note that we don't permute dimensions here, only reshape the data
    # see "Row-/Columns-major" in test/readwrite.jl for an example
    # why we don't need permutedims here
    return reshape(wrap(data), reverse(p.dims)...)
end

Base.size(vip::ValueInfoProto) = size(vip._type)
Base.size(tp::TypeProto) = size(tp.tensor_type)
Base.size(tp::TensorProto) = tp.dims
Base.size(tp_t::TypeProto_Tensor) = hasproperty(tp_t, :shape) ? size(tp_t.shape) : missing
Base.size(tsp::TensorShapeProto) = size.(Tuple(reverse(tsp.dim)))
Base.size(tsp_d::TensorShapeProto_Dimension) = hasproperty(tsp_d, :dim_value) ? tsp_d.dim_value : missing

"""
    attribute(p::AttributeProto)

Return attribute in `p` as a name => value pair.
"""
function attribute(p::AttributeProto)
    # Copy paste from ONNX.jl
    if (p._type != 0)
        field = [:f, :i, :s, :t, :g, :floats, :ints, :strings, :tensors, :graphs][p._type]
        if field === :s
            return Symbol(p.name) => String(getproperty(p, field))
        elseif  field === :strings
            return Symbol(p.name) => String.(getproperty(p, field))
        end
        return Symbol(p.name) => getproperty(p, field)
    end
end

Base.Dict(pa::AbstractVector{AttributeProto}) = Dict(attribute(p) for p in pa)


text:
using Documenter, Distributions
import Random: AbstractRNG, rand!

makedocs(
    format = :html,
    sitename = "Distributions.jl",
    modules = [Distributions],
    pages = [
        "index.md",
        "starting.md",
        "types.md",
        "univariate.md",
        "truncate.md",
        "multivariate.md",
        "matrix.md",
        "mixture.md",
        "fit.md",
        "extends.md",
    ]
)

deploydocs(
    repo = "github.com/JuliaStats/Distributions.jl.git",
    target = "build",
    julia  = "1.0",
    deps = nothing,
    make = nothing
)


text:
function get_casT1!(T1::Array{Float64,2}, Ccas::Array{Float64,1}, dets::Array{Determinant,1}, ref::Determinant, frozen::Int, ndocc::Int)

    for id in eachindex(dets)

        @inbounds D = dets[id]
        αexc = αexcitation_level(ref, D)
        βexc = βexcitation_level(ref, D)

        if αexc == 1 && βexc == 0

            i, = αexclusive(ref, D)
            a, = αexclusive(D, ref)

            p = phase(ref, D)

            # i is absolute. Take out the frozen orbitals to match the T arrays.
            # a is abolute. Take out the occupied orbitals to get relative index.
            @inbounds T1[i-frozen,a-ndocc] = Ccas[id]*p
        end
    end
end

function get_casT2!(T1::Array{Float64,2}, T2::Array{Float64,4}, Ccas::Array{Float64,1}, dets::Array{Determinant,1}, ref::Determinant, frozen::Int, ndocc::Int)

    for id in eachindex(dets)

        @inbounds D = dets[id]
        αexc = αexcitation_level(ref, D)
        βexc = βexcitation_level(ref, D)

        if αexc == 1 && βexc == 1

            i, = αexclusive(ref, D)
            j, = βexclusive(ref, D)
            a, = αexclusive(D, ref)
            b, = βexclusive(D, ref)

            p = phase(ref, D)

            # i and j are absolute. Take out the frozen orbitals to match the T arrays.
            # a and b are abolute. Take out the occupied orbitals to get relative index.
            @inbounds T2[i-frozen,j-frozen,a-ndocc,b-ndocc] = Ccas[id]*p - T1[i-frozen, a-ndocc]*T1[j-frozen,b-ndocc]
        end
    end
end

function get_casT3(idx::NTuple{6,Int}, Coef::Float64, frozen::Int, ndocc::Int, T1::Array{Float64,2}, T2::Array{Float64,4}, ref::Determinant)

    i,j,k,a,b,c = idx

    # Create the corresponding determinant
    p = 1.0
    _p, _det = annihilate(ref, i, 'α')
    p = _p*p
    _p, _det = annihilate(_det, j,  'β')
    p = _p*p
    _p, _det = annihilate(_det, k,  'α')
    p = _p*p

    _p, _det = create(_det, c,  'α')
    p = _p*p
    _p, _det = create(_det, b,  'β')
    p = _p*p
    _p, _det = create(_det, a, 'α')
    p = _p*p

    T3 = Coef*p

    # Shift absolute index, for cleanliness
    i, j, k = (i,j,k) .- frozen
    a, b, c = (a,b,c) .- ndocc

    @inbounds begin
        T3 -=  T1[i,a]*T1[j,b]*T1[k,c]
        T3 +=  T1[i,c]*T1[j,b]*T1[k,a]
        T3 -=  T1[i,a]*T2[k,j,c,b]
        T3 +=  T1[k,a]*T2[i,j,c,b]
        T3 +=  T1[i,c]*T2[k,j,a,b]
        T3 -=  T1[k,c]*T2[i,j,a,b]
        T3 +=  T1[j,b]*T2[k,i,a,c]
        T3 -=  T1[j,b]*T2[i,k,a,c]
    end
    return T3
end

function find_casT3(ridx::NTuple{6,Int}, C3coef::Array{Float64,1}, C3dets::Array{Determinant,1}, ref::Determinant, frozen::Int, ndocc::Int, T1::Array{Float64,2}, T2::Array{Float64,4})

    i,j,k,a,b,c = ridx
    # Change indexes to absolute
    i, j, k = (i,j,k) .+ frozen
    a, b, c = (a,b,c) .+ ndocc

    # Create the corresponding determinant
    p = 1.0
    _p, _det = annihilate(ref, i, 'α')
    p = _p*p
    _p, _det = annihilate(_det, j,  'β')
    p = _p*p
    _p, _det = annihilate(_det, k,  'α')
    p = _p*p

    _p, _det = create(_det, c,  'α')
    p = _p*p
    _p, _det = create(_det, b,  'β')
    p = _p*p
    _p, _det = create(_det, a, 'α')
    p = _p*p

    search = findall(d->d==_det, C3dets)

    if length(search) == 0
        return 0.0
    end

    Coef = C3coef[search[1]]
    T3 = Coef*p

    # Shift absolute index, for cleanliness
    i, j, k = (i,j,k) .- frozen
    a, b, c = (a,b,c) .- ndocc

    @inbounds begin
        T3 -=  T1[i,a]*T1[j,b]*T1[k,c]
        T3 +=  T1[i,c]*T1[j,b]*T1[k,a]
        T3 -=  T1[i,a]*T2[k,j,c,b]
        T3 +=  T1[k,a]*T2[i,j,c,b]
        T3 +=  T1[i,c]*T2[k,j,a,b]
        T3 -=  T1[k,c]*T2[i,j,a,b]
        T3 +=  T1[j,b]*T2[k,i,a,c]
        T3 -=  T1[j,b]*T2[i,k,a,c]
    end
    return T3
end

function get_casT4αβ(idx::NTuple{8,Int}, Coef::Float64, C3coef::Array{Float64,1}, C3dets::Array{Determinant,1},
                   ref::Determinant, frozen::Int, ndocc::Int, T1::Array{Float64,2}, T2::Array{Float64,4})

    i,j,k,l,a,b,c,d = idx

    # Create the corresponding determinant
    p = 1.0

    _p, _det = annihilate(ref,  i, 'α')
    p = _p*p
    _p, _det = annihilate(_det, j, 'β')
    p = _p*p
    _p, _det = annihilate(_det, k,  'α')
    p = _p*p
    _p, _det = annihilate(_det, l,  'β')
    p = _p*p

    _p, _det = create(_det, d, 'β')
    p = _p*p
    _p, _det = create(_det, c, 'α')
    p = _p*p
    _p, _det = create(_det, b, 'β')
    p = _p*p
    _p, _det = create(_det, a, 'α')
    p = _p*p

    T4 =  p*Coef

    # Shift absolute index, for cleanliness
    i, j, k, l = (i,j,k,l) .- frozen
    a, b, c, d = (a,b,c,d) .- ndocc

    @inbounds begin
        T4 -= T1[j,b]*T1[i,a]*T1[k,c]*T1[l,d]
        T4 += T1[l,b]*T1[i,a]*T1[k,c]*T1[j,d]
        T4 += T1[j,b]*T1[k,a]*T1[i,c]*T1[l,d]
        T4 -= T1[l,b]*T1[k,a]*T1[i,c]*T1[j,d]

        T4 -= T1[j,b]*T1[i,a]*T2[k,l,c,d]
        T4 += T1[l,b]*T1[i,a]*T2[k,j,c,d]
        T4 += T1[k,c]*T1[i,a]*T2[l,j,b,d]
        T4 -= T1[k,c]*T1[i,a]*T2[j,l,b,d]
        T4 += T1[j,d]*T1[i,a]*T2[k,l,c,b]
        T4 -= T1[l,d]*T1[i,a]*T2[k,j,c,b]
        T4 += T1[j,b]*T1[k,a]*T2[i,l,c,d]
        T4 -= T1[l,b]*T1[k,a]*T2[i,j,c,d]
        T4 -= T1[i,c]*T1[k,a]*T2[l,j,b,d]
        T4 += T1[i,c]*T1[k,a]*T2[j,l,b,d]
        T4 -= T1[j,d]*T1[k,a]*T2[i,l,c,b]
        T4 += T1[l,d]*T1[k,a]*T2[i,j,c,b]
        T4 += T1[i,c]*T1[j,b]*T2[k,l,a,d]
        T4 -= T1[k,c]*T1[j,b]*T2[i,l,a,d]
        T4 += T1[l,d]*T1[j,b]*T2[k,i,a,c]
        T4 -= T1[l,d]*T1[j,b]*T2[i,k,a,c]
        T4 -= T1[i,c]*T1[l,b]*T2[k,j,a,d]
        T4 += T1[k,c]*T1[l,b]*T2[i,j,a,d]
        T4 -= T1[j,d]*T1[l,b]*T2[k,i,a,c]
        T4 += T1[j,d]*T1[l,b]*T2[i,k,a,c]
        T4 -= T1[j,d]*T1[i,c]*T2[k,l,a,b]
        T4 += T1[l,d]*T1[i,c]*T2[k,j,a,b]
        T4 += T1[j,d]*T1[k,c]*T2[i,l,a,b]
        T4 -= T1[l,d]*T1[k,c]*T2[i,j,a,b]

        T4 -= T2[k,l,c,d]*T2[i,j,a,b]
        T4 += T2[k,j,c,d]*T2[i,l,a,b]
        T4 += T2[i,l,c,d]*T2[k,j,a,b]
        T4 -= T2[i,j,c,d]*T2[k,l,a,b]
        T4 -= T2[l,j,b,d]*T2[k,i,a,c]
        T4 += T2[l,j,b,d]*T2[i,k,a,c]
        T4 += T2[j,l,b,d]*T2[k,i,a,c]
        T4 -= T2[j,l,b,d]*T2[i,k,a,c]
        T4 += T2[k,l,c,b]*T2[i,j,a,d]
        T4 -= T2[k,j,c,b]*T2[i,l,a,d]
        T4 -= T2[i,l,c,b]*T2[k,j,a,d]
        T4 += T2[i,j,c,b]*T2[k,l,a,d]
    end

    # Rename the find_T3 function just for cleanliness
    T3(ridx) = find_casT3(ridx, C3coef, C3dets, ref, frozen, ndocc, T1, T2)

    @inbounds begin
        T4 -= T1[i,a]*T3((j,k,l,b,c,d))
        T4 += T1[k,a]*T3((j,i,l,b,c,d))
        T4 -= T1[j,b]*T3((i,l,k,a,d,c))
        T4 += T1[l,b]*T3((i,j,k,a,d,c))
        T4 += T1[i,c]*T3((j,k,l,b,a,d))
        T4 -= T1[k,c]*T3((j,i,l,b,a,d))
        T4 += T1[j,d]*T3((i,l,k,a,b,c))
        T4 -= T1[l,d]*T3((i,j,k,a,b,c))
    end
    return T4
end

function get_casT4αα(idx::NTuple{8,Int}, Coef::Float64, C3coef::Array{Float64,1}, C3dets::Array{Determinant,1},
                    ref::Determinant, frozen::Int, ndocc::Int, T1::Array{Float64,2}, T2::Array{Float64,4})

    i,j,k,l,a,b,c,d = idx

    # Create the corresponding determinant
    p = 1.0

    _p, _det = annihilate(ref,  i, 'α')
    p = _p*p
    _p, _det = annihilate(_det, j,  'β')
    p = _p*p
    _p, _det = annihilate(_det, k,  'α')
    p = _p*p
    _p, _det = annihilate(_det, l,  'α')
    p = _p*p

    _p, _det = create(_det, d,  'α')
    p = _p*p
    _p, _det = create(_det, c,  'α')
    p = _p*p
    _p, _det = create(_det, b, 'β')
    p = _p*p
    _p, _det = create(_det, a, 'α')
    p = _p*p

    T4 =  p*Coef

    # Shift absolute index, for cleanliness
    i, j, k, l = (i,j,k,l) .- frozen
    a, b, c, d = (a,b,c,d) .- ndocc

    @inbounds begin
        T4 -= T1[j,b]*T1[i,a]*T1[k,c]*T1[l,d]
        T4 += T1[j,b]*T1[i,a]*T1[l,c]*T1[k,d]
        T4 += T1[j,b]*T1[k,a]*T1[i,c]*T1[l,d]
        T4 -= T1[j,b]*T1[k,a]*T1[l,c]*T1[i,d]
        T4 -= T1[j,b]*T1[l,a]*T1[i,c]*T1[k,d]
        T4 += T1[j,b]*T1[l,a]*T1[k,c]*T1[i,d]

        T4 += T1[j,b]*T1[i,a]*T2[l,k,c,d]
        T4 -= T1[j,b]*T1[i,a]*T2[k,l,c,d]
        T4 -= T1[k,c]*T1[i,a]*T2[l,j,d,b]
        T4 += T1[l,c]*T1[i,a]*T2[k,j,d,b]
        T4 += T1[k,d]*T1[i,a]*T2[l,j,c,b]
        T4 -= T1[l,d]*T1[i,a]*T2[k,j,c,b]
        T4 -= T1[j,b]*T1[k,a]*T2[l,i,c,d]
        T4 += T1[j,b]*T1[k,a]*T2[i,l,c,d]
        T4 += T1[i,c]*T1[k,a]*T2[l,j,d,b]
        T4 -= T1[l,c]*T1[k,a]*T2[i,j,d,b]
        T4 -= T1[i,d]*T1[k,a]*T2[l,j,c,b]
        T4 += T1[l,d]*T1[k,a]*T2[i,j,c,b]
        T4 += T1[j,b]*T1[l,a]*T2[k,i,c,d]
        T4 -= T1[j,b]*T1[l,a]*T2[i,k,c,d]
        T4 -= T1[i,c]*T1[l,a]*T2[k,j,d,b]
        T4 += T1[k,c]*T1[l,a]*T2[i,j,d,b]
        T4 += T1[i,d]*T1[l,a]*T2[k,j,c,b]
        T4 -= T1[k,d]*T1[l,a]*T2[i,j,c,b]
        T4 -= T1[i,c]*T1[j,b]*T2[l,k,a,d]
        T4 += T1[i,c]*T1[j,b]*T2[k,l,a,d]
        T4 += T1[k,c]*T1[j,b]*T2[l,i,a,d]
        T4 -= T1[k,c]*T1[j,b]*T2[i,l,a,d]
        T4 -= T1[l,c]*T1[j,b]*T2[k,i,a,d]
        T4 += T1[l,c]*T1[j,b]*T2[i,k,a,d]
        T4 += T1[i,d]*T1[j,b]*T2[l,k,a,c]
        T4 -= T1[i,d]*T1[j,b]*T2[k,l,a,c]
        T4 -= T1[k,d]*T1[j,b]*T2[l,i,a,c]
        T4 += T1[k,d]*T1[j,b]*T2[i,l,a,c]
        T4 += T1[l,d]*T1[j,b]*T2[k,i,a,c]
        T4 -= T1[l,d]*T1[j,b]*T2[i,k,a,c]
        T4 -= T1[k,d]*T1[i,c]*T2[l,j,a,b]
        T4 += T1[l,d]*T1[i,c]*T2[k,j,a,b]
        T4 += T1[i,d]*T1[k,c]*T2[l,j,a,b]
        T4 -= T1[l,d]*T1[k,c]*T2[i,j,a,b]
        T4 -= T1[i,d]*T1[l,c]*T2[k,j,a,b]
        T4 += T1[k,d]*T1[l,c]*T2[i,j,a,b]

        T4 += T2[l,k,c,d]*T2[i,j,a,b]
        T4 -= T2[k,l,c,d]*T2[i,j,a,b]
        T4 -= T2[l,i,c,d]*T2[k,j,a,b]
        T4 += T2[i,l,c,d]*T2[k,j,a,b]
        T4 += T2[k,i,c,d]*T2[l,j,a,b]
        T4 -= T2[i,k,c,d]*T2[l,j,a,b]
        T4 += T2[l,j,d,b]*T2[k,i,a,c]
        T4 -= T2[l,j,d,b]*T2[i,k,a,c]
        T4 -= T2[k,j,d,b]*T2[l,i,a,c]
        T4 += T2[k,j,d,b]*T2[i,l,a,c]
        T4 += T2[i,j,d,b]*T2[l,k,a,c]
        T4 -= T2[i,j,d,b]*T2[k,l,a,c]
        T4 -= T2[l,j,c,b]*T2[k,i,a,d]
        T4 += T2[l,j,c,b]*T2[i,k,a,d]
        T4 += T2[k,j,c,b]*T2[l,i,a,d]
        T4 -= T2[k,j,c,b]*T2[i,l,a,d]
        T4 -= T2[i,j,c,b]*T2[l,k,a,d]
        T4 += T2[i,j,c,b]*T2[k,l,a,d]
    end

    # Rename the find_T3 function just for cleanliness
    T3(ridx) = find_casT3(ridx, C3coef, C3dets, ref, frozen, ndocc, T1, T2)

    @inbounds begin
        T4 -= T1[i,a]*T3((k,j,l,c,b,d))
        T4 += T1[k,a]*T3((i,j,l,c,b,d))
        T4 -= T1[l,a]*T3((i,j,k,c,b,d))
        T4 += T1[j,b]*T3((k,i,l,a,c,d))
        T4 -= T1[j,b]*T3((i,k,l,a,c,d))
        T4 += T1[j,b]*T3((i,l,k,a,c,d))
        T4 += T1[i,c]*T3((k,j,l,a,b,d))
        T4 -= T1[k,c]*T3((i,j,l,a,b,d))
        T4 += T1[l,c]*T3((i,j,k,a,b,d))
        T4 -= T1[i,d]*T3((k,j,l,a,b,c))
        T4 += T1[k,d]*T3((i,j,l,a,b,c))
        T4 -= T1[l,d]*T3((i,j,k,a,b,c))
    end
    return T4
end

function ec_T3onT1!(ridx::NTuple{6,Int}, T3::Float64, ecT1::Array{Float64,2}, Voovv::Array{Float64,4})

    i,j,k,a,b,c = ridx
    @inbounds begin
        ecT1[j,b] += 0.25*T3*(Voovv[i,k,a,c] - Voovv[k,i,a,c])
        ecT1[i,a] += T3*(1.5*Voovv[j,k,b,c] - 0.5*Voovv[k,j,b,c])
        ecT1[j,a] += 0.25*T3*(Voovv[k,i,b,c] - Voovv[i,k,b,c])
    end
end

function ec_T3onT2!(ridx::NTuple{6,Int}, T3::Float64, ecT2::Array{Float64,4}, T1::Array{Float64,2}, fov::Array{Float64,2}, Voovv::Array{Float64,4}, Vovvv::Array{Float64,4}, Vooov::Array{Float64,4})

    i,j,k,a,b,c = ridx

    # Simple
    @inbounds begin
        ecT2[j,i,b,a] += fov[k,c]*T3
        ecT2[i,j,a,b] += fov[k,c]*T3

        # Slice
        ecT2[j,i,:,a] += T3*Vovvv[k,:,c,b]
        ecT2[:,j,a,b] += 0.5*T3*Vooov[k,i,:,c]
        ecT2[:,j,a,b] += -0.5*T3*Vooov[i,k,:,c]
        ecT2[:,i,b,a] -= T3*Vooov[j,k,:,c]
        ecT2[j,:,b,a] += 0.5*T3*Vooov[k,i,:,c]
        ecT2[i,:,a,b] -= T3*Vooov[j,k,:,c]
        ecT2[j,:,b,a] += -0.5*T3*Vooov[i,k,:,c]
        ecT2[j,i,b,:] += -0.5*T3*Vovvv[k,:,a,c]
        ecT2[i,j,a,:] += T3*Vovvv[k,:,c,b]
        ecT2[j,i,b,:] += 0.5*T3*Vovvv[k,:,c,a]
        ecT2[i,j,:,b] += 0.5*T3*Vovvv[k,:,c,a]
        ecT2[i,j,:,b] += -0.5*T3*Vovvv[k,:,a,c]

        # Contraction
        sI,sA = size(T1)
        for I = 1:sI
            for A = 1:sA
                ecT2[j,i,A,a] -= T1[I,A]*T3*Voovv[k,I,c,b]
                ecT2[j,i,b,A] += 0.5*T1[I,A]*T3*Voovv[k,I,a,c]
                ecT2[j,i,b,A] += -0.5*T1[I,A]*T3*Voovv[I,k,a,c]
                ecT2[i,j,a,A] -= T1[I,A]*T3*Voovv[I,k,b,c]
                ecT2[i,I,a,b] -= T1[I,A]*T3*Voovv[j,k,A,c]
                ecT2[j,I,b,a] += 0.5*T1[I,A]*T3*Voovv[k,i,A,c]
                ecT2[j,I,b,a] += -0.5*T1[I,A]*T3*Voovv[i,k,A,c]
                ecT2[I,j,a,b] += 0.5*T1[I,A]*T3*Voovv[k,i,A,c]
                ecT2[I,j,a,b] += -0.5*T1[I,A]*T3*Voovv[i,k,A,c]
                ecT2[I,i,b,a] -= T1[I,A]*T3*Voovv[j,k,A,c]
                ecT2[i,j,a,b] -= T1[I,A]*T3*Voovv[k,I,A,c]
                ecT2[i,j,a,b] += 2.0*T1[I,A]*T3*Voovv[I,k,A,c]
                ecT2[j,i,b,a] -= T1[I,A]*T3*Voovv[k,I,A,c]
                ecT2[j,i,b,a] += 2.0*T1[I,A]*T3*Voovv[I,k,A,c]
                ecT2[i,j,A,b] += -0.5*T1[I,A]*T3*Voovv[I,k,a,c]
                ecT2[i,j,A,b] += 0.5*T1[I,A]*T3*Voovv[k,I,a,c]
            end
        end
    end
end

function ec_T4αβonT2!(ridx::NTuple{8,Int}, T4::Float64, ecT2::Array{Float64,4}, Voovv::Array{Float64,4})

    i,j,k,l,a,b,c,d = ridx

    ecT2[i,j,a,b] += T4*Voovv[l,k,d,c]
end

function ec_T4ααonT2!(ridx::NTuple{8,Int}, T4::Float64, ecT2::Array{Float64,4}, Voovv::Array{Float64,4})

    i,j,k,l,a,b,c,d = ridx

    ec = 0.25*T4*(Voovv[k,l,c,d] - Voovv[l,k,c,d])
    ecT2[i,j,a,b] += ec
    ecT2[j,i,b,a] += ec
end

function cas_decomposition(ref::Determinant, Casdata::Dict{String,Array}, ndocc::Int, frozen::Int, fov::Array{Float64,2}, Voovv::Array{Float64,4}, Vovvv::Array{Float64,4}, Vooov::Array{Float64,4})

    # Get T1 and T2
    T1 = zeros(size(fov))
    T2 = zeros(size(Voovv))
    @output "→ Getting T1..."
    get_casT1!(T1, Casdata["C1"], Casdata["C1dets"], ref, frozen, ndocc)
    get_casT1!(T1, Casdata["C1"], Casdata["C1dets"], ref, frozen, ndocc)
    @output " Done.\n"
    @output "→ Getting T2..."
    get_casT2!(T1, T2, Casdata["C2"], Casdata["C2dets"], ref, frozen, ndocc)
    @output " Done.\n"

    # Initialize arrays
    ecT1 = zeros(size(fov))
    ecT2 = zeros(size(Voovv))

    @output "→ Computing external correction from T3..."
    # Compute External Correction from T3
    for n in eachindex(Casdata["C3"])
        D = Casdata["C3dets"][n]
        C = Casdata["C3"][n]

        if abs(C) < 10^-10
            continue
        end

        αexc = αexcitation_level(ref, D)
        βexc = βexcitation_level(ref, D)

        if αexc != 2 || βexc != 1
            continue
        end

        k,i = αexclusive(ref, D)
        j,  = βexclusive(ref,D)
        c,a = αexclusive(D, ref)
        b,  = βexclusive(D, ref)

        T3 = get_casT3((i,j,k,a,b,c), C, frozen, ndocc, T1, T2, ref)

        i, j, k = (i,j,k) .- frozen
        a, b, c = (a,b,c) .- ndocc

        # Include all permutations of T3 into ecT1
        ec_T3onT1!((i,j,k,a,b,c), T3, ecT1, Voovv)
        ec_T3onT1!((k,j,i,a,b,c),-T3, ecT1, Voovv)
        ec_T3onT1!((i,j,k,c,b,a),-T3, ecT1, Voovv)
        ec_T3onT1!((k,j,i,c,b,a), T3, ecT1, Voovv)

        # Include all permutations of T3 into ecT2
        ec_T3onT2!((i,j,k,a,b,c), T3, ecT2, T1, fov, Voovv, Vovvv, Vooov)
        ec_T3onT2!((k,j,i,a,b,c),-T3, ecT2, T1, fov, Voovv, Vovvv, Vooov)
        ec_T3onT2!((i,j,k,c,b,a),-T3, ecT2, T1, fov, Voovv, Vovvv, Vooov)
        ec_T3onT2!((k,j,i,c,b,a), T3, ecT2, T1, fov, Voovv, Vovvv, Vooov)
    end
    @output "Done\n"

    # Compute External Correction from T4
    @output "→ Computing external correction from T4..."
    for n in eachindex(Casdata["C4"])
        D = Casdata["C4dets"][n]
        C = Casdata["C4"][n]

        if abs(C) < 10^-10
            continue
        end

        αexc = αexcitation_level(ref, D)
        βexc = βexcitation_level(ref, D)

        if αexc == 2 && βexc == 2

            k,i  = αexclusive(ref, D)
            l,j  = βexclusive(ref,D)
            c,a  = αexclusive(D, ref)
            d,b  = βexclusive(D, ref)

            T4 = get_casT4αβ((i,j,k,l,a,b,c,d), C, Casdata["C3"], Casdata["C3dets"], ref, frozen, ndocc, T1, T2)

            i, j, k, l = (i,j,k,l) .- frozen
            a, b, c, d = (a,b,c,d) .- ndocc

            # Include all permutations of T4αβ into ecT2
            ec_T4αβonT2!((i,j,k,l,a,b,c,d),  T4, ecT2, Voovv)
            ec_T4αβonT2!((i,j,k,l,a,d,c,b), -T4, ecT2, Voovv)
            ec_T4αβonT2!((i,j,k,l,c,d,a,b),  T4, ecT2, Voovv)
            ec_T4αβonT2!((i,j,k,l,c,b,a,d), -T4, ecT2, Voovv)
            ec_T4αβonT2!((k,j,i,l,a,d,c,b),  T4, ecT2, Voovv)
            ec_T4αβonT2!((k,j,i,l,a,b,c,d), -T4, ecT2, Voovv)
            ec_T4αβonT2!((k,j,i,l,c,d,a,b), -T4, ecT2, Voovv)
            ec_T4αβonT2!((k,j,i,l,c,b,a,d),  T4, ecT2, Voovv)
            ec_T4αβonT2!((k,l,i,j,a,d,c,b), -T4, ecT2, Voovv)
            ec_T4αβonT2!((k,l,i,j,a,b,c,d),  T4, ecT2, Voovv)
            ec_T4αβonT2!((k,l,i,j,c,d,a,b),  T4, ecT2, Voovv)
            ec_T4αβonT2!((k,l,i,j,c,b,a,d), -T4, ecT2, Voovv)
            ec_T4αβonT2!((i,l,k,j,a,d,c,b),  T4, ecT2, Voovv)
            ec_T4αβonT2!((i,l,k,j,a,b,c,d), -T4, ecT2, Voovv)
            ec_T4αβonT2!((i,l,k,j,c,d,a,b), -T4, ecT2, Voovv)
            ec_T4αβonT2!((i,l,k,j,c,b,a,d),  T4, ecT2, Voovv)

        elseif αexc == 3 && βexc == 1

            l,k,i  = αexclusive(ref, D)
            j,     = βexclusive(ref,D)
            d,c,a  = αexclusive(D, ref)
            b,     = βexclusive(D, ref)

            T4 = get_casT4αα((i,j,k,l,a,b,c,d), C, Casdata["C3"], Casdata["C3dets"], ref, frozen, ndocc, T1, T2)

            i, j, k, l = (i,j,k,l) .- frozen
            a, b, c, d = (a,b,c,d) .- ndocc

            # Include all permutations of T4αβ into ecT2
            ec_T4ααonT2!((i,j,k,l,a,b,c,d),  T4, ecT2, Voovv)
            ec_T4ααonT2!((i,j,k,l,c,b,d,a),  T4, ecT2, Voovv)
            ec_T4ααonT2!((i,j,k,l,a,b,d,c), -T4, ecT2, Voovv)
            ec_T4ααonT2!((i,j,k,l,c,b,a,d), -T4, ecT2, Voovv)
            ec_T4ααonT2!((i,j,k,l,d,b,a,c),  T4, ecT2, Voovv)
            ec_T4ααonT2!((i,j,k,l,d,b,c,a), -T4, ecT2, Voovv)
            ec_T4ααonT2!((k,j,i,l,c,b,d,a), -T4, ecT2, Voovv)
            ec_T4ααonT2!((k,j,i,l,a,b,d,c),  T4, ecT2, Voovv)
            ec_T4ααonT2!((k,j,i,l,c,b,a,d),  T4, ecT2, Voovv)
            ec_T4ααonT2!((k,j,i,l,a,b,c,d), -T4, ecT2, Voovv)
            ec_T4ααonT2!((k,j,i,l,d,b,a,c), -T4, ecT2, Voovv)
            ec_T4ααonT2!((k,j,i,l,d,b,c,a),  T4, ecT2, Voovv)
            ec_T4ααonT2!((k,j,l,i,c,b,d,a),  T4, ecT2, Voovv)
            ec_T4ααonT2!((k,j,l,i,a,b,d,c), -T4, ecT2, Voovv)
            ec_T4ααonT2!((k,j,l,i,c,b,a,d), -T4, ecT2, Voovv)
            ec_T4ααonT2!((k,j,l,i,a,b,c,d),  T4, ecT2, Voovv)
            ec_T4ααonT2!((k,j,l,i,d,b,a,c),  T4, ecT2, Voovv)
            ec_T4ααonT2!((k,j,l,i,d,b,c,a), -T4, ecT2, Voovv)
            ec_T4ααonT2!((l,j,k,i,c,b,d,a), -T4, ecT2, Voovv)
            ec_T4ααonT2!((l,j,k,i,a,b,d,c),  T4, ecT2, Voovv)
            ec_T4ααonT2!((l,j,k,i,c,b,a,d),  T4, ecT2, Voovv)
            ec_T4ααonT2!((l,j,k,i,a,b,c,d), -T4, ecT2, Voovv)
            ec_T4ααonT2!((l,j,k,i,d,b,a,c), -T4, ecT2, Voovv)
            ec_T4ααonT2!((l,j,k,i,d,b,c,a),  T4, ecT2, Voovv)
            ec_T4ααonT2!((l,j,i,k,c,b,d,a),  T4, ecT2, Voovv)
            ec_T4ααonT2!((l,j,i,k,a,b,d,c), -T4, ecT2, Voovv)
            ec_T4ααonT2!((l,j,i,k,c,b,a,d), -T4, ecT2, Voovv)
            ec_T4ααonT2!((l,j,i,k,a,b,c,d),  T4, ecT2, Voovv)
            ec_T4ααonT2!((l,j,i,k,d,b,a,c),  T4, ecT2, Voovv)
            ec_T4ααonT2!((l,j,i,k,d,b,c,a), -T4, ecT2, Voovv)
            ec_T4ααonT2!((i,j,l,k,c,b,d,a), -T4, ecT2, Voovv)
            ec_T4ααonT2!((i,j,l,k,a,b,d,c),  T4, ecT2, Voovv)
            ec_T4ααonT2!((i,j,l,k,c,b,a,d),  T4, ecT2, Voovv)
            ec_T4ααonT2!((i,j,l,k,a,b,c,d), -T4, ecT2, Voovv)
            ec_T4ααonT2!((i,j,l,k,d,b,a,c), -T4, ecT2, Voovv)
            ec_T4ααonT2!((i,j,l,k,d,b,c,a),  T4, ecT2, Voovv)
        end
    end
    @output "Done\n"
    return T1, T2, ecT1, ecT2
end


text:
@inline function calc_A_fifth_order(x::VInt32, ::Val{N}) where N
    p = 31; r = 3; q = 31; A = Int32(1518500249); C = Int32(-1871437695);
    E = Int32(2956927708 >> N)
    rounding = one(x) << (q - N - 1)
    x² = (x * x) >> N
    (rounding + A + (x² * (C + (x² * E) >> r) >> N) >> (p - q)) >> (q - N)
end

@inline function calc_B_fifth_order(x::VInt32, ::Val{N}) where N
    p = 32; r = 4; q = 30; B = Int32(-1192627308); D = Int32(1960919730);
    F = Int32(-3760127719 >> N);
    rounding = one(x) << (q - N - 1)
    x² = (x * x) >> N
    (rounding + x * (B + (x² * (D + (x² * F) >> r) >> N) >> (p - q)) >> N) >> (q - N)
end

@inline function calc_A(x::VInt32, bits::B) where B <: Union{Val{11}, Val{12}}
    calc_A_fifth_order(x, bits)
end
@inline function calc_B(x::VInt32, bits::B) where B <: Union{Val{11}, Val{12}}
    calc_B_fifth_order(x, bits)
end


text:
"""
    TransformedKernel(k::Kernel, t::Transform)

Kernel derived from `k` for which inputs are transformed via a [`Transform`](@ref) `t`.

The preferred way to create kernels with input transformations is to use the composition
operator [`∘`](@ref) or its alias `compose` instead of `TransformedKernel` directly since
this allows optimized implementations for specific kernels and transformations.

See also: [`∘`](@ref)
"""
struct TransformedKernel{Tk<:Kernel,Tr<:Transform} <: Kernel
    kernel::Tk
    transform::Tr
end

@functor TransformedKernel

(k::TransformedKernel)(x, y) = k.kernel(k.transform(x), k.transform(y))

# Optimizations for scale transforms of simple kernels to save allocations:
# Instead of a multiplying every element of the inputs before evaluating the metric,
# we perform a scalar multiplcation of the distance of the original inputs, if possible.
function (k::TransformedKernel{<:SimpleKernel,<:ScaleTransform})(
    x::AbstractVector{<:Real}, y::AbstractVector{<:Real}
)
    return kappa(k.kernel, _scale(k.transform, metric(k.kernel), x, y))
end

function _scale(t::ScaleTransform, metric::Euclidean, x, y)
    return first(t.s) * evaluate(metric, x, y)
end
function _scale(t::ScaleTransform, metric::Union{SqEuclidean,DotProduct}, x, y)
    return first(t.s)^2 * evaluate(metric, x, y)
end
_scale(t::ScaleTransform, metric, x, y) = evaluate(metric, t(x), t(y))

"""
    kernel ∘ transform
    ∘(kernel, transform)
    compose(kernel, transform)

Compose a `kernel` with a transformation `transform` of its inputs.

The prefix forms support chains of multiple transformations:
`∘(kernel, transform1, transform2) = kernel ∘ transform1 ∘ transform2`.

# Definition

For inputs ``x, x'``, the transformed kernel ``\\widetilde{k}`` derived from kernel ``k`` by
input transformation ``t`` is defined as
```math
\\widetilde{k}(x, x'; k, t) = k\\big(t(x), t(x')\\big).
```

# Examples

```jldoctest
julia> (SqExponentialKernel() ∘ ScaleTransform(0.5))(0, 2) == exp(-0.5)
true

julia> ∘(ExponentialKernel(), ScaleTransform(2), ScaleTransform(0.5))(1, 2) == exp(-1)
true
```

See also: [`TransformedKernel`](@ref)
"""
Base.:∘(k::Kernel, t::Transform) = TransformedKernel(k, t)
Base.:∘(k::TransformedKernel, t::Transform) = TransformedKernel(k.kernel, k.transform ∘ t)

# Simplify kernels with identity transformation of the inputs
Base.:∘(k::Kernel, ::IdentityTransform) = k
Base.:∘(k::TransformedKernel, ::IdentityTransform) = k

Base.show(io::IO, κ::TransformedKernel) = printshifted(io, κ, 0)

function printshifted(io::IO, κ::TransformedKernel, shift::Int)
    printshifted(io, κ.kernel, shift)
    return print(io, "\n" * ("\t"^(shift + 1)) * "- $(κ.transform)")
end

# Kernel matrix operations

function kernelmatrix_diag!(K::AbstractVector, κ::TransformedKernel, x::AbstractVector)
    return kernelmatrix_diag!(K, κ.kernel, _map(κ.transform, x))
end

function kernelmatrix_diag!(
    K::AbstractVector, κ::TransformedKernel, x::AbstractVector, y::AbstractVector
)
    return kernelmatrix_diag!(K, κ.kernel, _map(κ.transform, x), _map(κ.transform, y))
end

function kernelmatrix!(K::AbstractMatrix, κ::TransformedKernel, x::AbstractVector)
    return kernelmatrix!(K, κ.kernel, _map(κ.transform, x))
end

function kernelmatrix!(
    K::AbstractMatrix, κ::TransformedKernel, x::AbstractVector, y::AbstractVector
)
    return kernelmatrix!(K, κ.kernel, _map(κ.transform, x), _map(κ.transform, y))
end

function kernelmatrix_diag(κ::TransformedKernel, x::AbstractVector)
    return kernelmatrix_diag(κ.kernel, _map(κ.transform, x))
end

function kernelmatrix_diag(κ::TransformedKernel, x::AbstractVector, y::AbstractVector)
    return kernelmatrix_diag(κ.kernel, _map(κ.transform, x), _map(κ.transform, y))
end

function kernelmatrix(κ::TransformedKernel, x::AbstractVector)
    return kernelmatrix(κ.kernel, _map(κ.transform, x))
end

function kernelmatrix(κ::TransformedKernel, x::AbstractVector, y::AbstractVector)
    return kernelmatrix(κ.kernel, _map(κ.transform, x), _map(κ.transform, y))
end


text:
# v0.6

function insertionsort!(A::Array{T}) where T <: Number
    for i in 1:length(A)-1
        value = A[i+1]
        j = i
        while j > 0 && A[j] > value
            A[j+1] = A[j]
            j -= 1
        end
        A[j+1] = value
    end
    return A
end

x = randn(5)
@show x insertionsort!(x)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# light source parameters
n = 10
νm = [456812, 456808,456811, 456802]
Em = ones(length(νm))
σ = 20.0
γRate = 1.0e8

# detector parameters
deadtime = 10.0 #nanoseconds
resolution = 0.01 #nanoseconds
jitter = 0.015 #nanoseconds
efficiency = 0.9
darkcounts = 1.0e-8 #GHz

# Detector object construction
@test_throws AssertionError Detector(deadtime,resolution,jitter,1.1,darkcounts)

# create detector and light source objects
source = LightSource(n,Em,νm,σ,γRate)
detect = Detector(deadtime,resolution,jitter,efficiency,darkcounts)

readout(20,source,detect)


text:
using DifferentialEquations
using PyPlot
using DynamicalSystems
using Interpolations
using DelimitedFiles
using Printf
using LsqFit

################### simulation code for the mackey glass equation ##################
calc_dP(P, delayed_P, beta, n, gamma) = beta * delayed_P / (1 + delayed_P^n) - gamma * P

function mackey_glass_rhs(P, history, params, t)
  beta, n, gamma, tau = params
  delayed_P = history(params, t - tau)
  return calc_dP(P, delayed_P, beta, n, gamma)
end

const algorithm = MethodOfSteps(AutoTsit5(Rosenbrock23()))
const solve_kwargs = (abstol=1e-6, reltol=1e-6, maxiters=10^10)
const default_P0 = 0.5
const default_chaotic_params = (beta = 2, n = 9.65, gamma = 1, tau = 2)
const default_farmer_params = (beta = 0.2, n = 10.0, gamma = 0.1, tau = 16.8)
const farmer_taus = [17, 23, 23.8, 30]

function simulate_mackey_glass(beta, n, gamma, tau; tmax=100.0, t0=0.0, dt=0.1, P0=default_P0, dense=false)
    initial_history(p, t) = P0
    params = (beta, n, gamma, tau)
    prob = DDEProblem(mackey_glass_rhs, P0, initial_history, (t0, tmax), params; constant_lags=[tau])
    if dense
        return solve(prob, algorithm; solve_kwargs..., dense=true)
    else
        return solve(prob, algorithm; solve_kwargs..., saveat=dt)
    end
end

####################### some simple plots #######################
function plot_evolution()
    beta, n, gamma, _ = default_chaotic_params
    figure()
    subplot(2,1,1)
    tau = 1
    sol = simulate_mackey_glass(beta, n, gamma, tau)
    plot(sol.t, sol.u)
    xlabel("t")
    ylabel("P(t)")
    title("\$\\tau = $tau\$")

    subplot(2,1,2)
    tau = 2
    sol = simulate_mackey_glass(default_chaotic_params...)
    plot(sol.t, sol.u)
    xlabel("t")
    ylabel("P(t)")
    title("\$\\tau = $tau\$")

    suptitle("\$\\beta = $beta, n = $n, \\gamma = $gamma\$")
    tight_layout()
    savefig("machey_glass_evolution.pdf")
end

function plot_production_feedback_function()
    figure()
    x = range(0, 2, length=400)
    cmap = PyPlot.cm_get_cmap("viridis")
    n_list = [3, 5, 7, 10, 15, 20]
    for (i, n) in enumerate(n_list)
        plot(x, @.(x / (1 + x^n)), color=cmap((i - 1)/(length(n_list) - 1)),
             label="n = $n")
    end
    xlabel("x")
    ylabel("\$ x/(1 + x^n) \$")
    legend(ncol=2, framealpha=1)
    title("\$\\beta = 1\$")
    savefig("mackey_glass_feedback_function.pdf")
end

################## attractor of the mackey glass equation ######################
function mackey_glass_attractor(beta, n, gamma, tau;
                     embedding_delay=tau, do_3d=false,
                     periods=50, steps_per_period=100, transit_steps=0,
                     lw=1.5, do_plot=true, add_title=true)
    sol = simulate_mackey_glass(beta, n, gamma, tau;
                                   tmax=periods*tau,
                                   dt=embedding_delay / steps_per_period)
    if do_plot
        if do_3d
            plot3D(sol.u[1 + 2*steps_per_period + transit_steps:end],
                   sol.u[1 + steps_per_period + transit_steps:end-steps_per_period],
                   sol.u[1 + transit_steps:end-2*steps_per_period],
                   lw=lw)
        else
            plot(sol.u[1 + steps_per_period + transit_steps:end],
                 sol.u[1 + transit_steps:end-steps_per_period],
                 lw=lw)
        end
        xlabel("\$P(t)\$")
        ylabel("\$P(t - \\tau_\\mathrm{embedding})\$")
        if do_3d
            zlabel("\$P(t - 2\\tau_\\mathrm{embedding})\$")
        end
        if add_title
            title("\$\\beta = $beta, n = $n, \\gamma = $gamma, \\tau = $tau\$")
        end
    end
    return sol
end

########################### bifurcations of the macky glass ###########################
# n_list = [7, 7.75, 8.5, 8.79, 9.65, 9.696, 9.7056, 9.7451, 10, 20] # numbers from the scholarpedia page
# beta_list = [1.5, 1.7, 1.76, 1.782, 1.8, 2.0, 2.1]
function plot_tau_bifurcation()
    figure()
    rows, cols = 3, 2
    tau_list = [1.3, 1.4, 1.65, 1.71, 1.72, 1.74]
    beta, n, gamma, _ = default_chaotic_params
    for (i, tau) in enumerate(tau_list)
        subplot(rows, cols, i)
        mackey_glass_attractor(beta, n, gamma, tau;
                  add_title=false, lw=0.1, transit_steps=100*500, steps_per_period=100, periods=1000)
        title("\$\\tau = $tau\$")
    end
    suptitle("\$ \\beta = $beta, n = $n, \\gamma = $gamma, \\tau_\\mathrm{embedding} = \\tau\$")
    tight_layout()
    savefig("mackey_glass_tau_period_doubling_bifurcation.pdf")
end

############################### poincare section of the macky glass eq #############################
include("psos_from_fn.jl")

const farmer_delay = 10

# from the FARMER paper
function mackey_glass_poincare_section(beta, n, gamma, tau;
        crossing_val=0.85, delay=farmer_delay, direction=-1,
        time=10000, transit=100, bisections=time)
    tmax = delay * time
    @time sol = simulate_mackey_glass(beta, n, gamma, tau; tmax=tmax, dense=true)
    f(t) = SVector(sol(t), sol(t - delay), sol(t - 2*delay))
    @time psos = poincare_section_from_function(f, transit*delay, tmax, bisections,
            SVector(crossing_val, 0, 0), SVector(sign(direction), 0, 0))
    return psos
end

fractal_dimension(psos) = generalized_dim(psos; q=0.0)

psos_filename(tau) = "mackey_glass_psos_tau=$tau.dat"

function compute_poincare_sections()
    for tau in farmer_taus
        psos = mackey_glass_poincare_section(
            default_farmer_params.beta, default_farmer_params.n,
            default_farmer_params.gamma, tau;
            time=1_000_000, transit=500)
        writedlm(psos_filename(tau), psos)
    end
end

const psos_dim_filename = "mackey_glass_psos_fractal_dimension.dat"

function process_poincare_sections()
    Ds = Float64[]
    fig, axes = subplots(2, 2)
    for (i, tau) in enumerate(farmer_taus)
        psos = readdlm(psos_filename(tau))
        axes[i].plot(psos[:, 2], psos[:, 3], ".", ms=0.1)
        axes[i].set_xlabel("P(t - $farmer_delay)")
        axes[i].set_ylabel("(P(t - $(2*farmer_delay))")
        D = fractal_dimension(standardize(Dataset(@view(psos[:, 2:3]))))
        push!(Ds, D)
        axes[i].set_title(@sprintf("\$\\tau\$ = %.2f, D = %.2f", tau, D))
    end
    tight_layout()
    savefig("mackey_glass_psos_plot.png")
    writedlm(psos_dim_filename, Ds)
end

########################### lyapunov exponents of the macky glass eq ############################
# advance history vector by tau + dt
function mackey_glass_method_of_steps_map!(next_history_vector, initial_history_vector, params, _t)
    beta, n, gamma, tau = params
    dt = tau / (length(initial_history_vector) - 1)
    interp = LinearInterpolation(range(0, tau, length=length(initial_history_vector)), initial_history_vector)
    initial_history(p, t) = interp(t + tau)
    problem = DDEProblem(mackey_glass_rhs, initial_history_vector[end], initial_history, (0.0, tau + dt), params; constant_lags=[tau])
    solution = solve(problem, algorithm; solve_kwargs..., saveat=dt)
    next_history_vector[1:end] = @view(solution.u[2:end])
    return nothing
end

function mackey_glass_lyapunov_spectrum(beta, n, gamma, tau; N=80, transit=50, L=50)
    dds = DiscreteDynamicalSystem(mackey_glass_method_of_steps_map!, repeat([default_P0], N), (beta, n, gamma, tau))
    return lyapunovspectrum(dds, L; Ttr=transit, show_progress=true)
end

function lyapunov_dimension(spectrum)
    if spectrum[1] <= 0.0
        return NaN
    else
        j = findfirst(spectrum .< 0.0)
        @assert j < length(spectrum)
        return j + sum(spectrum[1:j]) / abs(spectrum[j + 1])
    end
end

const lyapunov_filename = "mackey_glass_lyapunov_spectra.dat"

function compute_lyapunov_spectra()
    beta, n, gamma, _ = default_farmer_params
    N = 30
    specs = Matrix{Float64}(undef, N, length(farmer_taus))
    for (i, tau) in enumerate(farmer_taus)
        specs[:, i] = mackey_glass_lyapunov_spectrum(beta, n, gamma, tau; N=N, transit=50, L=30)
    end
    writedlm(lyapunov_filename, specs)
end

const lyapunov_dim_filename = "mackey_glass_lyapunov_dim.dat"

function process_lyapunov_spectra()
    specs = readdlm(lyapunov_filename)
    figure()
    Ds = Float64[]
    for (i, tau) in enumerate(farmer_taus)
        spectrum = specs[:, i]
        D = lyapunov_dimension(spectrum)
        push!(Ds, D)
        step(1:length(spectrum), spectrum, where="mid", label="\$\\tau = $tau\$")
    end
    xlabel("i")
    ylabel("\$\\lambda_i\$")
    beta, n, gamma, _ = default_farmer_params
    title("\$\\beta = $beta, n = $n, \\gamma = $gamma\$")
    legend()
    savefig("mackey_glass_lyapunov_spectra_plot.pdf")
    writedlm(lyapunov_dim_filename, Ds)
end

function plot_dimensions()
    psos_dims = readdlm(psos_dim_filename)
    lyapunov_dims = readdlm(lyapunov_dim_filename)
    figure()
    plot(farmer_taus, psos_dims .+ 1, "x", label="Dimension from Poincare Section")
    plot(farmer_taus, lyapunov_dims, "+", label="Dimension from Lyapunov Spectrum")
    xlabel(raw"Delay, $\tau$")
    ylabel("Dimension, D")
    legend()
    savefig("mackey_glass_dim_plot.pdf")
end

############################################# main ###############################################
function main()
    plot_production_feedback_function()
    plot_evolution()
    plot_tau_bifurcation()
    compute_poincare_sections()
    process_poincare_sections()
    compute_lyapunov_spectra()
    process_lyapunov_spectra()
    plot_dimensions()
end

function test_fractal_dim(i; tol=0.25, do_standardize=false, k=20, w=1, z=-1)
    psos = Dataset(@view(readdlm(psos_filename(farmer_taus[i]))[:, 2:3]))
    if do_standardize
        psos = standardize(psos)
    end
    bs = estimate_boxsizes(psos; k=k, w=w, z=z)
    H = genentropy.(Ref(psos), bs; q=0.0)
    x = - log.(bs)
    idx, slops = linear_regions(x, H; tol=tol)
    for i in 1:length(idx)-1
        r = idx[i] + 1 : idx[i+1]
        plot(x[r], H[r], "x")
    end
    lr = findmax([idx[i+1] - idx[i] for i=1:length(idx)-1])[2]
    r = idx[lr] + 1 : idx[lr + 1]
    fit = curve_fit((x, p) -> @.(p[1]*x + p[2]), x[r], H[r], [1.0, 1.0])
    slope = fit.param[1]
    slope_error = sqrt(estimate_covar(fit)[1, 1])
    plot(x[r], @.(fit.param[2] + slope*x[r]), "-k")
    return slope, slope_error
end

max_stars_repo_name: cosmo-jana/numerics-physics-stuff

max_issues_repo_name: cosmo-jana/numerics-physics-stuff

max_forks_repo_name: cosmo-jana/numerics-physics-stuff