text:
# script to run inference for the OU SDEMEM model
using Pkg
using LinearAlgebra
using DataFrames

println("start run script for cpmmh")

# load functions
include(pwd()*"/src/SDEMEM OU neuron data/ou_sdemem.jl")
include(pwd()*"/src/SDEMEM OU neuron data/mcmc.jl")

M_subjects = 100

ρ = parse(Float64,ARGS[1])
nbr_particles = parse(Int,ARGS[2])*ones(Int64, M_subjects)
seed = parse(Int,ARGS[3])

y, prior_parameters_η, prior_parameters_σ_ϵ = set_up(M_subjects, seed)

# set number of obs for each data set
T_vec = zeros(Int64, M_subjects)
for i = 1:M_subjects; T_vec[i] = length(y[i].mV); end

# run MH-Gibbs

R = 100000 #15000 #10000
burn_in = 2000

job = string(M_subjects)*"_"*string(ρ)*"_"*string(nbr_particles[1])*"_"*"bridge"

startval_ϕ = ones(M_subjects,3)

for j = 1:3

    μ_0_j, M_0_j, α_j, β_j = prior_parameters_η[j,:]

    startval_ϕ[:,j,1] = μ_0_j*ones(M_subjects)

end

startval_η = zeros(6)

for j = 1:3

    # get paramaters
    μ_0_j, M_0_j, α_j, β_j = prior_parameters_η[j,:]

    startval_η[j] = μ_0_j # set start value to mean
    startval_η[j+3] = (α_j - 1)/β_j

end

startval_log_σ_ϵ = log(0.2)

# set time step
dt = diff(y[1].mSec)[1]

u_prop_old = [randn(nbr_particles[1], T_vec[1]+1)]
u_resample_old = [randn(T_vec[1], 2)]

for i in 2:M_subjects;
    append!(u_prop_old, [randn(nbr_particles[i],T_vec[i]+1)])
    append!(u_resample_old, [randn(T_vec[i],2)])
end

ll = @time bridge_pf(y, exp(startval_log_σ_ϵ), startval_ϕ, dt, u_prop_old, u_resample_old, nbr_particles, true)
sum(ll)

Σ_i_σ_ϵ = 1^2

Σ_i_ϕ = Matrix{Float64}[]
for i in 1:M_subjects; push!(Σ_i_ϕ, [0.004 0 0;0 0.04 0;0 0 0.004]); end

γ_ϕ_0 = 1.
γ_σ_ϵ_0 = 1.
μ_i_ϕ = startval_ϕ
μ_i_σ_ϵ = startval_log_σ_ϵ
α_star_ϕ = 0.25
α_star_σ_ϵ = 0.25
log_λ_i_ϕ = log.(2.4/sqrt(3)*ones(M_subjects))
log_λ_i_σ_ϵ = log(2.4)
update_interval = 1
α_power = 0.7
start_update = 100

# set random numbers
Random.seed!(seed)

run_time_cpmmh = @elapsed chain_ϕ_cpmmh, chain_log_σ_ϵ_cpmmh, chain_η_cpmmh, accept_vec_cpmmh, loglik_vec = gibbs_cpmmh(R,
                                                                                y,
                                                                                dt,
                                                                                Σ_i_σ_ϵ,
                                                                                Σ_i_ϕ,
                                                                                γ_ϕ_0,
                                                                                γ_σ_ϵ_0,
                                                                                μ_i_ϕ,
                                                                                μ_i_σ_ϵ,
                                                                                α_star_ϕ,
                                                                                α_star_σ_ϵ,
                                                                                log_λ_i_ϕ,
                                                                                log_λ_i_σ_ϵ,
                                                                                update_interval,
                                                                                start_update,
                                                                                α_power,
                                                                                startval_ϕ,
                                                                                startval_log_σ_ϵ,
                                                                                startval_η,
                                                                                prior_parameters_η,
                                                                                prior_parameters_σ_ϵ,
                                                                                nbr_particles,
                                                                                T_vec,
                                                                                ρ,
                                                                                bridge_pf);

println(run_time_cpmmh)
println(sum(accept_vec_cpmmh[1,:])/(M_subjects*R))
println(sum(accept_vec_cpmmh[2,:])/R)
println(sum(accept_vec_cpmmh[3,:])/(3*R))

sim_data = zeros(4,1)

sim_data[1] = run_time_cpmmh
sim_data[2] = sum(accept_vec_cpmmh[1,:])/(M_subjects*R)
sim_data[3] = sum(accept_vec_cpmmh[2,:])/R
sim_data[4] = sum(accept_vec_cpmmh[3,:])/(3*R)

# Save results

chain_ϕ_export = zeros(M_subjects*3, R)

idx = 0
for m = 1:M_subjects
    for j = 1:3
        global idx = idx + 1
        chain_ϕ_export[idx,:] = chain_ϕ_cpmmh[m,j,:]
    end
end

loglik_data = zeros(R,1)
loglik_data[:] = loglik_vec
DataFrame(loglik_data)

