text:
using DrWatson
@quickactivate "EPExperiments"

using EPExperiments
using Logging
using LoggingExtras
using JLD2
using FileIO
using Random
using AdvancedHMC
using CSV
using DataFrames
using StatsPlots
using ArgParse

import StatsFuns: logsumexp

const LOG_FILENAME = "out.log"

@with_kw struct RadonEstimateConfig{T<:AnISConfig}
    experiment_name::String = "test"
    num_samples::Int = 10
    num_betas::Int = 101
    anis_params::T = AnISHMCConfig{Turing.ReverseDiffAD{false}}(
        proposal=AdvancedHMC.StaticTrajectory(AdvancedHMC.Leapfrog(0.05), 10),
        num_steps=10,
        rejection_type=SimpleRejection()
    )
    seed::Int = 1234
    data_file::String = "radon.csv"
    radon_model::Expectation = radon_hierarchical
    cost_fn_name::Symbol = :radon_target_f
end

# Some old experimental results didn't include the cost_fn_name
# as a parameter in the config.
function old2new_config(old_config, cost_fn_name)
    return RadonEstimateConfig(
        experiment_name=old_config.experiment_name,
        num_samples=old_config.num_samples,
        anis_params=old_config.anis_params,
        num_betas=old_config.num_betas,
        data_file=old_config.data_file,
        radon_model=old_config.radon_model,
        seed=old_config.seed,
        cost_fn_name=cost_fn_name
    )
end

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--seed", "-s"
            help = "Random seed for the program."
            arg_type = Int
            required = true
    end

    return parse_args(s)
end

function load_radon_results(
    experiment_dirname
)
    result_folder = datadir("radon_estimate", experiment_dirname)
    results = load(joinpath(result_folder, "results.jld2"))

    cost_estimate = results["cost_estimate"]
    radon_tabi = results["radon_tabi"]
    experiment_config = results["experiment_config"]
    logger = TeeLogger(
        ConsoleLogger(),
        FileLogger(joinpath(result_folder, LOG_FILENAME))
    )

    return cost_estimate, radon_tabi, experiment_config, result_folder, logger
end

function reprocess_radon_results(
    experiment_dirname;
    old_experiment_config=false,
    cost_fn_name=nothing
)
    cost_estimate, radon_tabi, experiment_config, result_folder, logger = load_radon_results(
        experiment_dirname
    )

    if old_experiment_config
        @assert !isnothing(cost_fn_name)
        experiment_config = old2new_config(experiment_config, cost_fn_name)
    end

    process_radon_results(
        cost_estimate,
        radon_tabi,
        experiment_config,
        result_folder,
        logger
    )
end

function process_radon_results(
    cost_estimate,
    radon_tabi,
    experiment_config,
    result_folder,
    logger
)
    @unpack num_betas = experiment_config
    betas = (geomspace(1, 1001, num_betas) .- 1) ./ 1000

    @unpack cost_fn_name = experiment_config
    cost_fn = RADON_COST_FNS[cost_fn_name]

    @show size(radon_tabi[:Z2_info].info[:intermediate_log_weights])
    @show size(betas)
    savefig(plot_intermediate_weights(
        radon_tabi[:Z2_info].info[:intermediate_log_weights], betas),
        joinpath(result_folder, "Z2_intermediate_weights.png")
    )
    savefig(plot_intermediate_weights(
        radon_tabi[:Z1_positive_info].info[:intermediate_log_weights], betas),
        joinpath(result_folder, "Z1_intermediate_weights.png")
    )
    savefig(plot_intermediate_ess(
        radon_tabi[:Z2_info].info[:intermediate_log_weights], betas),
        joinpath(result_folder, "Z2_intermediate_ess.png")
    )
    savefig(plot_intermediate_ess(
        radon_tabi[:Z1_positive_info].info[:intermediate_log_weights], betas),
        joinpath(result_folder, "Z1_intermediate_ess.png")
    )
    params_to_plot = [:mu_alpha, :mu_beta]
    savefig(
        plot(radon_tabi[:Z2_info][:,params_to_plot,:]),
        joinpath(result_folder, "Z2_samples.png")
    )
    savefig(
        plot(radon_tabi[:Z1_positive_info][:,params_to_plot,:]),
        joinpath(result_folder, "Z1_samples.png")
    )

    alpha_names = namesingroup(radon_tabi[:Z1_positive_info], :alpha)
    savefig(
        plot(radon_tabi[:Z1_positive_info][:,alpha_names,:], colordim=:parameter),
        joinpath(result_folder, "Z1_counties.png")
    )
    savefig(
        plot(radon_tabi[:Z1_positive_info][:,[alpha_names[10]],:]),
        joinpath(result_folder, "Z1_county10.png")
    )
    savefig(
        plot(radon_tabi[:Z2_info][:,alpha_names,:], colordim=:parameter),
        joinpath(result_folder, "Z2_counties.png")
    )
    savefig(
        plot(radon_tabi[:Z2_info][:,[alpha_names[10]],:]),
        joinpath(result_folder, "Z2_county10.png")
    )

    Z2_chn = radon_tabi[:Z2_info]
    alpha = Array(Z2_chn[namesingroup(Z2_chn, :alpha)])
    beta = Array(Z2_chn[namesingroup(Z2_chn, :beta)])

    eps = 1
    radon_losses = map(1:size(alpha, 1)) do (i)
        cost_fn(
            alpha[i,:],
            beta[i,:],
            POST_PRED_COUNTY_IDX,
            POST_PRED_FLOOR
        )
    end

    Z2_lws = Array(Z2_chn[:log_weight])
    normalisation = exp(logsumexp(Z2_lws))
    anis_estimate = sum(exp.(Z2_lws) .* radon_losses) / normalisation

    Z1_lws = Array(radon_tabi[:Z1_positive_info][:log_weight])
    Z1_lws_retargeted = Z1_lws .- log.(radon_losses)
    Z2_lws_retargeted = Z2_lws .+ log.(radon_losses)

    ESS_retargeted = (
        Z1_positive_info = compute_ess(Z1_lws_retargeted),
        Z2_info = compute_ess(Z2_lws_retargeted)
    )

    # NOTE: Usually we have to divide each estimate by the number of samples used
    # but because those "normalisations" cancel in our downstream calculations we
    # avoid them.
    Z1_ws = exp.(Z1_lws)
    Z2_ws = exp.(Z2_lws)
    Z1s_unnormalised = cumsum(Z1_ws, dims=1)
    Z2s_unnormalised = cumsum(Z2_ws, dims=1)
    tabi_ests = Z1s_unnormalised ./ Z2s_unnormalised
    anis_ests = cumsum(vec(Z2_ws) .* radon_losses) ./ Z2s_unnormalised

    num_samples = size(tabi_ests, 1)
    ixs = 1:num_samples
    tabi_ixs = 2 * ixs
    convergence_plot = plot(tabi_ixs, tabi_ests[ixs], label="TAAnIS")
    plot!(
        convergence_plot,
        ixs,
        anis_ests[ixs],
        label="AnIS"
    )
    savefig(convergence_plot, joinpath(result_folder, "convergence_plot.png"))

    ess_plot = make_ess_plot(Z1_lws, Z2_lws, Z2_lws_retargeted)
    savefig(ess_plot, joinpath(result_folder, "ess_plot.png"))

    with_logger(logger) do
        @info "Expectation estimate TABI: $(cost_estimate)"
        @info "Expectation estimate AnIS: $(anis_estimate)"

        for (k, est_results) in pairs(radon_tabi)
            if isnothing(est_results)
                continue
            end
            msg = string([
                "$(string(k)):\n",
                "ESS: $(est_results.info[:ess])\n",
                "ESS retargeted: $(ESS_retargeted[k])\n",
                "Log evidence: $(est_results.logevidence)\n",
                "Log weights: $(get(est_results, :log_weight)[:log_weight])\n"
            ]...)

            @info "$(msg)"
        end
    end
end

function radon_estimate(experiment_config)
    @unpack seed, anis_params, num_betas, num_samples = experiment_config
    @unpack radon_model, cost_fn_name, data_file = experiment_config
    Random.seed!(seed)

    result_folder = mkpath(joinpath(
        datadir("radon_estimate"),
        savename(experiment_config)
    ))

    logger = TeeLogger(
        ConsoleLogger(),
        FileLogger(joinpath(result_folder, LOG_FILENAME))
    )

    radon_df = DataFrame(CSV.File(datadir("exp_raw", data_file)))
    log_radon, county_idx, floor, num_counties = preprocess_radon_df(
        radon_df, data_file
    )

    radon_exp = radon_model(
        log_radon,
        county_idx,
        floor,
        num_counties;
        cost_fn=RADON_COST_FNS[cost_fn_name]
    )

    # Default:
    betas = (geomspace(1, 1001, num_betas) .- 1) ./ 1000
    # 200_400 ratio:
    # low_beta = 3 * 1e-3
    # low_betas = ((geomspace(1, 1001, 200) .- 1) ./ 1000) .* low_beta
    # high_betas = ((geomspace(1, 1001, 202) .- 1) ./ 1000) .* (1 - low_beta) .+ low_beta
    # betas = vcat(low_betas, high_betas[2:end])
    anis_alg = get_anis_alg(anis_params, betas)
    tabi = TABI(
        TuringAlgorithm(anis_alg, num_samples),
        TuringAlgorithm(anis_alg, 0),
        TuringAlgorithm(anis_alg, num_samples)
    )

    @time begin
    cost_estimate, radon_tabi = estimate_expectation(
        radon_exp,
        tabi;
        store_intermediate_samples=true,
        progress=true
    )
    end

    @tagsave(joinpath(result_folder, "results.jld2"), Dict(
        "cost_estimate" => cost_estimate,
        "radon_tabi" => radon_tabi,
        "experiment_config" => experiment_config
    ))

    process_radon_results(
        cost_estimate,
        radon_tabi,
        experiment_config,
        result_folder,
        logger
    )
end

radon_estimate(RadonEstimateConfig(
    experiment_name="eps_super_simple_model",
    num_samples=1000,
    anis_params=AnISHMCConfig{Turing.ReverseDiffAD{false}}(
        proposal=AdvancedHMC.NUTS{MultinomialTS,GeneralisedNoUTurn}(
            AdvancedHMC.Leapfrog(0.044)
        ),
        num_steps=1,
        rejection_type=SimpleRejection()
    ),
    num_betas=201,
    data_file="small_radon_num_counties=20.csv",
    radon_model=radon_hierarchical,
    seed=600,
    cost_fn_name=:radon_target_f
))


text:
module LogLikeNumbers
using TropicalNumbers
using LogarithmicNumbers
using NiLang, NiLang.AD

import NiLang: i_sum
using ..Reversible: TropicalTypes

export Tropical

include("base.jl")
include("instructs.jl")
include("blas.jl")
include("einsum.jl")

end # module

ext: jl


text:
if basename(pwd()) == "aoc"
    cd("2017/11")
end

function hexdistance(path::AbstractString)
    position = [0, 0]
    for step in split(path, ',')
        position .+= if step == "n"
            [0, 2]
        elseif step == "ne"
            [1, 1]
        elseif step == "se"
            [1, -1]
        elseif step == "nw"
            [-1, 1]
        elseif step == "sw"
            [-1, -1]
        else
            [0, -2]
        end
    end
    sum(abs.(position)) ÷ 2
end

@assert hexdistance("ne,ne,ne") == 3
@assert hexdistance("ne,ne,sw,sw") == 0
@assert hexdistance("ne,ne,s,s") == 2
@assert hexdistance("se,sw,se,sw,sw") == 3

hexdistance(readline("input.txt"))

