text:
# Digital down conversion

function read_seg_samples_ddc(boardnum,numblocks,n::Integer,seg_len,window,v_offset=0f0,v_conv=0.350f0)
    # read from the card and downconvert + average segments, and return as an array of complex IQ points
    seg_IQ = read_seg_waveforms_ddc(boardnum,numblocks,n,seg_len,v_offset,v_conv)
    # return windowed segment averages
    return average_IQ_seg(seg_IQ, window)
end

function read_seg_waveforms_ddc(boardnum,numblocks,n::Integer,seg_len,v_offset=0f0,v_conv=0.350f0)
    # read from the card, downconvert segments, and return as an array of complex IQ points
    # n is the ratio of the sampling frequency to the IF, it must be an integer multiple of 4
    # first, trigger the board and download data
    setup_acquire(boardnum,numblocks)
    if iszero(first(get_averager(boardnum)))
    signal = get_volts_12(boardnum,numblocks,v_offset,v_conv)
    else # averaging is enabled
        signal = get_volts_16(boardnum,numblocks,v_offset,v_conv)
    end
    # calc the length of the actual signal, given the segment length
    sig_len = (length(signal)÷seg_len) * seg_len
    if n > 4 && (n÷4)*4 == n # if needed, decimate
        signal = decim_fir(signal, n÷4, sig_len)
        seg_len ÷= (n÷4)
        resize!(signal, (length(signal)÷2)*2) # resize to even length
    elseif n == 4 # resize vector to contain whole number of segments
        resize!(signal, sig_len)
    else
        throw(ArgumentError("Fs/IF ratio must be an integer multiple of 4 that is also ≥ 4"))
    end
    # downmix with efficient method (x2 decimation since signal is split into I & Q components)
    baseband = ddc4!(signal)
    # return IQ segments as complex array
    return reshape(reinterpret(Complex{Float32},baseband), seg_len÷2, :)
end

function decim_fir(signal::AbstractVector, rate::Integer, sig_len::Integer=length(signal))
    # create LPF FIR decimator object
    lpf = FIRDecimator(rate)
    # resize signal vector to right length and decimate with FIR filter
    req_zeros = inputlength(lpf, sig_len÷rate) - sig_len
    return decimate(lpf, resize_signal!(signal, sig_len, req_zeros))
end

function ddc4!(signal::Vector)
    # Downmix the signal to dc, assuming that the signal is sampled at 4x the
    # signal frequency. This method avoids computing sine & cosine values.
    # Returns an array with points [I[1],Q[1],I[2],Q[2],I[3],Q[3],...]
    @inbounds for i in 1:4:length(signal)
        # Multiply every other sample (I2 & Q2) by -1 to downconvert
        signal[i+2] = -signal[i+2]
        signal[i+3] = -signal[i+3]
    end
    return signal
end

function ddcn(signal::Vector{T},n::Integer) where {T<:Real}
	# Downconvert signal assuming that the sampling frequency is an integer
	# multiple of the IF (if n==4, use the ddc4! function)
	# Precompute sine and cos
	X = Matrix{T}(undef,2,n)
	for i in 1:n
		X[1,i] = cospi(T(2*(i-1))/n)
		X[2,i] = sinpi(T(2*(i-1))/n)
	end
	# Multiply every sample to compute I & Q
    D = Matrix{T}(undef,2,length(signal))
    @inbounds for i in 1:length(signal)
        D[1,i] = X[1,mod1(i,n)]*signal[i]
        D[2,i] = X[2,mod1(i,n)]*signal[i]
    end
	return D
end

function average_IQ_seg(seg_IQ::AbstractMatrix{T}) where {T<:Complex}
    # Average over segments (the first dimension)
    IQz = mean(seg_IQ,dims=1)
    # Return complex IQ value vector
    return vec(IQz)
end

function average_IQ_seg(seg_IQ::AbstractMatrix{T},window) where {T<:Complex}
    # Trim each segment as specified by "window"
    seg_len = size(seg_IQ,1)
    cut_i = floor(Int,window[1]/sum(window)*seg_len) + 1
    cut_f = ceil(Int,(window[1]+window[2])/sum(window)*seg_len)
    # Average each windowed segment (the first dimension)
    IQz = mean(view(seg_IQ,cut_i:cut_f,:),dims=1)
    # Return complex IQ values
    return vec(IQz)
end

function resize_signal!(signal::Vector{T}, sig_len::Integer, zero_pad::Integer=0) where T
	# Resize the signal to (sig_len+zero_pad). The indices past sig_len are zeroed-out
    tot_len = sig_len + zero_pad
    resize!(signal, tot_len)
    @simd for i in (sig_len+1) : tot_len
        @inbounds signal[i] = zero(T)
    end
    return signal
end

function resize_signal!(signal::AbstractVector{T}, sig_len::Integer, zero_pad::Integer=0) where T
	# Resize the signal to (sig_len+zero_pad). The indices past sig_len are zeroed-out
    return vcat(signal[1:sig_len], zeros(T,zero_pad))
end


text:
println("Compile on main node ...")
include("imports.jl")
using Distributed

rmprocs(workers())
addprocs(20)  # 20
println("Compile on workers ...")
@everywhere include("imports.jl")

# Simulations to run.
sims = dict_list(Dict(:skew => [true, false], :tdist => [true, false], :K => collect(1:5)))

# Run simulations in parallel.
@time res = pmap(run, sims, on_error=identity)

# Send results to S3
Util.s3sync(from="$(Info.resultsdir_compare)/$(simname)",
            to="$(Info.awsbucket_compare)/$(simname)",
            tags=`--exclude '*.nfs'`)

# Postproessing
include("postprocess.jl")

ext: jl


text:
module NMF_MC

include("../Misc.jl")
include("BackwardKernel.jl")
using .Misc, .BackwardKernel

using Distributions, SpecialFunctions

export Particle, smc_weight, particle_mcmc

mutable struct Particle{ℜ <: Real}
    S_K::Array{ℜ,1}
    S_IK::Array{ℜ,2}
    S_JK::Array{ℜ,2}
    function Particle(I::Int,J::Int,K::Int; γ::ℜ=0.1) where {ℜ<:Real}
        a::ℜ = I*J*γ
        return new{ℜ}(fill(a/K,K), fill(a/(K*I),I,K), fill(a/(K*J),J,K))
    end
    function Particle(S_K::Array{ℜ,1}, S_IK::Array{ℜ,2}, S_JK::Array{ℜ,2}) where {ℜ<:Real}
        return new{ℜ}(S_K, S_IK, S_JK)
    end
end

function smc_weight(X::Array{ℜ,2}, K::Ƶ, N::Ƶ=1; μ::ℜ=nanmean(X), γ::ℜ=0.1, adaptive::Bool=true) where {ℜ<:Real, Ƶ<:Int}
    I::Ƶ, J::Ƶ = size(X)
    T::Ƶ = Ƶ(sum(X))

    ESS::Array{ℜ} = zeros(ℜ,T)
    ESS_min::Array{ℜ} = adaptive ? fill(N/2.0,T) : fill(N+1.0,T)

    a::ℜ, b::ℜ = I*J*γ, γ/μ
    log_Zᵗ::ℜ = a*log(b) - (a+T)*log(b + 1.0) - sum(lgamma.(X .+ 1.0))

    P = [Particle(I,J,K; γ=γ) for n=1:N]
    P_temp = [Particle(I,J,K; γ=γ) for n=1:N]

    log_wᵗ::Array{ℜ}, log_Wᵗ::Array{ℜ} = fill(log_Zᵗ,N), fill(-log(N),N)
    Wᵗ::Array{ℜ}, cum_Wᵗ::Array{ℜ} = fill(1.0/N,N), zeros(ℜ,N)

    log_νᵗ::ℜ, log_q_k::Array{ℜ}, q_k::Array{ℜ}= 0.0, zeros(ℜ,K), zeros(ℜ,K)
    iᵗ::Ƶ, jᵗ::Ƶ, kᵗ::Ƶ, uᵗ::ℜ  = 0, 0, 0, 0.0

    for (t,(iᵗ, jᵗ)) in enumerate(EventQueue(X))
        for (n,p) in enumerate(P)
            log_q_k .= log.(p.S_JK[jᵗ,:]) .+ log.(p.S_IK[iᵗ,:]) .- log.(p.S_K)
            log_νᵗ = logsumexp(log_q_k)
            log_q_k .-= log_νᵗ
            q_k .= exp.(log_q_k)
            kᵗ = rand(Categorical(q_k))

            p.S_K[kᵗ] += 1.0
            p.S_JK[jᵗ,kᵗ] += 1.0
            p.S_IK[iᵗ,kᵗ] += 1.0

            log_wᵗ[n] += log_νᵗ
        end

        log_Zᵗ = logsumexp(log_wᵗ)
        log_Wᵗ .= log_wᵗ .- log_Zᵗ
        log_Zᵗ -= log(N)
        Wᵗ .= exp.(log_Wᵗ)

        ESS[t] = 1.0/sum(Wᵗ .* Wᵗ)

        if ESS[t] < ESS_min[t]
            for n=1:N
                P_temp[n].S_K .= P[n].S_K
                P_temp[n].S_JK .= P[n].S_JK
                P_temp[n].S_IK .= P[n].S_IK
            end

            cum_Wᵗ .= cumsum(Wᵗ)
            uᵗ = rand()/N
            p_id = 1

            for n=1:N # systematic resampling
                while cum_Wᵗ[p_id] < uᵗ
                    p_id += 1
                end
                P[n].S_K .= P_temp[p_id].S_K
                P[n].S_JK .= P_temp[p_id].S_JK
                P[n].S_IK .= P_temp[p_id].S_IK
                uᵗ += 1.0/N
            end
            log_wᵗ .= log_Zᵗ
        end
    end
    p_id = rand(Categorical(Wᵗ))
    P[p_id].S_K .-= a/K
    P[p_id].S_JK .-= a/(K*J)
    P[p_id].S_IK .-= a/(I*K)
    return log_Zᵗ, ESS, P[p_id]
end

function particle_mcmc(X::Array{ℜ,2}, K::Ƶ, N::Ƶ=1; μ::ℜ=nanmean(X), γ::ℜ=0.1, EPOCHS::Ƶ=1, adaptive::Bool=false) where {ℜ<:Real, Ƶ<:Int}
    I::Ƶ, J::Ƶ = size(X)
    a::ℜ, b::ℜ = I*J*γ, γ/μ

    log_PX::Array{ℜ}, stays::Ƶ = zeros(ℜ, EPOCHS), 0
    log_PX[1], ESS::Array{ℜ}, part::Particle{ℜ} = smc_weight(X, K, N; μ=μ, γ=γ, adaptive=adaptive)
    log_PX_prop::ℜ, ESS_prop::Array{ℜ}, part_prop::Particle{ℜ} = smc_weight(X, K, N; μ=μ, γ=γ, adaptive=adaptive)

    for eph=2:EPOCHS
        log_PX_prop, ESS_prop, part_prop = smc_weight(X, K, N; μ=μ, γ=γ, adaptive=adaptive)

        if log(rand()) < log_PX_prop - log_PX[eph-1]
            log_PX[eph] = log_PX_prop
            ESS .= ESS_prop
            part.S_K .= part_prop.S_K
            part.S_IK .= part_prop.S_IK
            part.S_JK .= part_prop.S_JK
        else
            log_PX[eph] = log_PX[eph-1]
            stays += 1
        end
    end
    log_Z_τ::ℜ = logsumexp(log_PX) - log(EPOCHS)
    return log_Z_τ, log_PX, 1.0 - stays./EPOCHS, part
end

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
#!/usr/bin/env julia

# # setup
# import Pkg
# Pkg.pkg"add Makie#master AbstractPlotting#master GLMakie#master
#         Observables#master VideoIO"
import Pkg
Pkg.activate(".")
using Makie, Dates, Observables, VideoIO

# load the video

f = VideoIO.openvideo("test.mp4")

# Determine the width and height of the Scene
pixelaspectratio = VideoIO.aspect_ratio(f)
h = f.height
w = round(typeof(h), f.width * pixelaspectratio)

# To flip or not to flip?
flipx = false
flipy = false

rsc() = Scene(;camera = campixel!, raw = true, backgroundcolor = :black)

# define the initial Scene
scene = Scene(resolution = (w, h), backgroundcolor = :black)

# define some Observables
time = Node(0.0)
buf = Node(read(f))

# plot the image.  We're plotting an Observable, so on its update it will update
# the image as well.
img = image!(scene, 1:h, 1:w, buf, show_axis = false, scale_plot = false)[end]
scene

# Rotate and scale the scene
Makie.rotate!(scene, -0.5π)
Makie.scale!(scene, flipx ? -1 : 1,  flipy ? -1 : 1)
# Update all aspects of the Scene.
update_limits!(scene)
update_cam!(scene)
update!(scene)

subbuf = buf[]
# when time is updated, lazy-load the video.
lift(time) do t
    seek(f, t)
    read!(f, subbuf)
    buf[] = subbuf
end
# update the scene to make sure
update!(scene)
# display it
scene

# set the time to some value - check all this works
time[] = 1.0

# print in minutes
function print_time(i)
    return string(Time(0) + Microsecond(round(Int, 1e6*i)))#Microsecond(round(Int, 1000000 * i)))
end
# create a slider to control time (a time machine⁉)
timeslider = slider!(rsc(), 0:1/f.framerate:10, valueprinter = print_time,
            textcolor = :white)
# When the value of the slider is changed, change the time.
lift(timeslider[end][:value]) do value
    time[] = value