CSV.write("data/SDEMEM OU neuron data/cpmmh/sim_data_"*string(seed)*"_"*job*".csv", DataFrame(sim_data))
CSV.write("data/SDEMEM OU neuron data/cpmmh/chain_sigma_epsilon_"*string(seed)*"_"*job*".csv", DataFrame(chain_log_σ_ϵ_cpmmh'))
CSV.write("data/SDEMEM OU neuron data/cpmmh/chain_eta_"*string(seed)*"_"*job*".csv", DataFrame(chain_η_cpmmh'))
CSV.write("data/SDEMEM OU neuron data/cpmmh/chain_phi_"*string(seed)*"_"*job*".csv", DataFrame(chain_ϕ_export'))
CSV.write("data/SDEMEM OU neuron data/cpmmh/loglik_data_"*string(seed)*"_"*job*".csv", DataFrame(loglik_data))

max_stars_repo_path: src/SDEMEM OU neuron data/run_script_cpmmh_bridge.jl


text:
function savefig(figname, obj, pdir=mktempdir())
  println(pdir)
  fexts = [".png", ".pdf", ".svg", ".tex"]
  for ext in fexts
    isfile(figname * ext) && rm(figname * ext)
  end
  pgfsave( joinpath(pdir, figname * ".png"), obj)
  pgfsave( joinpath(pdir, figname * ".pdf"), obj)
  run(`pdf2svg $(joinpath(pdir, figname * ".pdf")) $(joinpath(pdir, figname * ".svg"))`)
  pgfsave(joinpath(pdir, figname * ".tex"), obj);
  return nothing
end

ext: jl


text:
using BenchmarkTools, FinancialMonteCarlo, DualNumbers

S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01;
D = 90.0;

Nsim = 10000;
Nstep = 30;
sigma = dual(0.2, 1.0);
mc = MonteCarloConfiguration(Nsim, Nstep);
toll = 1e-3;

rfCurve = ZeroRate(r);

FwdData = Forward(T * 2.0)
EUData = EuropeanOption(T * 1.1, K)
AMData = AmericanOption(T * 4.0, K)
BarrierData = BarrierOptionDownOut(T * 0.5, K, D)
AsianFloatingStrikeData = AsianFloatingStrikeOption(T * 3.3)
AsianFixedStrikeData = AsianFixedStrikeOption(T * 1.2, K)
Model = BlackScholesProcess(sigma, Underlying(S0, d));

optionDatas = [FwdData, EUData, AMData, BarrierData, AsianFloatingStrikeData, AsianFixedStrikeData]

@btime (FwdPrice, EuPrice, AMPrice, BarrierPrice, AsianPrice1, AsianPrice2) = pricer(Model, rfCurve, mc, optionDatas);


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct PolicySetDefinitionsApi <: SwaggerApi
    client::Swagger.Client
end

"""
Creates or updates a policy set definition.
This operation creates or updates a policy set definition in the given subscription with the given name.
Param: policySetDefinitionName::String (required)
Param: parameters::PolicySetDefinition (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PolicySetDefinition
"""
function _swaggerinternal_policySetDefinitionsCreateOrUpdate(_api::PolicySetDefinitionsApi, policySetDefinitionName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", PolicySetDefinition, "/subscriptions/{subscriptionId}/providers/Microsoft.Authorization/policySetDefinitions/{policySetDefinitionName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "policySetDefinitionName", policySetDefinitionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policySetDefinitionsCreateOrUpdate(_api::PolicySetDefinitionsApi, policySetDefinitionName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsCreateOrUpdate(_api, policySetDefinitionName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policySetDefinitionsCreateOrUpdate(_api::PolicySetDefinitionsApi, response_stream::Channel, policySetDefinitionName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsCreateOrUpdate(_api, policySetDefinitionName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Creates or updates a policy set definition.
This operation creates or updates a policy set definition in the given management group with the given name.
Param: policySetDefinitionName::String (required)
Param: parameters::PolicySetDefinition (required)
Param: api_version::String (required)
Param: managementGroupId::String (required)
Return: PolicySetDefinition
"""
function _swaggerinternal_policySetDefinitionsCreateOrUpdateAtManagementGroup(_api::PolicySetDefinitionsApi, policySetDefinitionName::String, parameters, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", PolicySetDefinition, "/providers/Microsoft.Management/managementgroups/{managementGroupId}/providers/Microsoft.Authorization/policySetDefinitions/{policySetDefinitionName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "policySetDefinitionName", policySetDefinitionName)  # type String
    Swagger.set_param(_ctx.path, "managementGroupId", managementGroupId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policySetDefinitionsCreateOrUpdateAtManagementGroup(_api::PolicySetDefinitionsApi, policySetDefinitionName::String, parameters, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsCreateOrUpdateAtManagementGroup(_api, policySetDefinitionName, parameters, api_version, managementGroupId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policySetDefinitionsCreateOrUpdateAtManagementGroup(_api::PolicySetDefinitionsApi, response_stream::Channel, policySetDefinitionName::String, parameters, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsCreateOrUpdateAtManagementGroup(_api, policySetDefinitionName, parameters, api_version, managementGroupId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Deletes a policy set definition.
This operation deletes the policy set definition in the given subscription with the given name.
Param: policySetDefinitionName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_policySetDefinitionsDelete(_api::PolicySetDefinitionsApi, policySetDefinitionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/providers/Microsoft.Authorization/policySetDefinitions/{policySetDefinitionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "policySetDefinitionName", policySetDefinitionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policySetDefinitionsDelete(_api::PolicySetDefinitionsApi, policySetDefinitionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsDelete(_api, policySetDefinitionName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policySetDefinitionsDelete(_api::PolicySetDefinitionsApi, response_stream::Channel, policySetDefinitionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsDelete(_api, policySetDefinitionName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Deletes a policy set definition.
This operation deletes the policy set definition in the given management group with the given name.
Param: policySetDefinitionName::String (required)
Param: api_version::String (required)
Param: managementGroupId::String (required)
Return: Nothing
"""
function _swaggerinternal_policySetDefinitionsDeleteAtManagementGroup(_api::PolicySetDefinitionsApi, policySetDefinitionName::String, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/providers/Microsoft.Management/managementgroups/{managementGroupId}/providers/Microsoft.Authorization/policySetDefinitions/{policySetDefinitionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "policySetDefinitionName", policySetDefinitionName)  # type String
    Swagger.set_param(_ctx.path, "managementGroupId", managementGroupId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policySetDefinitionsDeleteAtManagementGroup(_api::PolicySetDefinitionsApi, policySetDefinitionName::String, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsDeleteAtManagementGroup(_api, policySetDefinitionName, api_version, managementGroupId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policySetDefinitionsDeleteAtManagementGroup(_api::PolicySetDefinitionsApi, response_stream::Channel, policySetDefinitionName::String, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsDeleteAtManagementGroup(_api, policySetDefinitionName, api_version, managementGroupId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Retrieves a policy set definition.
This operation retrieves the policy set definition in the given subscription with the given name.
Param: policySetDefinitionName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PolicySetDefinition
"""
function _swaggerinternal_policySetDefinitionsGet(_api::PolicySetDefinitionsApi, policySetDefinitionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PolicySetDefinition, "/subscriptions/{subscriptionId}/providers/Microsoft.Authorization/policySetDefinitions/{policySetDefinitionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "policySetDefinitionName", policySetDefinitionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policySetDefinitionsGet(_api::PolicySetDefinitionsApi, policySetDefinitionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsGet(_api, policySetDefinitionName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policySetDefinitionsGet(_api::PolicySetDefinitionsApi, response_stream::Channel, policySetDefinitionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsGet(_api, policySetDefinitionName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Retrieves a policy set definition.
This operation retrieves the policy set definition in the given management group with the given name.
Param: policySetDefinitionName::String (required)
Param: api_version::String (required)
Param: managementGroupId::String (required)
Return: PolicySetDefinition
"""
function _swaggerinternal_policySetDefinitionsGetAtManagementGroup(_api::PolicySetDefinitionsApi, policySetDefinitionName::String, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PolicySetDefinition, "/providers/Microsoft.Management/managementgroups/{managementGroupId}/providers/Microsoft.Authorization/policySetDefinitions/{policySetDefinitionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "policySetDefinitionName", policySetDefinitionName)  # type String
    Swagger.set_param(_ctx.path, "managementGroupId", managementGroupId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policySetDefinitionsGetAtManagementGroup(_api::PolicySetDefinitionsApi, policySetDefinitionName::String, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsGetAtManagementGroup(_api, policySetDefinitionName, api_version, managementGroupId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policySetDefinitionsGetAtManagementGroup(_api::PolicySetDefinitionsApi, response_stream::Channel, policySetDefinitionName::String, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsGetAtManagementGroup(_api, policySetDefinitionName, api_version, managementGroupId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Retrieves a built in policy set definition.
This operation retrieves the built-in policy set definition with the given name.
Param: policySetDefinitionName::String (required)
Param: api_version::String (required)
Return: PolicySetDefinition
"""
function _swaggerinternal_policySetDefinitionsGetBuiltIn(_api::PolicySetDefinitionsApi, policySetDefinitionName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PolicySetDefinition, "/providers/Microsoft.Authorization/policySetDefinitions/{policySetDefinitionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "policySetDefinitionName", policySetDefinitionName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policySetDefinitionsGetBuiltIn(_api::PolicySetDefinitionsApi, policySetDefinitionName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsGetBuiltIn(_api, policySetDefinitionName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policySetDefinitionsGetBuiltIn(_api::PolicySetDefinitionsApi, response_stream::Channel, policySetDefinitionName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsGetBuiltIn(_api, policySetDefinitionName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Retrieves the policy set definitions for a subscription.
This operation retrieves a list of all the policy set definitions in the given subscription.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PolicySetDefinitionListResult
"""
function _swaggerinternal_policySetDefinitionsList(_api::PolicySetDefinitionsApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PolicySetDefinitionListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Authorization/policySetDefinitions", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policySetDefinitionsList(_api::PolicySetDefinitionsApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsList(_api, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policySetDefinitionsList(_api::PolicySetDefinitionsApi, response_stream::Channel, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsList(_api, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Retrieves built-in policy set definitions.
This operation retrieves a list of all the built-in policy set definitions.
Param: api_version::String (required)
Return: PolicySetDefinitionListResult
"""
function _swaggerinternal_policySetDefinitionsListBuiltIn(_api::PolicySetDefinitionsApi, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PolicySetDefinitionListResult, "/providers/Microsoft.Authorization/policySetDefinitions", ["azure_auth"])
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policySetDefinitionsListBuiltIn(_api::PolicySetDefinitionsApi, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsListBuiltIn(_api, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policySetDefinitionsListBuiltIn(_api::PolicySetDefinitionsApi, response_stream::Channel, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsListBuiltIn(_api, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Retrieves all policy set definitions in management group.
This operation retrieves a list of all the a policy set definition in the given management group.
Param: api_version::String (required)
Param: managementGroupId::String (required)
Return: PolicySetDefinitionListResult
"""
function _swaggerinternal_policySetDefinitionsListByManagementGroup(_api::PolicySetDefinitionsApi, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PolicySetDefinitionListResult, "/providers/Microsoft.Management/managementgroups/{managementGroupId}/providers/Microsoft.Authorization/policySetDefinitions", ["azure_auth"])
    Swagger.set_param(_ctx.path, "managementGroupId", managementGroupId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policySetDefinitionsListByManagementGroup(_api::PolicySetDefinitionsApi, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsListByManagementGroup(_api, api_version, managementGroupId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policySetDefinitionsListByManagementGroup(_api::PolicySetDefinitionsApi, response_stream::Channel, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policySetDefinitionsListByManagementGroup(_api, api_version, managementGroupId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export policySetDefinitionsCreateOrUpdate, policySetDefinitionsCreateOrUpdateAtManagementGroup, policySetDefinitionsDelete, policySetDefinitionsDeleteAtManagementGroup, policySetDefinitionsGet, policySetDefinitionsGetAtManagementGroup, policySetDefinitionsGetBuiltIn, policySetDefinitionsList, policySetDefinitionsListBuiltIn, policySetDefinitionsListByManagementGroup


text:
# Defines the Column type and associated methods. Consistent
# with the rest of this package a column is a sequence of values
# from Pascal's triangle (n, k) where n ranges from k upwards.

abstract type AbstractColumn{V} <: AbstractVector{V} end
Base.IteratorEltype(::Type{<:AbstractColumn}) = Base.HasEltype()
Base.IndexStyle(::Type{<:AbstractColumn}) = IndexLinear()
isatleft(c::AbstractColumn) = isfirst(c)

"""
    Column(colnum, data)
    Column(colnum, datasize)
    Column(V::Type, colnum, datasize)
    Column(c::Column)

Represents a column of Pascal's triangle. That is, the top-most
`datasize` elements of `binomial(n,colnum)` where `n ≥ colnum`,
but this is calculated more efficiently. As with all types in this
package, columns are numbered from 0.

# Examples
```jldoctest
julia> Column(2, 3)
3-element Column{Integer}:
 1
 3
 6

julia> Column(Float64, 2, 3)
3-element Column{Float64}:
 1.0
 3.0
 6.0
```
"""
mutable struct Column{V} <: AbstractColumn{V}
    colnum::Int
    data::Array{V,1}
    function Column(colnum::Integer, data::Array)
        colnum ≥ 0 || throw(DomainError("colnum must be nonnegative"))
        return new{eltype(data)}(colnum, data)
    end
end
Column(c::Column) = Column(c.colnum, deepcopy(c.data))
function Column{V}(colnum, datasize) where {V <: Real}
    datasize ≥ 0 || throw(DomainError("datasize must be nonnegative"))
    arr = ones(V,datasize)
    entry = Entry(colnum,colnum,one(V))
    for i ∈ 2:datasize
        down!(entry)
        arr[i] = entry.val
    end
    return Column(colnum, arr)
end
Column(colnum::Integer, datasize::Integer) = Column{Int}(colnum, datasize)

colnumber(c::Column) = c.colnum
value(c::Column) = c.data
Base.values(c::Column) = value(c)

Base.length(c::Column) = length(c.data)
Base.size(c::Column) = (length(c),)
Base.getindex(c::Column, i::Int) = c.data[i-c.colnum+1]
Base.firstindex(c::Column) = c.colnum
Base.lastindex(c::Column) = c.colnum + length(c.data) - 1
Base.eltype(c::Column) = eltype(c.data)
Base.axes(c::Column) = (firstindex(c):lastindex(c),)

isfirst(c::Column) = c.colnum == 0
isvalid(c::Column) = c.colnum ≥ 0 && c.data == Column(c.colnum, length(c.data)).data
toarray(c::Column) = Entry.(firstindex(c):lastindex(c), c.colnum, c.data)

mutable struct LazyColumn{V} <: AbstractColumn{V}
    colnum::Int
    data::Dict{Int, V}
    function LazyColumn(colnum::Integer, data::Dict)
        colnum ≥ 0 || throw(DomainError("colnum must be nonnegative"))
        return new{valtype(data)}(colnum, data)
    end
end
LazyColumn(c::LazyColumn) = LazyColumn(c.colnum, copy(c.data))
LazyColumn(c::Column) = LazyColumn(c.colnum, 1:length(c) .=> c.data)
LazyColumn{V}(colnum) where {V <: Real} = LazyColumn(colnum, Dict{Int,V}(1 => 1))
LazyColumn(colnum::Integer) = LazyColumn{Int}(colnum)

function Column(c::LazyColumn)
    data = eltype(c)
    i = c.colnum
    while haskey(c.data, i)
        push!(data, c.data[i])
        i += 1
    end
    return Column(c.colnum, data)
end

colnumber(c::LazyColumn) = c.colnum
function Base.getindex(c::LazyColumn, i::Int)
    c.colnum ≤ i || throw(BoundsError(c, i))
    c.colnum == 0 && return one(valtype(c.data))
    j = i - c.colnum + 1
    if haskey(c.data, j)
        return c.data[j]
    end
    e = Entry{valtype(c.data)}(i, c.colnum)
    uprange = (j-1):-1:max(1,j-PRECALC_NUMBER)
    a = Entry(e)
    for k ∈ uprange
        up!(a)
        if !haskey(c.data, k)
            c.data[k] = a.val
        end
    end
    downrange = (j+1):(j+PRECALC_NUMBER)
    a = Entry(e)
    for k ∈ downrange
        down!(a)
        if !haskey(c.data, k)
            c.data[k] = a.val
        end
    end
    return e.val
end
Base.getindex(c::LazyColumn, I) = [c[i] for i ∈ I]
Base.firstindex(c::LazyColumn) = c.colnum
Base.IteratorSize(::LazyColumn) = Base.IsInfinite()
Base.eltype(::LazyColumn{V}) where {V <: Real} = V
Base.keys(c::LazyColumn) = (c.colnum - 1) .+ keys(c.data)

function Base.show(io::IO, ::MIME"text/plain", c::LazyColumn)
    print(io, "$(typeof(c).name){$(valtype(c.data).name)}")
    print(io, "($(c.colnum))\n")
    print(io, "($(length(c.data)) elements precalculated, ")
    print(io, "up to row $(maximum(keys(c.data))+c.colnum-1))")
end
function Base.show(io::IO, c::LazyColumn)
    print(io, "$(typeof(c).name)")
    print(io, "($(c.colnum))")
end

Base.:(==)(a::LazyColumn, b::LazyColumn) = a.colnum == b.colnum &&
    isvalid(a) && isvalid(b)

isfirst(c::LazyColumn) = c.colnum == 0
function isvalid(c::LazyColumn)
    c.colnum < 0 && return false
    for (r, val) ∈ c.data
        if binomial(big(c.colnum + r - 1), big(c.colnum)) ≠ val
            return false
        end
    end
    return true
end

toarray(c::LazyColumn) = Entry.(keys(c), c.colnum, values(c.data))


text:
######### Stan program example  ###########

using StanSample, MCMCChains

eightschools ="
data {
  int<lower=0> J; // number of schools
  real y[J]; // estimated treatment effects
  real<lower=0> sigma[J]; // s.e. of effect estimates
}
parameters {
  real mu;
  real<lower=0> tau;
  real eta[J];
}
transformed parameters {
  real theta[J];
  for (j in 1:J)
    theta[j] <- mu + tau * eta[j];
}
model {
  eta ~ normal(0, 1);
  y ~ normal(theta, sigma);
}
"

schools8data = Dict("J" => 8,
  "y" => [28,  8, -3,  7, -1,  1, 18, 12],
  "sigma" => [15, 10, 16, 11,  9, 11, 10, 18],
  "tau" => 25
)

sm = SampleModel("schools8", eightschools)

rc = stan_sample(sm, data=schools8data)

if success(rc)
  chns = read_samples(sm, :mcmcchains)
end


text:
#
# This file contains functions that deal with physical ingredients
# in the model.
#

# Some coding conventions:
#
#  model : a solution of the system of equations.
#

using Unitful, UnitfulRecipes, UnitfulAstro
using PhysicalConstants.CODATA2018
G = NewtonianConstantOfGravitation

#
# Utility functions. These calculate quantities in units of some scale.
# The advantage is that the code is agnostic to whether the user provided
# units for these quantities or now, but of course bad errors can happen as
# no error checks are made.
#
"Temperature in units of 10^4 K"
T4(T) = ustrip(T/(1e4*u"K"))
"The number of ionizing photons in units of 10^{49} s^{-1}"
Q0p49(Q0) = ustrip(Q0/(1e49*u"1/s"))
"The mean energy of the ionizing photons in units of 18 eV"
hnui18(hnui) = ustrip(hnui/(18.0*u"eV"))
"The dust cross section in units of 10^{-21} cm^{2}"
σ21(σ) = ustrip(σ/(1e-21*u"1/cm^2"))

"The Case B recombination coefficient from Draine 2011"
alpha_B(T) = 2.56e-13*T4(T)^(-0.83)*u"cm^3/s"

#
# Scaling parameters of the model
#

"The characteristic density used in defining u = n0/n - Eq 4 in D11"
n0(;T=1e4*u"K", Q0=1e49*u"1/s", hnui=18*u"eV", σ=nothing) = 4.54e5*(T4(T)^4.66/Q0p49(Q0))*hnui18(hnui)^(-3)*u"1/cm^3"

"The characteristic length scale used to define y = r/λ0 - Eq 5 in D11"
λ0(;T=1e4*u"K", Q0=1e49*u"1/s", hnui=18*u"eV", σ=nothing) = 2.47e16*(Q0p49(Q0)/T4(T)^2.83)*hnui18(hnui)^2*u"cm"

"The γ parameter used in the equation definition"
gamma(;T=1e4*u"K", Q0=nothing, hnui=18*u"eV", σ=1e-21*u"1/cm^2") = 11.2*T4(T)^1.83*σ21(σ)/hnui18(hnui)

"The characteristic mass"
M0(T; kwargs...) = 4π*λ0(;T=T4(T), kwargs...)^3*n0(;T=T4(T), kwargs...)*ProtonMass/3
"The characteristic gravitational energy"
Eg(T; kwargs...) = G*ProtonMass*M0(T4(T); kwargs...)/λ0(;T=T4(T), kwargs...)
"The characteristic thermal energy"
Et(T) = 2*BoltzmannConstant*T4(T)*1e4u"K"/3

#
# To-from (n, r) <-> (u, y)
#
"Convert density to u"
uval(n; kwargs...) = n0(; kwargs...)/n
"Calculate radius from y"
radius(y; kwargs...) = λ0(; kwargs...)*y
"Calculate y from radius"
yval(r; kwargs...) = r/λ0(; kwargs...)
"Calculate density from the solution"
density(y, model::AbstractDustyHII) = n0(;kwargs...)/ufunc(y, model)

#
# Physical properties of the solution
#

"RMS density given a solution - eq 12 in D11"
function n_rms(model; kwargs...)
    integrand_rms(y) = y^2/ufunc(y, model)^2
    int, dint = quadgk(integrand_rms, model.ymin, model.ymax)
    return n0(;kwargs...)*sqrt(3.0*int/model.ymax^3)
end

"Mean density given a solution - eq 11 in D11"
function n_mean(model; kwargs...)
    integrand_mean(y) = y^2/ufunc(y, model)
    int, dint = quadgk(integrand_mean, model.ymin, model.ymax)
    return 3.0*n0(;kwargs...)*int/model.ymax^3
end

"The Stromgren radius for n(RMS)"
function Rs(sol; Q0=1e49u"1/s", T=1e4u"K", kwargs...)
    nrms_3 = n_rms(sol; Q0=Q0, T=T, kwargs...)/(1e3u"1/cm^3")
    return 2.10e18*Q0p49(Q0)^(1/3.)*T4(T)^0.25/nrms_3^(2/3.)*u"cm"
end

# RMS optical depth
taud0(sol; σ=1e-21*u"1/cm^2", kwargs...) = n_rms(sol; kwargs...)*Rs(sol; kwargs...)*σ

"Pressure on the edge of the model"
p_edge(model::AbstractDustyHII; T=1e4u"K", kwargs...) = 2.0*n0(;kwargs...)*Unitful.k*T/ufunc(ymax(model), model)

"The radius of the ionised zone"
Rion(sol;kwargs...) = ymax(sol)*λ0(;kwargs...)

"Fraction of absorbed E>13.6 eV photons absorbed by H"
fion(sol;kwargs...) = (Rion(sol; kwargs...)/Rs(sol; kwargs...))^2


text:
module Day13
export day13

TESTINPUT = """6,10
0,14
9,10
0,3
10,4
4,11
6,0
6,12
4,1
0,13
10,12
3,4
3,0
8,4
1,10
2,14
8,10
9,0

fold along y=7
fold along x=5"""

function readInput(inputfile="inputs/input13.txt")
    lines = Iterators.Stateful(eachline(inputfile))

    # parse dots
    inds = Vector{Int64}[]
    for line in lines
        length(line) == 0 && break # empty line separates coordinates from folds
        push!(inds, parse.(Int64, split(line, ",")))
    end

    # parse folds
    folds = []
    for line in lines
        axis = line[12] == 'x' ? 1 : 2
        values = parse(Int64, line[14:end])
        push!(folds, (axis, values)) # +1 to account for Julia arrays
    end
    return inds, folds
end

function day13((indices, folds) = readInput("inputs/input13.txt"); output="inputs/output13.txt")
    part1 = length(_fold!(indices, folds[1]))
    [_fold!(indices, fold) for fold in folds[2:end]]

    # visualize
    maxX = maximum(x->x[1], indices)
    maxY = maximum(x->x[2], indices)
    arr = fill('.', maxX+1, maxY+1)

    arr[[CartesianIndex(ind...) + CartesianIndex(1,1) for ind in indices]] .= '#'
    open(output, "w") do out
        println.(Ref(out), String.(eachcol(arr)))
    end

    return (part1, "See $output")
end

function _fold!(inds, (axis, coordinate))
    for ind in inds
        if ind[axis] > coordinate
            ind[axis] = 2*coordinate - ind[axis]
        end
    end
    unique!(inds)
end

end #module


text:
#=== 2d ===# # add moment tensor source #
function addmt!(wf::Union{elwf2d,acwf2d}, sou::MTsource, it::Int64)

    wf.vx[sou.BDnloc[1],sou.BDnloc[2]] = wf.vx[sou.BDnloc[1],sou.BDnloc[2]] + (sou.mt[1] + sou.mt[3])*sou.waveform[it]
    wf.vx[sou.BDnloc[1],sou.BDnloc[2]-1] = wf.vx[sou.BDnloc[1],sou.BDnloc[2]-1] - sou.mt[1]*sou.waveform[it]
    wf.vx[sou.BDnloc[1]-1,sou.BDnloc[2]] = wf.vx[sou.BDnloc[1]-1,sou.BDnloc[2]] - sou.mt[3]*sou.waveform[it]
    wf.vz[sou.BDnloc[1],sou.BDnloc[2]] = wf.vz[sou.BDnloc[1],sou.BDnloc[2]] + (sou.mt[2] + sou.mt[3])*sou.waveform[it]
    wf.vz[sou.BDnloc[1]-1,sou.BDnloc[2]] = wf.vz[sou.BDnloc[1]-1,sou.BDnloc[2]] - sou.mt[2]*sou.waveform[it]
    wf.vz[sou.BDnloc[1],sou.BDnloc[2]-1] = wf.vz[sou.BDnloc[1],sou.BDnloc[2]-1] + sou.mt[3]*sou.waveform[it]

end

#=== 2d ===# # add single force source #
function addsf!(wf::Union{elwf2d,acwf2d}, sou::SFsource, it::Int64)

    wf.vx[sou.BDnloc[1],sou.BDnloc[2]] = wf.vx[sou.BDnloc[1],sou.BDnloc[2]] + sou.coeff[2]*sou.waveform[it]
    wf.vz[sou.BDnloc[1],sou.BDnloc[2]] = wf.vz[sou.BDnloc[1],sou.BDnloc[2]] + sou.coeff[1]*sou.waveform[it]

end

#=== 3d ===#
function addmt!(wf::Union{elwf3d,acwf3d}, sou::MTsource, it::Int64)

    wf.vx[sou.BDnloc[1],sou.BDnloc[2],sou.BDnloc[3]] = wf.vx[sou.BDnloc[1],sou.BDnloc[2],sou.BDnloc[3]] + (sou.mt[1] + sou.mt[4] + sou.mt[5])*sou.waveform[it]
    wf.vx[sou.BDnloc[1],sou.BDnloc[2]-1,sou.BDnloc[3]] = wf.vx[sou.BDnloc[1],sou.BDnloc[2]-1,sou.BDnloc[3]] - sou.mt[1]*sou.waveform[it]
    wf.vx[sou.BDnloc[1],sou.BDnloc[2],sou.BDnloc[3]-1] = wf.vx[sou.BDnloc[1],sou.BDnloc[2],sou.BDnloc[3]-1] - sou.mt[4]*sou.waveform[it]
    wf.vx[sou.BDnloc[1]-1,sou.BDnloc[2],sou.BDnloc[3]] = wf.vx[sou.BDnloc[1]-1,sou.BDnloc[2],sou.BDnloc[3]] - sou.mt[5]*sou.waveform[it]

    wf.vy[sou.BDnloc[1],sou.BDnloc[2],sou.BDnloc[3]] = wf.vy[sou.BDnloc[1],sou.BDnloc[2],sou.BDnloc[3]] + (sou.mt[2] + sou.mt[4] + sou.mt[6])*sou.waveform[it]
    wf.vy[sou.BDnloc[1],sou.BDnloc[2],sou.BDnloc[3]-1] = wf.vy[sou.BDnloc[1],sou.BDnloc[2],sou.BDnloc[3]-1] - sou.mt[2]*sou.waveform[it]
    wf.vy[sou.BDnloc[1],sou.BDnloc[2]-1,sou.BDnloc[3]] = wf.vy[sou.BDnloc[1],sou.BDnloc[2]-1,sou.BDnloc[3]] - sou.mt[4]*sou.waveform[it]
    wf.vy[sou.BDnloc[1]-1,sou.BDnloc[2],sou.BDnloc[3]] = wf.vy[sou.BDnloc[1]-1,sou.BDnloc[2],sou.BDnloc[3]] - sou.mt[6]*sou.waveform[it]

    wf.vz[sou.BDnloc[1],sou.BDnloc[2],sou.BDnloc[3]] = wf.vz[sou.BDnloc[1],sou.BDnloc[2],sou.BDnloc[3]] + (sou.mt[3] + sou.mt[5] + sou.mt[6])*sou.waveform[it]
    wf.vz[sou.BDnloc[1]-1,sou.BDnloc[2],sou.BDnloc[3]] = wf.vz[sou.BDnloc[1]-1,sou.BDnloc[2],sou.BDnloc[3]] - sou.mt[2]*sou.waveform[it]
    wf.vz[sou.BDnloc[1],sou.BDnloc[2]-1,sou.BDnloc[3]] = wf.vz[sou.BDnloc[1],sou.BDnloc[2]-1,sou.BDnloc[3]] - sou.mt[5]*sou.waveform[it]
    wf.vz[sou.BDnloc[1],sou.BDnloc[2],sou.BDnloc[3]-1] = wf.vz[sou.BDnloc[1],sou.BDnloc[2],sou.BDnloc[3]-1] - sou.mt[6]*sou.waveform[it]

end

#=== 3 ===# # add single force source #
function addsf!(wf::Union{elwf3d,acwf3d}, sou::SFsource, it::Int64)

    wf.vx[sou.BDnloc[1],sou.BDnloc[2],sou.BDnloc[3]] = wf.vx[sou.BDnloc[1],sou.BDnloc[2],sou.BDnloc[3]] + sou.coeff[2]*sou.waveform[it] #

    wf.vy[sou.BDnloc[1],sou.BDnloc[2],sou.BDnloc[3]] =  wf.vy[sou.BDnloc[1],sou.BDnloc[2],sou.BDnloc[3]] + sou.coeff[3]*sou.waveform[it] #

    wf.vz[sou.BDnloc[1],sou.BDnloc[2],sou.BDnloc[3]] = wf.vz[sou.BDnloc[1],sou.BDnloc[2],sou.BDnloc[3]] + sou.coeff[1]*sou.waveform[it] #

end


text:
summarize(A) = vec(mean(mapslices(x -> sum(x) , A, dims = 3), dims = 1))
summarize_diff(A) = vec(mean(diff(mapslices(sum, A, dims = 3)[:,:,1], dims = 2), dims = 1))

function cleanup!(A)
    for i in 1:size(A,1)
        for j in 2:size(A,2)
            if all(x -> ismissing(x), A[i,j,:])
                A[i,j,:] = A[i,j-1,:]
            end
        end
    end
    A
end

function ma(x::Vector{T}, k::Int) where T<:AbstractFloat
    n = length(x)
    y = fill(zero(T), n)
    for i in 1:n
        start = max(1, i - k)
        finish = min(n, i)
        y[i] = mean(x[start:finish])
    end
    return y[k:end]
end

function refresh_identifiers(n = 5)
    global pal =  Colors.distinguishable_colors(n, parse.(Colorant, ["dodgerblue", "green"]))
    global labels = ["Bulk crafting, current system",  "Playset-by-playset crafting, current system",
    "Playset-by-playset crafting, duplicate prevention", "Bulk crafting, duplicate prevention",
    "3 then 4 playset crafting, current system"]
end

function cumulativeplot(A)
    col = "#" * hex( popfirst!(pal) )
    label = popfirst!(labels)

    plot(summarize(A), label = label, color = col)
end

function diffplot(D, yline, w = 30)
    mD = ma(D,w)
    col = "#" * hex(popfirst!(pal))
    label = popfirst!(labels)
    plot(w:(w+length(mD)-1), mD, color = col, label = label)
    axhline(y = yline, linestyle = "dashed", color = col)
end

function ff1_dom(x)
    x["set"] == :DOM && x["rarity"] == 3 &&
    !in(x["name"], ["Chandra's Outburst", "Firesong and Sunspeaker", "Niambi, Faithful Healer"])
end

function ff1_grn(x)
    x["set"] == :GRN && x["rarity"] == 3 &&
    !in(x["name"], ["Ral's Dispersal", "Vraska's Stoneglare"])
end

function ff1_xln(x)
    x["set"] == :XLN && x["rarity"] == 3 &&
    !in(x["name"], ["Sun-Blessed Mount", "Grasping Current"])
end

function ff2_xln(x)
    in(x["name"],
    [
    "Search for Azcanta // Azcanta, the Sunken Ruin",
    "Vraska's Contempt",
    "Settle the Wreckage",
    "Drowned Catacomb",
    "Growing Rites of Itlimoc // Itlimoc, Cradle of the Sun",
    "Legion's Landing // Adanto, the First Fort",
    "Sunpetal Grove",
    "Vanquisher's Banner",
    "Glacial Fortress",
    "Dragonskull Summit",
    "Treasure Map // Treasure Cove",
    "Sorcerous Spyglass",
    "Primal Amulet // Primal Wellspring",
    "Hostage Taker",
    "Rootbound Crag",
    "Thaumatic Compass // Spires of Orazca",
    "Ripjaw Raptor",
    "Shapers' Sanctuary",
    "Arguel's Blood Fast // Temple of Aclazotz",
    "Deathgorge Scavenger",
    "Sunbird's Invocation",
    "Regisaur Alpha",
    "River's Rebuke"
    ]
    )
end

function ff2_dom(x)
    in(x["name"],
    [
    "Sulfur Falls",
    "Woodland Cemetery",
    "Goblin Chainwhirler",
    "Clifftop Retreat",
    "Shalai, Voice of Plenty",
    "Steel Leaf Champion",
    "Isolated Chapel",
    "Helm of the Host",
    "Hinterland Harbor",
    "Gilded Lotus",
    "Benalish Marshal",
    "Tempest Djinn"
    ]
    )
end

function ff2_grn(x)
    in(x["name"],
    [
    "Assassin's Trophy",
    "Steam Vents",
    "Risk Factor",
    "Watery Grave",
    "Sacred Foundry",
    "Overgrown Tomb",
    "Temple Garden",
    "Runaway Steam-Kin",
    "Knight of Autumn",
    "Ionize",
    "Experimental Frenzy",
    "Mission Briefing",
    "Chromatic Lantern",
    "Legion Warboss",
    "Pelt Collector",
    "Expansion // Explosion",
    "Thief of Sanity",
    "Niv-Mizzet, Parun",
    "Deafening Clarion",
    "Ritual of Soot",
    "Find // Finality",
    "Venerated Loxodon",
    "Beast Whisperer",
    "Tajic, Legion's Edge"
    ]
    )
end

function case(filterfun, setfun, sets; nsim = 200)
    pars = SimParameters(
    bonus_packs = Dict{Symbol,Int}(),
    track_collection_progress = true,
    max_track_progress_packs = 600,
    nextset = setfun,
    )

    par2 = deepcopy(pars)
    par2.prevent_duplicates = true

    rares = filter(filterfun, card_db)
    ix3 = findall(filterfun, card_db)

    deck = [map(x -> (x["name"], 4), rares)]
    deck2 = map(x -> [x], deck[1])
    deck3 = vcat(
        collect( [(n,3)] for (n,a) in deck[1] ),
        collect( [(n,4)] for (n,a) in deck[1] )
    );

    deck = deckinfo.(deck)
    deck3 = deckinfo.(deck3)
    deck2 = deckinfo.(deck2)

    output = Vector{ArenaSim.SimOutput}()

    push!(output, simulate(nsim, deck; parameters=pars, sets = sets))
    push!(output, simulate(nsim, deck2; parameters=pars, sets = sets))
    push!(output, simulate(nsim, deck2; parameters=par2, sets = sets))
    push!(output, simulate(nsim, deck; parameters=par2, sets = sets))
    push!(output, simulate(nsim, deck3; parameters=pars, sets = sets))

    cumulative = map(x -> cleanup!(x.collection_progress)[:,:,ix3], output)
    difference = map(x -> summarize_diff(x), cumulative)

    height = 6
    figure(figsize = (height*634/475,height))
    refresh_identifiers()
    for a in cumulative
        cumulativeplot(a)
    end
    legend(fontsize = 8)
    ymax = maximum(map(x -> maximum(summarize(x)), cumulative))
    yticks(0:4:ymax, fontsize = 7)
    grid()
    xlabel("Packs opened")
    ylabel("Average desired rares acquired")
    tight_layout()

    lines = map(x -> summarize(x), cumulative)
    lines = map(x -> last(x) / length(x), lines)

    figure(figsize = (height*634/475,height))
    refresh_identifiers()
    for (i,d) in enumerate(difference)
        diffplot(d, lines[i])
    end
    legend(fontsize=7)
    xlabel("Packs opened")
    ylabel("Average desired rares per pack in previous 30 packs")
    tight_layout()
end

max_stars_repo_name: asieev/VaultAnalysis

max_issues_repo_name: asieev/VaultAnalysis

max_forks_repo_name: asieev/VaultAnalysis


text:
# Handers for execute_request and related messages, which are
# the core of the Jupyter protocol: execution of Julia code and
# returning results.

import Base.Libc: flush_cstdio
import Pkg

# global variable so that display can be done in the correct Msg context
execute_msg = Msg(["julia"], Dict("username"=>"jlkernel", "session"=>uuid4()), Dict())
# global variable tracking the number of bytes written in the current execution
# request
const stdio_bytes = Ref(0)

import REPL: helpmode

# use a global array to accumulate "payloads" for the execute_reply message
const execute_payloads = Dict[]

# An array of macros to run on the contents of each cell
const cell_macros = []

function run_cell_code(code)
    do_auto_macros = true

    # Check for cell macros
    cell_macro_calls = []
    line_no = 1
    if startswith(code, "@@")
        lines = split(code, '\n')
        for line in lines
            line = strip(line)
            if length(line) <= 0
                line_no += 1
                continue
            end
            startswith(line, "@@") || break
            mac_call = Meta.parse(line[nextind("abc", 0, 2):end])
            @assert Meta.isexpr(mac_call, :macrocall)
            if mac_call.args[1] == Symbol("@noauto")
                do_auto_macros = false
            else
                push!(cell_macro_calls, mac_call)
            end
            line_no += 1
        end
        code = join((@view lines[line_no:end]), '\n')
    end

    # Apply macros to ast
    mod = current_module[]
    file = "In[$n]"

    ast = parse_file(code, file, false)
    # More compatible?
    #ast = Meta.parse("begin\n$(code)\nend", raise=false)
    ## Make block top level
    #if Meta.isexpr(ast, :block)
    #    ast.head = :toplevel
    #    line_node = LineNumberNode(line_no, file)
    #    push!(ast.args, 1, line_node)
    #end

    for mac_call in Iterators.reverse(cell_macro_calls)
        push!(mac_call.args, ast)
        ast = macroexpand(mod, mac_call, recursive=false)
    end
    if do_auto_macros
        if SOFTSCOPE[]
            ast = _apply_macro(var"@softscope", ast, mod)
        end
        for mac in Iterators.reverse(cell_macros)
            ast = _apply_macro(mac, ast, mod)
        end
    end

    # Run
    Base.eval(mod, ast)
end

function _apply_macro(mac, ast, mod)
    macro_expr = :(@macro_placehoder $ast)
    @assert Meta.isexpr(macro_expr, :macrocall)
    macro_expr.args[1] = mac
    macroexpand(mod, macro_expr, recursive=false)
end

"""
    parse_file(content, fname; raise=true, depwarn=true)

Parse the entire string as a file, reading multiple expressions.  Equivalent to
`Meta.parse()` but for more than one expression.
"""
function parse_file(content::AbstractString, fname::AbstractString,
                    raise::Bool=true, depwarn::Bool=true)
    # returns (expr, end_pos). expr is () in case of parse error.
    bcontent = String(content)
    bfname = String(fname)
    # For now, assume all parser warnings are depwarns
    ex = Meta.with_logger(depwarn ? Meta.current_logger() : Meta.NullLogger()) do
        ccall(:jl_parse_all, Any,
              (Ptr{UInt8}, Csize_t, Ptr{UInt8}, Csize_t),
              bcontent, sizeof(bcontent), fname, sizeof(fname))
    end
    if raise && Meta.isexpr(ex, :error)
        throw(ParseError(ex.args[1]))
    end
    ex
end

function execute_request(socket, msg)
    code = msg.content["code"]
    @vprintln("EXECUTING ", code)
    global execute_msg = msg
    global n, In, Out, ans
    stdio_bytes[] = 0
    silent = msg.content["silent"]
    store_history = get(msg.content, "store_history", !silent)
    empty!(execute_payloads)

    if !silent
        n += 1
        send_ipython(publish[],
                     msg_pub(msg, "execute_input",
                             Dict("execution_count" => n,
                                          "code" => code)))
    end

    silent = silent || REPL.ends_with_semicolon(code)
    if store_history
        In[n] = code
    end

    # "; ..." cells are interpreted as shell commands for run
    code = replace(code, r"^\s*;.*$" =>
                   m -> string(replace(m, r"^\s*;" => "Base.repl_cmd(`"),
                               "`, stdout)"))

    # "] ..." cells are interpreted as pkg shell commands
    if occursin(r"^\].*$", code)
        code = "IJulia.Pkg.REPLMode.do_cmd(IJulia.minirepl[], \"" *
            escape_string(code[2:end]) * "\"; do_rethrow=true)"
    end

    # a cell beginning with "? ..." is interpreted as a help request
    hcode = replace(code, r"^\s*\?" => "")

    try
        for hook in preexecute_hooks
            invokelatest(hook)
        end

        ans = result = if hcode != code # help request
            Core.eval(Main, helpmode(hcode))
        else
            #run the code!
            occursin(magics_regex, code) && match(magics_regex, code).offset == 1 ? magics_help(code) :
                run_cell_code(code)
        end

        if silent
            result = nothing
        elseif result !== nothing
            if store_history
                Out[n] = result
            end
        end

        user_expressions = Dict()
        for (v,ex) in msg.content["user_expressions"]
            try
                value = include_string(current_module[], ex)
                # Like the IPython reference implementation, we return
                # something that looks like a `display_data` but also has a
                # `status` field:
                # https://github.com/ipython/ipython/blob/master/IPython/core/interactiveshell.py#L2609-L2614
                user_expressions[v] = Dict("status" => "ok",
                                           "data" => display_dict(value),
                                           "metadata" => metadata(value))
            catch e
                # The format of user_expressions[v] is like `error` except that
                # it also has a `status` field:
                # https://jupyter-client.readthedocs.io/en/stable/messaging.html#execution-errors
                user_expressions[v] = Dict("status" => "error",
                                           error_content(e, catch_backtrace())...)
            end
        end

        for hook in postexecute_hooks
            invokelatest(hook)
        end

        # flush pending stdio
        flush_all()

        undisplay(result) # dequeue if needed, since we display result in pyout
        invokelatest(display) # flush pending display requests

        if result !== nothing
            result_metadata = invokelatest(metadata, result)
            result_data = invokelatest(display_dict, result)
            send_ipython(publish[],
                         msg_pub(msg, "execute_result",
                                 Dict("execution_count" => n,
                                              "metadata" => result_metadata,
                                              "data" => result_data)))

        end
        send_ipython(requests[],
                     msg_reply(msg, "execute_reply",
                               Dict("status" => "ok",
                                            "payload" => execute_payloads,
                                            "execution_count" => n,
                                            "user_expressions" => user_expressions)))
        empty!(execute_payloads)
    catch e
        bt = catch_backtrace()
        try
            # flush pending stdio
            flush_all()
            for hook in posterror_hooks
                invokelatest(hook)
            end
        catch
        end
        empty!(displayqueue) # discard pending display requests on an error
        content = error_content(e,bt)
        send_ipython(publish[], msg_pub(msg, "error", content))
        content["status"] = "error"
        content["execution_count"] = n
        send_ipython(requests[], msg_reply(msg, "execute_reply", content))
    end
end


text:
using DiffEqBase, OrdinaryDiffEq
using Plots; plotly()
############### Parameters

#=
cₕ = 15.0
cₖ = 10.0
κₕ = .2
κₖ = 7.0
dₕ = 1.0
dₖ = 1.0
nₕ = 2
nₖ = 2
m = 2
=#

cₕ = 3
cₖ = 5
κₕ = 1
κₖ = 1
dₕ = 0.4
dₖ = 0.4
nₕ = 2
nₖ = 2
m = 2
aₕ = 0.05
aₖ = 0.05

#=
const cₕ = 3
const cₖ = 5
const κₕ = 1
const κₖ = 1
const dₕ = 0.4
const dₖ = 0.4
const nₕ = 2
const nₖ = 2
const m = 2
const aₕ = 0.05
const aₖ = 0.05
=#

################# Setup the problem

type HoxKrox
  RAin::Float64
end
function (p::HoxKrox)(du,u,p,t)
  gₕ,gₖ = u
  du[1] .= (cₕ.*gₕ.^nₕ .+ (κₕ*p.RAin).^m) ./ (1 .+ cₕ.*gₕ.^nₕ .+ cₖ.*gₖ.^nₖ .+ (κₕ*p.RAin).^m) - dₕ.*gₕ
  du[2] .= (cₖ.*gₖ.^nₖ .+ (κₖ*p.RAin).^m) ./ (1 .+ cₕ.*gₕ.^nₕ .+ cₖ.*gₖ.^nₖ .+ (κₖ*p.RAin).^m) - dₖ.*gₖ
end
u0 = [0.5,2.0]
hk = HoxKrox(2.0)
tspan = (0.0,500.0)
prob = ODEProblem(hk,u0,tspan)

####################### Setup Bif values

SShox = Vector{Float64}()
SSkrox = Vector{Float64}()
SShox2 = Vector{Float64}()
SSkrox2 = Vector{Float64}()
####################### Get steady states at different RAin

RAins = 0:0.01:2.5

for RAin_scalar in RAins
  prob.f.RAin = RAin_scalar
  sol = solve(prob,Tsit5(),save_everystep=false)
  push!(SShox,sol[end][1])
  push!(SSkrox,sol[end][2])
end

prob.u0 = [2.0,0.5]

for RAin_scalar in RAins
  prob.f.RAin = RAin_scalar
  sol = solve(prob,Tsit5(),save_everystep=false)
  push!(SShox2,sol[end][1])
  push!(SSkrox2,sol[end][2])
end

p1 = plot(RAins,[SShox,SSkrox],title="Steady State Krox Start",xaxis="RAin",
     yaxis="Steady State Value",label=["Hox" "Krox"],
     lw = 3);

p2 = plot(RAins,[SShox2,SSkrox2],title="Steady State Hox Start",xaxis="RAin",
     yaxis="Steady State Value",label=["Hox" "Krox"],
     lw = 3);
display(plot(p1,p2))

#=
####################### Get steady states at different RAin

RAin = ra_sol[end][:,:,2][:,21:61]
x = X[5,21:61]
RAins = RAin[5,:]
SShox = Vector{Float64}()
SSkrox = Vector{Float64}()
SShox2 = Vector{Float64}()
SSkrox2 = Vector{Float64}()

prob.u0 = [0.5,2.0]

for RAin_scalar in RAins
  prob.f.RAin = RAin_scalar
  sol = solve(prob,Tsit5(),save_everystep=false)
  push!(SShox,sol[end][1])
  push!(SSkrox,sol[end][2])
end

p21 = plot(x,[SShox,SSkrox],title="Steady State Hox/Krox vs RAin",xaxis="x (Anterior -> Posterior)",
     yaxis="Steady State Value",label=["Hox" "Krox"],
     lw = 3)

prob.u0 = [2.0,0.5]

for RAin_scalar in RAins
  prob.f.RAin = RAin_scalar
  sol = solve(prob,Tsit5(),save_everystep=false)
  push!(SShox2,sol[end][1])
  push!(SSkrox2,sol[end][2])
end

p22= plot(x,[SShox2,SSkrox2],title="Steady State Hox/Krox vs RAin",xaxis="x (Anterior -> Posterior)",
     yaxis="Steady State Value",label=["Hox" "Krox"],
     lw = 3)

display(plot(p21,p22))
=#

max_stars_repo_name: ChrisRackauckas/RAMeanVar

max_issues_repo_name: ChrisRackauckas/RAMeanVar

max_forks_repo_name: ChrisRackauckas/RAMeanVar


text:
mutable struct CustomObjective <: AbstractObjective
    name::String
    sense::String
    fun::Function
    args::NamedTuple # ex: (a = [0, 0, 0], b = "hello")
    CustomObjective(name, sense, fun, args) = new(name, sense, fun, args)
    CustomObjective() =
        new("nothing fun", "max", () -> nothing, (default_arg_1 = 0, default_arg_2 = 0))
end


text:
using Distributions, Random, Test

## Helper: check a decision by simulation
function check_by_sim(d :: ExtremeValueDecision{T},  value_ixM :: Array{T,2},
    prob_ixM :: Array{T,2},  eVal_iV :: Array{T,1})    where T <: AbstractFloat

    rng = MersenneTwister(123);
    nTypes, n = size(value_ixM);
    success = true;

    for iType = 1 : nTypes
        success1 = check_one_by_sim(d, vec(value_ixM[iType,:]),
            vec(prob_ixM[iType,:]),  eVal_iV[iType],  rng);
        success = success && success1;
    end

    return success
end

## Check one type by simulation
# Returns a Bool (success or failure)
function check_one_by_sim(d :: ExtremeValueDecision{F},  valueV :: Vector{F},
    probV :: Vector{F},  eVal :: F,  rng :: AbstractRNG)  where  F <: AbstractFloat

    # Random vars for all alternatives
    n = length(valueV);
    nSim = Int64(5e5);
    simProbV, simVal = sim_extreme_value(d, valueV, nSim, rng);

    success = true;
    if !isapprox(simProbV, probV, atol = 5e-3)
       println("True probs:  ", probV);
       println("Sim probs:   ", simProbV);
       success = false;
    end

    # Check expected value
    if abs(simVal ./ eVal - 1.0) > 5e-3
       @warn "Discrepancy in values (demeaned = $(demeaned(d))): $simVal,  $eVal"
       success = false;
    end

    return success
end

function sim_extreme_value(d :: ExtremeValueDecision{F},  valueV :: Vector{F},
    nSim :: Integer, rng :: AbstractRNG) where F

    # Random vars for all alternatives
    n = length(valueV);
    maxV = zeros(F, nSim);
    countV = zeros(Int, n);
    for i1 = 1 : nSim
        randV = draw_gumbel_shocks(rng, pref_scale(d), n; demeaned = demeaned(d));
        # Shocks must be added, not subtracted.
        maxV[i1], maxIdx = findmax(valueV .+ randV);
        countV[maxIdx] += 1;
    end
    simProbV = countV ./ nSim;
    simVal = mean(maxV);
    return simProbV, simVal
end

function one_option_test()
    @testset "One option" begin
        prefScale = 1.4;
        nTypes = 3;
        value_iV = collect(range(32.0, 48.0, length = nTypes));

        # Demeaned
        d = ExtremeValueDecision(prefScale, true, true);
        prob_iV, eVal_iV = extreme_value_decision(d, value_iV);
        @test prob_iV == ones(Float64, nTypes)
        @test isapprox(eVal_iV, value_iV, rtol = 0.02)

        # Not demeaned
        d2 = ExtremeValueDecision(prefScale, false, true);
        prob_iV, eVal2_iV = extreme_value_decision(d2, value_iV);
        @test prob_iV == ones(Float64, nTypes)
        # If Gumbel shocks are subtracted (they have positive mean), demeaned expected values should be higher
        @test isapprox(eVal2_iV .- prefScale * EconLH.EulerConst, eVal_iV)
    end
end

function one_type_test()
    @testset "One type" begin
        prefScale = 0.8;
        nx = 5;
        valueV = (-0.2, 0.3, 0.9, -1.8, 2.0);
        @assert length(valueV) == nx;
        probV, eVal = extreme_value_decision_one(valueV, prefScale; demeaned = true);
        prob2V, eVal2 = extreme_value_decision_one([valueV...], prefScale;
            demeaned = true);
        @test isapprox(probV, prob2V);
        @test isapprox(eVal, eVal2);
    end
end

function many_types_test()
    @testset "Many types" begin
        prefScale = 0.8;
        nx = 3;
        nTypes = 5;
        value_ixV = [collect(LinRange(1.0, 2.0, nx) .+ 0.1 * j)  for j = 1 : nTypes];

        # Tuple input
        valueTuples = (value_ixV..., );
        prob1_ixV, eVal1_jV = extreme_value_decision(valueTuples, prefScale);

        # Matrix input. Columns are alternative
        value_ixM = zeros(nTypes, nx);
        for j = 1 : nTypes
            value_ixM[j,:] = value_ixV[j];
        end
        prob_ixM, eVal_jV = extreme_value_decision(value_ixM, prefScale);

        for j = 1 : nTypes
            @test isapprox(prob1_ixV[j], prob_ixM[j,:]);
        end
        @test isapprox(eVal_jV, eVal1_jV);
    end
end

function check_by_sim_test()
    @testset "Check by simulation" begin
        for demeaned in [true, false]
            n = 4;
            prefScale = 1.4;
            d = ExtremeValueDecision(prefScale, demeaned, true);

            # Row vector of values (as matrix)
            value_xV = collect(range(50.0, 55.0, length = n)');

            for nTypes in [1, 3]
                if nTypes == 1
                    value_ixM = value_xV;
                else
                    value_ixM = collect(range(1.0, 0.9, length = nTypes)) *
                        collect(range(50.0, 55.0, length = n)');
                end

                prob_ixM, eVal_iV = EconLH.extreme_value_decision(d, value_ixM);
                # This makes sure that there is something meaningful to simulate
                @test all(prob_ixM .< 0.95)
                @test check_by_sim(d, value_ixM, prob_ixM, eVal_iV);
            end
        end
    end
end

function draw_gumbel_test()
    @testset "Draw gumbel" begin
        rng = MersenneTwister(123);
        prefScale = 0.4;
        nSim = Int(1e3);

        shockV = draw_gumbel_shocks(rng, prefScale, nSim; demeaned = true);
        @test isapprox(mean(shockV), 0.0, atol = 0.05)
        @test size(shockV) == (nSim,)

        shockM = draw_gumbel_shocks(rng, prefScale, (nSim, 2); demeaned = false);
        @test isapprox(mean(shockM), prefScale * EconLH.EulerConst, atol = 0.05)
        @test size(shockM) == (nSim, 2)
    end
end

@testset "ExtremeValueDecision" begin
    one_option_test();
    one_type_test();
    many_types_test();
    check_by_sim_test();
    draw_gumbel_test();
end

# ---------------


text:
""""
Data Structure for the economical parameters of renewable generation technologies.
    The data structure can be called calling all the fields directly or using named fields.
    All the limits are defined by NamedTuples and some fields can take ```nothing```

    ## Examples

"""
struct EconRenewable <: TechnicalParams
    curtailpenalty::Float64 # [$/MWh]
    variablecost::Union{Float64,Nothing} # [$/MWh]
end

EconRenewable(; curtailcost = 0.0, variablecost = nothing) = EconRenewable(curtailcost, variablecost)

""""
Data Structure for the economical parameters of thermal generation technologies.
    The data structure can be called calling all the fields directly or using named fields.
    All the limits are defined by NamedTuples and some fields can take ```nothing```

    ## Examples

"""
struct EconThermal <: TechnicalParams
    capacity::Float64 # [MW]
    variablecost::Union{Function,Array{Tuple{Float64, Float64}}} # [$/MWh]
    fixedcost::Float64            # [$/h]
    startupcost::Float64          # [$]
    shutdncost::Float64           # [$]
    annualcapacityfactor::Union{Float64,Nothing}  # [0-1]
end

EconThermal(;   capacity = 0.0,
            variablecost = [(0.0,1.0)],
            fixedcost = 0.0,
            startupcost = 0.0,
            shutdncost = 0.0,
            annualcapacityfactor = nothing
        ) = EconThermal(capacity, variablecost, fixedcost, startupcost, shutdncost, annualcapacityfactor)


text:
#!/usr/bin/env julia

using Revise,SEM
using FastGaussQuadrature
using Plots,LinearAlgebra

n = 16
x,w = gausslobatto(n)
xo  = linspace(-1,1,20*n)

J = interpMat(xo,x)
D = derivMat(x)

# neumann operators
B = diagm(w)
A = D'*B *D
C(v) = B * diagm(v) * D

## CASE SETUP
function utrue(x,v,t)
    xx = @. x - v*t
    return @. sin(2pi*xx)
end
f = @. 0+0*x
vl= 1.0
v = @. vl+0*x
ν = 0e-0
u = utrue(x,vl,0.0)
id= Matrix(I,n,n)

#R = id[2:end-1,:] # double dirichlet
R = id[1:end-0,:] # periodic
Q = semq(1,n,true)

CFL = 0.1
dx = minimum(diff(x))
dt = dx * CFL / vl
println("CFL=$CFL, dt=$dt")

## global system (full rank)
xx = R*x
AA = R*A*R' .* ν
BB = R*B
CC = R*C(v)*R'
bb = BB*f
AC = AA + CC

## local system (rank deficient)
b  = R'*bb
Al = R'*AA*R
Bl = R'*BB
Cl = R'*CC*R
Hl = Al + Bl

#println("er: ",norm(u - ut,Inf))
p=plot(xo,J*u,width=3,ylims=(-1.5,1.5))
display(p)

uh = [zero(u) for i in 1:3] # histories
gh = [zero(u) for i in 1:3] # explicit term

T = 1.0
t = zeros(4,1)
nstep = Int(floor(T/dt))

for istep=1:nstep
    global u,rhs,Hl,bdfA,bdfB

    updateHist!(t)
    updateHist!(u,uh)

    t[1] += dt
    bdfA,bdfB = bdfExtK(t)

    gh[1] = -Cl*uh[1] # convection

    rhs = Bl * f
    for i=1:3
        gh[i] = -Cl*uh[i]                 # convection
        rhs  += -bdfB[1+i] .* (Bl *uh[i]) # histories
        rhs  +=  bdfA[i]   .* gh[i]       # explicit term
    end

    Hl = Al + bdfB[1]*Bl

    Hl  = Q*Q'*Hl
    rhs = Q*Q'*rhs
    u = pcg(rhs,Hl)

    ut = utrue(x,vl,t[1])
    er = norm(u .- ut,Inf)
    println("Pointwise error: $er")

    plt = plot(xo,J*u,width=3)
    plt = plot!(title="Step $istep, Time $(t[1])")
    plt = plot!(ylims=(-1.5,1.5))
    display(plt)
end
#----------------------------------------------------------------------#
nothing


text:
using WingBiology

using LinearAlgebra

using Test
using BenchmarkTools

include("Beam.jl")
include("Straight.jl")
include("WingBody.jl")
include("Interpolation.jl")
include("SimpleWing.jl")
include("Links.jl")

max_stars_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/WingBiology


text:
using ConservationLawsParticles
using SpecialFunctions, RecursiveArrayTools, DifferentialEquations, Plots
using Random
Random.seed!(42)

# model
V(t, x) = 0.0 # 2 + 0.5sin(2x) + 0.3sin(2sqrt(2)*x) + 0.2cos(2sqrt(7)*x)
W(t, x) = 0.0 # x^2/2 - abs(x)
W′(t, x) = 0.0 # x - sign(x)
mob(ρ) = 0.0 # max((ρ-1)^2 * (2ρ+1), 0)
smodel = DiffusiveSampledModel((V,), ((W′,),), (mob,), (Diffusion(1),))
imodel = DiffusiveIntegratedModel((V,), ((W,),), (mob,), (Diffusion(1),))

# initial condition
x0 = ArrayPartition(sort(randn(50)))

boxsol(t, x; a, b) = (erf((x-a)/2sqrt(t)) - erf((x-b)/2sqrt(t))) / (2*(b-a))
refsol(t, x) = sum(boxsol(t, x; a=x0[i-1], b=x0[i])
    for i in eachindex(x0)[2:end]) / (length(x0)-1)

# time span
tspan = (0., 0.1)

# ODE system for the particles
sprob = ODEProblem(velocities_diff!, x0, tspan, smodel)
# iprob = ODEProblem(velocities_diff!, x0, tspan, imodel)

abstol = reltol = 2e-8

@time ssol = solve(sprob, BS5(); abstol=abstol, reltol=reltol);
# @time isol = solve(iprob, BS5(); abstol=abstol, reltol=reltol);

# plot the particle trajectories
plot(legend=false, title="Trajectories", xlabel="time", ylabel="position")
plot!(ssol, vars=eachindex(x0)[1:5:end]; color=:blue)
savefig("diffusion5.png")

# plot an animation of the density
anim = @animate for t in range(sqrt.(tspan)...; length=101).^2
    plot(title="Density", xlabel="position", ylabel="density",
        xlims=(-5,3), ylims=(0,.8))
    plot!(x -> refsol(t, x); color=:blue, label="analytic")
    densityplot!(ssol(t); color=:red, label="particles")
end
gif(anim, "diffusion5.gif")

# plot an animation of the density
tspan = (0., .001)
anim = @animate for t in range(sqrt.(tspan)...; length=101).^2
    plot(title="Density", xlabel="position", ylabel="density",
        xlims=(-4,2), ylims=(0,1), legend=:topleft)
    plot!(x -> refsol(t, x); color=:blue, label="analytic")
    densityplot!(ssol(t); color=:red, label="particles")
end
gif(anim, "diffusion5.1.gif")


text:
ENV["EDITOR"] = "vim"
# ENV["PYTHON"] = "python"

function activate_default_packages()
    @eval using Revise
    @eval using OhMyREPL
    @eval colorscheme!("GruvboxDark")
    @eval using AbbreviatedStackTraces
    @eval ENV["JULIA_STACKTRACE_MINIMAL"] = true
    @eval using BenchmarkTools
    @eval using Logging: global_logger
    @eval using TerminalLoggers
    @eval global_logger(TerminalLogger())
    @eval using Eyeball
end

try
    activate_default_packages()
catch err
    @error "Could not load startup packages. To install, use `_setup()`" err
end

function _setup()
    @eval using Pkg
    Pkg.add("Revise")
    Pkg.add("OhMyREPL")
    Pkg.add("BenchmarkTools")
    Pkg.add("TerminalLoggers")
    Pkg.add(url="https://github.com/BioTurboNick/AbbreviatedStackTraces.jl")
    Pkg.add("Eyeball")
    activate_default_packages()
end

max_stars_repo_name: nathbo/dotfiles

max_issues_repo_name: nathbo/dotfiles

max_forks_repo_name: nathbo/dotfiles


text:
# Generate documentation with this command:
# (cd docs && julia --color=yes make.jl)

push!(LOAD_PATH, "..")

using Documenter
using SSHT

makedocs(; sitename="SSHT", format=Documenter.HTML(), modules=[SSHT])

deploydocs(; repo="github.com/eschnett/SSHT.jl.git", devbranch="main", push_preview=true)


text:
struct RandomStrategy{G <: Game} <: Strategy end

function BoardGames.getmove(board, s::RandomStrategy)
    rand(getmoves(board))
end

function BoardGames.getvarsnames(s::RandomStrategy)
    return String[]
end

function BoardGames.getvalues(s::RandomStrategy)
    return ()
end

BoardGames.name(::RandomStrategy) = "Random Selection"

function Base.copy(s::RandomStrategy)
    return s
end