function maxhexdistance(path::AbstractString)
    position = [0, 0]
    maxdist = 0
    for step in split(path, ',')
        position .+= if step == "n"
            [0, 2]
        elseif step == "ne"
            [1, 1]
        elseif step == "se"
            [1, -1]
        elseif step == "nw"
            [-1, 1]
        elseif step == "sw"
            [-1, -1]
        else
            [0, -2]
        end
        maxdist = max(maxdist, sum(abs.(position)) ÷ 2)
    end
    maxdist
end

maxhexdistance(readline("input.txt"))

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
"""
    KL(p::Beta, q::Beta)

See [KL Beta](https://en.wikipedia.org/wiki/Beta_distribution#Quantities_of_information_(entropy))
"""
function KL(p::Beta, q::Beta)
    αp, βp = params(p)
    αq, βq = params(q)
    return logbeta(αq, βq) - logbeta(αp, βp) + (αp - αq) * digamma(αp) +
        (βp - βq) * digamma(βp) + (αq - αp + βq - βp) * digamma(αp + βp)
end

"""
    KL(p::Exponential, q::Exponential)

See [KL Exponential](https://en.wikipedia.org/wiki/Exponential_distribution#Kullback%E2%80%93Leibler_divergence)
"""
function KL(p::Exponential, q::Exponential)
    λp = scale(p)
    λq = scale(q)
    return log(λp) - log(λq) + λq / λp - 1
end

"""
    KL(p::Gamma, q::Gamma)

See [KL Gamma](https://en.wikipedia.org/wiki/Gamma_distribution#Kullback%E2%80%93Leibler_divergence)
"""
function KL(p::Gamma, q::Gamma)
    # We use the parametrization with the rate β
    αp, αq = shape.((p, q))
    βp, βq = rate.((p, q))
    return (αp - αq) * digamma(αp) - loggamma(αp) + loggamma(αq) +
        αq * (log(βp) - log(βq)) + αp * (βq - βp) / βp
end

"""
    KL(p::InverseGamma, q::InverseGamma)

See [KL Inverse-Gamma](https://en.wikipedia.org/wiki/Inverse-gamma_distribution#Properties)
"""
function KL(p::InverseGamma, q::InverseGamma)
    # We can reuse the implementation of Gamma
    return KL(Gamma(shape(p), rate(p)), Gamma(shape(q), rate(q)))
end

"""
    KL(p::Normal, q::Normal)

See [KL Gaussian](https://en.wikipedia.org/wiki/Normal_distribution#Other_properties)
"""
function KL(p::Normal, q::Normal)
    μp, σp = params(p)
    μq, σq = params(q)
    return 0.5 * (abs2(σp / σq) + abs2((μp - μq) / σq) - 1 + 2 * (log(σq) - log(σp)))
end

"""
    KL(p::Poisson, q::Poisson)

See [KL Poisson](https://en.wikipedia.org/wiki/Poisson_distribution#Other_properties)
"""
function KL(p::Poisson, q::Poisson)
    λp, λq = rate.((p, q))
    return λq - λp + λp * (log(λp) - log(λq))
end


text:
baremodule _Compiler

import Base
using Core.Compiler: SimpleVector
for n in Base.names(Core.Compiler; all = true, imported = true)
    n === :typeinf_code && continue
    Base.startswith(Base.string(n), "#") && continue
    if Base.isdefined(Core.Compiler, n)
        Base.@eval _Compiler using Core.Compiler: $n
    end
end

#! format: off

function typeinf_code(interp::AbstractInterpreter, method::Method, @nospecialize(atypes), sparams::SimpleVector, run_optimizer::Bool)
    mi = specialize_method(method, atypes, sparams)::MethodInstance
    # ccall(:jl_typeinf_begin, Cvoid, ())
    result = InferenceResult(mi)
    frame = InferenceState(result, false, interp)
    frame === nothing && return (nothing, Any)
    if typeinf(interp, frame) && run_optimizer
        opt_params = OptimizationParams(interp)
        opt = OptimizationState(frame, opt_params, interp)
        optimize(opt, opt_params, result.result)
        opt.src.inferred = true
    end
    # ccall(:jl_typeinf_end, Cvoid, ())
    frame.inferred || return (nothing, Any)
    return (frame.src, widenconst(result.result))
end

end


text:
function hfun_bar(vname)
  val = Meta.parse(vname[1])
  return round(sqrt(val), digits=2)
end

function hfun_m1fill(vname)
  var = vname[1]
  return pagevar("index", var)
end

function lx_baz(com, _)
  # keep this first line
  brace_content = Franklin.content(com.braces[1]) # input string
  # do whatever you want here
  return uppercase(brace_content)
end

####
function lx_it(com, _)
  # keep this first line
  brace_content = Franklin.content(com.braces[1]) # input string
  # do whatever you want here
  return "*$(brace_content)*"
end

function lx_bf(com, _)
  # keep this first line
  brace_content = Franklin.content(com.braces[1]) # input string
  # do whatever you want here
  return "**$(brace_content)**"
end

function lx_pdf(com, _)
  # keep this first line
  brace_content = Franklin.content(com.braces[1]) # input string
  # do whatever you want here
  return "[[pdf]]($(brace_content))"
end

function lx_doi(com, _)
  # keep this first line
  brace_content = Franklin.content(com.braces[1]) # input string
  # do whatever you want here
  return "[[doi]](https://doi.org/$(brace_content))"
end

function lx_arXiv(com, _)
  # keep this first line
  brace_content = Franklin.content(com.braces[1]) # input string
  # do whatever you want here
  return "[[arXiv]](https://arXiv.org/abs/$(brace_content))"
end


text:
using StructTypes
using FilePaths: AbstractPath

const MIME_TYPES = Dict(
    "acc" => "audio/acc",
    "abw" => "application/x-abiword",
    "arc" => "application/x-freearc",
    "avi" => "video/x-msvideo",
    "azw" => "application/vnd.amazon.ebook",
    "bin" => "application/octet-stream",
    "bmp" => "image/bmp",
    "bz"  => "application/x-bzip",
    "bz2" => "application/x-bzip2",
    "csh" => "application/x-csh",
    "css" => "text/css",
    "csv" => "text/csv",
    "doc" => "application/msword",
    "docx" => "application/vnd.openxmlformats-officedocument.wordprocessingml.document",
    "eot" => "application/vnd.ms-fontobject",
    "epub" => "application/epub+zip",
    "gz" => "application/gzip",
    "gif" => "image/gif",
    "htm" => "text/html",
    "html" => "text/html",
    "ico" =>  "image/vnd.microsoft.icon",
    "ics" =>  "text/calendar",
    "jar" => "application/java-archive",
    "jpeg" => "image/jpeg",
    "jpg" => "image/jpeg",
    "js" => "text/javascript",
    "json" => "application/json",
    "jsonld" => "application/ld+json",
    "mid" => "audio/midi",
    "mjs" => "text/javascript",
    "mp3" => "audio/mpeg",
    "mpeg" => "video/mpeg",
    "mpkg" => "application/vnd.apple.installer+xml",
    "odp" => "application/vnd.oasis.opendocument.presentation",
    "ods" => "application/vnd.oasis.opendocument.spreadsheet",
    "odt" => "application/vnd.oasis.opendocument.text",
    "oga" => "audio/ogg",
    "ogv" => "video/ogg",
    "ogx" => "application/ogg",
    "opus" => "audio/opus",
    "otf" => "font/otf",
    "png" => "image/png",
    "pdf" => "application/pdf",
    "php" => "application/x-httpd-php",
    "ppt" => "application/vnd.ms-powerpoint",
    "pptx" => "application/vnd.openxmlformats-officedocument.presentationml.presentation",
    "rar" => "application/vnd.rar",
    "rtf" => "application/rtf",
    "sh" => "application/x-sh",
    "svg" => "	image/svg+xml",
    "swf" => "application/x-shockwave-flash",
    "tar" => "application/x-tar",
    "tif" => "image/tiff",
    "tiff" => "image/tiff",
    "ts" => "video/mp2t",
    "ttf" => "font/ttf",
    "txt" => "text/plain",
    "vsd"=> "application/vnd.visio",
    "wav" => "audio/wav",
    "weba" => "audio/webm",
    "webm" => "video/webm",
    "webp" => "image/webp",
    "woff" => "font/woff",
    "woff2" => "font/woff2",
    "xhtml" => "application/xhtml+xml",
    "xls" => "application/vnd.ms-excel",
    "xlsx" => "application/vnd.openxmlformats-officedocument.spreadsheetml.sheet",
    "xml" => "application/xml",
    "xul" => "application/vnd.mozilla.xul+xml",
    "zip" => "application/zip",
    "7z" => "application/x-7z-compressed"
)

function mime_type(p::AbstractPath)
    ext = extension(p)
    get(MIME_TYPES, ext, "text/plain")
end

function mime_type(t::DataType)
   if StructTypes.StructType(t) != StructTypes.NoStructType()
        return "application/json"
   elseif t <: NamedTuple
        return "application/json"
   elseif t <: AbstractString
        return "text/plain"
   else
        return nothing
   end
end

# This is needed to allow type inference in io.jl
mime_type(::Any) = nothing

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
# In Julia, do: include("onTheFlyGH.jl"); compare();
using FastGaussQuadrature
using Printf

function compare()

	m = 6
	exa = sqrt(pi)/1^((m+1)/2)*prod((2:2:m).-1)/2^(m/2) # equals \int_{-\infty}^\infty p^m exp(-[k=1]*p^2) dp for even m
	ns = 2 .^(8:2:18)
	print(exa, ", ns=", ns, "\n")
	rests = NaN*ones(4,length(ns))
	s1 = "Time (s) ASY  "
	s2 = "Time (s) cref{Aherm}??  " #Add and check \cref in .tex
	s3 = "Error ASY  "
	s4 = "Error cref{Aherm}??  "
	for ni = 1:length(ns)
		n = ns[ni]
		rests[1,ni] = @elapsed igh = gh(2*n)
		s1 = string(s1, " & ", @sprintf("%7.1e", rests[1,ni] ))
		rests[2,ni] = @elapsed ig = onTheFly(n, ceil(Int64, 17*sqrt(n*1.0)) )
		s2 = string(s2, " & ", @sprintf("%7.1e", rests[2,ni] ))
		@time ig = onTheFly(n, ceil(Int64,17*sqrt(n*1.0)) )
		@time igh = gh(2*n)
		@test (ig-exa)/exa < 1e-10
		@test (igh-exa)/exa < 1e-10
		# print(n, ", ", (ig-exa)/exa, "=ig, igh=", (igh-exa)/exa, "\n")
		rests[3,ni] = (igh-exa)/exa
		s3 = string(s3, " & ", @sprintf("%7.1e", rests[3,ni] ))
		rests[4,ni] = (ig-exa)/exa
		s4 = string(s4, " & ", @sprintf("%7.1e", rests[4,ni] ))
	end
	print(rests, "\n\n asdf \n", s1, "\n", s2, "\n", s3, "\n", s4, "\n") # Prints the table
end