end

# Go back and forward by 1 frame each
# backbutton = button!(rsc(), "<", textcolor = :white)
fwdbutton  = button!(rsc(), ">", textcolor = :white)
# tie the slider value to the buttons as well
# lift(backbutton[end][:clicks]) do clicks
    # timeslider[end][:value][] -= 1/f.framerate
# end
lift(fwdbutton[end][:clicks]) do clicks
    timeslider[end][:value][] += 1/f.framerate
end

hbox(vbox(timeslider, fwdbutton), scene)

# Lo and behold!  A minimal video "player"!

##################################################################################

# Play button not implemented yet.
#=
playing = Node(true)
playlabel = lift(playing) do p
    if p
        "→"
    else
        "||"
    end
end
playbutton = button!(Scene(t), playlabel, textcolor = :white)
lift(playbutton[end][:clicks]) do clicks
    playing[] = !playing[]
end
playing[] = true
playlabel[]
=#


text:
# # Notebook Template

# Use this when using the end-to-end bash scripts to run ClimateMachine.jl and VizCLIMA.jl. These scripts use Literate.jl to automatically convert this .jl script to a Jupyter notebook in a pre-specified location.

# ## Load modules
using Plots
using NCDatasets
#using Statistics: mean
#using DelimitedFiles
#using PrettyTables
#using PaddedViews
#using Dates
#using Printf

# ## Setup run-time enviromnent
ENV["GKSwstype"] = "100"

# ## Extract data

# Specify location of .nc files
CLIMA_NETCDF = "../netcdf/";

# Get all netcdf file names in the CLIMA_NETCDF directory and select the last one
fnames = filter(x -> occursin(".nc", x), readdir( CLIMA_NETCDF ) );
filename = "$CLIMA_NETCDF"fnames[end]

# Open file and print its info
ds = NCDataset(filename, "r");
println(ds)
close(ds)

# ## View/modify your notebook on a local machine
# on remote terminal (e.g. login2): jupyter notebook --no-browser --port=9999
# on local terminal: ssh -N -f -L 9999:localhost:9999 <yourusername>@login2.hpc.caltech.edu
# in local browser: localhost:9999
# see https://github.com/CliMA/ClimateMachine.jl/wiki/Visualization


text:
# Parameters -
# baseAssetPrice::Float64
# volatility::Float64
# timeToExercise::Float64
# numberOfLevels::Float64
# strikePrice::Float64
# riskFreeRate::Float64
# dividendRate::Float64

"""
    delta(assetSet::Set{PSAbstractAsset}, parameters::PSBinaryLatticeModel, underlyingAssetPrice::Float64;
        earlyExercise::Bool = false)::PSResult

Compute the change in the option price for a 1-dollar increase in the underlying stock price

# Arguments
- `assetSet::Set{PSAbstractAsset}`: A set containing the put and call models involved in this trade.
- `parameters::PSBinaryLatticeModel`: A PSBinaryLatticeModel object containing the parameters for the lattice
- `underlyingAssetPrice::Float64`: Underlying stock price at the time of purchase of the contract
- `earlyExercise::Bool = false`: Can this option contract be excercised early (true for American options)
"""
function delta(assetSet::Set{PSAbstractAsset}, parameters::PSBinaryLatticeModel, underlyingAssetPrice::Float64;
    earlyExercise::Bool = false)::PSResult

    # TODO: checks ...
    # check - asset price > 0
    if (underlyingAssetPrice <= zero(underlyingAssetPrice))
        return (PSResult(ArgumentError("Underlying asset price must be positive")))
    end

    # base - we are looking at how the options price changes with $1 increase in the underlying -
    result = option_contract_price(assetSet,parameters,underlyingAssetPrice; earlyExercise=earlyExercise)
    if (isa(result.value,Exception) == true)
        return result
    end
    results_tuple = result.value;
    base_option_price = first(results_tuple.cost_calculation_result.option_contract_price_array);

    # perturbed - what is the price if the underlying increases by $1 -
    result = option_contract_price(assetSet,parameters, (underlyingAssetPrice + 1.0); earlyExercise=earlyExercise)
    if (isa(result.value,Exception) == true)
        return result
    end
    results_tuple = result.value;
    perturbed_option_price = first(results_tuple.cost_calculation_result.option_contract_price_array);

    # compute delta -
    delta = perturbed_option_price - base_option_price

    # return -
    return PSResult(delta)

end

"""
    theta(assetSet::Set{PSAbstractAsset}, parameters::PSBinaryLatticeModel, underlyingAssetPrice::Float64;
        earlyExercise::Bool = false)::PSResult

Compute the change in the price of the option contract for a one day decrease in the number of days left until expiration

# Arguments
- `assetSet::Set{PSAbstractAsset}`: A set containing the put and call models involved in this trade.
- `parameters::PSBinaryLatticeModel`: A PSBinaryLatticeModel object containing the parameters for the lattice
- `underlyingAssetPrice::Float64`: Underlying stock price at the time of purchase of the contract
- `earlyExercise::Bool = false`: Can this option contract be excercised early (true for American options)
"""
function theta(assetSet::Set{PSAbstractAsset}, parameters::PSBinaryLatticeModel, underlyingAssetPrice::Float64;
    earlyExercise::Bool = false)::PSResult

    # TODO: checks ...
    # check - asset price > 0
    if (underlyingAssetPrice <= zero(underlyingAssetPrice))
        return (PSResult(ArgumentError("Underlying asset price must be positive")))
    end

    # base - we are looking at how the options price changes with $1 increase in the underlying -
    result = option_contract_price(assetSet,parameters,underlyingAssetPrice; earlyExercise=earlyExercise)
    if (isa(result.value,Exception) == true)
        return result
    end
    results_tuple = result.value;
    base_option_price = first(results_tuple.cost_calculation_result.option_contract_price_array);

    # create a new parameters w/one less day -
    perturbedParameters = deepcopy(parameters)
    perturbedParameters.timeToExercise = (parameters.timeToExercise - 1.0)

    # perturbed - we are looking at how the options price changes with $1 increase in the underlying -
    result = option_contract_price(assetSet, perturbedParameters, underlyingAssetPrice; earlyExercise=earlyExercise)
    if (isa(result.value,Exception) == true)
        return result
    end
    results_tuple = result.value;
    perturbed_option_price = first(results_tuple.cost_calculation_result.option_contract_price_array);

    # compute theta -
    theta = perturbed_option_price - base_option_price

    # return -
    return PSResult(theta)
end

"""
    gamma(assetSet::Set{PSAbstractAsset}, parameters::PSBinaryLatticeModel, underlyingAssetPrice::Float64;
        earlyExercise::Bool = false)::PSResult

Compute the rate of change of the delta parameter

# Arguments
- `assetSet::Set{PSAbstractAsset}`: A set containing the put and call models involved in this trade.
- `parameters::PSBinaryLatticeModel`: A PSBinaryLatticeModel object containing the parameters for the lattice
- `underlyingAssetPrice::Float64`: Underlying stock price at the time of purchase of the contract
- `earlyExercise::Bool = false`: Can this option contract be excercised early (true for American options)
"""
function gamma(assetSet::Set{PSAbstractAsset}, parameters::PSBinaryLatticeModel, underlyingAssetPrice::Float64;
    earlyExercise::Bool = false)::PSResult

    # TODO: checks -

    # check - asset price > 0
    if (underlyingAssetPrice <= zero(underlyingAssetPrice))
        return (PSResult(ArgumentError("Underlying asset price must be positive")))
    end

    # compute a base delta -
    result = delta(assetSet,parameters, underlyingAssetPrice; earlyExercise=earlyExercise)
    if (isa(result.value,Exception) == true)
        return result
    end
    base_delta = result.value;

    # compute an updated delta -
    result = delta(assetSet, parameters, (underlyingAssetPrice + 1.0); earlyExercise=earlyExercise)
    if (isa(result.value,Exception) == true)
        return result
    end
    perturbed_delta = result.value;

    # diff -
    gamma = (perturbed_delta - base_delta)

    # compute the difference and return -
    return PSResult(gamma)
end

"""

    vega(assetSet::Set{PSAbstractAsset}, parameters::PSBinaryLatticeModel, underlyingAssetPrice::Float64;
        earlyExercise::Bool = false)::PSResult

Compute the change in the price of an option contract for a 1-percent increase in the implied volatility

# Arguments
- `assetSet::Set{PSAbstractAsset}`: A set containing the put and call models involved in this trade.
- `parameters::PSBinaryLatticeModel`: A PSBinaryLatticeModel object containing the parameters for the lattice
- `underlyingAssetPrice::Float64`: Underlying stock price at the time of purchase of the contract
- `earlyExercise::Bool = false`: Can this option contract be excercised early (true for American options)
"""
function vega(assetSet::Set{PSAbstractAsset}, parameters::PSBinaryLatticeModel, underlyingAssetPrice::Float64;
    earlyExercise::Bool = false)::PSResult

    # TODO: checks ...
    if (underlyingAssetPrice <= zero(underlyingAssetPrice))
        return (PSResult(ArgumentError("Underlying asset price must be positive")))
    end

    # base - we are looking at how the options price changes with $1 increase in the underlying -
    result = option_contract_price(assetSet,parameters,underlyingAssetPrice; earlyExercise=earlyExercise)
    if (isa(result.value,Exception) == true)
        return result
    end
    results_tuple = result.value;
    base_option_price = first(results_tuple.cost_calculation_result.option_contract_price_array);

    # create a new parameters w/increased volatility -
    epsilon = 0.01
    baseVolatility = parameters.volatility
    perturbedParameters = deepcopy(parameters)
    perturbedParameters.volatility = baseVolatility*(1.0+epsilon)

    # perturbed - we are looking at how the options price changes with $1 increase in the underlying -
    result = option_contract_price(assetSet, perturbedParameters, underlyingAssetPrice; earlyExercise=earlyExercise)
    if (isa(result.value,Exception) == true)
        return result
    end
    results_tuple = result.value;
    perturbed_option_price = first(results_tuple.cost_calculation_result.option_contract_price_array);

    # compute theta -
    vega = perturbed_option_price - base_option_price

    # return -
    return PSResult(vega)
end

"""

    rho(assetSet::Set{PSAbstractAsset}, parameters::PSBinaryLatticeModel, underlyingAssetPrice::Float64;
        earlyExercise::Bool = false)::PSResult

Compute the change in the price of an option contract for a 1% increase in the risk free rate

# Arguments
- `assetSet::Set{PSAbstractAsset}`: A set containing the put and call models involved in this trade.
- `parameters::PSBinaryLatticeModel`: A PSBinaryLatticeModel object containing the parameters for the lattice
- `underlyingAssetPrice::Float64`: Underlying stock price at the time of purchase of the contract
- `earlyExercise::Bool = false`: Can this option contract be excercised early (true for American options)
"""
function rho(assetSet::Set{PSAbstractAsset}, parameters::PSBinaryLatticeModel, underlyingAssetPrice::Float64;
    earlyExercise::Bool = false)::PSResult

    # TODO: checks ...
    # check - asset price > 0
    if (underlyingAssetPrice <= zero(underlyingAssetPrice))
        return (PSResult(ArgumentError("Underlying asset price must be positive")))
    end

    # base - we are looking at how the options price changes with $1 increase in the underlying -
    result = option_contract_price(assetSet,parameters,underlyingAssetPrice; earlyExercise=earlyExercise)
    if (isa(result.value,Exception) == true)
        return result
    end
    results_tuple = result.value;
    base_option_price = first(results_tuple.cost_calculation_result.option_contract_price_array);

    # create a new parameters w/one less day -
    epsilon = 0.01
    baseRiskFreeRate = parameters.riskFreeRate
    perturbedParameters = deepcopy(parameters)
    perturbedParameters.riskFreeRate = baseRiskFreeRate*(1.0+epsilon)

    # perturbed - we are looking at how the options price changes with $1 increase in the underlying -
    result = option_contract_price(assetSet, perturbedParameters, underlyingAssetPrice; earlyExercise=earlyExercise)
    if (isa(result.value,Exception) == true)
        return result
    end
    results_tuple = result.value;
    perturbed_option_price = first(results_tuple.cost_calculation_result.option_contract_price_array);

    # compute -
    rho = perturbed_option_price - base_option_price

    # return -
    return PSResult(rho)
end


text:
T_FWHM = 0.2ps
T_T0 = T_FWHM/(2*log(1+sqrt(2)))
β₂ = -0.013ps^2/m
γnl = 0.000032W/m
α = 1.75e-5/m
coupling = 1.75e-5
P_in= 55.6e-3W
FSR = 0.0182THz
L = 0.0119m
detuning = 0.0012
comb_ikeda = ToyModel(;betacoeff=[0.0/m, 0.0ps/m, β₂],
                    nonlinearcoeff=γnl,
                    linearloss=α,
                    power_in=P_in,
                    pulsetime=T_FWHM,
                    length=L,
                    FSR=FSR,
                    detuning=detuning,
                    coupling=coupling)
prob_ikeda = build_problem(comb_ikeda, DynamicIkeda();
                           time_window=25.0ps,
                           time_pts=2^12)
sol_ikeda = solve(prob_ikeda)