function onTheFly(n::Int64, mi::Int64)

	alpha = -0.5 # Expect Julia to inline this into the following formulae
	d = 1/(4*n+2*alpha+2) # Could specify all variables to be d::Float64
	s = 0.0
	csn = ceil(Int64,sqrt(n*1.0))
	for k=min(n,mi):-1:(csn+1)
		p = (4*n-4*k+3.0)/(4*n+1)
		t = pi^2*(p-1)^2/16
		for i = 1:6 # Newton iterations for t
			t = t-(p*pi+2*sqrt(t-t^2) -acos(2*t-1))/2*sqrt(t/(1-t))
		end
		x = 0.0
		w = 0.0
		# We aim for extremely high n so only four terms suffice
		if n < 1265
			x = x -d^5/181440 *(1-t)^2/t^2 *(10797500*(1-t)^(-10) - 43122800*(1-t)^(-9) + 66424575*(1-t)^(-8) -48469876*(1-t)^(-7) + 193536*alpha^6 + 16131880*(1-t)^(-6) + 80*(315*alpha^4 - 630*alpha^2 -221)*(1-t)^(-4) - 1727136*(1-t)^(-5) - 967680*alpha^4 - 320*(63*alpha^4 - 126*alpha^2 +43)*(1-t)^(-3)  + 384*(945*alpha^6 - 4620*alpha^4 + 6405*alpha^2 - 1346)*(1-t)^(-2) +1354752*alpha^2   - 23040*(21*alpha^6 - 105*alpha^4 + 147*alpha^2 - 31)*(1-t)^(-1) -285696)
			w = w - (1-t)^3/90720/t^3 *d^6 *(43190000*(1-t)^(-12) -204917300*(1-t)^(-11) + 393326325*(1-t)^(-10)  - 386872990*(1-t)^(-9) + 201908326*(1-t)^(-8) +80*(315*alpha^4 - 630*alpha^2 + 53752)*(1-t)^(-6)  - 50986344*(1-t)^(-7) - 320*(189*alpha^4 -378*alpha^2 - 89)*(1-t)^(-5) + 480*(63*alpha^4 - 126*alpha^2 + 43)*(1-t)^(-4)  -384*(315*alpha^6 - 1470*alpha^4 + 1995*alpha^2 - 416)*(1-t)^(-3) + 2304*(21*alpha^6 -105*alpha^4 + 147*alpha^2 - 31)*(1-t)^(-2) )
		end
		if n < 22026
			x = x +d^3*(1-t)/t/720*(1600*(1-t)^(-6) - 3815*(1-t)^(-5) + 480*alpha^4 +2814*(1-t)^(-4) - 576*(1-t)^(-3)   - 960*alpha^2 - 48*(15*alpha^4 - 30*alpha^2 + 7)*(1-t)^(-1) -16*(1-t)^(-2) + 224)
			w = w +(1-t)^2/720/t^2*d^4*(8000*(1-t)^(-8) - 24860*(1-t)^(-7) + 27517*(1-t)^(-6) - 12408*(1-t)^(-5) + 1712*(1-t)^(-4) +16*(15*alpha^4 - 30*alpha^2 + 7)*(1-t)^(-2) + 32*(1-t)^(-3))
		end
		x = sqrt(t*(4*n+1) -1/(3*(4*n+1))*(5/4*(1-t)^(-2) - (1-t)^(-1) - 1/4) + x)
		s = s+((-x)^6 + (x)^6)/x*exp(-x^2)*pi*sqrt(t/(1-t))*(1 + (2*t + 3)/6/(4*n+1)^2/(t-1)^3 +w)
	end
	for k=csn:-1:1
		jak = pi*(k-1/2)
		x = 0.0
		w = 0.0
		if n < 259
			x= x + (10644*jak^8 + 60*(887*alpha^2 - 2879)*jak^6 + (125671*alpha^4 -729422*alpha^2 + 1456807)*jak^4 + 3*(63299*alpha^6 - 507801*alpha^4 + 1678761*alpha^2 - 2201939)*jak^2 + 2*(107959*alpha^8 - 1146220*alpha^6 + 5095482*alpha^4 -10087180*alpha^2 + 6029959) )*d^8/42525
	   		w = w + (215918*alpha^8 + 53220*jak^8 + 240*(887*alpha^2 - 2879)*jak^6 -2292440*alpha^6 + 3*(125671*alpha^4 - 729422*alpha^2 + 1456807)*jak^4 + 10190964*alpha^4  + 6*(63299*alpha^6 - 507801*alpha^4 + 1678761*alpha^2 -2201939)*jak^2 - 20174360*alpha^2 + 12059918)/42525*d^8
		end
		if n < 1265
			x = x + (657*jak^6 +36*jak^4*(73*alpha^2-181) +2*jak^2*(2459*alpha^4 -10750*alpha^2 +14051) + 4*(1493*alpha^6 -9303*alpha^4 +19887*alpha^2 - 12077) )*d^6/2835
			w = w + (1493*alpha^6 + 657*jak^6 + 27*(73*alpha^2 - 181)*jak^4 - 9303*alpha^4  + (2459*alpha^4 -10750*alpha^2 + 14051)*jak^2 + 19887*alpha^2 - 12077)*4/2835*d^6
		end
		if n < 22026
			x = x + (11*jak^4 +3*jak^2*(11*alpha^2-19) +46*alpha^4 -140*alpha^2 +94)*d^4/45
			w = w + (46*alpha^4 + 33*jak^4 +6*jak^2*(11*alpha^2 -19) -140*alpha^2 +94)/45*d^4
		end
		x = pi*(k-1/2)*sqrt(1/(4*n+1) +(pi^2*(k-1/2)^2 -3/2)/(3*(4*n+1)^3) + x*d )
		s = s+((-x)^6 + (x)^6)*exp(-x^2)/(x*(n+1/4)*(besselj(-3/2, pi*(k-1/2)) )^2)*(1 + (2*pi^2*(k-1/2)^2 -3/2)/(3*(4*n+1)^2) +w)/2
	end
	return s
end

function gh(n)
	(x,w) = gausshermite(2n)
	sum(w .* x.^6)
end


text:
function filter_records(records, gCols, eCols)
	#Helper function to transform and integrate the measurement records

	#Extract the mean ground and excited traces
	meanGround = squeeze(mean(records[:,gCols,:], [2,3]), [2,3]);
	meanExcited = squeeze(mean(records[:,eCols,:], [2,3]), [2,3]);

	#Take the total mean as the shift operator
	centres = 0.5*(meanGround+meanExcited);

	#Take the difference to get the weighting and angle for the optimum quadrature
	D = meanExcited-meanGround;
	rotAngles = angle(D);

	#We square the weights to normalize to +/- 0.5
	weights = abs(D)/(sum(abs(D).^2));

	#Now we can unwind
	unwoundRecords = (records .- centres) .* exp(-1im*rotAngles);

	#Take the weighted sum and look at the good quadrature
	intRecords = squeeze(sum(weights.*real(unwoundRecords),1),1);

	return intRecords

end

function comp_to_paulis(records)
	#Helper function to convert from computational basis calibrations to pauli basis II, IZ, ZI, ZZ etc.

	#Matrix to convert from computation basis to Pauli basis diagonals
	A = 0.5*[1 1 1 1; 1 -1 1 -1; 1 1 -1 -1; 1 -1 -1 1]

	#Look at mean and variance transformed
	meanSig = mean(records, 2)
	varSig = cov(transpose(records))

	betas = A*meanSig;
	varBetas = A*varSig*A';

	#And the standard error of the mean
	betasErrors = sqrt(diag(varBetas)/size(records,2));

	return betas, betasErrors

end


text:
using ExprParsers

EP.@exprparser struct MySymbol
  symbol = EP.anything
end
@test MySymbol <: EP.ExprParser
@test MySymbol <: EP.ExprParserWithParsed
@test MySymbol_Parsed <: EP.ExprParsed

@test EP.ExprParsed(MySymbol) === MySymbol_Parsed

@test parse_expr(MySymbol(), symbol = 4) == MySymbol_Parsed(symbol = 4)


text:
struct Filters
    options::Vector{Option}
    Filters(options::AbstractVector{<:Option}) = new(convert(Vector{Option}, options))
end

function Filters(df)
    cols = Tables.columns(df)
    selectors = Option[]
    for name in Tables.columnnames(cols)
        v = Tables.getcolumn(cols, name)
        selector = if iscontinuous(v)
            range = histrange(v, 50)
            Option(string(name), RangeSelector(range), Observable(false))
        else
            unique = uniquesorted(v) # TODO: optimize for strings
            options = [Option(string(val), val, Observable(true)) for val in unique]
            Option(string(name), Checkboxes(options), Observable(false))
        end
        push!(selectors, selector)
    end
    return Filters(selectors)
end

function filter_column!(previous::AbstractVector{Bool}, v::AbstractVector, cb::Checkboxes)
    accepted_options = Set(option.value for option in cb.options if option.selected[])
    previous .&= in(accepted_options).(v)
    return previous
end

function filter_column!(previous::AbstractVector{Bool}, v::AbstractVector, rg::RangeSelector)
    min, max = map(getindex, rg.selected)
    @. previous &= min ≤ v ≤ max
    return previous
end

filter_column!(::Nothing, v::AbstractVector, rg) = filter_column!(fill(true, length(v)), v,  rg)

function selected_data(f::Filters, table)
    data, options = Tables.columns(table), f.options

    selected = foldl(options, init=nothing) do acc, option
        colname = Symbol(option.key)
        wdg = option.value
        column = Tables.getcolumn(data, colname)
        return filter_column!(acc, column, wdg)
    end
    return mapcols!(t -> t[selected], SimpleTable(data))
end

jsrender(session::Session, filters::Filters) = jsrender(session, Togglers(filters.options))


text:
using Tar

PNG_SUITE_DIR = "PngSuite"
PNG_SUITE_PATH = joinpath(PNG_TEST_PATH, PNG_SUITE_DIR)
PNG_SUITE_FILE = joinpath(PNG_TEST_PATH, "PngSuite.tgz")

# See https://github.com/JuliaIO/PNGFiles.jl/pull/5 for discussion of these errors
const imdiff_tolerance = Dict(
    "basi2c16.png" => 0.010046128738460332,
    "basi6a08.png" => 0.021200980392156808,
    "basn2c16.png" => 0.010046128738460332,
    "basn6a08.png" => 0.021200980392156808,
    "bgan6a08.png" => 0.021200980392156808,
    "bgwn6a08.png" => 0.021200980392156808,
    "g03n2c08.png" => 0.05052849264705868,
    "g05n2c08.png" => 0.041383272058823424,
    "g07n2c08.png" => 0.09764476102941164,
    "g25n2c08.png" => 0.203565410539216,
    "oi1n2c16.png" => 0.010046128738460332,
    "oi2n2c16.png" => 0.010046128738460332,
    "oi4n2c16.png" => 0.010046128738460332,
    "oi9n2c16.png" => 0.010046128738460332,
    "pp0n2c16.png" => 0.010046128738460332,
    "ps1n2c16.png" => 0.010046128738460332,
    "ps2n2c16.png" => 0.010046128738460332,
    "tp0n2c08.png" => 0.012859987745098002,
    "tp0n3p08.png" => 0.012331495098039179,
)

if !isdir(PNG_SUITE_PATH)
    mkdir(PNG_SUITE_PATH)
    download("https://github.com/JuliaIO/PNG.jl/releases/download/PngSuite-2017jul19/PngSuite-2017jul19.tar.gz", PNG_SUITE_FILE)
    Tar.extract(PNG_SUITE_FILE, PNG_SUITE_PATH)
    rm(PNG_SUITE_FILE)
end

struct _Palleted; end
const pngsuite_colormap = Dict("0g" => Gray, "2c" => RGB, "3p" => _Palleted, "4a" => GrayA, "6a" => RGBA)

function parse_pngsuite(x::AbstractString)
    code = splitext(x)[1]
    return (
        case=code[1:end-5],
        is_interlaced=code[end-4]=='i',
        color_type=pngsuite_colormap[code[end-3:end-2]],
        bit_depth=parse(Int, code[end-1:end])
    )
end
parse_pngsuite(x::Symbol) = parse_pngsuite(String(x))

@testset "PngSuite" begin
    for fpath in glob(joinpath("./**/$(PNG_SUITE_DIR)", "[!x]*[!_new][!_new_im].png"))
        case = last(splitpath(fpath))
        case_info = parse_pngsuite(case)
        C = case_info.color_type
        b = case_info.bit_depth
        is_gray = C in (Gray, GrayA)

        @testset "$(case)" begin
            global read_in_pngf = PNGFiles.load(fpath, gamma = is_gray ? 1.0 : nothing)
            @test typeof(read_in_pngf) <: AbstractMatrix

            path, ext = splitext(fpath)
            newpath = path * "_new" * ext

            open(io->PNGFiles.save(io, read_in_pngf), newpath, "w") #test IO method
            @test PNGFiles.save(newpath, read_in_pngf) == 0
            global read_in_immag = _standardize_grayness(ImageMagick.load(fpath))

            @testset "$(case): PngSuite/ImageMagick read type equality" begin
                if case_info.case in ("tbb", "tbg", "tbr", "tbw")  # transaprency adds an alpha layer
                    if C === RGB
                        C = RGBA
                    elseif C === Gray
                        C = GrayA
                    end
                end
                if C === _Palleted
                    # _Palleted images could have an alpha channel, but its not evident from "case"
                    @test eltype(read_in_pngf) == eltype(read_in_immag)
                else
                    basictype = b > 8 ? Normed{UInt16,16} : Normed{UInt8,8}
                    @test eltype(read_in_pngf) == C{basictype}
                end
            end
            if b >= 8 # ImageMagick.jl does not read in sub 8 bit images correctly
                @testset "$(case): ImageMagick read values equality" begin
                    imdiff_val = imdiff(collect(read_in_pngf), read_in_immag)
                    onfail(@test imdiff_val <= get(imdiff_tolerance, case, 0.01)) do
                        PNGFiles._inspect_png_read(fpath)
                        _add_debugging_entry(fpath, case, imdiff_val)
                    end
                end
            end
            @testset "$(case): IO is idempotent" begin
                @test all(read_in_pngf .≈ PNGFiles.load(newpath))
                @test all(read_in_pngf .≈ open(io->PNGFiles.load(io), newpath))
            end
        end
    end

    ## TODO: Malformed pngs that should error. This throws `signal (6): Aborted` since we
    ## don't work with `png_jmpbuf` properly.
    # for fpath in glob(joinpath("./**/$(PNG_SUITE_DIR)", "[x]*.png"))
    #     case = splitpath(fpath)[end]
    #     @info case
    #     @testset "$(case)" begin
    #         @test_throws ErrorException PNGFiles.load(fpath)
    #     end
    # end
end


text:
# assumes that the plane in which the CoP lies is the x-y plane of the frame in which the wrench is expressed
function center_of_pressure(wrench::Wrench)
    T = eltype(wrench)
    normal = StaticArrays.SVector(zero(T), zero(T), one(T))
    torque = wrench.angular
    force = wrench.linear
    Point3D(wrench.frame, normal × torque / (normal ⋅ force))
end

max_stars_repo_name: tkoolen/VariableHeightMomentumControl

max_issues_repo_name: tkoolen/VariableHeightMomentumControl

max_forks_repo_name: tkoolen/VariableHeightMomentumControl


text:
struct Options
    nx::Int
    parallel::Union{Bool, Symbol}
    cudatype::Type
    cudabatchsize::Int
end
function Options(;nx::Int, parallel=false, cudatype::Type=Float64, cudabatchsize::Int=1)
    if !(parallel in (false, :threads, :cuda))
        throw(ArgumentError("parallel should be one of false, :threads, :cuda"))
    end
    Options(nx, parallel, cudatype, cudabatchsize)
end

# log unnormalized density
fk(m::Softcore, x, xbefore) = exp(logfk(m, x, xbefore))
function logfk(m::Softcore, x, xbefore)
  d((x1, y1), (x2, y2)) = sqrt((x1-x2)^2 + (y1-y2)^2)
  -sum((m.kernel(d(x,y)) for y in xbefore))
end

function quadps(window, nx)
  (x1,x2), (y1,y2) = window.x, window.y
  ny = ceil(Int, nx*(y2-y1)/(x2-x1))
  qx = quadps(x1, x2, nx)
  qy = quadps(y1, y2, ny)

  qx, qy, step(qx)*step(qy)
end
function quadps(x1, x2, n)
  s = (x2-x1)/(2n)
  LinRange(x1+s, x2-s, n)
end

function compute_integral_y!(m::Softcore, Ix, x, xs, qy)
  fill!(Ix, 0.0)
  for y in qy
    s1 = 0.0 #Base.TwicePrecision(0.0)
    for i in eachindex(xs)
      @inbounds x2, y2 = xs[i]
      d = sqrt((x-x2)^2 + (y-y2)^2)
      s1 += m.kernel_integral(d)
      @inbounds Ix[i] += exp(-Float64(s1))
    end
  end
  Ix
end

# normconstants(xs, f, window, nx::Int) = normconstants(xs, f, window, (nx=nx,))
function normconstants(m::Softcore, xs, window, intpar)
  intpar.nx == 0 && return ones(length(xs))
  qx, qy, w = quadps(window, intpar.nx)
  if intpar.parallel === false
    I = compute_integral(m, qx, qy, xs, (nx=intpar.nx,))
  elseif intpar.parallel === :threads
    I = compute_integral(m, qx, qy, xs, (nx=intpar.nx, threads=true))
  elseif intpar.parallel === :cuda
    if !cuda_available
      throw(ArgumentError("CUDA not loaded."))
    end
    I = compute_integral(m, qx, qy, xs, (nx=intpar.nx, type=intpar.cudatype,
      batchsize=intpar.cudabatchsize))
  else
    throw(ArgumentError("Invalid Options given."))
  end
  I .*= w
end
# Using separate accumulator for each grid line prodives a little better accuracy
# at the cost of doubling the memory requirement.
function compute_integral(m::Softcore, qx, qy, xs, ::NamedTuple{(:nx,)})
  I = zeros(length(xs))
  Ix = similar(I)
  for x in qx
    compute_integral_y!(m, Ix, x, xs, qy)
    I .+= Ix
  end
  I
end

function compute_integral(m::Softcore, qx, qy, xs, ::NamedTuple{(:nx, :threads)})
  I = [zeros(length(xs)) for _ in 1:Threads.nthreads()]
  Ix = [similar(I[1]) for _ in 1:Threads.nthreads()]
  Threads.@threads for x in qx
    compute_integral_y!(m, Ix[Threads.threadid()], x, xs, qy)
    I[Threads.threadid()] .+= Ix[Threads.threadid()]
  end
  sum(I)
end

max_stars_repo_name: mikkoku/SequentialPointProcesses

max_issues_repo_name: mikkoku/SequentialPointProcesses

max_forks_repo_name: mikkoku/SequentialPointProcesses


text:
struct LinearCombination{T, As<:Tuple{Vararg{LinearMap}}, Ts<:Tuple} <: LinearMap{T}
    maps::As
    coeffs::Ts
    function LinearCombination{T, As, Ts}(maps::As, coeffs::Ts) where {T, As, Ts}
        N = length(maps)
        N == length(coeffs) || error("Number of coefficients doesn't match number of terms")
        sz = size(maps[1])
        for n = 1:N
            size(maps[n]) == sz || throw(DimensionMismatch("LinearCombination"))
            promote_type(T, eltype(maps[n]), typeof(coeffs[n])) == T || throw(InexactError())
        end
        new{T, As, Ts}(maps, coeffs)
    end
end

LinearCombination{T}(maps::As, coeffs::Ts) where {T, As, Ts} = LinearCombination{T, As, Ts}(maps, coeffs)

# basic methods
Base.size(A::LinearCombination) = size(A.maps[1])
LinearAlgebra.issymmetric(A::LinearCombination) = all(issymmetric, A.maps) # sufficient but not necessary
LinearAlgebra.ishermitian(A::LinearCombination) = all(ishermitian, A.maps) && all(isreal, A.coeffs) # sufficient but not necessary
LinearAlgebra.isposdef(A::LinearCombination) = all(isposdef, A.maps) && all(isposdef, A.coeffs) # sufficient but not necessary

# adding linear maps
function Base.:(+)(A1::LinearCombination, A2::LinearCombination)
    size(A1) == size(A2) || throw(DimensionMismatch("+"))
    T = promote_type(eltype(A1), eltype(A2))
    return LinearCombination{T}(tuple(A1.maps..., A2.maps...), tuple(A1.coeffs..., A2.coeffs...))
end
function Base.:(+)(A1::LinearMap, A2::LinearCombination)
    size(A1) == size(A2) || throw(DimensionMismatch("+"))
    T = promote_type(eltype(A1), eltype(A2))
    return LinearCombination{T}(tuple(A1, A2.maps...), tuple(one(T), A2.coeffs...))
end
Base.:(+)(A1::LinearCombination, A2::LinearMap) = +(A2, A1)
function Base.:(+)(A1::LinearMap, A2::LinearMap)
    size(A1)==size(A2) || throw(DimensionMismatch("+"))
    T = promote_type(eltype(A1), eltype(A2))
    return LinearCombination{T}(tuple(A1, A2), tuple(one(T), one(T)))
end
function Base.:(-)(A1::LinearCombination, A2::LinearCombination)
    size(A1) == size(A2) || throw(DimensionMismatch("-"))
    T = promote_type(eltype(A1), eltype(A2))
    return LinearCombination{T}(tuple(A1.maps..., A2.maps...), tuple(A1.coeffs..., map(-, A2.coeffs)...))
end
function Base.:(-)(A1::LinearMap, A2::LinearCombination)
    size(A1) == size(A2) || throw(DimensionMismatch("-"))
    T = promote_type(eltype(A1), eltype(A2))
    return LinearCombination{T}(tuple(A1, A2.maps...), tuple(one(T), map(-, A2.coeffs)...))
end
function Base.:(-)(A1::LinearCombination, A2::LinearMap)
    size(A1) == size(A2) || throw(DimensionMismatch("-"))
    T = promote_type(eltype(A1), eltype(A2))
    return LinearCombination{T}(tuple(A1.maps..., A2), tuple(A1.coeffs..., -one(T)))
end
function Base.:(-)(A1::LinearMap, A2::LinearMap)
    size(A1) == size(A2) || throw(DimensionMismatch("-"))
    T = promote_type(eltype(A1), eltype(A2))
    return LinearCombination{T}(tuple(A1, A2), tuple(one(T), -one(T)))
end

# scalar multiplication
Base.:(-)(A::LinearMap) = LinearCombination{eltype(A)}(tuple(A), tuple(-one(eltype(A))))
Base.:(-)(A::LinearCombination) = LinearCombination{eltype(A)}(A.maps, map(-, A.coeffs))

function Base.:(*)(α::Number, A::LinearMap)
    T = promote_type(eltype(α), eltype(A))
    return LinearCombination{T}(tuple(A), tuple(α))
end
Base.:(*)(A::LinearMap, α::Number) = *(α, A)
function Base.:(*)(α::Number, A::LinearCombination)
    T = promote_type(eltype(α), eltype(A))
    return LinearCombination{T}(A.maps, map(x->α*x, A.coeffs))
end
Base.:(*)(A::LinearCombination, α::Number) = *(α, A)

function Base.:(\)(α::Number, A::LinearMap)
    T = promote_type(eltype(α), eltype(A))
    return LinearCombination{T}(tuple(A), tuple(1/α))