text:
function debug_callback(msg_severity, msg_type, callback_data_ptr, user_data_ptr)
    callback_data_ptr == C_NULL && return UInt32(0)
    callback_data = unsafe_load(callback_data_ptr)
    msg = unsafe_string(callback_data.pMessage)

    # ignore messages about available device extensions
    if !startswith(msg, "Device Extension: VK")
        id_name = unsafe_string(callback_data.pMessageIdName)
        # out of date as of 1.2.190
        id_name == "VUID-VkShaderModuleCreateInfo-pCode-04147" && contains(msg, "SPV_EXT_physical_storage_buffer") && return UInt32(0)

        msg_type = @match msg_type begin
            &Vk.DEBUG_UTILS_MESSAGE_TYPE_GENERAL_BIT_EXT => "General"
            &Vk.DEBUG_UTILS_MESSAGE_TYPE_VALIDATION_BIT_EXT => "Validation"
            &Vk.DEBUG_UTILS_MESSAGE_TYPE_PERFORMANCE_BIT_EXT => "Performance"
            _ => error("Unknown message type $msg_type")
        end

        log = string("$msg_type ($id_name): $msg")

        # defer logging to avoid problems when the callback is run inside a finalizer
        @async @match msg_severity begin
            &Vk.DEBUG_UTILS_MESSAGE_SEVERITY_VERBOSE_BIT_EXT => @debug(log)
            &Vk.DEBUG_UTILS_MESSAGE_SEVERITY_INFO_BIT_EXT => @info(log)
            &Vk.DEBUG_UTILS_MESSAGE_SEVERITY_WARNING_BIT_EXT => @warn(log)
            &Vk.DEBUG_UTILS_MESSAGE_SEVERITY_ERROR_BIT_EXT => @error(log)
            _ => error("Unknown message severity $msg_severity")
        end
    end
    UInt32(0)
end

struct Instruction
    name::Symbol
    args::Vector{Any}
    kwargs::Vector{NamedTuple}
end

function Base.show(io::IO, inst::Instruction)
    printstyled(io, inst.name; color = :cyan)
    print(io, '(', join(inst.args, ", "))
    if !isempty(inst.kwargs)
        print(io, "; ", join(inst.kwargs, ", "))
    end
    print(io, ')')
end

struct SnoopCommandBuffer <: CommandBuffer
    records::Vector{Instruction}
end

SnoopCommandBuffer() = SnoopCommandBuffer([])

function Base.show(io::IO, ::MIME"text/plain", snoop::SnoopCommandBuffer)
    println(io, "SnoopCommandBuffer($(length(snoop.records)) commands recorded):")
    foreach(snoop.records) do record
        println(io, "  ", record)
    end
end

name(ex) = @match ex begin
    ::Symbol => ex
    :($x::$_) => x
end

macro snoopdef(ex)
    (args, kwargs) = @match ex begin
        :($cmd(command_buffer, $(args...); $(kwargs...))::$_) => (args, kwargs)
        :($cmd(command_buffer, $(args...))::$_) => (args, [])
    end
    @match ex begin
        :($cmd(command_buffer, $(_...))::$rtype) || :($cmd(command_buffer, $(_...); $(_...))::$rtype) => begin
            (rtype_annotation, last_stmt) = @match rtype begin
                :(ResultTypes.Result{Result,VulkanError}) => (:(Result{Vk.Result,Vk.VulkanError}), :(Vk.SUCCESS))
                :Cvoid => (:Cvoid, :nothing)
            end
            push_ex = :(push!(snoop.records, Instruction($(QuoteNode(cmd)), collect([$(name.(args)...)]), collect(kwargs))))
            esc(:(Vk.$cmd(snoop::SnoopCommandBuffer, $(args...); kwargs...)::$rtype_annotation = ($push_ex; $last_stmt)))
        end
        _ => error("Malformed expression: $ex")
    end
end

@snoopdef begin_command_buffer(command_buffer, info::Vk.CommandBufferBeginInfo)::ResultTypes.Result{Result,VulkanError}
@snoopdef end_command_buffer(command_buffer)::ResultTypes.Result{Result,VulkanError}
@snoopdef cmd_bind_pipeline(command_buffer, pipeline_bind_point::Vk.PipelineBindPoint, pipeline)::Cvoid
@snoopdef cmd_set_viewport(command_buffer, viewports::AbstractArray)::Cvoid
@snoopdef cmd_set_scissor(command_buffer, scissors::AbstractArray)::Cvoid
@snoopdef cmd_set_line_width(command_buffer, line_width::Real)::Cvoid
@snoopdef cmd_set_depth_bias(command_buffer, depth_bias_constant_factor::Real, depth_bias_clamp::Real, depth_bias_slope_factor::Real)::Cvoid
@snoopdef cmd_set_blend_constants(command_buffer, blend_constants::NTuple{4, Float32})::Cvoid
@snoopdef cmd_set_depth_bounds(command_buffer, min_depth_bounds::Real, max_depth_bounds::Real)::Cvoid
@snoopdef cmd_set_stencil_compare_mask(command_buffer, face_mask::Vk.StencilFaceFlag, compare_mask::Integer)::Cvoid
@snoopdef cmd_set_stencil_write_mask(command_buffer, face_mask::Vk.StencilFaceFlag, write_mask::Integer)::Cvoid
@snoopdef cmd_set_stencil_reference(command_buffer, face_mask::Vk.StencilFaceFlag, reference::Integer)::Cvoid
@snoopdef cmd_bind_descriptor_sets(command_buffer, pipeline_bind_point::Vk.PipelineBindPoint, layout, first_set::Integer, descriptor_sets::AbstractArray, dynamic_offsets::AbstractArray)::Cvoid
@snoopdef cmd_bind_index_buffer(command_buffer, buffer, offset::Integer, index_type::Vk.IndexType)::Cvoid
@snoopdef cmd_bind_vertex_buffers(command_buffer, buffers::AbstractArray, offsets::AbstractArray)::Cvoid
@snoopdef cmd_draw(command_buffer, vertex_count::Integer, instance_count::Integer, first_vertex::Integer, first_instance::Integer)::Cvoid
@snoopdef cmd_draw_indexed(command_buffer, index_count::Integer, instance_count::Integer, first_index::Integer, vertex_offset::Integer, first_instance::Integer)::Cvoid
@snoopdef cmd_draw_indirect(command_buffer, buffer, offset::Integer, draw_count::Integer, stride::Integer)::Cvoid
@snoopdef cmd_draw_indexed_indirect(command_buffer, buffer, offset::Integer, draw_count::Integer, stride::Integer)::Cvoid
@snoopdef cmd_dispatch(command_buffer, group_count_x::Integer, group_count_y::Integer, group_count_z::Integer)::Cvoid
@snoopdef cmd_dispatch_indirect(command_buffer, buffer, offset::Integer)::Cvoid
@snoopdef cmd_copy_buffer(command_buffer, src_buffer, dst_buffer, regions::AbstractArray)::Cvoid
@snoopdef cmd_copy_image(command_buffer, src_image, src_image_layout::Vk.ImageLayout, dst_image, dst_image_layout::Vk.ImageLayout, regions::AbstractArray)::Cvoid
@snoopdef cmd_blit_image(command_buffer, src_image, src_image_layout::Vk.ImageLayout, dst_image, dst_image_layout::Vk.ImageLayout, regions::AbstractArray, filter::Vk.Filter)::Cvoid
@snoopdef cmd_copy_buffer_to_image(command_buffer, src_buffer, dst_image, dst_image_layout::Vk.ImageLayout, regions::AbstractArray)::Cvoid
@snoopdef cmd_copy_image_to_buffer(command_buffer, src_image, src_image_layout::Vk.ImageLayout, dst_buffer, regions::AbstractArray)::Cvoid
@snoopdef cmd_update_buffer(command_buffer, dst_buffer, dst_offset::Integer, data_size::Integer, data::Ptr{Cvoid})::Cvoid
@snoopdef cmd_fill_buffer(command_buffer, dst_buffer, dst_offset::Integer, size::Integer, data::Integer)::Cvoid
@snoopdef cmd_clear_color_image(command_buffer, image, image_layout::Vk.ImageLayout, color::Vk.ClearColorValue, ranges::AbstractArray)::Cvoid
@snoopdef cmd_clear_depth_stencil_image(command_buffer, image, image_layout::Vk.ImageLayout, depth_stencil::Vk.ClearDepthStencilValue, ranges::AbstractArray)::Cvoid
@snoopdef cmd_clear_attachments(command_buffer, attachments::AbstractArray, rects::AbstractArray)::Cvoid
@snoopdef cmd_resolve_image(command_buffer, src_image, src_image_layout::Vk.ImageLayout, dst_image, dst_image_layout::Vk.ImageLayout, regions::AbstractArray)::Cvoid
@snoopdef cmd_set_event(command_buffer, event, stage_mask::Vk.PipelineStageFlag)::Cvoid
@snoopdef cmd_reset_event(command_buffer, event, stage_mask::Vk.PipelineStageFlag)::Cvoid
@snoopdef cmd_wait_events(command_buffer, events::AbstractArray, memory_barriers::AbstractArray, buffer_memory_barriers::AbstractArray, image_memory_barriers::AbstractArray; src_stage_mask = 0, dst_stage_mask = 0)::Cvoid
@snoopdef cmd_pipeline_barrier(command_buffer, src_stage_mask::Vk.PipelineStageFlag, dst_stage_mask::Vk.PipelineStageFlag, memory_barriers::AbstractArray, buffer_memory_barriers::AbstractArray, image_memory_barriers::AbstractArray; dependency_flags = 0)::Cvoid
@snoopdef cmd_begin_query(command_buffer, query_pool, query::Integer; flags = 0)::Cvoid
@snoopdef cmd_end_query(command_buffer, query_pool, query::Integer)::Cvoid
@snoopdef cmd_begin_conditional_rendering_ext(command_buffer, conditional_rendering_begin::Vk.ConditionalRenderingBeginInfoEXT)::Cvoid
@snoopdef cmd_end_conditional_rendering_ext(command_buffer)::Cvoid
@snoopdef cmd_reset_query_pool(command_buffer, query_pool, first_query::Integer, query_count::Integer)::Cvoid
@snoopdef cmd_write_timestamp(command_buffer, pipeline_stage::Vk.PipelineStageFlag, query_pool, query::Integer)::Cvoid
@snoopdef cmd_copy_query_pool_results(command_buffer, query_pool, first_query::Integer, query_count::Integer, dst_buffer, dst_offset::Integer, stride::Integer; flags = 0)::Cvoid
@snoopdef cmd_push_constants(command_buffer, layout, stage_flags::Vk.ShaderStageFlag, offset::Integer, size::Integer, values::Ptr{Cvoid})::Cvoid
@snoopdef cmd_begin_render_pass(command_buffer, render_pass_begin::Vk.RenderPassBeginInfo, contents::Vk.SubpassContents)::Cvoid
@snoopdef cmd_next_subpass(command_buffer, contents::Vk.SubpassContents)::Cvoid
@snoopdef cmd_end_render_pass(command_buffer)::Cvoid
@snoopdef cmd_execute_commands(command_buffer, command_buffers::AbstractArray)::Cvoid
@snoopdef cmd_debug_marker_begin_ext(command_buffer, marker_info::Vk.DebugMarkerMarkerInfoEXT)::Cvoid
@snoopdef cmd_debug_marker_end_ext(command_buffer)::Cvoid
@snoopdef cmd_debug_marker_insert_ext(command_buffer, marker_info::Vk.DebugMarkerMarkerInfoEXT)::Cvoid
@snoopdef cmd_execute_generated_commands_nv(command_buffer, is_preprocessed::Bool, generated_commands_info::Vk.GeneratedCommandsInfoNV)::Cvoid
@snoopdef cmd_preprocess_generated_commands_nv(command_buffer, generated_commands_info::Vk.GeneratedCommandsInfoNV)::Cvoid
@snoopdef cmd_bind_pipeline_shader_group_nv(command_buffer, pipeline_bind_point::Vk.PipelineBindPoint, pipeline, group_index::Integer)::Cvoid
@snoopdef cmd_push_descriptor_set_khr(command_buffer, pipeline_bind_point::Vk.PipelineBindPoint, layout, set::Integer, descriptor_writes::AbstractArray)::Cvoid
@snoopdef cmd_set_device_mask(command_buffer, device_mask::Integer)::Cvoid
@snoopdef cmd_dispatch_base(command_buffer, base_group_x::Integer, base_group_y::Integer, base_group_z::Integer, group_count_x::Integer, group_count_y::Integer, group_count_z::Integer)::Cvoid
@snoopdef cmd_push_descriptor_set_with_template_khr(command_buffer, descriptor_update_template, layout, set::Integer, data::Ptr{Cvoid})::Cvoid
@snoopdef cmd_set_viewport_w_scaling_nv(command_buffer, viewport_w_scalings::AbstractArray)::Cvoid
@snoopdef cmd_set_discard_rectangle_ext(command_buffer, discard_rectangles::AbstractArray)::Cvoid
@snoopdef cmd_set_sample_locations_ext(command_buffer, sample_locations_info::Vk.SampleLocationsInfoEXT)::Cvoid
@snoopdef cmd_begin_debug_utils_label_ext(command_buffer, label_info::Vk.DebugUtilsLabelEXT)::Cvoid
@snoopdef cmd_end_debug_utils_label_ext(command_buffer)::Cvoid
@snoopdef cmd_insert_debug_utils_label_ext(command_buffer, label_info::Vk.DebugUtilsLabelEXT)::Cvoid
@snoopdef cmd_write_buffer_marker_amd(command_buffer, pipeline_stage::Vk.PipelineStageFlag, dst_buffer, dst_offset::Integer, marker::Integer)::Cvoid
@snoopdef cmd_begin_render_pass_2(command_buffer, render_pass_begin::Vk.RenderPassBeginInfo, subpass_begin_info::Vk.SubpassBeginInfo)::Cvoid
@snoopdef cmd_next_subpass_2(command_buffer, subpass_begin_info::Vk.SubpassBeginInfo, subpass_end_info::Vk.SubpassEndInfo)::Cvoid
@snoopdef cmd_end_render_pass_2(command_buffer, subpass_end_info::Vk.SubpassEndInfo)::Cvoid
@snoopdef cmd_draw_indirect_count(command_buffer, buffer, offset::Integer, count_buffer, count_buffer_offset::Integer, max_draw_count::Integer, stride::Integer)::Cvoid
@snoopdef cmd_draw_indexed_indirect_count(command_buffer, buffer, offset::Integer, count_buffer, count_buffer_offset::Integer, max_draw_count::Integer, stride::Integer)::Cvoid
@snoopdef cmd_set_checkpoint_nv(command_buffer, checkpoint_marker::Ptr{Cvoid})::Cvoid
@snoopdef cmd_bind_transform_feedback_buffers_ext(command_buffer, buffers::AbstractArray, offsets::AbstractArray; sizes = C_NULL)::Cvoid
@snoopdef cmd_begin_transform_feedback_ext(command_buffer, counter_buffers::AbstractArray; counter_buffer_offsets = C_NULL)::Cvoid
@snoopdef cmd_end_transform_feedback_ext(command_buffer, counter_buffers::AbstractArray; counter_buffer_offsets = C_NULL)::Cvoid
@snoopdef cmd_begin_query_indexed_ext(command_buffer, query_pool, query::Integer, index::Integer; flags = 0)::Cvoid
@snoopdef cmd_end_query_indexed_ext(command_buffer, query_pool, query::Integer, index::Integer)::Cvoid
@snoopdef cmd_draw_indirect_byte_count_ext(command_buffer, instance_count::Integer, first_instance::Integer, counter_buffer, counter_buffer_offset::Integer, counter_offset::Integer, vertex_stride::Integer)::Cvoid
@snoopdef cmd_set_exclusive_scissor_nv(command_buffer, exclusive_scissors::AbstractArray)::Cvoid
@snoopdef cmd_bind_shading_rate_image_nv(command_buffer, image_layout::Vk.ImageLayout; image_view = C_NULL)::Cvoid
@snoopdef cmd_set_viewport_shading_rate_palette_nv(command_buffer, shading_rate_palettes::AbstractArray)::Cvoid
@snoopdef cmd_set_coarse_sample_order_nv(command_buffer, sample_order_type::Vk.CoarseSampleOrderTypeNV, custom_sample_orders::AbstractArray)::Cvoid
@snoopdef cmd_draw_mesh_tasks_nv(command_buffer, task_count::Integer, first_task::Integer)::Cvoid
@snoopdef cmd_draw_mesh_tasks_indirect_nv(command_buffer, buffer, offset::Integer, draw_count::Integer, stride::Integer)::Cvoid
@snoopdef cmd_draw_mesh_tasks_indirect_count_nv(command_buffer, buffer, offset::Integer, count_buffer, count_buffer_offset::Integer, max_draw_count::Integer, stride::Integer)::Cvoid
@snoopdef cmd_copy_acceleration_structure_nv(command_buffer, dst, src, mode::Vk.CopyAccelerationStructureModeKHR)::Cvoid
@snoopdef cmd_copy_acceleration_structure_khr(command_buffer, info::Vk.CopyAccelerationStructureInfoKHR)::Cvoid
@snoopdef cmd_copy_acceleration_structure_to_memory_khr(command_buffer, info::Vk.CopyAccelerationStructureToMemoryInfoKHR)::Cvoid
@snoopdef cmd_copy_memory_to_acceleration_structure_khr(command_buffer, info::Vk.CopyMemoryToAccelerationStructureInfoKHR)::Cvoid
@snoopdef cmd_write_acceleration_structures_properties_khr(command_buffer, acceleration_structures::AbstractArray, query_type::Vk.QueryType, query_pool, first_query::Integer)::Cvoid
@snoopdef cmd_write_acceleration_structures_properties_nv(command_buffer, acceleration_structures::AbstractArray, query_type::Vk.QueryType, query_pool, first_query::Integer)::Cvoid
@snoopdef cmd_build_acceleration_structure_nv(command_buffer, info::Vk.AccelerationStructureInfoNV, instance_offset::Integer, update::Bool, dst, scratch, scratch_offset::Integer; instance_data = C_NULL, src = C_NULL)::Cvoid
@snoopdef cmd_trace_rays_khr(command_buffer, raygen_shader_binding_table::Vk.StridedDeviceAddressRegionKHR, miss_shader_binding_table::Vk.StridedDeviceAddressRegionKHR, hit_shader_binding_table::Vk.StridedDeviceAddressRegionKHR, callable_shader_binding_table::Vk.StridedDeviceAddressRegionKHR, width::Integer, height::Integer, depth::Integer)::Cvoid
@snoopdef cmd_trace_rays_nv(command_buffer, raygen_shader_binding_table_buffer, raygen_shader_binding_offset::Integer, miss_shader_binding_offset::Integer, miss_shader_binding_stride::Integer, hit_shader_binding_offset::Integer, hit_shader_binding_stride::Integer, callable_shader_binding_offset::Integer, callable_shader_binding_stride::Integer, width::Integer, height::Integer, depth::Integer; miss_shader_binding_table_buffer = C_NULL, hit_shader_binding_table_buffer = C_NULL, callable_shader_binding_table_buffer = C_NULL)::Cvoid
@snoopdef cmd_trace_rays_indirect_khr(command_buffer, raygen_shader_binding_table::Vk.StridedDeviceAddressRegionKHR, miss_shader_binding_table::Vk.StridedDeviceAddressRegionKHR, hit_shader_binding_table::Vk.StridedDeviceAddressRegionKHR, callable_shader_binding_table::Vk.StridedDeviceAddressRegionKHR, indirect_device_address::Integer)::Cvoid
@snoopdef cmd_set_ray_tracing_pipeline_stack_size_khr(command_buffer, pipeline_stack_size::Integer)::Cvoid
@snoopdef cmd_set_performance_marker_intel(command_buffer, marker_info::Vk.PerformanceMarkerInfoINTEL)::ResultTypes.Result{Result, VulkanError}
@snoopdef cmd_set_performance_stream_marker_intel(command_buffer, marker_info::Vk.PerformanceStreamMarkerInfoINTEL)::ResultTypes.Result{Result, VulkanError}
@snoopdef cmd_set_performance_override_intel(command_buffer, override_info::Vk.PerformanceOverrideInfoINTEL)::ResultTypes.Result{Result, VulkanError}
@snoopdef cmd_set_line_stipple_ext(command_buffer, line_stipple_factor::Integer, line_stipple_pattern::Integer)::Cvoid
@snoopdef cmd_build_acceleration_structures_khr(command_buffer, infos::AbstractArray, build_range_infos::AbstractArray)::Cvoid
@snoopdef cmd_build_acceleration_structures_indirect_khr(command_buffer, infos::AbstractArray, indirect_device_addresses::AbstractArray, indirect_strides::AbstractArray, max_primitive_counts::AbstractArray)::Cvoid
@snoopdef cmd_set_cull_mode_ext(command_buffer; cull_mode = 0)::Cvoid
@snoopdef cmd_set_front_face_ext(command_buffer, front_face::Vk.FrontFace)::Cvoid
@snoopdef cmd_set_primitive_topology_ext(command_buffer, primitive_topology::Vk.PrimitiveTopology)::Cvoid
@snoopdef cmd_set_viewport_with_count_ext(command_buffer, viewports::AbstractArray)::Cvoid
@snoopdef cmd_set_scissor_with_count_ext(command_buffer, scissors::AbstractArray)::Cvoid
@snoopdef cmd_bind_vertex_buffers_2_ext(command_buffer, buffers::AbstractArray, offsets::AbstractArray; sizes = C_NULL, strides = C_NULL)::Cvoid
@snoopdef cmd_set_depth_test_enable_ext(command_buffer, depth_test_enable::Bool)::Cvoid
@snoopdef cmd_set_depth_write_enable_ext(command_buffer, depth_write_enable::Bool)::Cvoid
@snoopdef cmd_set_depth_compare_op_ext(command_buffer, depth_compare_op::Vk.CompareOp)::Cvoid
@snoopdef cmd_set_depth_bounds_test_enable_ext(command_buffer, depth_bounds_test_enable::Bool)::Cvoid
@snoopdef cmd_set_stencil_test_enable_ext(command_buffer, stencil_test_enable::Bool)::Cvoid
@snoopdef cmd_set_stencil_op_ext(command_buffer, face_mask::Vk.StencilFaceFlag, fail_op::Vk.StencilOp, pass_op::Vk.StencilOp, depth_fail_op::Vk.StencilOp, compare_op::Vk.CompareOp)::Cvoid
@snoopdef cmd_set_patch_control_points_ext(command_buffer, patch_control_points::Integer)::Cvoid
@snoopdef cmd_set_rasterizer_discard_enable_ext(command_buffer, rasterizer_discard_enable::Bool)::Cvoid
@snoopdef cmd_set_depth_bias_enable_ext(command_buffer, depth_bias_enable::Bool)::Cvoid
@snoopdef cmd_set_logic_op_ext(command_buffer, logic_op::Vk.LogicOp)::Cvoid
@snoopdef cmd_set_primitive_restart_enable_ext(command_buffer, primitive_restart_enable::Bool)::Cvoid
@snoopdef cmd_copy_buffer_2_khr(command_buffer, copy_buffer_info::Vk.CopyBufferInfo2KHR)::Cvoid
@snoopdef cmd_copy_image_2_khr(command_buffer, copy_image_info::Vk.CopyImageInfo2KHR)::Cvoid
@snoopdef cmd_blit_image_2_khr(command_buffer, blit_image_info::Vk.BlitImageInfo2KHR)::Cvoid
@snoopdef cmd_copy_buffer_to_image_2_khr(command_buffer, copy_buffer_to_image_info::Vk.CopyBufferToImageInfo2KHR)::Cvoid
@snoopdef cmd_copy_image_to_buffer_2_khr(command_buffer, copy_image_to_buffer_info::Vk.CopyImageToBufferInfo2KHR)::Cvoid
@snoopdef cmd_resolve_image_2_khr(command_buffer, resolve_image_info::Vk.ResolveImageInfo2KHR)::Cvoid
@snoopdef cmd_set_fragment_shading_rate_khr(command_buffer, fragment_size::Vk.Extent2D, combiner_ops::NTuple{2, Vk.FragmentShadingRateCombinerOpKHR})::Cvoid
@snoopdef cmd_set_fragment_shading_rate_enum_nv(command_buffer, shading_rate::Vk.FragmentShadingRateNV, combiner_ops::NTuple{2, Vk.FragmentShadingRateCombinerOpKHR})::Cvoid
@snoopdef cmd_set_vertex_input_ext(command_buffer, vertex_binding_descriptions::AbstractArray, vertex_attribute_descriptions::AbstractArray)::Cvoid
@snoopdef cmd_set_color_write_enable_ext(command_buffer, color_write_enables::AbstractArray)::Cvoid
@snoopdef cmd_set_event_2_khr(command_buffer, event, dependency_info::Vk.DependencyInfoKHR)::Cvoid
@snoopdef cmd_reset_event_2_khr(command_buffer, event, stage_mask::Integer)::Cvoid
@snoopdef cmd_wait_events_2_khr(command_buffer, events::AbstractArray, dependency_infos::AbstractArray)::Cvoid
@snoopdef cmd_pipeline_barrier_2_khr(command_buffer, dependency_info::Vk.DependencyInfoKHR)::Cvoid
@snoopdef cmd_write_timestamp_2_khr(command_buffer, stage::Integer, query_pool, query::Integer)::Cvoid
@snoopdef cmd_write_buffer_marker_2_amd(command_buffer, stage::Integer, dst_buffer, dst_offset::Integer, marker::Integer)::Cvoid


text:
struct CutoffApprox
    cutoff::Float64
end
struct DropArgMinApprox
    N::Int
end
struct SORApprox
    N::Int
end

function approximate!(a::CutoffApprox, f::BayesFilter)
    idxtodrop = sort!(findall(x -> x < a.cutoff, f.probabilities), rev = true)
    for i in idxtodrop
        splice!(f.probabilities, i)
        splice!(f.posteriors, i)
    end
    f.probabilities ./= sum(f.probabilities)
end
function approximate!(a::DropArgMinApprox, f::BayesFilter)
    M = length(f.probabilities)
    if M > a.N
        i = argmin(f.probabilities)
        splice!(f.probabilities, i)
        splice!(f.posteriors, i)
    end
    f.probabilities ./= sum(f.probabilities)
end
function approximate!(a::SORApprox, f::BayesFilter)
    M = length(f.probabilities)
    if M > a.N
        dropidx = SOR!(f.probabilities, a.N, f.rng)
        for i in sort!(dropidx, rev = true)
            splice!(f.probabilities, i)
            splice!(f.posteriors, i)
        end
    end
    f.probabilities ./= sum(f.probabilities)
end
function SOR!(w, N, rng)
    c, keepabove, perm = findc(w, N)
    dropidx = Int[]
    u = c*rand(rng)
    for i in reverse(perm[1:keepabove])
        u -= w[i]
        if u <= 0
            w[i] = c  #This is in Fearnhead & Liu 2007
            u += c
        else
            push!(dropidx, i)
            w[i] = 0.
        end
    end
    dropidx