end
Base.:(/)(A::LinearMap, α::Number) = \(α, A)
function Base.:(\)(α::Number, A::LinearCombination)
    T = promote_type(eltype(α), eltype(A))
    return LinearCombination{T}(A.maps, map(x->α\x, A.coeffs))
end
Base.:(/)(A::LinearCombination, α::Number) = \(α, A)

# comparison of LinearCombination objects
Base.:(==)(A::LinearCombination, B::LinearCombination) = (eltype(A)==eltype(B) && A.maps==B.maps && A.coeffs==B.coeffs)

# special transposition behavior
LinearAlgebra.transpose(A::LinearCombination) = LinearCombination{eltype(A)}(map(transpose, A.maps), A.coeffs)
LinearAlgebra.adjoint(A::LinearCombination)   = LinearCombination{eltype(A)}(map(adjoint, A.maps), map(conj, A.coeffs))

# multiplication with vectors
function A_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector)
    # no size checking, will be done by individual maps
    A_mul_B!(y, A.maps[1], x)
    A.coeffs[1] == 1 || lmul!(A.coeffs[1], y)
    z = similar(y)
    for n=2:length(A.maps)
        A_mul_B!(z, A.maps[n], x)
        axpy!(A.coeffs[n], z, y)
    end
    return y
end
function At_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector)
    # no size checking, will be done by individual maps
    At_mul_B!(y, A.maps[1], x)
    A.coeffs[1] == 1 || lmul!(A.coeffs[1], y)
    z = similar(y)
    for n = 2:length(A.maps)
        At_mul_B!(z, A.maps[n], x)
        axpy!(A.coeffs[n], z, y)
    end
    return y