end
function findc(w, N)
    M = length(w)
    idx = sortperm(w)
    A = M - N + 1
    Σ = sum(w[idx[i]] for i in 1:A)
    c = 0.
    while true
        c = Σ/(N + A - M)
        (A == M || w[idx[A]] < c < w[idx[A+1]]) && break
        A += 1
        Σ += w[idx[A]]
    end
    (c = c, keepabove = A, perm = idx)
end


text:
#
# fft-tests.jl -
#
# Tests for FFT and circulant convolution operators.
#
module TestingLazyAlgebraFFT

using Test
using LazyAlgebra
using LazyAlgebra.Foundations
using AbstractFFTs, FFTW

floats = (Float32, Float64)
types = (Float32, #Float64, Complex{Float32},
         Complex{Float64})
alphas = (0, 1, -1,  2.71, π)
betas = (0, 1, -1, -1.33, Base.MathConstants.φ)

lpad_args(n::Integer, x...) = lpad(string(x...), n)
lpad_with_parentheses(n::Integer, x...) = lpad(string("(", x..., ")"), n)

@testset "FFT utilities                          " begin
    dims1 = (1, 2, 3, 4, 5, 7, 9, 287, 511)
    dims2 = (1, 2, 3, 4, 5, 8, 9, 288, 512)
    dims3 = (1, 2, 3, 4, 5, 8, 9, 288, 512)
    @test LazyAlgebra.FFTs.goodfftdims(dims1) == dims2
    @test LazyAlgebra.FFTs.goodfftdims(map(Int16, dims1)) == dims2
    @test LazyAlgebra.FFTs.goodfftdims(dims1...) == dims2
    @test LazyAlgebra.FFTs.rfftdims(1,2,3,4,5) == (1,2,3,4,5)
    @test LazyAlgebra.FFTs.rfftdims(2,3,4,5) == (2,3,4,5)
    @test LazyAlgebra.FFTs.rfftdims(3,4,5) == (2,4,5)
    @test LazyAlgebra.FFTs.rfftdims(4,5) == (3,5)
    @test LazyAlgebra.FFTs.rfftdims(5) == (3,)
    @test LazyAlgebra.FFTs.fftfreq(1) == [0]
    @test LazyAlgebra.FFTs.fftfreq(2) == [0,-1]
    @test LazyAlgebra.FFTs.fftfreq(3) == [0,1,-1]
    @test LazyAlgebra.FFTs.fftfreq(4) == [0,1,-2,-1]
    @test LazyAlgebra.FFTs.fftfreq(5) == [0,1,2,-2,-1]
end # testset

@testset "FFT operator $(lpad_with_parentheses(26, T))" for T in types
    R = real(T)
    ϵ = sqrt(eps(R)) # relative tolerance, can certainly be much tighter
    for dims in ((45,), (20,), (33,12), (30,20), (4,5,6))
        # Create a FFT operator given its argument.
        x = rand(T, dims)
        n = length(x)
        xsav = vcopy(x)
        F = FFTOperator(x)
        @test x == xsav # check that input has been preserved
        @test MorphismType(F) == (T<:Complex ? Endomorphism() : Morphism())
        @test input_size(F) == dims
        @test input_size(F) == ntuple(i->input_size(F,i), ndims(x))
        @test output_size(F) == (T<:Complex ? dims :
                                 Tuple(AbstractFFTs.rfft_output_size(x, 1:ndims(x))))
        @test output_size(F) == ntuple(i->output_size(F,i), ndims(x))
        @test input_eltype(F) == T
        @test output_eltype(F) == typeof(complex(zero(R)))
        @test input_ndims(F) == ndims(x)
        @test output_ndims(F) == ndims(x)
        if T<:Complex
            @test LazyAlgebra.FFTs.destroys_input(F.forward) == true
            @test LazyAlgebra.FFTs.destroys_input(F.backward) == true
        else
            @test LazyAlgebra.FFTs.preserves_input(F.forward) == true
            @test LazyAlgebra.FFTs.preserves_input(F.backward) == false
        end
        @test LazyAlgebra.identical(F, F) == true
        io = IOBuffer()
        show(io, F)
        @test String(take!(io)) == "FFT"
        @test typeof(F') === Adjoint{typeof(F)}
        @test typeof(inv(F)) === Inverse{typeof(F)}
        @test typeof(inv(F)') === InverseAdjoint{typeof(F)}
        @test typeof(inv(F')) === InverseAdjoint{typeof(F)}
        @test F'*F == length(x)*Identity()
        @test F*F' == length(x)*Identity()
        @test inv(F)*F == Identity()
        @test F*inv(F) == Identity()
        @test inv(F')*F' == Identity()
        @test F'*inv(F') == Identity()
        @test inv(F')*inv(F) == (1//length(x))*Identity()
        @test inv(F)*inv(F') == (1//length(x))*Identity()

        # Create operators which should be considered as the same as F.
        F1 =  FFTOperator(T, dims...)
        @test LazyAlgebra.identical(F1, F) == true
        F2 =  FFTOperator(T, map(Int16, dims))
        @test LazyAlgebra.identical(F2, F) == true

        # Check applying operator.
        xbad = rand(T, ntuple(i -> (i == 1 ? dims[i]+1 : dims[i]), length(dims)))
        @test_throws DimensionMismatch F*xbad
        y = rand(Complex{R}, (T<:Complex ? input_size(F) :
                              output_size(F)))
        z = (T<:Complex ? fft(x) : rfft(x))
        @test x == xsav # check that input has been preserved
        ysav = vcopy(y)
        w = (T<:Complex ? ifft(y) : irfft(y, dims[1]))
        @test F*x ≈ z atol=0 rtol=ϵ norm=vnorm2
        @test F\y ≈ w atol=0 rtol=ϵ norm=vnorm2
        @test F(x) == F*x
        @test F'(y) == F'*y
        @test inv(F)(y) == F\y
        @test inv(F')(x) == F'\x
        @test x == xsav # check that input has been preserved
        @test y == ysav # check that input has been preserved
        @test y == ysav # check that input has been preserved
        for α in alphas,
            β in betas,
            scratch in (false, true)
            @test apply!(α, Direct, F, x, scratch, β, vcopy(y)) ≈
                R(α)*z + R(β)*y atol=0 rtol=ϵ
            if scratch
                vcopy!(x, xsav)
            else
                @test x == xsav # check that input has been preserved
            end
            @test apply!(α, Adjoint, F, y, scratch, β, vcopy(x)) ≈
                R(n*α)*w + R(β)*x atol=0 rtol=ϵ
            if scratch
                vcopy!(y, ysav)
            else
                @test y == ysav # check that input has been preserved
            end
            @test apply!(α, Inverse, F, y, scratch, β, vcopy(x)) ≈
                R(α)*w + R(β)*x atol=0 rtol=ϵ
            if scratch
                vcopy!(y, ysav)
            else
                @test y == ysav # check that input has been preserved
            end
            @test apply!(α, InverseAdjoint, F, x, scratch, β, vcopy(y)) ≈
                R(α/n)*z + R(β)*y atol=0 rtol=ϵ
            if scratch
                vcopy!(x, xsav)
            else
                @test x == xsav # check that input has been preserved
            end
        end
    end
end # testset

@testset "Circular convolution $(lpad_with_parentheses(18, T))" for T in types
    R = real(T)
    ϵ = sqrt(eps(R)) # relative tolerance, can certainly be much tighter
    n1, n2, n3 = 18, 12, 4
    for dims in ((n1,), (n1,n2), (n1,n2,n3))
        # Basic methods.
        x = rand(T, dims)
        n = length(x)
        h = rand(T, dims)
        H = CirculantConvolution(h, shift=false, flags=FFTW.ESTIMATE)
        @test MorphismType(H) == Endomorphism()
        @test input_size(H) == dims
        @test input_size(H) == ntuple(i->input_size(H,i), ndims(x))
        @test output_size(H) == dims
        @test output_size(H) == ntuple(i->output_size(H,i), ndims(x))
        @test input_eltype(H) == T
        @test output_eltype(H) == T
        @test input_ndims(H) == ndims(x)
        @test output_ndims(H) == ndims(x)
        @test eltype(H) == T
        @test size(H) == (dims..., dims...)
        @test ndims(H) == 2*length(dims)
        @test (size(H)..., 1) == ntuple(i->size(H, i), ndims(H)+1)

        # Test apply! method.
        F = FFTOperator(x)
        G = F\Diag(F*h)*F
        y = rand(T, dims)
        xsav = vcopy(x)
        ysav = vcopy(y)
        y1 = H*x
        @test H(x) == y1
        y2 = G*x
        y3 = (T<:Real ? real(ifft(fft(h).*fft(x)))
              :              ifft(fft(h).*fft(x)))
        @test y1 ≈ y2 atol=0 rtol=ϵ
        @test y1 ≈ y3 atol=0 rtol=ϵ
        if T<:Real
            y4 = irfft(rfft(h).*rfft(x), n1)
            @test y1 ≈ y4 atol=0 rtol=ϵ
        end
        z1 = H'*y
        z2 = G'*y
        z3 = (T<:Real ? real(ifft(conj.(fft(h)).*fft(y)))
              :              ifft(conj.(fft(h)).*fft(y)))
        @test z1 ≈ z2 atol=0 rtol=ϵ
        @test z1 ≈ z3 atol=0 rtol=ϵ
        if T<:Real
            z4 = irfft(conj.(rfft(h)).*rfft(y), n1)
            @test z1 ≈ z4 atol=0 rtol=ϵ
        end
        for α in alphas,
            β in betas,
            scratch in (false, true)
            @test apply!(α, Direct, H, x, scratch, β, vcopy(y)) ≈
                R(α)*y1 + R(β)*y atol=0 rtol=ϵ
            if scratch
                vcopy!(x, xsav)
            else
                @test x == xsav # check that input has been preserved
            end
            @test apply!(α, Adjoint, H, y, scratch, β, vcopy(x)) ≈
                R(α)*z1 + R(β)*x atol=0 rtol=ϵ
            if scratch
                vcopy!(y, ysav)
            else
                @test y == ysav # check that input has been preserved
            end
        end
    end
end # testset
nothing

end # module


text:
function with_prefix(ls::AbstractVector{<:AbstractString}, prefix)
    ls[startswith.(ls, prefix)]
end


text:
using Revise
using LaTeXStrings
using PyPlot
using Infiltrator
using Bem2d

"""
    plotfields_contours(els, xobs, yobs, idx, field, title)

Contouring of a single field for six-panel plotfields.
"""
function plotfields_contours(els, xobs, yobs, idx, field, titlestring)
    ncontours = 40
    xlim = [minimum(xobs) maximum(xobs)]
    ylim = [minimum(yobs) maximum(yobs)]
    subplot(2, 3, idx)
    scale = 1.0
    fieldmax = maximum(@.abs(field))
    contourf(xobs, yobs, reshape(field, size(xobs)), ncontours,
        vmin = -scale * fieldmax, vmax = scale * fieldmax, cmap = rycroftcmap())
    clim(-scale * fieldmax, scale * fieldmax)
    colorbar(fraction = 0.020, pad = 0.05, extend = "both")
    contour(xobs, yobs, reshape(field, size(xobs)), ncontours,
        vmin = -scale * fieldmax, vmax = scale * fieldmax, linewidths = 0.25, colors = "k")
    title(titlestring)
    stylesubplots(xlim, ylim)
    plotelements(els)
    return nothing
end

"""
    stylesubplots(xlim, ylim)

Common styling for subplots generated by plotfields.
"""
function stylesubplots(xlim, ylim)
    gca().set_aspect("equal")
    gca().set_xlim([xlim[1], xlim[2]])
    gca().set_ylim([ylim[1], ylim[2]])
    gca().set_xticks([xlim[1], 0, xlim[2]])
    gca().set_yticks([ylim[1], 0, ylim[2]])
    return nothing
end

"""
    plotfieldslocal(els, xobs, yobs, disp, stress, suptitlestring)

Set up six panel plots for displacement and stress fields.
"""
function plotfieldslocal(els, xobs, yobs, disp, stress, suptitlestring)
    figure(figsize = (15, 5))

    subplot(2, 3, 1)
    npts = 100
    uxreduced = reshape(disp[:, 1], npts, npts);
    uyreduced = reshape(disp[:, 2], npts, npts);
    uxreduced = uxreduced[1:10:end, 1:10:end];
    uyreduced = uyreduced[1:10:end, 1:10:end];
    xobsreduced = reshape(xobs, npts, npts);
    yobsreduced = reshape(yobs, npts, npts);
    xobsreduced = xobsreduced[1:10:end, 1:10:end];
    yobsreduced = yobsreduced[1:10:end, 1:10:end];
    quiver(xobsreduced[:], yobsreduced[:], uxreduced[:], uyreduced[:], units = "width", color = "b")
    stylesubplots([minimum(xobs) maximum(xobs)], [minimum(yobs) maximum(yobs)])
    plotelements(els)
    ylabel("y (m)")
    title(L"\mathbf{u}")

    plotfields_contours(els, xobs, yobs, 2, disp[:, 1], L"u_x")
    plotfields_contours(els, xobs, yobs, 3, disp[:, 2], L"u_y")
    plotfields_contours(els, xobs, yobs, 4, stress[:, 1], L"\sigma_{xx}")
    xlabel("x (m)")
    ylabel("y (m)")
    plotfields_contours(els, xobs, yobs, 5, stress[:, 2], L"\sigma_{yy}")
    xlabel("x (m)")
    plotfields_contours(els, xobs, yobs, 6, stress[:, 3], L"\sigma_{xy}")
    xlabel("x (m)")
    tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
    return nothing
end

"""
    basicelements()

Compare farfield displacements and stresses from constant displacement
and quadratic displacement elements.  The point of this is to highlight
the fact the quadratic elements can be used to create linear slip elements
that can accurately model a linear displacement gradient with no spurious
singularities.
"""
function basicelements()
    mu = 3e10
    nu = 0.25

    # Create a flat fault
    els = Elements(Int(1e5))
    nels = 1
    x1, y1, x2, y2 = discretizedline(-10e3, 0, 10e3, 0, nels)
    for i in 1:length(x1)
        els.x1[els.endidx + i] = x1[i]
        els.y1[els.endidx + i] = y1[i]
        els.x2[els.endidx + i] = x2[i]
        els.y2[els.endidx + i] = y2[i]
        els.name[els.endidx + i] = "fault"
    end
    standardize_elements!(els)

    # Convenience dictionary for element names
    idx = getidxdict(els)

    # Set fault slip
    xcenters = els.xcenter[1:els.endidx]
    ycenters = els.ycenter[1:els.endidx]
    xnodes = sort(els.xnodes[1:els.endidx, :][:])
    ynodes = sort(els.ynodes[1:els.endidx, :][:])

    # Constant x-slip only
    constxslip = ones(nels)
    constyslip = zeros(nels)
    quadxslip = [0.25 0.50 -0.15]
    quadyslip = [0.00 0.00 0.00]

    # Observation coordinates for far-field calculation
    npts = 100
    obswidth = 20e3
    xobs, yobs = obsgrid(-obswidth, -obswidth, obswidth, obswidth, npts)

    # Displacements and stresses
    dispconst, stressconst = constdispstress(slip2dispstress, xobs, yobs, els, idx["fault"], constxslip, constyslip, mu, nu)
    dispquad, stressquad = quaddispstress(slip2dispstress, xobs, yobs, els, idx["fault"], quadxslip, quadyslip, mu, nu)

    # Plot
    close("all")
    plotfieldslocal(els, reshape(xobs, npts, npts), reshape(yobs, npts, npts),
               dispconst, stressconst, "CS element")
    plotfieldslocal(els, reshape(xobs, npts, npts), reshape(yobs, npts, npts),
               dispquad, stressquad, "3QN quadratic slip element")
    return nothing
end
basicelements()


text:
"""
    Verlet(∇V!, M, Δt)

Set up the Verlet integrator.

### Fields

* ∇V!   - In place gradient of the potential
* M     - Mass (either scalar or vector)
* Δt    - Time step
"""
struct Verlet{TGV, TF<:AbstractFloat, TM} <: SecondOrderNonMetropolisSampler
    ∇V!::TGV
    M::TM
    Δt::TF
end

mutable struct VerletState{Tq, Tx<:AbstractVector{Tq}} <:SecondOrderNonMetropolisSamplerState
    x::Tx
    ∇V::Tq
    p_mid::Tq
end

function InitState!(x₀, sampler::Verlet)
    ∇V = similar(x₀[1])
    sampler.∇V!(∇V , x₀[1]);
    return VerletState(x₀, copy(∇V) , similar(x₀[1]));
end

function InitState(x₀, sampler::Verlet)
    ∇V = similar(x₀[1])
    sampler.∇V!(∇V , x₀[1]);
    return VerletState(deepcopy(x₀), copy(∇V) , similar(x₀[1]));
end

function UpdateState!(state::VerletState, sampler::Verlet)
    @. state.p_mid = state.x[2] - 0.5 * sampler.Δt * state.∇V;
    @. state.x[1] = state.x[1] + sampler.Δt * state.p_mid/sampler.M;
    sampler.∇V!(state.∇V, state.x[1]);
    @. state.x[2] = state.p_mid - 0.5 * sampler.Δt * state.∇V;
    state
end

max_stars_repo_name: gideonsimpson/JuBasicMD


text:
using Statistics

# # After implementing, I think I saw there is a Iterators.product that does this?
function combine(dict_arr, key, vals)
    new_dict_arr = []
    for dict in dict_arr
        for val in vals
            new_dict = deepcopy(Dict(dict))
            new_dict[key] = val
            push!(new_dict_arr, new_dict)
        end
    end
    return new_dict_arr
end

function split_algo(algo_ic, swept_params)
    # println(swept_params)
    dicts = [Dict()]
    for k in swept_params
        list_of_params = diff(algo_ic)[k]
        dicts = combine(dicts, k, list_of_params)
    end
    return dicts
end

function get_best(ic, sweep_params, metric, percentage=(0, 1))
    splits = split_algo(ic,sweep_params)
    errors = ones(length(splits)) * Inf
    for (ind, split) in enumerate(splits)
        res = load_results(search(ic, split), metric)
        if length(size(res)) == 3
            num_log = size(res)[2]
            start_index = floor(Int, num_log * percentage[1])
            if (start_index == 0)
                start_index = 1
            end
            end_index = ceil(Int, num_log * percentage[2])

            # println("start_index $(start_index), end_index $(end_index)")
            # println(size(res))
            res = res[:, start_index:end_index, :]
        end

        error = mean(res)
        errors[ind] = error
    end
    low_err, low_err_ind = findmin(errors)

    return search(ic, splits[low_err_ind])
end

function print_params(ic, algo_params, sweep_params)
    println("For IC: ")
    for p in algo_params
        println(p,"  ", ic[1].parsed_args[p])
    end
    for p in sweep_params
        println(p,"  ", ic[1].parsed_args[p])
    end
    println()
end

function load_results(ic, logger_key; return_type = "tensor")
    num_results = length(ic)
    results = []
    for itm in ic.items
        file_path = joinpath(itm.folder_str, "results.jld2")
        if (isfile(file_path))
            try
                data = FileIO.load(file_path)["results"]
                push!(results,data[logger_key])
            catch e
                println("Error for file $(file_path)")
                println(e)
                continue
            end

        end
    end

    if (size(results)[1] == 0)
        println("empty collection")
        return [2^31]
    end

    if return_type == "tensor"
        return cat(results..., dims = 3)
    elseif return_type == "array"
        return results
    end
end

function single_onehot(data, num_gvfs)
    onehot_enc = zeros(num_gvfs)
    for gvf_i in 1:num_gvfs
        onehot_enc[gvf_i] += (data .== gvf_i)
    end
    return onehot_enc
end

function get_single_goal_percentage(visits; step_size=10)
    num_gvfs = 4
    m_length = length(visits)
    goal_visits = zeros(num_gvfs, ceil(Integer, m_length / step_size))

    for episode in 1:length(visits)
        index = ceil(Integer, episode / step_size)
        goal_visits[:, index] += single_onehot(visits[episode], num_gvfs)
    end
    goal_visit_percentage = goal_visits / step_size
end

function smooth_lines(data, step)
    # really simple smooth line function
    means = []
    for start in 1:step:size(data)[1]
        chunk = data[start:min(start + step - 1, size(data)[1])]
        chunk_mean = mean(chunk)
        push!(means, chunk_mean)
    end
    return means
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
"PDDL state description."
abstract type State end

"Returns an integer index for specified `state`. Defaults to hashing."
stateindex(domain::Domain, state::State) = hash(state)

"Returns an iterator over objects in the `state`."
get_objects(state::State) = error("Not implemented.")

"Returns an iterator over objects in the `state` with a particular `type`."
get_objects(state::State, type::Symbol) = error("Not implemented.")

"Returns a map from state objects to their types."
get_objtypes(state::State) = error("Not implemented.")

"Returns the type of an `object` in a `state`."
get_objtype(state::State, object) = get_objtypes[object]

"Returns an iterator over true Boolean predicates in a `state`."
get_facts(state::State) = error("Not implemented.")

"Gets the value of a fluent specified by `term`."
get_fluent(state::State, term::Term) = error("Not implemented.")

"Sets the value of a fluent specified by `term` to `val`."
set_fluent!(state::State, val, term::Term) = error("Not implemented.")

"Returns a map from fluent names to values (false predicates may be omitted)."
get_fluents(state::State) = error("Not implemented.")

"Returns the names of fluents in a state (false predicates may be omitted)."
get_fluent_names(state::State) = (k for (k, v) in get_fluents(state))

"Returns the values of fluents in a state (false predicates may be omitted)."
get_fluent_values(state::State) = (v for (k, v) in get_fluents(state))

Base.getindex(state::State, term::Term) = get_fluent(state, term)

Base.setindex!(state::State, val, term::Term) = set_fluent!(state, val, term)

Base.pairs(state::State) = get_fluents(state)

Base.keys(state::State) = get_fluent_names(state)

Base.values(state::State) = get_fluent_values(state)


text:
module PersistentArrays

# package code goes here
include("node.jl")
include("array.jl")

export PersistentArray, version, lookup, update!, commit!, increment, Node
end # module


text:
using Push
using Test

cfg_path = joinpath(dirname(@__FILE__), "configuration/float.cfg")
cfg = Push.load_configuration(cfg_path)

# FLOAT.*
s = Push.run("(1.0 FLOAT.*)", cfg)
@test s.float == [1.0]
s = Push.run("(2.5 3.0 50.0 0.1 FLOAT.*)", cfg)
@test s.float == [2.5, 3.0, 5.0]
s = Push.run("(2.5 3.0 -50.0 10.0 FLOAT.*)", cfg)
@test s.float == [2.5, 3.0, -500.0]

# FLOAT.+
s = Push.run("(1.0 FLOAT.+)", cfg)
@test s.float == [1.0]
s = Push.run("(7.0 3.0 2.5 10.25 FLOAT.+)", cfg)
@test s.float == [7.0, 3.0, 12.75]
s = Push.run("(7.0 3.0 2.75 -10.25 FLOAT.+)", cfg)
@test s.float == [7.0, 3.0, -7.5]

# FLOAT.-
s = Push.run("(1.0 FLOAT.-)", cfg)
@test s.float == [1.0]
s = Push.run("(7.0 3.0 10.75 2.25 FLOAT.-)", cfg)
@test s.float == [7.0, 3.0, 8.5]
s = Push.run("(7.0 3.0 10.0 25.0 FLOAT.-)", cfg)
@test s.float == [7.0, 3.0, -15.0]

# FLOAT./
s = Push.run("(1.0 FLOAT./)", cfg)
@test s.float == [1.0]
s = Push.run("(7.0 3.0 10.0 2.0 FLOAT./)", cfg)
@test s.float == [7.0, 3.0, 5.0]
s = Push.run("(7.0 3.0 -10.0 4.0 FLOAT./)", cfg)
@test s.float == [7.0, 3.0, -2.5]
s = Push.run("(7.0 3.0 -10.0 0.0 FLOAT./)", cfg)
@test s.float == [7.0, 3.0, -10.0, 0.0]

# FLOAT.<
s = Push.run("(1.0 FLOAT.<)", cfg)
@test s.float == [1.0] && isempty(s.boolean)
s = Push.run("(100.0 900.0 2.0 3.0 FLOAT.<)", cfg)
@test s.float == [100.0, 900.0] && s.boolean == [true]
s = Push.run("(100.0 900.0 5.0 -10.0 FLOAT.<)", cfg)
@test s.float == [100.0, 900.0] && s.boolean == [false]
s = Push.run("(100.0 900.0 5.0 5.0 FLOAT.<)", cfg)
@test s.float == [100.0, 900.0] && s.boolean == [false]

# FLOAT.=
s = Push.run("(1.0 FLOAT.=)", cfg)
@test s.float == [1.0] && isempty(s.boolean)
s = Push.run("(100.0 900.0 -10.0 90.0 FLOAT.=)", cfg)
@test s.float == [100.0, 900.0] && s.boolean == [false]
s = Push.run("(100.0 900.0 -10.0 -10.0 FLOAT.=)", cfg)
@test s.float == [100.0, 900.0] && s.boolean == [true]

# FLOAT.>
s = Push.run("(1.0 FLOAT.>)", cfg)
@test s.float == [1.0] && isempty(s.boolean)
s = Push.run("(100.0 900.0 2.0 3.0 FLOAT.>)", cfg)
@test s.float == [100.0, 900.0] && s.boolean == [false]
s = Push.run("(100.0 900.0 5.0 -10.0 FLOAT.>)", cfg)
@test s.float == [100.0, 900.0] && s.boolean == [true]
s = Push.run("(100.0 900.0 5.0 5.0 FLOAT.>)", cfg)
@test s.float == [100.0, 900.0] && s.boolean == [false]

# FLOAT.DUP
s = Push.run("(10.0 20.0 30.0 FLOAT.DUP)", cfg)
@test s.float == [10.0, 20.0, 30.0, 30.0]

# FLOAT.FLUSH
s = Push.run("(10.0 20.0 30.0 FLOAT.FLUSH)", cfg)
@test s.float == []

# FLOAT.FROMBOOLEAN
s = Push.run("(10.0 20.0 TRUE FALSE FLOAT.FROMBOOLEAN)", cfg)
@test s.float == [10.0, 20.0, 0.0] && s.boolean == [true]
s = Push.run("(10.0 20.0 FALSE TRUE FLOAT.FROMBOOLEAN)", cfg)
@test s.float == [10.0, 20.0, 1.0] && s.boolean == [false]

# FLOAT.FROMINTEGER
s = Push.run("(89 FLOAT.FROMINTEGER)", cfg)
@test s.integer == [] && s.float == [89.0]

# FLOAT.MAX
s = Push.run("(10.0 10.0 FLOAT.MAX)", cfg)
@test s.float == [10.0]
s = Push.run("(20.0 10.0 FLOAT.MAX)", cfg)
@test s.float == [20.0]
s = Push.run("(10.0 20.0 FLOAT.MAX)", cfg)
@test s.float == [20.0]

# FLOAT.MIN
s = Push.run("(10.0 10.0 FLOAT.MIN)", cfg)
@test s.float == [10.0]
s = Push.run("(20.0 10.0 FLOAT.MIN)", cfg)
@test s.float == [10.0]
s = Push.run("(10.0 20.0 FLOAT.MIN)", cfg)
@test s.float == [10.0]

# FLOAT.POP
s = Push.run("(10.0 20.0 30.0 FLOAT.POP)", cfg)
@test s.float == [10.0, 20.0]

# FLOAT.ROT
s = Push.run("(1.0 FLOAT.ROT)", cfg)
@test s.float == [1.0]
s = Push.run("(0.0 10.0 20.0 30.0 FLOAT.ROT)", cfg)
@test s.float == [0.0, 30.0, 10.0, 20.0]

# FLOAT.SHOVE
s = Push.run("(1.0 FLOAT.SHOVE)", cfg)
@test s.float == [1.0]
s = Push.run("(2.0 0 FLOAT.SHOVE)", cfg)
@test s.float == [2.0]
s = Push.run("(2.0 9000 FLOAT.SHOVE)", cfg)
@test s.float == [2.0]
s = Push.run("(1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 900.0 -20 FLOAT.SHOVE)", cfg)
@test s.float == [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,900.0]
s = Push.run("(1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 900.0 0 FLOAT.SHOVE)", cfg)
@test s.float == [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 900.0]
s = Push.run("(1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 900.0 1 FLOAT.SHOVE)", cfg)
@test s.float == [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 900.0, 8.0]
s = Push.run("(1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 900.0 2 FLOAT.SHOVE)", cfg)
@test s.float == [1.0,2.0,3.0,4.0,5.0,6.0,900.0,7.0,8.0]
s = Push.run("(1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 900.0 9 FLOAT.SHOVE)", cfg)
@test s.float == [900.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0]
s = Push.run("(1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 900.0 1000 FLOAT.SHOVE)", cfg)
@test s.float == [900.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0]

# FLOAT.STACKDEPTH
s = Push.run("(FLOAT.STACKDEPTH)", cfg)
@test s.integer == [0]
s = Push.run("(1 20.0 30.0 40.0 FLOAT.STACKDEPTH)", cfg)
@test s.integer == [1, 3] && s.float == [20.0, 30.0, 40.0]

# FLOAT.SWAP
s = Push.run("(1.0 FLOAT.SWAP)", cfg)
@test s.float == [1.0]
s = Push.run("(1.0 2.0 3.0 4.0 5.0 6.0 FLOAT.SWAP)", cfg)
@test s.float == [1.0,2.0,3.0,4.0,6.0,5.0]

# FLOAT.YANK
s = Push.run("(1.0 FLOAT.YANK)", cfg)
@test s.float == [1.0]
s = Push.run("(10.0 0 FLOAT.YANK)", cfg)
@test s.float == [10.0]
s = Push.run("(10.0 20.0 30.0 40.0 50.0 -90 FLOAT.YANK)", cfg)
@test s.float == [10.0, 20.0, 30.0, 40.0, 50.0]
s = Push.run("(10.0 20.0 30.0 40.0 50.0 0 FLOAT.YANK)", cfg)
@test s.float == [10.0, 20.0, 30.0, 40.0, 50.0]
s = Push.run("(10.0 20.0 30.0 40.0 50.0 1 FLOAT.YANK)", cfg)
@test s.float == [10.0, 20.0, 30.0, 50.0, 40.0]
s = Push.run("(10.0 20.0 30.0 40.0 50.0 2 FLOAT.YANK)", cfg)
@test s.float == [10.0, 20.0, 40.0, 50.0, 30.0]
s = Push.run("(10.0 20.0 30.0 40.0 50.0 5 FLOAT.YANK)", cfg)
@test s.float == [20.0, 30.0, 40.0, 50.0, 10.0]
s = Push.run("(10.0 20.0 30.0 40.0 50.0 987 FLOAT.YANK)", cfg)
@test s.float == [20.0, 30.0, 40.0, 50.0, 10.0]

# FLOAT.YANKDUP
s = Push.run("(1.0 FLOAT.YANKDUP)", cfg)
@test s.float == [1.0]
s = Push.run("(10.0 20.0 30.0 40.0 50.0 0 FLOAT.YANKDUP)", cfg)
@test s.float == [10.0, 20.0, 30.0, 40.0, 50.0, 50.0]
s = Push.run("(10.0 20.0 30.0 40.0 50.0 1 FLOAT.YANKDUP)", cfg)
@test s.float == [10.0, 20.0, 30.0, 40.0, 50.0, 40.0]
s = Push.run("(10.0 20.0 30.0 40.0 50.0 2 FLOAT.YANKDUP)", cfg)
@test s.float == [10.0, 20.0, 30.0, 40.0, 50.0, 30.0]
s = Push.run("(10.0 20.0 30.0 40.0 50.0 3 FLOAT.YANKDUP)", cfg)
@test s.float == [10.0, 20.0, 30.0, 40.0, 50.0, 20.0]
s = Push.run("(10.0 20.0 30.0 40.0 50.0 4 FLOAT.YANKDUP)", cfg)
@test s.float == [10.0, 20.0, 30.0, 40.0, 50.0, 10.0]
s = Push.run("(10.0 20.0 30.0 40.0 50.0 987 FLOAT.YANKDUP)", cfg)
@test s.float == [10.0, 20.0, 30.0, 40.0, 50.0, 10.0]
s = Push.run("(10.0 20.0 30.0 40.0 50.0 -987 FLOAT.YANKDUP)", cfg)
@test s.float == [10.0, 20.0, 30.0, 40.0, 50.0, 50.0]

# FLOAT.TAN
s = Push.run("(90.0 FLOAT.TAN)", cfg)
@test s.float == Float32[tan(90.0)]

# FLOAT.COS
s = Push.run("(90.0 FLOAT.COS)", cfg)
@test s.float == Float32[cos(90.0)]

# FLOAT.SIN
s = Push.run("(90.0 FLOAT.SIN)", cfg)
@test s.float == Float32[sin(90.0)]

# FLOAT.DEFINE
s = Push.run("(X FLOAT.DEFINE X)", cfg)
@test isempty(s.float) && s.name == [:X, :X]
s = Push.run("(3.0 X FLOAT.DEFINE X)", cfg)
@test isempty(s.name) && s.float == [3.0]

# FLOAT.RAND
s = Push.run("(FLOAT.RAND)", cfg)
@test length(s.float) == 1

# FLOAT.%


text:
module Prim

using ..Omega
using ..IDS: ID, uid
using ..NonDet: RandVar, URandVar, MaybeRV, isconstant, ppapl, apl, reify, elemtype
import ..NonDet: name, ppapl, ciid

import ..Omega:Ω , params, lift

using ..Util
using Spec
import Statistics: quantile
import Distributions
const Djl = Distributions
import Random
using DocStringExtensions: SIGNATURES

export  bernoulli,
        betarv,
        β,
        categorical,
        # dirichlet,
        exponential,
        gammarv,
        Γ,
        invgamma,
        kumaraswamy,
        logistic,
        # mvnormal,
        normal,
        poisson,
        rademacher,
        uniform

"Primitive random variable of known distribution"
abstract type PrimRandVar <: RandVar end

name(t::T) where {T <: PrimRandVar} = T.name.name

"Parameters of `rv`"
@generated function params(rv::PrimRandVar)
  fields = [Expr(:., :rv, QuoteNode(f)) for f in fieldnames(rv) if f !== :id]
  Expr(:tuple, fields...)
end

ppapl(rv::PrimRandVar, ωπ) = rvtransform(rv)(ωπ, reify(ωπ, params(rv))...)

# Helper for having primitives with multivariate parameters
@generated function anysize(args::Union{<:AbstractArray, Real}...)
  isarr = (arg -> arg <: AbstractArray).([args...])
  firstarr = findfirst(isarr)
  if isempty(isarr)
    :(())
  else
    :(size(args[$firstarr]))
  end
end
@spec same(size.(filter(a -> a isa AbstractArray, args)))

include("univariate.jl")      # Univariate Distributions
include("multivariate.jl")    # Multivariate Distributions
include("randx.jl")           # Interfacing with normal julia code

end


text:
mutable struct LocalGIFunctionApproximator{G<:AbstractGrid} <: LocalFunctionApproximator
	grid::G
	gvalues::Vector{Float64}
end
LocalGIFunctionApproximator(grid::G) where {G<:AbstractGrid} = LocalGIFunctionApproximator(grid,zeros(length(grid)))

################ INTERFACE FUNCTIONS ################
function n_interpolating_points(gifa::LocalGIFunctionApproximator)
    return length(gifa.grid)
end

function get_all_interpolating_points(gifa::LocalGIFunctionApproximator)
	return vertices(gifa.grid)
end

function get_all_interpolating_values(gifa::LocalGIFunctionApproximator)
	return gifa.gvalues
end

function get_interpolating_nbrs_idxs_wts(gifa::LocalGIFunctionApproximator, v::AbstractVector{Float64})
    return interpolants(gifa.grid, v)
end

function compute_value(gifa::LocalGIFunctionApproximator, v::AbstractVector{Float64})
    return interpolate(gifa.grid, gifa.gvalues, v)
end

function compute_value(gifa::LocalGIFunctionApproximator, v_list::AbstractVector{V}) where V <: AbstractVector{Float64}
    @assert length(v_list) > 0
    vals = [compute_value(gifa, pt) for pt in v_list]
    return vals
end

function set_all_interpolating_values(gifa::LocalGIFunctionApproximator, gvalues::AbstractVector{Float64})
    gifa.gvalues = copy(gvalues)
end

function finite_horizon_extension(lfa::LocalFunctionApproximator, hor::StepRange{Int64,Int64})
    cut_points = lfa.grid.cutPoints
    extended_grid = RectangleGrid(cut_points..., hor)

    return LocalGIFunctionApproximator(extended_grid)
end


text:
a = 355 // 113

Float64(a)


text:
using DSP
using VerilogFIR

# Create a low pass filter using DPS
coeffs = digitalfilter(Lowpass(5; fs=50), FIRWindow(hamming(8)));
# Create the Filter
generate_fir("filter.v",coeffs,w_input=16,w_output=16)


text:
"""
# ExcelReaders

## Synopsis
ExcelReaders allows you to read data from Excel files.

## Exported functions
* ``openxl``: Open an Excel file for repeated reads.
* ``readxl``: Read data from a specified range in an Excel file.
* ``readxlsheet``: Read data from a specified sheet in an Excel file.

## Exported types
* ``ExcelFile``: A handle to an Excel file that was opened with ``openxl``.
* ``ExcelErrorCell``: Value returned for cells with an Excel error code.

## Further information
The ``ExcelReaders.tutorial`` help topic provides a tutorial to the package.

The package homepage is [https://github.com/davidanthoff/ExcelReaders.jl](https://github.com/davidanthoff/ExcelReaders.jl). Please
report any issues with the package there.
"""
ExcelReaders

"""
# ExcelReaders.tutorial

## Basic usage

The most basic usage is this:

````julia
using ExcelReaders
data = readxl("Filename.xlsx", "Sheet1!A1:C4")
````

This will return an array with all the data in the cell range A1 to
C4 on Sheet1 in the Excel file Filename.xlsx.

If you expect to read multiple ranges from the same Excel file you can get much
better performance by opening the Excel file only once:

````julia
using ExcelReaders
f = openxl("Filename.xlsx")
data1 = readxl(f, "Sheet1!A1:C4")
data2 = readxl(f, "Sheet2!B4:F10")
````

## Reading a whole sheet

The ``readxlsheet`` function reads complete Excel sheets, without a need to
specify precise range information. The most basic usage is

````julia
using ExcelReaders
data = readxlsheet("Filename.xlsx", "Sheet1")
````

This will read all content on Sheet1 in the file Filename.xlsx. Eventual blank
rows and columns at the top and left are skipped. ``readxlsheet`` takes a number
of optional keyword arguments:

- ``skipstartrows`` accepts either ``:blanks`` (default) or a positive integer.
With ``:blank`` any empty initial rows are skipped. An integer skips as many rows
as specified.
- ``skipstartcols`` accepts either ``:blanks`` (default) or a positive integer.
With ``:blank`` any empty initial columns are skipped. An integer skips as many
columns as specified.
- ``nrows`` accepts either ``:all`` (default) or a positive integer. With
``:all``, all rows (except skipped ones) are read. An integer specifies the exact
number of rows to be read.
- ``ncols`` accepts either ``:all`` (default) or a postiive integer. With ``:all``,
all columns (except skipped ones) are read. An integer specifies the exact number of columns to be read.

``readxlsheet`` also accepts an ``ExcelFile`` (as obtained from ``openxl``) as its
first argument.
"""
tutorial = nothing


text:
# Define delay differential equation for solute concentrations in a
# stream analogue mesocosm (SAM)
macro sam_model(reactionKinetics, shape)
    if shape == "powerLaw"
        if reactionKinetics == "firstOrder"
        return esc(:(
            function sam_model!(du,u,h,p,t)
                τ₋, τ₊, α, k, qPerUnitHZ, sHZ, sCH = p
                if(t < τ₋)
                    integral = h(p, t-τ₋) * quadgk(τ -> E(τ, α, τ₋, τ₊)*firstOrder(τ, k),τ₋, τ₊)[1]
                elseif(τ₋<t<τ₊)
                    integral = h(p, t-τ₊) * quadgk(τ -> E(τ, α, τ₋, τ₊)*firstOrder(τ, k), t, τ₊)[1] +
                    quadgk(τ -> E(τ, α, τ₋, τ₊)* h(p,t-τ)[1]*firstOrderRFP(τ, k),τ₋, t)[1]
                else
                    integral = quadgk(τ -> E(τ, α, τ₋, τ₊)* h(p,t-τ)[1]*firstOrderRFP(τ, k),τ₋, τ₊)[1]
                end
                du[1] = qPerUnitHZ*(sHZ/sCH)*(integral - u[1])
              end
                ))
        elseif reactionKinetics == "firstOrderProductionAndConsumption"
        return esc(:(
            function sam_model!(du,u,h,p,t)
                τ₋, τ₊, α, qPerUnitHZ, sHZ, sCH, kprod, kcons, substrate0 = p
                if(t < τ₋)
                    integral = quadgk(τ -> E(τ, α, τ₋, τ₊)* firstOrderProductionAndConsumption(τ, kprod, kcons, substrate0, h(p,t-τ₋)[1]),τ₋, τ₊)[1]
                elseif(τ₋<t<τ₊)
                    integral = quadgk(τ -> E(τ, α, τ₋, τ₊)* firstOrderProductionAndConsumption(τ, kprod, kcons, substrate0, h(p,t-τ₊)[1]), t, τ₊)[1] +
                    quadgk(τ -> E(τ, α, τ₋, τ₊)* firstOrderProductionAndConsumption(τ, kprod, kcons, substrate0, h(p,t-τ)[1]),τ₋, t)[1]
                else
                    integral = quadgk(τ -> E(τ, α, τ₋, τ₊)* firstOrderProductionAndConsumption(τ, kprod, kcons, substrate0, h(p,t-τ)[1]),τ₋, τ₊)[1]
                end
                du[1] = qPerUnitHZ*(sHZ/sCH)*(integral - u[1])
              end
                ))
        end
    end
  end

max_stars_repo_name: elizabethjmohr/hyporheicBGC

max_issues_repo_name: elizabethjmohr/hyporheicBGC

max_forks_repo_name: elizabethjmohr/hyporheicBGC


text:
# Autogenerated wrapper script for Libical_jll for x86_64-linux-gnu-cxx11
export libical, libical_cxx, libicalss, libicalss_cxx, libicalvcal

using BerkeleyDB_jll
using Glib_jll
using ICU_jll
using XML2_jll
JLLWrappers.@generate_wrapper_header("Libical")
JLLWrappers.@declare_library_product(libical, "libical.so.3")
JLLWrappers.@declare_library_product(libical_cxx, "libical_cxx.so.3")
JLLWrappers.@declare_library_product(libicalss, "libicalss.so.3")
JLLWrappers.@declare_library_product(libicalss_cxx, "libicalss_cxx.so.3")
JLLWrappers.@declare_library_product(libicalvcal, "libicalvcal.so.3")
function __init__()
    JLLWrappers.@generate_init_header(BerkeleyDB_jll, Glib_jll, ICU_jll, XML2_jll)
    JLLWrappers.@init_library_product(
        libical,
        "lib/libical.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libical_cxx,
        "lib/libical_cxx.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libicalss,
        "lib/libicalss.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libicalss_cxx,
        "lib/libicalss_cxx.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libicalvcal,
        "lib/libicalvcal.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#=

  Magic sequence problem in Julia ConstraintSolver.jl

  http://www.dcs.st-and.ac.uk/~ianm/CSPLib/prob/prob019/spec.html
  """
  A magic sequence of length n is a sequence of integers x0 . . xn-1 between
  0 and n-1, such that for all i in 0 to n-1, the number i occurs exactly xi
  times in the sequence. For instance, 6,2,1,0,0,0,1,0,0,0 is a magic sequence
  since 0 occurs 6 times in it, 1 occurs twice, ...
  """

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

# Algorithmic approach
function magic_sequence_alg(n::Int64)
    s = zeros(Int64,n)
    s[1:3] .= n-4,2,1
    s[end-3] = 1
    return s
end

function magic_sequence(n=4,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            # "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    @variable(model, 0 <= x[1:n] <= n-1, Int)

    @constraint(model, sum(x) == n)
    scalar_product(model,collect(0:n-1),x,n)

    for i in 0:n-1
        count_ctr(model,x,:(==),i,x[i+1])
    end

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        # println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                # println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                println("x:$x_val")

            end
        end
    else
        println("status:$status")
    end

    return status
end

for n in 2:30
    print("n: $n  ")
    @time magic_sequence(n)
end

println("\nAlgorithmic approach")
@time magic_sequence_alg(100_000_000)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
import CompScienceMeshes.coordtype

export numfunctions, coordtype, scalartype, assemblydata, geometry, refspace, valuetype

abstract type RefSpace{T,D} end
abstract type AbstractSpace end
abstract type Space{T} <: AbstractSpace end

scalartype{T}(s::RefSpace{T}) = T
scalartype{T}(s::Space{T}) = T

"""
    refspace(basis)

Returns the ReferenceSpace of local shape functions on which the basis is built.
"""
function refspace end

"""
    coordtype(s)

The scalar field over which the argument to a basis function or operator's integration
kernel are defined. This is always a salar data type, even if the function or kernel
is defined over a multi-dimensional space.
"""
coordtype{T}(::RefSpace{T}) = T
coordtype{T}(::Space{T}) = T

"""
    numfunctions(r)

Return the number of functions in a `Space` or `RefSpace`.
"""
numfunctions{T,D}(rs::RefSpace{T,D}) = D

"""
    scalartype(x)

The scalar field over which the values of a global or local basis function, or an
operator are defined. This should always be a scalar type, even if the basis or
operator takes on values in a vector or tensor space. This data type is used to
determine the `eltype` of assembled discrete operators.
"""
function scalartype end
scalartype(x1, xs...) = Base.promote_type(scalartype(x1), scalartype(xs...))

"""
    geometry(basis)

Returns an iterable collection of geometric elements on which the functions
in `basis` are defined. The order the elements are encountered needs correspond
to the element indices used in the data structure returned by `assemblydata`.
"""
geometry(s::Space) = s.geo
basisfunction(s::Space, i) = s.fns[i]
numfunctions(space::Space) = length(space.fns)

type DirectProductSpace{T} <: AbstractSpace
    factors::Vector{Space{T}}
end

import Base.cross
cross{T}(a::Space{T}, b::Space{T}) = DirectProductSpace(Space{T}[a,b])
cross{T}(a::DirectProductSpace{T}, b::Space{T}) = DirectProductSpace(Space{T}[a.factors; b])
numfunctions(S::DirectProductSpace) = sum([numfunctions(s) for s in S.factors])
scalartype{T}(s::DirectProductSpace{T}) = T
geometry(x::DirectProductSpace) = weld(x.geo...)

type Shape{T}
  cellid::Int
  refid::Int
  coeff::T
end

immutable AssemblyData{T}
    data::Array{Tuple{Int,T},3}
end

Base.getindex(ad::AssemblyData, c, r) = ADIterator(c,r,size(ad.data,1),ad)

immutable ADIterator{T}
    c::Int
    r::Int
    I::Int
    ad::AssemblyData{T}
end

Base.start(it::ADIterator) = 1
Base.next(it::ADIterator, i) = (it.ad.data[i,it.r,it.c], i+1)
Base.done(it::ADIterator, i) = (it.I < i || it.ad.data[i,it.r,it.c][1] < 1)

function add!{T}(bf::Vector{Shape{T}}, cellid, refid, coeff)
    push!(bf, Shape(cellid, refid, coeff))
end

"""
    assemblydata(basis)

Given a Basis this function returns a data structure containing the information
required for matrix assemble. More precise the following expressions are valid
for the returned object `ad`:

```
ad[c,r,i].globalindex
ad[c,r,i].coefficient
```

Here, `c` and `r` are indices in the iterable set of geometric elements and the
set of local shape functions on each element. `i` ranges from 1 to the maximum
number of basis functions local shape function `r` on element `r` contributes
to.

For a triplet `(c,r,i)`, `globalindex` is the index in the Basis of the
`i`-th basis function that has a contribution from local shape function `r` on
element `r`. `coefficient` is the coefficient of that contribution in the
linear combination defining that basis function in terms of local shape
function.

*Note*: the indices `c` correspond to the position of the corresponding
element whilst iterating over `geometry(basis)`.
"""
function assemblydata(basis::Space)

    @assert numfunctions(basis) != 0

    T = coordtype(basis)

    geo = geometry(basis)
    num_cells = numcells(geo)

    num_bfs  = numfunctions(basis)
    num_refs = numfunctions(refspace(basis))

    # determine the maximum number of function defined over a given cell
    celltonum = zeros(Int, num_cells, num_refs)
    for b in 1 : num_bfs
        for shape in basisfunction(basis, b)
            c = shape.cellid
            r = shape.refid
            celltonum[c,r] += 1
        end
    end

    # mark the active elements, i.e. elements over which
    # at least one function is defined.
    active = falses(num_cells)
    index_among_actives = fill(0, num_cells)
    num_active_cells = 0
    for i in 1:num_cells
        if maximum(@view celltonum[i,:]) > 0
            num_active_cells += 1
            active[i] = true
            index_among_actives[i] = num_active_cells
        end
    end

    @assert num_active_cells != 0
    E = typeof(chart(geo, first(cells(geo))))
    elements = Vector{E}(num_active_cells)
    j = 1
    for (i,cell) in enumerate(cells(geo))
        active[i] || continue
        elements[j] = chart(geo, cell)
        j += 1
    end
    #elements = map(x->chart(geo,x[2]), filter(x->x[1], zip(active, cells(geo))))

    max_celltonum = maximum(celltonum)
    fill!(celltonum, 0)
    data = fill((0,zero(T)), max_celltonum, num_refs, num_active_cells)
    for b in 1 : num_bfs
        for shape in basisfunction(basis, b)
            c = shape.cellid
            l = index_among_actives[c]
            @assert 0 < l <= num_active_cells
            r = shape.refid
            w = shape.coeff
            k = (celltonum[c,r] += 1)
            data[k,r,l] = (b,w)
        end
    end

    return elements, AssemblyData(data)
end


text:
# K. Roosa and G. Chowell.
# Assessing parameter identifiability in compartmental dynamic models using a computational approach: application to infectious disease transmission models
# https://doi.org/10.1186/s12976-018-0097-6
# Model 2
using Logging

using StructuralIdentifiability

logger = Logging.SimpleLogger(stdout, Logging.Debug)
global_logger(logger)

ode = @ODEmodel(
    S'(t) = -b * S(t) * (I(t) + J(t) + q * A(t)) * Ninv,
    E'(t) = b * S(t) * (I(t) + J(t) + q * A(t)) * Ninv - k * E(t),
    A'(t) = k * (1 - r) * E(t) - g1 * A(t),
    I'(t) = k * r * E(t) - (alpha + g1) * I(t),
    J'(t) = alpha * I(t) - g2 * J(t),
    C'(t) = alpha * I(t),
    y(t) = C(t),
    y2(t) = Ninv
)

@time println(assess_global_identifiability(ode))


text:
using LinearAlgebra

g(x, y) = x + sin(x)*y^2

function g(x::Vector, y::Vector)
    z = similar(x)
    for i in 1:length(x)
        z[i] = g(x[i], y[i])
    end
    return z
end

x = randn(10); y = randn(10);

@time g(x,y)

@time g(x,y)

# Use another input vector of type Float32
z = randn(Float32, 10);

@time g(z,y)

## Memory Layout
using BenchmarkTools

function sum_by_row(x::AbstractMatrix)
    res = zero(eltype(x))
    for i in 1:size(x, 1)
        res += sum(x[i,:])
    end
    return res
end

function sum_by_col(x::AbstractMatrix)
    res = zero(eltype(x))
    for i in 1:size(x, 2)
        res += sum(x[:,i])
    end
    return res
end

x = randn(Float32, 10_000, 10_000);

@btime sum_by_row($x) # 713.33 ms (20000 allocations: 382.23 MiB)
@btime sum_by_col($x) # 99.842 ms (20000 allocations: 382.23 MiB)

## SIMD
function sum_vectors!(x, y, z)
    n = length(x)
    for i in 1:n
        x[i] = y[i] + z[i]
    end
end

function sum_vectors_simd!(x, y, z)
    n = length(x)
    @inbounds @simd for i in 1:n
        x[i] = y[i] + z[i]
    end
end

a = zeros(Float32, 100_000);
b = randn(Float32, 100_000);
c = randn(Float32, 100_000);

@btime sum_vectors!($a,$b,$c) # 79.932 μs (0 allocations: 0 bytes)
@btime sum_vectors_simd!($a,$b,$c) # 18.960 μs (0 allocations: 0 bytes)

## For the interested
@code_llvm sum_vectors_simd!(a,b,c)

## Multithreading
using Base.Threads

function sum_vectors_threading!(x, y, z)
    n = length(x)
    @inbounds Threads.@threads for i in 1:n
        x[i] = y[i] + z[i]
    end
end

@btime sum_vectors_threading!($a,$b,$c) # 132.019 μs (6 allocations: 560 bytes)