end
function Ac_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector)
    # no size checking, will be done by individual maps
    Ac_mul_B!(y, A.maps[1], x)
    A.coeffs[1] == 1 || lmul!(conj(A.coeffs[1]), y)
    z = similar(y)
    for n=2:length(A.maps)
        Ac_mul_B!(z, A.maps[n], x)
        axpy!(conj(A.coeffs[n]), z, y)
    end
    return y
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -140.37215752831116 <= q <= 140.41841564400983)

                     add_NL_constraint(m, :(log(1 + exp(0.9038900968817085 + 0.1284043831226005*log(1 + exp(0.06165362792774953 + 0.6866075626761576*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + -0.8737233458006566*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + -0.3942911833193632*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + 0.6872447034250162*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.5711071882051262*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + 0.26042276084039084*log(1 + exp(0.8058455478649482 + 0.7751639240502244*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + -0.13474005943220124*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + 0.6437669970537465*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + 0.5932192098608322*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.5325023046383324*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + 0.6458808355119907*log(1 + exp(-0.1540966309548315 + 0.690406062481447*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + 0.4718956794952942*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + 0.8703243890856593*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + -0.8079540211479861*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.683716372184894*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + -0.6567358944011805*log(1 + exp(-0.2486694730930381 + -0.9006272851005916*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + -0.9841514005441985*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + 0.7501520395123897*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + -0.43486617405885397*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + -0.4797454164387962*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + -0.41502277707459934*log(1 + exp(-0.9928729682863509 + 0.9588446333785003*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + 0.7154937103302039*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + -0.17921833371514895*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + -0.7009330672719543*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.15311755902026114*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + log(1 + exp(-0.8820041740138964 + 0.8813585352462416*log(1 + exp(0.06165362792774953 + 0.6866075626761576*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + -0.8737233458006566*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + -0.3942911833193632*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + 0.6872447034250162*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.5711071882051262*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + -0.30334290832683886*log(1 + exp(0.8058455478649482 + 0.7751639240502244*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + -0.13474005943220124*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + 0.6437669970537465*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + 0.5932192098608322*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.5325023046383324*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + 0.6303754843616343*log(1 + exp(-0.1540966309548315 + 0.690406062481447*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + 0.4718956794952942*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + 0.8703243890856593*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + -0.8079540211479861*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.683716372184894*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + 0.9344814932534002*log(1 + exp(-0.2486694730930381 + -0.9006272851005916*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + -0.9841514005441985*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + 0.7501520395123897*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + -0.43486617405885397*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + -0.4797454164387962*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + -0.8925753470795712*log(1 + exp(-0.9928729682863509 + 0.9588446333785003*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + 0.7154937103302039*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + -0.17921833371514895*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + -0.7009330672719543*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.15311755902026114*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + log(1 + exp(-0.6413973695460733 + -0.5291254814289021*log(1 + exp(0.06165362792774953 + 0.6866075626761576*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + -0.8737233458006566*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + -0.3942911833193632*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + 0.6872447034250162*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.5711071882051262*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + -0.3315514300310589*log(1 + exp(0.8058455478649482 + 0.7751639240502244*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + -0.13474005943220124*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + 0.6437669970537465*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + 0.5932192098608322*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.5325023046383324*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + 0.9076257425527015*log(1 + exp(-0.1540966309548315 + 0.690406062481447*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + 0.4718956794952942*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + 0.8703243890856593*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + -0.8079540211479861*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.683716372184894*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + -0.9763381142518175*log(1 + exp(-0.2486694730930381 + -0.9006272851005916*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + -0.9841514005441985*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + 0.7501520395123897*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + -0.43486617405885397*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + -0.4797454164387962*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + -0.06362970408750801*log(1 + exp(-0.9928729682863509 + 0.9588446333785003*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + 0.7154937103302039*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + -0.17921833371514895*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + -0.7009330672719543*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.15311755902026114*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + log(1 + exp(0.9122346082558659 + 0.2864216745322903*log(1 + exp(0.06165362792774953 + 0.6866075626761576*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + -0.8737233458006566*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + -0.3942911833193632*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + 0.6872447034250162*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.5711071882051262*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + 0.7014729451110155*log(1 + exp(0.8058455478649482 + 0.7751639240502244*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + -0.13474005943220124*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + 0.6437669970537465*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + 0.5932192098608322*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.5325023046383324*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + -0.27913045286088156*log(1 + exp(-0.1540966309548315 + 0.690406062481447*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + 0.4718956794952942*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + 0.8703243890856593*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + -0.8079540211479861*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.683716372184894*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + 0.014698309183167169*log(1 + exp(-0.2486694730930381 + -0.9006272851005916*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + -0.9841514005441985*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + 0.7501520395123897*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + -0.43486617405885397*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + -0.4797454164387962*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + -0.9872323021033105*log(1 + exp(-0.9928729682863509 + 0.9588446333785003*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + 0.7154937103302039*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + -0.17921833371514895*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + -0.7009330672719543*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.15311755902026114*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + log(1 + exp(0.6835137249354726 + 0.06184589975771315*log(1 + exp(0.06165362792774953 + 0.6866075626761576*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + -0.8737233458006566*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + -0.3942911833193632*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + 0.6872447034250162*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.5711071882051262*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + 0.6514470819403977*log(1 + exp(0.8058455478649482 + 0.7751639240502244*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + -0.13474005943220124*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + 0.6437669970537465*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + 0.5932192098608322*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.5325023046383324*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + 0.1907919388573478*log(1 + exp(-0.1540966309548315 + 0.690406062481447*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + 0.4718956794952942*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + 0.8703243890856593*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + -0.8079540211479861*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.683716372184894*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + 0.3162639033084287*log(1 + exp(-0.2486694730930381 + -0.9006272851005916*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + -0.9841514005441985*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + 0.7501520395123897*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + -0.43486617405885397*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + -0.4797454164387962*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))) + 0.2338844135686573*log(1 + exp(-0.9928729682863509 + 0.9588446333785003*log(1 + exp(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5]))) + 0.7154937103302039*log(1 + exp(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5]))) + -0.17921833371514895*log(1 + exp(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5]))) + -0.7009330672719543*log(1 + exp(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5]))) + 0.15311755902026114*log(1 + exp(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using GLMakie
using PolaronPathIntegrals

# Physical constants
const ħ = 1.05457162825e-34; # Reduced Planck's constant (kg m^2 s^{-1})
const eV = 1.602176487e-19; # Electron Volt (kg m^2 s^{-2})
const m_e = 9.10938188e-31; # Electron Mass (kg)
const k_B = 1.3806504e-23; # Boltzmann's constant (kg m^2 K^{-1} s^2)
const ϵ_0 = 8.854e-12 # Dielectric constant (C^2 N^{-1} m^{-2})
const c = 2.99792458e8 # Speed of light (m s^{-1})

fig = Figure(resolution = (1800, 900), backgroundcolor = RGBf0(0.98, 0.98, 0.98))

energy_axis = fig[1, 3] = Axis(fig, xlabel = "T [K]", ylabel = "F(T) [meV]", title = "Free Energy")
vary_axis = fig[2, 3] = Axis(fig, xlabel = "T [K]", ylabel = "v/w(T) [Hz]", title = "Variational Parameters")
mass_axis = fig[1, 4] = Axis(fig, xlabel = "T [K]", ylabel = "M(T) [m_e kg]", title = "Fictitious Mass")
spring_axis = fig[2, 4] = Axis(fig, xlabel = "T [K]", ylabel = "κ(T) [me kg/s^2]", title = "Fictitious Spring Constant")
mobility_axis_temp = fig[1, 2] = Axis(fig, xlabel = "T [K]", ylabel = "μ(T) [cm^2/Vs]", title = "Mobility vs Temperature")
absorption_axis_temp = fig[2, 2] = Axis(fig, xlabel = "T [K]", ylabel = "Γ(T) [cm^-1]", title = "Optical Absorption vs Temperature")
mobility_axis_freq = fig[1, 5] = Axis(fig, xlabel = "Ω [Hz]", ylabel = "μ(Ω) [cm^2/Vs]", title = "Mobility vs Frequency")
absorption_axis_freq = fig[2, 5] = Axis(fig, xlabel = "Ω [Hz]", ylabel = "Γ(Ω) [cm^-1]", title = "Optical Absorption vs Frequency")

lsgrid = labelslidergrid!(
    fig,
    ["ϵ_0", "ϵ_∞", "m", "ω"],
    [1:0.01:100, 1:0.01:100, 0.1:0.01:1, 0.1:0.01:10],
    formats = [x -> "$(round(x, digits = 3)) $s" for s in ["F/m", "F/m", "m_e", "THz"]],
    width = 1000
)

set_close_to!(lsgrid.sliders[1], 24.1)
set_close_to!(lsgrid.sliders[2], 4.5)
set_close_to!(lsgrid.sliders[3], 0.12)
set_close_to!(lsgrid.sliders[4], 2.25)

fig[3, 1:5] = lsgrid.layout

Ω_slider = Slider(fig[1:2, 1], range = 0.001:0.1:10.001, startvalue = 0.001, horizontal = false, tellwidth = true, height = nothing, width = Auto())
freq_sliders = [s.value for s in vcat(lsgrid.sliders , [Ω_slider])]
Ω_slider_label = Label(fig[1:2, 0], lift(s -> "Ω: $(round(s[], digits = 2)) Hz", Ω_slider.value), rotation = pi/2)

polaron_temp = lift(freq_sliders...) do slvalues...
    make_polaron(slvalues[2], slvalues[1], slvalues[4] * 1e12, slvalues[3]; temp = 100:10:400, efield_freq = slvalues[5], verbose = false)
end

T = lift(p -> p.T, polaron_temp)
v = lift(p -> Vector{Float64}(p.v), polaron_temp)
w = lift(p -> Vector{Float64}(p.w), polaron_temp)
κ = lift(p -> Vector{Float64}(p.κ), polaron_temp)
M = lift(p -> Vector{Float64}(p.M), polaron_temp)
F = lift(p -> Vector{Float64}(p.F), polaron_temp)
μ_T = lift(p -> vcat.(p.μ...)[1], polaron_temp)
Γ_T = lift(p -> vcat.(p.Γ...)[1], polaron_temp)

T_slider = Slider(fig[1:2, 7], range = 100:1:400, startvalue = 300, horizontal = false, tellwidth = true, height = nothing, width = Auto())
T_slider_label = Label(fig[1:2, 8], lift(s -> "T: $(s[]) K", T_slider.value), rotation = -pi/2)

temp_sliders = [s.value for s in vcat(lsgrid.sliders , [T_slider])]

polaron_freq = lift(temp_sliders...) do slvalues...
    make_polaron(slvalues[2], slvalues[1], slvalues[4] * 1e12, slvalues[3]; temp = slvalues[5], efield_freq = 0.001:0.5:10.001, verbose = false)
end

Ω = lift(p -> p.Ω, polaron_freq)
μ_Ω = lift(p -> Vector{Float64}(p.μ[1]), polaron_freq)
Γ_Ω = lift(p -> Vector{Float64}(p.Γ[1]), polaron_freq)

scatterlines!(mobility_axis_temp, T, μ_T)
scatterlines!(absorption_axis_temp, T, Γ_T)
scatterlines!(energy_axis, T, F)
v_plot = scatterlines!(vary_axis, T, v, label = "v")
w_plot = scatterlines!(vary_axis, T, w, label = "w", marker = :diamond)
scatterlines!(spring_axis, T, κ)
scatterlines!(mass_axis, T, M)

axislegend(vary_axis, [v_plot, w_plot], ["v", "w"], position = :rb)

on(polaron_temp) do x
    autolimits!(mobility_axis_temp)
    autolimits!(absorption_axis_temp)
    autolimits!(energy_axis)
    autolimits!(vary_axis)
    autolimits!(spring_axis)
    autolimits!(mass_axis)
end

scatterlines!(mobility_axis_freq, Ω, μ_Ω)
scatterlines!(absorption_axis_freq, Ω, Γ_Ω)

on(polaron_freq) do x
    autolimits!(mobility_axis_freq)
    autolimits!(absorption_axis_freq)
end

fig

# μ_Ω = Vector{Any}()
# Γ_Ω = Vector{Any}()
# for i in 1:Int(floor(length(T[])/(N - 1))):length(T[])
#     push!(μ_Ω, lift(p -> Vector{Float64}(p.μ[i]), polaron))
#     push!(Γ_Ω, lift(p -> Vector{Float64}(p.Γ[i]), polaron))
# end
# μ_T = Vector{Any}()
# Γ_T = Vector{Any}()
# for i in 1:Int(floor(length(Ω[])/(N - 1))):length(Ω[])
#     push!(μ_T, lift(p -> vcat.(p.μ...)[i], polaron))
#     push!(Γ_T, lift(p -> vcat.(p.Γ...)[i], polaron))
# end
# for (μ, Γ) in zip(μ_T, Γ_T)
#     scatterlines!(mobility_axis, T, μ, colormap = :viridis)
#     scatterlines!(absorption_axis, T, Γ, color = :blue)
# end


text:
using Juno
import Zygote: Params, gradient

"""
    update!(x, x̄)

Update the array `x` according to `x .-= x̄`.
"""
function update!(x::AbstractArray, x̄)
  x .-= x̄
end

"""
    update!(opt, p, g)
    update!(opt, ps::Params, gs)

Perform an update step of the parameters `ps` (or the single parameter `p`)
according to optimizer `opt`  and the gradients `gs` (the gradient `g`).

As a result, the parameters are mutated and the optimizer's internal state may change.
The gradient could be mutated as well.
"""
function update!(opt, x, x̄)
  x̄r = ArrayInterface.restructure(x, x̄) # address some cases where Zygote's
                                          # output are not mutable, see #1510
  x .-= apply!(opt, x, x̄r)
end

function update!(opt, xs::Params, gs)
  for x in xs
    isnothing(gs[x]) && continue
    update!(opt, x, gs[x])
  end
end

# Callback niceties
call(f, xs...) = f(xs...)
runall(f) = f
runall(fs::AbstractVector) = () -> foreach(call, fs)

struct SkipException <: Exception end

"""
    skip()

Call `Flux.skip()` in a callback to indicate when a callback condition is met.
This will trigger the train loop to skip the current data point and not update with the calculated gradient.

# Examples
```julia
cb = function ()
  loss() > 1e7 && Flux.skip()
end
```
"""
function skip()
  throw(SkipException())
end

struct StopException <: Exception end

"""
    stop()

Call `Flux.stop()` in a callback to indicate when a callback condition is met.
This will trigger the train loop to stop and exit.

# Examples
```julia
cb = function ()
  accuracy() > 0.9 && Flux.stop()
end
```
"""
function stop()
  throw(StopException())
end

batchmemaybe(x) = tuple(x)
batchmemaybe(x::Tuple) = x

"""
    train!(loss, params, data, opt; cb)

For each datapoint `d` in `data`, compute the gradient of  `loss` with
respect to `params` through backpropagation and call the optimizer `opt`.

If `d` is a tuple of arguments to `loss` call `loss(d...)`, else call `loss(d)`.

A callback is given with the keyword argument `cb`. For example, this will print
"training" every 10 seconds (using [`Flux.throttle`](@ref)):
    train!(loss, params, data, opt, cb = throttle(() -> println("training"), 10))

The callback can call [`Flux.stop`](@ref) to interrupt the training loop.

Multiple optimisers and callbacks can be passed to `opt` and `cb` as arrays.
"""
function train!(loss, ps, data, opt; cb = () -> ())
  ps = Params(ps)
  cb = runall(cb)
  @progress for d in data
    try
      gs = gradient(ps) do
        loss(batchmemaybe(d)...)
      end
      update!(opt, ps, gs)
      cb()
    catch ex
      if ex isa StopException
        break
      elseif ex isa SkipException
        continue
      else
        rethrow(ex)
      end
    end
  end
end

"""
    @epochs N body

Run `body` `N` times. Mainly useful for quickly doing multiple epochs of
training in a REPL.

# Examples
```jldoctest
julia> Flux.@epochs 2 println("hello")
[ Info: Epoch 1
hello
[ Info: Epoch 2
hello
```
"""
macro epochs(n, ex)
  :(@progress for i = 1:$(esc(n))
      @info "Epoch $i"
      $(esc(ex))
    end)
end


text:
module Church
using Distributions

export value, resample, gc_church, background_sampler, n_samples, pdf, RV, @lift, Det, nocond

#Types
type NoCond
end
const nocond = NoCond()
type NoSampler
end
abstract RV
const nosampler = NoSampler()
type Sample{V, S<:Union(NoSampler, Function)} <: RV
    det::Any
    dist::Distribution
    value::V
    logp::Float64
    deps::Vector
    sampler::S
end
type Det{F, T <: Tuple} <: RV
    f::F
    args::T
    deps::Vector
end
type GetIndex{S, T <: Tuple} <: RV
    struct::S
    args::T
    deps::Vector
    #Value only used in GC phase.
    value::Any
end
WRV = Union(WeakRef, RV)

#Store declared samples/branches
const samples = Array(Union(Sample, WeakRef), 0)
const dets = Array(Union(Det, WeakRef), 0)
const getindexes = Array(Union(GetIndex, WeakRef), 0)
const conditions = Array(Union(Sample, WeakRef), 0)

#Show functions
import Base.show
show(io::IO, s::RV) = print(io, "RV(", value(s), ")")

#Constructors
fallback_sampler(d::Distribution, v) = rand(d)
import Distributions.sample
sample(det, cond::NoCond, sampler::Union(NoSampler, Function)) = begin
    dist = value(det)
    val = rand(dist)
    lp = logp(dist, val)
    s = Sample(det, dist, val, lp, Array(WRV, 0), sampler)
    push!(samples, s)
    add_dep(s, det)
    s
end
sample(det, val, sampler::Union(NoSampler, Function)) = begin
    dist = value(det)
    lp = logp(dist, val)
    s = Sample(det, dist, val, lp, Array(WRV, 0), sampler)
    push!(conditions, s)
    add_dep(s, det)
    s
end
Det(f, args) = begin
    d = Det(f, args, Array(WRV, 0))
    push!(dets, d)
    add_dep(d, f)
    for arg in args
        add_dep(d, arg)
    end
    d
end

GetIndex(struct, args...) = begin
    g = GetIndex(struct, args, Array(WRV, 0), nothing)
    push!(getindexes, g)
    add_dep(g, struct[map(value, args)...])
    for arg in args
        add_dep(g, arg)
    end
    g
end
#Allow indexing into multidimensional samples.
getindex(s::Sample, args...) = Det((s_, args_...) -> s_[args_...], tuple(s, args...))

#Lifting functions.
a(i::Int) = symbol("a$i")
type_a(i::Int, b::Bool) =
    b ? :($(a(i))::RV) : a(i)
lift(f::Symbol, issdg::Vector{Bool}) =
    Expr(:(=), Expr(:call, f, {type_a(i,  issdg[i]) for i = 1:length(issdg)}...),
               Expr(:call, :(Church.Det), f, Expr(:tuple, map(a, 1:length(issdg))...)))
lift_gi(typ, issdg::Vector{Bool}) =
    Expr(:(=), Expr(:call, :getindex, :(s::$typ),{type_a(i,  issdg[i]) for i = 1:length(issdg)}...),
               Expr(:call, :GetIndex, :s, map(a, 1:length(issdg))...))
lift(f, n::Int, inner_func::Function) = begin
    defs = Expr[]
    for nsamps in n:(-1):1
        for c in combinations(1:n, nsamps)
            issdg = fill(false, n)
            for i in c
                issdg[i] = true
            end
            push!(defs, inner_func(f, issdg))
        end
    end
    Expr(:block, defs...)
end
macro lift_gi(typ, n)
    esc(lift(typ, n, lift_gi))
end
macro lift(f, n)
    esc(quote
        if Main == Base.function_module($f)
            $(lift(f, n, lift))
        else
            eval(parse(string("import ", Base.function_module($f), ".", $f)))
            $(lift(f, n, lift))
        end
    end)
end

add_dep(s::RV, arg::RV) =
    if !in(s, arg.deps)
        push!(arg.deps, s)
    end
add_dep(s::RV, arg) = nothing

#Get the value of an expression.
value(s::Sample) = s.value
value(det::Det) =
    value(det.f)(map(value, det.args)...)
value(g::GetIndex) = begin
    #Should stop previous arg from being dependent.
    res = g.struct[map(value, g.args)...]
    #Add g as a dependent of struct [arg]
    add_dep(g, res)
    value(res)
end
value(s) = s

#Don't care whether cdf or pdf is needed.
logp(d::ContinuousDistribution, x) = logpdf(d, x)
logp(d::DiscreteDistribution, x) = logpdf(d, x)

#Find dependents
deps_inner(_, s::Sample, deps::Vector{Sample}) = begin
    if !in(s, deps)
        push!(deps, s)
    end
    nothing
end
deps_inner(_, s::Det, deps::Vector{Sample}) =
    deps_recurse(s, deps)
deps_inner(origin, s::GetIndex, deps::Vector{Sample}) =
    if in(origin, s.args) || (origin == s.struct[map(value, s.args)...])
        deps_recurse(s, deps)
    end
deps_recurse(s, deps::Vector{Sample}) = begin
    for dep in s.deps
        deps_inner(s, dep, deps)
    end
    nothing
end

#Remove redundant dependents.
remove_deps(origin::RV) =
    filter!(x -> isdependent(origin, x), origin.deps)
remove_deps(as::Vector) = begin
    for a in as
        remove_deps(a)
    end
    nothing
end

#Propose, accept/reject.
resample_inner(s::Sample) = begin
    old_val  = s.value
    deps = Sample[]
    deps_recurse(s, deps)
    (s.value, old_logp, new_logp) = propose(s.dist, old_val, s.sampler)
    old_logp += length(deps) == 0 ? 0. : mapreduce(dep->dep.logp, +, deps)
    new_dists = map(s -> value(s.det), deps)
    new_logps = map((d, s) -> logp(d, s.value), new_dists, deps)
    new_logp += length(new_logps) == 0 ? 0. : sum(new_logps)
    if !(exp(new_logp - old_logp) > rand())
        #Reject change
        s.value = old_val
    else
        #Accept changes
        for i = 1:length(deps)
            deps[i].dist = new_dists[i]
            deps[i].logp = new_logps[i]
        end
    end
    nothing
end
propose(d::Distribution, _, s::NoSampler) =
    (rand(d), 0., 0.)
propose(d::Distribution, old_val, s::Function) = begin
    prop_dist_forward = s(d, old_val)
    prop_val = rand(prop_dist_forward)
    prop_dist_backward = s(d, prop_val)
    (prop_val,
     logp(d, old_val )-logp(prop_dist_backward, old_val ),
     logp(d, prop_val)-logp(prop_dist_forward , prop_val))
end
resample() =
    if length(samples) != 0
        index = rand(1:length(samples))
        resample_inner(samples[index])
    end
resample(iter::Int) =
    for i = 1:iter
        resample()
    end

#Run sampling + garbage collection in background.
background_sampler(gc=true) =
    @async while true
        tic()
        resample()
        sleep(1E-6)
        if gc && rand() < toq()
            gc_church()
        end
    end

#Define distributions.
for dist in {:MvNormal, :MvTDist, filter!(isleaftype, subtypes(Distribution))...}
    ssym = match(r"[a-zA-Z]+", string(dist)).match
    tsym = symbol(ssym)
    fsym = symbol(lowercase(ssym))
    eval(quote
        $fsym(args...; condition=nocond, sampler=nosampler) = begin
            det = any(issdg, args) ? Det($tsym, args) : $tsym(args...)
            sample(det, condition, sampler)
        end
    end)
    eval(parse("export $(string(fsym))"))
end
issdg(x::RV) = true
issdg(x) = false

#Overload operators.
for op in [+, -, *, /, \, .*, ./, .\]
    op(a::RV, b::RV) = Det(op, (a, b))
    op(a::RV, b) = Det(op, (a, b))
    op(a, b::RV) = Det(op, (a, b))
    op(a::RV) = Det(op, (a,))
end
include("gc.jl")
include("datastructures.jl")
end


text:
using Test, ABG

function test_my_styled_printing()
    for i in 1:256
        printstyled(' ', lpad("$i", 3, '0'), ' ', color=i)
        (i % 16 == 0) && println()
    end
    println()

    println("\nAn average messages using color default\n")
    ABG.title("This is a title using bold, color cyan and 154")
    ABG.subtitle("This is a subtitle using color cyan and 154")
    ABG.item("Item 1 using color 74")
    ABG.item("Item 2")
    ABG.message("A description message using color :light_magenta")
    ABG.warning("This is a warning using color 229", "Another warning message")
    ABG.errormsg("One error", "Another error")
    ABG.println()
    ABG.done("Message done, Using color 40.")
    ABG.message("All message tests should be OK")
end

@test 1==1
test_my_styled_printing()


text:
#New for 0.7, splat changed syntax.
#TODO: Merge with newplot method below
newplot(itr::Base.Iterators.Pairs, args...; kwargs...) = newplot(collect(itr), args...; kwargs...)

#Common code across all plots
function newplot(kwargs::Vector{T}; ec_charttype::Union{String, Nothing} = nothing) where T

    #Create new chart
    ec = deepcopy(EChart(ec_charttype = ec_charttype))

    #Process keyword args after defined functionality
	kwargs!(ec, kwargs)

    return ec

end

#Separate arrays as inputs to array of dicts
function arrayofdicts(; kwargs...)

    #Determine number of values in array, check that all are equal length
    lengths = [length(x[2]) for x in kwargs]
    all(x -> x == lengths[1], lengths) ? nothing : error("All arrays need to have same length.")

    #Parse from arrays to dicts
    k = [x[1] for x in kwargs]
    v = [x[2] for x in kwargs]
    datafmt = []

    for j in 1:lengths[1]
        temp = Dict()
        for i in 1:length(k)
            temp[string(k[i])] = v[i][j]
        end
        push!(datafmt, temp)
    end

    return datafmt
end

#Combine arrays into array of arrays by each value
arrayofarray(x::AbstractVector,y::AbstractVector) = [[x,y] for (x,y) in zip(x,y)]
arrayofarray(x::AbstractVector,y::AbstractVector,z::AbstractVector) = [[x,y,z] for (x,y,z) in zip(x,y,z)]

#Common kwargs code for all plots
#For convenience, let color be specified as a string, even though it's always an array in echarts.js
function kwargs!(ec::EChart, kwargs::Vector{T}) where T

	for (k, v) in kwargs
        k == :color && typeof(v) in [String, JSFunction] ? setfield!(ec, k, [v]) : setfield!(ec, k, v)
    end

end

# Fill area inside areaStyle
function fill!(ec::EChart, cols::Int, fill::Union{Bool, Vector})

    typeof(fill) == Bool ? fill = [fill for i in 1:cols] : fill = fill

    for i in 1:cols
        fill[i] ? ec.series[i].areaStyle = ItemStyle() : nothing
    end

end

function boxplotstat(data::AbstractVector{<:Union{Missing, Int, AbstractFloat, Rational}})

    data_ = collect(skipmissing(data))
    #Calculate stats
    ss = summarystats(data_)
    iqr15 = 1.5 * (ss.q75 - ss.q25)
    upperbound = ss.q75 + iqr15
    lowerbound = ss.q25 - iqr15

    #Calculate outliers for scatterplot
    outliers = filter(x -> (x >= upperbound) || (x <= lowerbound), data_)

    return BoxPlotStats([lowerbound, ss.q25, ss.median, ss.q75, upperbound], outliers)

end


text:
module MathLinkExtras

using MathLink

MLTypes=Union{MathLink.WExpr,MathLink.WSymbol}

include("Operators.jl")
include("Wprint.jl")

end # module

max_stars_repo_name: fremling/MathLinkExtras

max_issues_repo_name: fremling/MathLinkExtras


text:
using ABM4OSN

p_1 = ABM4OSN.Post(0.5, 0.2, 1, 1)
p_2 = ABM4OSN.Post(0.5, 0.1, 1, 1)

@test p_1 > p_2


text:
function area_triangle(a,b,c)
a,b,c=reverse(sort([a,b,c]))
area = .25*sqrt((a+(b+c))*(c-(a-b))*(c+(a-b))*(a+(b-c)))
return area
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 5ca81b3a-7926-11eb-38e4-37053a873b78
begin
	using Markdown
	using InteractiveUtils
	using Images
	using Plots
	using PlutoUI
end

# ╔═╡ 8d041392-7926-11eb-308b-55a47a287cea
using JImages

# ╔═╡ f590466a-7926-11eb-18d7-3109f09d4f85

# ╔═╡ cab9623c-7926-11eb-024e-0fe16e3170d6
begin
	path = string(datadir, "/", "lena.h5")
	println("path to Lena is $path")
	lena = load_h5(path, "lena")
end

# ╔═╡ 02f670b8-7927-11eb-3d2e-8bb684cc0747
kg1 = Kernel.gaussian(1)

# ╔═╡ 6ae9cd0a-7927-11eb-256d-adfecbd5f27e
@time lena_g1 = imfilter(lena, reflect(centered(kg1)), Fill(0))

# ╔═╡ 85b1ec0a-7927-11eb-0d1e-8f132863920d
Gray.(lena)

# ╔═╡ aea0934e-7927-11eb-1025-ef5ebadbfb4a
iter = 100

# ╔═╡ b27aed28-7927-11eb-3e8a-17e5edee2584
md" Running RL with $iter iterations"

# ╔═╡ 7da73da2-7928-11eb-08cd-39c737d26622
@time lena_d, DFR = richardson_lucy(lena_g1, kg1,
                                    iterations=iter,
                                    view_window =1,
                                    fom_window=64,
                                    debug_level=quiet)

# ╔═╡ b2a9fe70-7928-11eb-1e88-d9749038d9ab
DFR

# ╔═╡ 3808c506-7933-11eb-3b76-1bcc9c18b823
plot(DFR[3:end])

# ╔═╡ a11f9b24-7928-11eb-0259-53dc3facfaae
md"DFR shows the variation of the relative difference. About 50 iterations are enough"

# ╔═╡ c23a1db0-7933-11eb-19bb-e7670600f387
Gray.(lena_d)

# ╔═╡ ccc3275e-7933-11eb-3b6d-8f88a26c5390
md"Running the fast version"

# ╔═╡ fef6fd5c-7933-11eb-2c9c-977ec2ea4a55
@time lena_d2 = frichardson_lucy(lena_g1, kg1, iterations=iter)

# ╔═╡ 962d09b2-7933-11eb-0cc0-f7f5cdf38116
md"## An alternative deconvolution"

# ╔═╡ a58bf18e-7d0c-11eb-2b5b-f52987f7c47a
md"### Gaussian kernel"

# ╔═╡ f60fab46-7d0c-11eb-1fc3-05956cd1a30b
function pprint(x)
	with_terminal() do
    	println(x)
	end
end

# ╔═╡ b7b4b8d0-7d0c-11eb-257e-8d2f1c78c9bd
begin
	G(x,y, σ=0.5) = 1 / sqrt(2πσ^2) * e^(-(x^2 + y^2) / (2σ^2) )
	G(I::CartesianIndex{2}, σ=0.5) = 1 / sqrt(2πσ^2) * e^(-(x^2 + y^2) / (2σ^2) )
end

# ╔═╡ 8c12e61c-7d0d-11eb-0955-f35ae7f18e89
zz = [[10., 30.] [20.,40.] ]

# ╔═╡ 12906938-7d1a-11eb-1560-9bb1872423e5
begin
	t(x,y) = x + y
	TT(I::CartesianIndex{2}) = t(I.I...)
end

# ╔═╡ 55d11c48-7d0f-11eb-22a3-d129c2a7e5e1
begin
	for i in CartesianIndices(zz)
          println(" cartesian index i = ", i ," value = ", zz[i])
		  println(" i.I = ", i.I... ," value = ", zz[i])
			println(" i.I = ", i.I... ," t(i.I...) = ", t(i.I...))
	end
end

# ╔═╡ 23619db8-7d1a-11eb-2144-c9766bda5ae5
begin
	for i in CartesianIndices(zz)
          println(" cartesian index i = ", i ," value = ", zz[i])
			println(" TT($i) = ", TT(i))
	end
end

# ╔═╡ cff0e784-7d27-11eb-1c8a-c73c55c769cd
l = 4 * ceil(Int, 0.5) + 1

# ╔═╡ 02791abe-7d28-11eb-00f0-a90f55eee923
w = 5÷2

# ╔═╡ 1bb6d316-7d28-11eb-18e7-fb4a5fcc031e
CartesianIndices((-w:w, -w:w))

# ╔═╡ 30c12eb4-7d28-11eb-35e0-83ae2b97ea51
bar(a,x...) = (a,b,x)

# ╔═╡ d5ba26e4-7d27-11eb-2b71-81b1385785b3

# ╔═╡ Cell order:
# ╠═5ca81b3a-7926-11eb-38e4-37053a873b78
# ╠═8d041392-7926-11eb-308b-55a47a287cea
# ╠═f590466a-7926-11eb-18d7-3109f09d4f85
# ╠═cab9623c-7926-11eb-024e-0fe16e3170d6
# ╠═02f670b8-7927-11eb-3d2e-8bb684cc0747
# ╠═6ae9cd0a-7927-11eb-256d-adfecbd5f27e
# ╠═85b1ec0a-7927-11eb-0d1e-8f132863920d
# ╠═aea0934e-7927-11eb-1025-ef5ebadbfb4a
# ╠═b27aed28-7927-11eb-3e8a-17e5edee2584
# ╠═7da73da2-7928-11eb-08cd-39c737d26622
# ╠═b2a9fe70-7928-11eb-1e88-d9749038d9ab
# ╠═3808c506-7933-11eb-3b76-1bcc9c18b823
# ╠═a11f9b24-7928-11eb-0259-53dc3facfaae
# ╠═c23a1db0-7933-11eb-19bb-e7670600f387
# ╠═ccc3275e-7933-11eb-3b6d-8f88a26c5390
# ╠═fef6fd5c-7933-11eb-2c9c-977ec2ea4a55
# ╠═962d09b2-7933-11eb-0cc0-f7f5cdf38116
# ╠═a58bf18e-7d0c-11eb-2b5b-f52987f7c47a
# ╠═f60fab46-7d0c-11eb-1fc3-05956cd1a30b
# ╠═b7b4b8d0-7d0c-11eb-257e-8d2f1c78c9bd
# ╠═8c12e61c-7d0d-11eb-0955-f35ae7f18e89
# ╠═12906938-7d1a-11eb-1560-9bb1872423e5
# ╠═55d11c48-7d0f-11eb-22a3-d129c2a7e5e1
# ╠═23619db8-7d1a-11eb-2144-c9766bda5ae5
# ╠═cff0e784-7d27-11eb-1c8a-c73c55c769cd
# ╠═02791abe-7d28-11eb-00f0-a90f55eee923
# ╠═1bb6d316-7d28-11eb-18e7-fb4a5fcc031e
# ╠═30c12eb4-7d28-11eb-35e0-83ae2b97ea51
# ╠═d5ba26e4-7d27-11eb-2b71-81b1385785b3

max_stars_repo_name: jjgomezcadenas/JImages

max_issues_repo_name: jjgomezcadenas/JImages

max_forks_repo_name: jjgomezcadenas/JImages


text:
# ---
# title: 1603. Design Parking System
# id: problem1603
# author: Qling
# date: 2021-03-25
# difficulty: Easy
# categories: Design
# link: <https://leetcode.com/problems/design-parking-system/description/>
# hidden: true
# ---
#
# Design a parking system for a parking lot. The parking lot has three kinds of
# parking spaces: big, medium, and small, with a fixed number of slots for each
# size.
#
# Implement the `ParkingSystem` class:
#
#   * `ParkingSystem(int big, int medium, int small)` Initializes object of the `ParkingSystem` class. The number of slots for each parking space are given as part of the constructor.
#   * `bool addCar(int carType)` Checks whether there is a parking space of `carType` for the car that wants to get into the parking lot. `carType` can be of three kinds: big, medium, or small, which are represented by `1`, `2`, and `3` respectively. **A car can only park in a parking space of its**`carType`. If there is no space available, return `false`, else park the car in that size space and return `true`.
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["ParkingSystem", "addCar", "addCar", "addCar", "addCar"]
#     [[1, 1, 0], [1], [2], [3], [1]]
#     **Output**
#     [null, true, true, false, false]
#
#     **Explanation**
#     ParkingSystem parkingSystem = new ParkingSystem(1, 1, 0);
#     parkingSystem.addCar(1); // return true because there is 1 available slot for a big car
#     parkingSystem.addCar(2); // return true because there is 1 available slot for a medium car
#     parkingSystem.addCar(3); // return false because there is no available slot for a small car
#     parkingSystem.addCar(1); // return false because there is no available slot for a big car. It is already occupied.
#
#
#
#
# **Constraints:**
#
#   * `0 <= big, medium, small <= 1000`
#   * `carType` is `1`, `2`, or `3`
#   * At most `1000` calls will be made to `addCar`
#
#
## @lc code=start
using LeetCode

struct ParkingSystem
    A::Vector{Int}
    ParkingSystem(big, median, small) = new([big, median, small])
end

function add_car(parking::ParkingSystem, car_type::Int)::Bool
    parking.A[car_type] -= 1
    return parking.A[car_type] >= 0
end
## @lc code=end


text:
module EvolutionaryStrategies

using Cambrian
using Random
using LinearAlgebra

abstract type ESState end

include("individual.jl")
include("exponential_nes.jl")
include("separable_nes.jl")

end


text:
let
	# if maxlayers is redefined, must remake the calibration cache or run without it (set usecache=false)
	global const maxlayers = 50			# 15-100 (model gets unstable if too high relative to time step, see stabilitycheck())

	secondsperyear	= 365*24*60*60		# [s/year]
	k				= 1.5/1e4			# thermal diffusivity of deep ocean [m^2/s]  (default 1.5 cm^2/s)
	w				= 4/secondsperyear	# mean ocean upwelling rate [m/s]  (default 4 m/yr)
	polarwaterratio	= 0.2				# ratio of Δtemp of downwelling polar water to mean Δtemp of surface ocean [unitless]
	mixeddepth		= 70.0				# depth of mixed layer [m]
	rho				= 1026				# density of seawater [kg/m3]
	c				= 3996				# specific heat capacity of seawater [J/kgK]
	oceanarea		= 0.71				# fraction of earth's surface covered by ocean [unitless]
	a_L_O_exch		= 0.31
	b_L_O_exch		= 1.59
	marinewarming 	= 1.3				# marine surface air warming enhancement [unitless]
	oceandepth		= 3688				# average depth of the ocean [m]

	zstep			= (oceandepth - mixeddepth)/maxlayers	# vertical length of each discrete layer below surface layer [m]
	midlayers 		= 2:maxlayers-1

	function derivative(T, step, position)
		if position == 1
			# This one-sided three point approximation is more accurate than (T[2] - T[1])/step.
			# It is O(step^2) instead of O(step).
			(-3*T[1] + 4*T[2] - T[3]) / (2*step)
		elseif position == length(T)
			# ditto
			(3*T[end] - 4*T[end-1] + T[end-2]) / (2*step)
		else
			# This central difference is actually also a three point approximation, and also O(step^2).
			(T[position+1] - T[position-1]) / (2*step)
		end
	end

	function secondderivative(T, step, position)
		# The position check is quite expensive, so let's comment it out and make sure we only call this in the interior.
		# If at the extremes, calculate from position one step inward.
		#position = clamp.(position, 2, length(T)-1)

		# This central three point approximation is also correct for the extreme left or right (if shifted inward one step).
		# It is O(step^2) in the interior and O(step) at the extremes.
		(T[position-1] - 2*T[position] + T[position+1]) / step^2
	end

	# Calculate ocean temperature using an upwelling-diffusion energy balance model.
	# Solve the advection-diffusion differential equation using explicit Euler finite differences.
	# The discretization is forward-in-time, central-in-space (FTCS).
	# The function does only one iteration of the time step but calculates temperature in all ocean layers (vectorized).
	global function temperatures!(s::ClimateState, p::ClimateParams)
		@unpack TotalRadiativeForcing, Temp, Temp_land, Temp_global = s
		@unpack timestep, lambda = p

		lambda_L = lambda
		lambda_O = lambda/marinewarming
		k_landocean = 1 # + b_L_O_exch - a_L_O_exch*lambda_L		#  land-ocean heat exchange coefficient [(W/m^2)/K]

		# We need some spatial derivatives of the temperature profile
		dz_Temp_1 = derivative(Temp, zstep, 1)
		dz_Temp = derivative.(Ref(Temp), zstep, midlayers)
		dz2_Temp = secondderivative.(Ref(Temp), zstep, midlayers)

		# First boundary condition (type Dirichlet):
		# Temp_mixed [°C]: temperature (change since preindustrial times) of mixed layer (i.e. surface ocean)
		Temp_mixed = Temp[1]

		# Temp_land [°C]: land temperature (change since preindustrial times)
		Temp_land = (TotalRadiativeForcing + k_landocean/(1-oceanarea)*marinewarming*Temp_mixed) /
						(1/lambda_L + k_landocean/(1-oceanarea))

		# heat transfer from land to ocean [W/m^2]
		Heatflux_land2ocean = k_landocean * (Temp_land - marinewarming*Temp_mixed)

		# heat transfer from mixed layer to deep ocean [W/m^2]
		Heatflux_mixed2deep = -rho*c*k * dz_Temp_1 - rho*c*w * (1-polarwaterratio) * Temp_mixed
		#Heatflux_mixed2deep = -rho*c*k * dz_Temp_1 - rho*c*w * (Temp[2] - polarwaterratio*Temp_mixed)

		# time derivative of temperature [°C/year]
		# timestep unit is [per year], so all time derivatives are also [per year]
		dt_Temp = Vector{Float64}(undef, maxlayers)

		# Energy balance for the mixed layer
		dt_Temp[1] = secondsperyear / (rho*mixeddepth*c) *
			(TotalRadiativeForcing - Temp_mixed/lambda_O - Heatflux_mixed2deep + Heatflux_land2ocean/oceanarea)

		# Main advection-diffusion differential equation:
		# We could speed things up a tiny bit by assigning Temp[midlayers+1] and Temp[midlayers-1] to variables, but meh.
		dt_Temp[midlayers] = secondsperyear * (k * dz2_Temp + w * dz_Temp)

		# Temp [°C]: temperature change since preindustrial times in each ocean layer
		Temp[1:maxlayers-1] += timestep * dt_Temp[1:maxlayers-1]

		# Second boundary condition (type Robin): approximate the derivative to get an equation for Temp[max]
		dz_Temp_max = w / k * (polarwaterratio*Temp_mixed - Temp[maxlayers])
		Temp[maxlayers] = Temp[maxlayers-1] + zstep * dz_Temp_max
		# more accurate approximation below, but accuracy is not important at ocean bottom
		#Temp[maxlayers] = -4*Temp[maxlayers-1] + Temp[maxlayers-2] + 2*3*zstep * dz_Temp_maxlayers

		# Temp_global [°C]: average land & ocean temperature change since preindustrial times
		Temp_global = marinewarming*oceanarea*Temp[1] + (1-oceanarea)*Temp_land

		@pack! s = Temp, Temp_land, Temp_global
	end

	global function init_temperatures!(s::ClimateState)
		Temp = zeros(maxlayers)
		Temp_global = 0.0
		Temp_land = 0.0
		@pack! s = Temp, Temp_global, Temp_land
	end

	global function stabilitycheck(timestep)
		# Courant number, see p56 of Kajishima (2017), Computational Fluid Dynamics
		courant = w * timestep * secondsperyear / zstep

		# cell Reynolds number, see Kajishima p59-60
		reynolds = w * zstep / k

		# stability relation
		@assert courant <= min(reynolds/2, 2/reynolds) "Instable discretization of temperature differential equation, " *
			"try reducing timestep or increasing maxlayers."
	end
end

max_stars_repo_name: niclasmattsson/ClimateCalculator

max_issues_repo_name: niclasmattsson/ClimateCalculator

max_forks_repo_name: niclasmattsson/ClimateCalculator