text:
module CausalityTools

    using DynamicalSystems
    using StaticArrays

    export Dataset, trajectory, DiscreteDynamicalSystem, ContinuousDynamicalSystem
    export SVector, @SVector, MVector, @MVector
    using Requires

    using Reexport
    @reexport using CausalityToolsBase
    @reexport using StateSpaceReconstruction
    @reexport using PerronFrobenius
    @reexport using TransferEntropy
    @reexport using CrossMappings

    # Modules defined in CausalityTools
    include("JointDistanceDistribution/JointDistanceDistribution.jl")

    display_update = true
    version = "v0.9.0"
    update_name = "update_$version"

    if display_update
    if !isfile(joinpath(@__DIR__, update_name))
    printstyled(stdout,
    """
    \nUpdate message: CausalityTools $version
    ----------------
    New functionality
    ----------------
    #### TransferEntropy.jl

    - Added [`NearestNeighbourMI`](@ref) transfer entropy estimator.
    - Added [`transferentropy(::Any, ::TEVars, ::NearestNeighbourMI)`](@ref) method.
    - Transfer entropy estimators now contain a field `b` which gives the base of the logarithm
        used during transfer entropy computations, and hence dictates the unit of the transfer
        entropy. By default, `b = 2`, which gives the transfer entropy in bits. The keyword `b`
        is thus obsolete in all transfer entropy methods that used it before.

    #### CausalityTools.jl

    - Added [`NearestNeighbourMITest`](@ref) transfer entropy test.
    - Added [`NormalisedPredictiveAsymmetryTest`](@ref), for which the predictive asymmetry
        is normalised to some fraction of the mean of the raw values of the predictive statistic.
    - All subtypes of `CausalityTest` are now mutable. This allows adjusting test
        parameters during sensitivity analyses without creating a new test every time.

    #### TimeseriesSurrogates.jl

    - Plotting functionality is behind a Requires-block. To use the plotting functionality, you
        now need to do `using Plots` beforehand.

    ### Documentation

    - Fixed small error in documentation for `DiscreteSystemSetup`.
    - Created high-level pages for transfer entropy. Re-did documentation.
    - Created high-level pages for cross mappings.
    """; color = :light_magenta)
    touch(joinpath(@__DIR__, update_name))
    end
    end

    # Example systems
    include("systems/Systems.jl")

    # SystemModels module
    include("SystemModels/SystemModels.jl")

    # Wrappers of the different methods.
    include("method_wrappers/highlevel_methods.jl")
    include("method_wrappers/transferentropy_from_timeseries.jl")

    # Various algorithsm that are implemented here and not in subpackages
    include("algorithms/smeasure.jl")

    # High-level estimators
    include("CausalityTests/CausalityTests.jl")

    ################################################################
    # Integration with DynamicalSystems.jl (for all causality tests)
    ################################################################
    include("integration_dynamicalsystems/IntegrationDynamicalSystems.jl")

    ################################################################
    # Integration with UncertainData.jl (for all causality tests)
    ################################################################
    include("integration_uncertaindata/IntegrationUncertainData.jl")

    function __init__()

        # Surrogate wrappers for embeddings and Datasets
        @require TimeseriesSurrogates="c804724b-8c18-5caa-8579-6025a0767c70" @eval include("surrogates/surrogates.jl")
    end

    # Plot recipes, also for all sub-packages
    #include("plot_recipes/CausalityToolsPlotRecipes.jl")

end # module


text:
function MixtureMultivariateNormalDistribution(arg0::List)
    return MixtureMultivariateNormalDistribution((List,), arg0)
end

function MixtureMultivariateNormalDistribution(arg0::RandomGenerator, arg1::List)
    return MixtureMultivariateNormalDistribution((RandomGenerator, List), arg0, arg1)
end

function MixtureMultivariateNormalDistribution(arg0::Vector{jdouble}, arg1::Vector{Vector{jdouble}}, arg2::Vector{Vector{jdouble}})
    return MixtureMultivariateNormalDistribution((Vector{jdouble}, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}), arg0, arg1, arg2)
end

function density(obj::MixtureMultivariateRealDistribution, arg0::Vector{jdouble})
    return jcall(obj, "density", jdouble, (Vector{jdouble},), arg0)
end

function get_components(obj::MixtureMultivariateRealDistribution)
    return jcall(obj, "getComponents", List, ())
end

function reseed_random_generator(obj::MixtureMultivariateRealDistribution, arg0::jlong)
    return jcall(obj, "reseedRandomGenerator", void, (jlong,), arg0)
end

function sample(obj::MixtureMultivariateRealDistribution)
    return jcall(obj, "sample", Vector{jdouble}, ())
end

ext: jl


text:
using Turing
using Distributions
using StatsPlots
using StatsBase
using Random

# Setup Turing
Turing.turnprogress(false)
Turing.setadbackend(:forward_diff)

# Binning Function
bin(data, edges) = fit(Histogram, data, edges);

# Sample Signal/Background
sample_signal(n; μ=1, σ=0) = rand(Normal(μ, σ), n)
sample_background(n; θ=1) = rand(Exponential(θ), n)

# Sample Generator
function sample_full(n; θ=1, μ=1, σ=0, r=1)
    return (sample_signal(round(Int, n * r); μ=μ, σ=σ),
            sample_background(n; θ=θ))
end

# Signal/Background Estimate Generator
function expected(n, binning; θ=1, μ=1, σ=0, r=1)
    s, b = sample_full(n; θ=θ, μ=μ, σ=σ, r=r)
    sₕ = bin(s, binning)
    bₕ = bin(b, binning)
    return sₕ, bₕ, merge(bₕ, sₕ)
end

# Fake Data Generator
function observed(α, n, binning; θ=1, μ=1, σ=0, r=1)
    s, b = sample_full(n; θ=θ, μ=μ, σ=σ, r=r)
    sₕ = bin(s, binning)
    sₕ.weights = map(w -> round(Int, α * w), sₕ.weights)
    return merge(sₕ, bin(b, binning))
end

# Fake Data Generator
function observed(α, sₕ, bₕ)
    dsₕ = deepcopy(sₕ)
    dsₕ.weights = map(w -> round(Int, α * w), dsₕ.weights)
    return merge(dsₕ, bₕ)
end

# True Signal Strength
true_α = 0.25

# Make Expected Counts and Fake Data
Random.seed!(1)
N = 4000
binning = 0:0.1:2.5
θ = 0.56
μ = 1.12
σ = 0.09
r = 0.1
signal, background, total = expected(N, binning; θ=θ, μ=μ, σ=σ, r=r)
data = observed(true_α, N, binning; θ=θ, μ=μ, σ=σ, r=r)

# "Bump Hunt" Plot
plot(total, title="Bump Hunt", label="expected signal")
plot!(background, label="expected background")
scatter!(data, label="data", color=:black)

# Counting Experiment Model
@model counting_experiment(data, signal, background) = begin
    # Signal Strength Prior
    α ~ Flat()

    # Draw a Poisson Process in each Bin
    for i in eachindex(data)
        data[i] ~ Poisson(α * signal[i] + background[i])
    end
end

iterations = 1000

# HMC:
ϵ = 0.001
τ = 1
hmc_sampler = HMC(iterations, ϵ, τ);

# MH:
mh_sampler = MH(iterations);

# NUTS:
adapts = 200
δ = 0.65
nuts_sampler = NUTS(iterations, adapts, δ)

# Sample Many Markov Chains
sample_many(n, args...) = mapreduce(_ -> sample(args...), chainscat, 1:n)

# Compute the Posterior
chains = 12
algorithm = mh_sampler
model_instance = counting_experiment(data.weights,
                                     signal.weights,
                                     background.weights)
observed_posterior = sample_many(chains, model_instance, algorithm)

println("true_α = $true_α")
println("median(α) = $(median(observed_posterior[:α].value))")
histogram(observed_posterior[:α], bins=70)

# 95% Confidence Limit for Signal Strength
limit(posterior; q=0.95) = quantile(posterior, q=[q])[1, 2]

# Observed Limit
observed_limit = limit(observed_posterior[:α])

chains = round(Int, chains / 3)
limit_samples = 100
expected_limits = Vector{Float64}(undef, limit_samples)
for i in 1:limit_samples
    sb = bin(sample_background(N; θ=θ), binning)
    toy_model = counting_experiment(sb.weights,
                                    signal.weights,
                                    background.weights)
    expected_limits[i] = limit(sample_many(chains, toy_model, algorithm))
end

histogram(expected_limits, bins=50)

expected_bands = quantile(expected_limits, [0.025, 0.16, 0.5, 0.84, 0.975])

plus_2σ = expected_bands[5]
excess_found = plus_2σ < observed_limit
println(excess_found ? "excess:" : "upper limit:")
println("  +2σ=", plus_2σ, " $(excess_found ? '<' : '>') obs=", observed_limit)

plus_5σ = quantile(expected_limits, 0.999999426697)

@model counting_experiment(data, control) = begin
    N = length(data)
  	signal = Vector(undef, N)
    background = Vector(undef, N)

 	  # Signal Strength Prior
    α ~ Flat()

    # Background in Control Region
    τ ~ Flat()

    # Fit Signal to Normal(μ, σ²)
    σ² ~ InverseGamma(2, 3)
  	μ ~ Normal(0, √σ²)

  	# Fit Background to Exponential(θ)
    θ ~ Flat()

    # Draw a Poisson Process in each Bin
    for i in 1:N
        # Fit Signal and Background
      	signal[i] ~ Normal(μ, σ²)
      	background[i] ~ Exponential(θ)
        # Estimate Expected Counts
        data[i] ~ Poisson(α * signal[i] + background[i])
    	  control[i] ~ Poisson(τ * background[i])
    end
end

max_stars_repo_name: JuliaTagBot/notebooks

max_issues_repo_name: JuliaTagBot/notebooks

max_forks_repo_name: JuliaTagBot/notebooks


text:
# Tests for Profit DEA Models

# Struct for testing DEA model without monetary option
struct wrongProfitDEAmodel <: AbstractProfitDEAModel
    n::Int64
    eff::Vector
end

@testset "ProfitDEAModel" begin

    ## Test Profit DEA Model
    # Test with Zofio, Pastor and Aparicio (2013) data
    X = [1 1; 1 1; 0.75 1.5; 0.5 2; 0.5 2; 2 2; 2.75 3.5; 1.375 1.75]
    Y = [1 11; 5 3; 5 5; 2 9; 4 5; 4 2; 3 3; 4.5 3.5]
    P = [2 1; 2 1; 2 1; 2 1; 2 1; 2 1; 2 1; 2 1]
    W = [2 1; 2 1; 2 1; 2 1; 2 1; 2 1; 2 1; 2 1]

    GxGydollar = 1 ./ (sum(P, dims = 2) + sum(W, dims = 2))
    GxGydollar = repeat(GxGydollar, 1, 2)

    deaprofitdollar = deaprofit(X, Y, W, P, Gx = GxGydollar, Gy = GxGydollar)

    @test typeof(deaprofitdollar) == ProfitDEAModel

    @test nobs(deaprofitdollar) == 8
    @test ninputs(deaprofitdollar) == 2
    @test noutputs(deaprofitdollar) == 2
    @test efficiency(deaprofitdollar, :Economic)   ≈ [2; 2; 0; 2; 2; 8; 12; 4] atol = 1e-3
    @test efficiency(deaprofitdollar, :Technical)  ≈ [0; 0; 0; 0; 0; 6; 12; 3] atol = 1e-3
    @test efficiency(deaprofitdollar, :Allocative) ≈ [2; 2; 0; 2; 2; 2; 0; 1] atol = 1e-3

    @test efficiency(deaprofit(X, Y, W, P, Gx = :Monetary, Gy = :Monetary)) == efficiency(deaprofitdollar)

    @test efficiency(deaprofit(targets(deaprofitdollar, :X), targets(deaprofitdollar, :Y), W, P, Gx = :Monetary, Gy = :Monetary)) ≈ zeros(8)

    @test normfactor(deaprofitdollar) == ones(8);

    @test peersmatrix(deaprofitdollar) == deaprofitdollar.lambda

    # Check directions checking technical efficiency
    @test efficiency(deaprofit(X, Y, W, P, Gx = :Zeros, Gy = :Ones), :Technical) == efficiency(deaddf(X, Y, Gx = :Zeros, Gy = :Ones, rts = :VRS))
    @test efficiency(deaprofit(X, Y, W, P, Gx = :Ones, Gy = :Zeros), :Technical) == efficiency(deaddf(X, Y, Gx = :Ones, Gy = :Zeros, rts = :VRS))
    @test efficiency(deaprofit(X, Y, W, P, Gx = :Observed, Gy = :Observed), :Technical) == efficiency(deaddf(X, Y, Gx = :Observed, Gy = :Observed, rts = :VRS))
    @test efficiency(deaprofit(X, Y, W, P, Gx = :Mean, Gy = :Mean), :Technical) == efficiency(deaddf(X, Y, Gx = :Mean, Gy = :Mean, rts = :VRS))

    # Print
    show(IOBuffer(), deaprofitdollar)

    # Check normalization factor with different direction
    @test normfactor(deaprofit(X, Y, W, P, Gx = :Ones, Gy = :Ones)) == vec(sum(P .* ones(size(Y)), dims = 2) .+ sum(W .* ones(size(X)), dims = 2))

    # Check monetary option
    @test ismonetary(deaprofitdollar) == false

    deaprofitmonetary = deaprofit(X, Y, W, P, Gx = :Ones, Gy = :Ones, monetary = true)
    @test efficiency(deaprofitmonetary, :Economic)   ≈ [2; 2; 0; 2; 2; 8; 12; 4] atol = 1e-3
    @test efficiency(deaprofitmonetary, :Technical)  ≈ [0; 0; 0; 0; 0; 6; 12; 3] atol = 1e-3
    @test efficiency(deaprofitmonetary, :Allocative) ≈ [2; 2; 0; 2; 2; 2; 0; 1] atol = 1e-3
    @test ismonetary(deaprofitmonetary) == true

    # Test errors
    @test_throws DimensionMismatch deaprofit([1; 2 ; 3], [4 ; 5], [1; 1; 1], [4; 5], Gx = [1; 2 ; 3], Gy = [4 ; 5]) #  Different number of observations
    @test_throws DimensionMismatch deaprofit([1; 2; 3], [4; 5; 6], [1; 2; 3; 4], [4; 5; 6], Gx = [1; 2; 3], Gy = [4; 5; 6]) # Different number of observation in input prices
    @test_throws DimensionMismatch deaprofit([1; 2; 3], [4; 5; 6], [1; 2; 3], [4; 5; 6; 7], Gx = [1; 2; 3], Gy = [4; 5; 6]) # Different number of observation in output prices
    @test_throws DimensionMismatch deaprofit([1 1; 2 2; 3 3], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3], [4; 5; 6], Gx = [1 1; 2 2; 3 3], Gy = [4; 5; 6]) # Different number of input prices and inputs
    @test_throws DimensionMismatch deaprofit([1; 2; 3], [4 4; 5 5; 6 6], [1; 2; 3], [4 4 4; 5 5 5; 6 6 6], Gx = [1; 2; 3], Gy = [4 4; 5 5; 6 6]) # Different number of oputput prices and outputs
    @test_throws DimensionMismatch deaprofit([1 1; 2 2; 3 3], [4; 5; 6], [1 1; 2 2; 3 3], [4; 5; 6], Gx = [1 1 1; 2 2 2; 3 3 3], Gy = [4; 5; 6]) # Different size of inputs direction
    @test_throws DimensionMismatch deaprofit([1; 2; 3], [4 4; 5 5; 6 6], [1; 2; 3], [4 4; 5 5; 6 6], Gx = [1; 2; 3], Gy = [4 4 4; 5 5 5; 6 6 6]) # Different size of outputs direction
    @test_throws ArgumentError deaprofit([1; 2; 3], [1; 2; 3], [1; 1; 1], [1; 1; 1], Gx = :Error, Gy = :Ones) # Invalid inputs direction
    @test_throws ArgumentError deaprofit([1; 2; 3], [1; 2; 3], [1; 1; 1], [1; 1; 1], Gx = :Ones, Gy = :Error) # Invalid outputs direction
    @test_throws ArgumentError targets(deaprofitdollar, :error) # Invalid target

    # Test struct defaults and errors
    @test_throws ArgumentError ismonetary(wrongProfitDEAmodel(3, [1; 2; 3])) # Model does not have info on peers

    # ------------------
    # Test Vector and Matrix inputs and outputs
    # ------------------
    # Tests against results in R

    # Inputs is Matrix, Outputs is Vector
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6	8]
    Y = [1; 1; 1; 1; 1; 1; 1; 1]
    W = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1]
    P = [1; 1; 1; 1; 1; 1; 1; 1]

    @test efficiency(deaprofit(X, Y, W, P, Gx = X, Gy = Y)) ≈ [0; 0.1666666667; 0.1666666667; 0.375; 0.5454545455; 0.375; 0.375; 0.5283018868]

    # Inputs is Vector, Output is Matrix
    X = [1; 1; 1; 1; 1; 1; 1; 1]
    Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5]
    W = [1; 1; 1; 1; 1; 1; 1; 1]
    P = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1]

    @test efficiency(deaprofit(X, Y, W, P, Gx = X, Gy = Y)) ≈ [0; 0.1538461538; 0.1538461538; 0.6666666667; 1.142857143; 0.3636363636; 0.3636363636; 1]

    # Inputs is Vector, Output is Vector
    X = [2; 4; 8; 12; 6; 14; 14; 9.412]
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353]
    W = [1; 1; 1; 1; 1; 1; 1; 1]
    P = [1; 1; 1; 1; 1; 1; 1; 1]

    @test efficiency(deaprofit(X, Y, W, P, Gx = X, Gy = Y)) ≈ [0.6666666667; 0; 0.0625; 0.1904761905; 0.4444444444; 0.3809523810; 0.2608695652; 0.6849978751]

end


text:
# ==================================================================================================== #
#                                                                                #
#                                          Helper Functions                                            #
#                                                                                #
# ==================================================================================================== #
function readConfig(file)
    @info "Reading Inputs..."

    tml = TOML.parsefile(file)
    sim = tml["Simulation"]
    χfill = nothing
    rr = r"^fixed:(?P<start>\N+):(?P<stop>\N+)"

    if tml["Simulation"]["chi_unusable_fill_value"] == "0"
        χfill = zero_χ_fill
    elseif tml["Simulation"]["chi_unusable_fill_value"] == "chi_lambda"
        χfill = lambda_χ_fill
    elseif tml["Simulation"]["chi_unusable_fill_value"] == "chi"
        χfill = χ_fill
    else
        error("could not parse chi fill value")
    end
    tc_type_f = Symbol(lowercase(tml["Simulation"]["fermionic_tail_correction"]))
    if !(tc_type_f in [:nothing, :richardson, :shanks])
        error("Unrecognized tail correction type \"$(tc_type_f)\"")
    end
    tc_type_b = Symbol(lowercase(tml["Simulation"]["bosonic_tail_correction"]))
    if !(tc_type_b in [:nothing, :richardson, :shanks, :coeffs])
        error("Unrecognized tail correction type \"$(tc_type_b)\"")
    end
    λc_type = Symbol(lowercase(tml["Simulation"]["lambda_correction"]))
    if !(λc_type in [:nothing, :sp, :sp_ch])
        error("Unrecognized tail correction type \"$(λc_type)\"")
    end
    ωsum_inp = lowercase(tml["Simulation"]["bosonic_sum_range"])
    m = match(rr, ωsum_inp)
    ωsum_type = if m !== nothing
        tuple(parse.(Int, m.captures)...)
    elseif ωsum_inp in ["common", "individual", "full"]
        Symbol(ωsum_inp)
    else
        error("Could not parse bosonic_sum. Allowed input is: common, individual, full, fixed:N:M")
    end
    smoothing = Symbol(lowercase(tml["Simulation"]["omega_smoothing"]))
    if !(smoothing in [:nothing, :range, :full])
        error("Unrecognized smoothing type \"$(smoothing)\"")
    end

    λrhs_type = Symbol(lowercase(tml["Simulation"]["rhs"]))
    !(λrhs_type in [:native, :fixed, :error_comp]) && error("Could not parse rhs type for lambda correction. Options are native, fixed, error_comp.")

    dbg_full_eom_omega = (haskey(tml["Debug"], "full_EoM_omega") && tml["Debug"]["full_EoM_omega"]) ? true : false

    env = EnvironmentVars(   tml["Environment"]["inputDataType"],
                             tml["Environment"]["writeFortran"],
                             tml["Environment"]["loadAsymptotics"],
                             tml["Environment"]["inputDir"],
                             tml["Environment"]["freqFile"],
                             tml["Environment"]["inputVars"],
                             tml["Environment"]["asymptVars"],
                             tml["Environment"]["cast_to_real"],
                             String([(i == 1) ? uppercase(c) : lowercase(c)
                                     for (i, c) in enumerate(tml["Environment"]["loglevel"])]),
                             lowercase(tml["Environment"]["logfile"]),
                             tml["Environment"]["progressbar"]
                            )
    mP = ModelParameters(tml["Model"]["U"],
                            tml["Model"]["mu"],
                            tml["Model"]["beta"],
                            tml["Model"]["nden"])

    freqFilePath = ""
    if !isfile(env.freqFile)
        m = match(r"b(?<bf>\d+)f(?<ff>\d+)",env.freqFile)
        nBose = parse(Int, m[:bf])
        nFermi = parse(Int, m[:ff])
        freqFilePath = "$(@__DIR__)/../scripts/freqList.jld2"
        gen_mesh(nBose, nFermi, 0, freqFilePath)
    else
        freqFilePath = env.freqFile
    end
    JLD2.@load freqFilePath freqRed_map freqList freqList_min parents ops nFermi nBose shift base offset

    sP = SimulationParameters(nBose,nFermi,shift,
                               tc_type_f,
                               tc_type_b,
                               λc_type,
                               ωsum_type,
                               λrhs_type,
                               tml["Simulation"]["force_full_bosonic_chi"],
                               χfill,
                               tml["Simulation"]["bosonic_tail_coeffs"],
                               tml["Simulation"]["fermionic_tail_coeffs"],
                               tml["Simulation"]["usable_prct_reduction"],
                               smoothing,
                               dbg_full_eom_omega
    )
    kGrids = []
    for Nk in tml["Simulation"]["Nk"]
        kGrid = gen_kGrid(tml["Model"]["kGrid"], Nk)
        push!(kGrids, kGrid)
    end
    kGridLoc = gen_kGrid(tml["Model"]["kGrid"], 1)
    return mP, sP, env, kGrids, kGridLoc, freqRed_map, freqList, freqList_min, parents, ops, nFermi, nBose, shift, base, offset
end

function convertGF!(GF, storedInverse, storeFull)
    if storedInverse
        GF = 1 ./ GF
    end
    if storeFull
        NH = size(GF, 1)
        tmp = copy(GF)
        GF = Array{eltype(GF)}(undef, 2*NH - 1)
        GF[1:(NH-1)] = reverse(conj.(tmp[2:NH]))
        GF[NH:end] = tmp[1:NH]
    end
end

function print_chi_bubble(qList, res, simParams)
    for j in 1:size(res,1)
        print(" ========== ω = $(j-(simParams.n_iω + 1)) =============== \n")
        for k in 1:size(res,2)
            print(" ---------- ν = $(k-1) -------------- \n")
            for (qi,q) in enumerate(qList)
                @printf("   q = (%.2f,%.2f): %.2f + %.2fi\n", q[1],q[2], real(res[j, k, qi]), imag(res[j, k, qi]))
            end
        end
    end
end

function reduce_3Freq(inp, freqBox, simParams)
    n_iω = simParams.n_iω
    n_iν = simParams.n_iν
    start_ωₙ =  -(freqBox[1,1] + n_iω) + 1
    start_νₙ =  -(freqBox[2,1] + n_iν) + 2

    if (length(start_ωₙ:(start_ωₙ+2*n_iω)) != size(inp,1)) ||
       (length(start_νₙ:(start_νₙ + 2*n_iν - 1)) != size(inp,2))
        @error "Input frequency range does not match input data. This can lead to wrong results! Input size: ", size(inp),
        " got ", length(start_ωₙ:(start_ωₙ+2*n_iω)), ",", length(start_νₙ:(start_νₙ + 2*n_iν - 1))
    end
    # This selects the range of [-ωₙ, ωₙ] and [-νₙ, νₙ-1] from gamma
    return inp[start_ωₙ:(start_ωₙ+2*n_iω), start_νₙ:(start_νₙ + 2*n_iν - 1), start_νₙ:(start_νₙ +2*n_iν - 1)]
end

reduce_3Freq!(inp, freqBox, simParams) = inp = reduce_3Freq(inp, freqBox, simParams)

function convert_from_fortran(simParams, env, loadFromBak=false)
    @info "Reading Fortran Input, this can take several minutes."
        g0 = readFortranSymmGF(simParams.n_iν+2*simParams.n_iω, env.inputDir*"/g0mand", storedInverse=true)
        gImp = readFortranSymmGF(simParams.n_iν+2*simParams.n_iω, env.inputDir*"/gm_wim", storedInverse=false)
        freqBox, Γch, Γsp = readFortranΓ(env.inputDir*"/gamma_dir")
        @info "Done Reading Gamma"
        χDMFTch, χDMFTsp = readFortranχDMFT(env.inputDir*"/chi_dir")
        @info "Done Reading chi"
    Γch = -1.0 .* reduce_3Freq(Γch, freqBox, simParams)
    Γsp = -1.0 .* reduce_3Freq(Γsp, freqBox, simParams)
    χDMFTch = reduce_3Freq(χDMFTch, freqBox, simParams)
    χDMFTsp = reduce_3Freq(χDMFTsp, freqBox, simParams)
    if env.writeFortran
        @info "Writing HDF5 (vars.jdl) and Fortran (fortran_out/) output."
        writeFortranΓ("fortran_out", "gamma", simParams, Γch, Γsp)
        writeFortranΓ("fortran_out", "chi", simParams, 0.5 .* (χDMFTch .+ χDMFTsp), 0.5 .* (χDMFTch .- χDMFTsp))
    end
    if length(env.inputVars) > 0
        JLD2.@save env.inputVars Γch Γsp χDMFTch χDMFTsp gImp g0
    end
end

# ==================================================================================================== #
#                                                                                #
#                                             JLD2 Input                                               #
#                                                                                #
# ==================================================================================================== #
function readEDAsymptotics_julia(env)
    @error("Asymptotics from jld2 not implemented yet")
end

# ==================================================================================================== #
#                                                                                #
#                                            Parquet Input                                             #
#                                                                                #
# ==================================================================================================== #
function readFortranSymmGF_pq(nFreq, filename; storedInverse, storeFull=false)
    GF = redirect_stdout(open("/dev/null","w")) do
        load(filename) |> @orderby(_.__index_level_0__) |> @take(nFreq)|> @map((_.Re) + (_.Im)*1im) |> collect
    end
    convertGF!(GF, storedInverse, storeFull)
    println(GF)
    return GF
end

function readFortranΓ_pq(fileName::String)
    #in = redirect_stdout(open("/dev/null","w")) do
        #load(filename) |> @orderby(_.__index_level_0__, _.__index_level_1__)
    #end
    #ωₙ, freqBox, Γcharge0, Γspin0
    Γcharge = Array{Complex{Float64}}(undef, length(files), size(Γspin0,1), size(Γspin0,2))
    Γspin   = Array{Complex{Float64}}(undef, length(files), size(Γspin0,1), size(Γspin0,2))
    Γcharge[1,:,:] = Γcharge0
    Γspin[1,:,:]   = Γspin0
    ω_min = ωₙ
    ω_max = ωₙ

    for (i,file) in enumerate(files[2:end])
        ωₙ, _, Γcharge_new, Γspin_new = readFortran3FreqFile(dirName * "/" * file, sign=-1.0)
        ω_min = if ωₙ < ω_min ωₙ else ω_min end
        ω_max = if ωₙ > ω_max ωₙ else ω_max end
        Γcharge[i,:,:] = Γcharge_new
        Γspin[i,:,:] = Γspin_new
    end
    freqBox = [ω_min ω_max; freqBox[1,1] freqBox[1,2]; freqBox[2,1] freqBox[2,2]]
    #Γcharge = permutedims(Γcharge, [3,1,2])
    #Γspin   = permutedims(Γspin, [3,1,2])
    return freqBox, Γcharge, Γspin
end

function readFortranχDMFT_pq(filename::String)
    in = redirect_stdout(open("/dev/null","w")) do
        load(filename) |> collect
    end
    #_, _, χup, χdown = readFortran3FreqFile(dirName * "/" * files[1], sign=1.0, freqInteger=false)

    for file in files[2:end]
        _, _, χup_new, χdown_new = readFortran3FreqFile(dirName * "/" * file, sign=1.0, freqInteger=false)
        χup = cat(χup, χup_new, dims=3)
        χdown = cat(χdown, χdown_new, dims=3)
    end
    χup = permutedims(χup, [3,1,2])
    χdown   = permutedims(χdown, [3,1,2])
    χCharge = χup .+ χdown
    χSpin   = χup .- χdown
    return χCharge, χSpin
end

function convert_from_fortran_pq(simParams, env)
    println("Reading Fortran Parquet-File Input, this can take several minutes.")
    g0 = readFortranSymmGF_pq(simParams.n_iν+simParams.n_iω, env.inputDir*"/g0mand.parquet", storedInverse=true)
    gImp = readFortranSymmGF_pq(simParams.n_iν+simParams.n_iω, env.inputDir*"/gm_wim.parquet", storedInverse=false)
    χDMFTCharge, χDMFTSpin = readFortranχDMFT_pq(env.inputDir*"/vert_chi.parquet")
    freqBox, Γcharge, Γspin = readFortranΓ_pq(env.inputDir*"/GAMMA_DM_FULLRANGE.parquet")
    Γcharge = -1.0 .* reduce_3Freq(Γcharge, freqBox, simParams)
    Γspin = -1.0 .* reduce_3Freq(Γspin, freqBox, simParams)
    χDMFTCharge = reduce_3Freq(χDMFTCharge, freqBox, simParams)
    χDMFTSpin = reduce_3Freq(χDMFTSpin, freqBox, simParams)
    println("Writing HDF5 (vars.jdl) and Fortran (fortran_out/) output.")
    writeFortranΓ("fortran_out", "gamma", simParams, Γcharge, Γspin)
    writeFortranΓ("fortran_out", "chi", simParams, 0.5 .* (χDMFTCharge .+ χDMFTSpin), 0.5 .* (χDMFTCharge .- χDMFTSpin))
    JLD2.@save env.inputVars Γch Γsp χDMFTch χDMFTsp gImp g0
end

# ==================================================================================================== #
#                                                                                #
#                                            Text Input                                                #
#                                                                                #
# ==================================================================================================== #
function readFortranqωFile(filename, nDims; readq = false, data_cols = 2)
    InString = open(filename, "r") do f
        readlines(f)
    end
    start = if readq 1 else nDims+1  end
    InArr = parse.(Float64,hcat(split.(InString)...)[start:end,:])
    if readq
        qVecArr = InArr[1:3,:]
    else
        qVecArr = []
        nDims = 0
    end
    data_λ = InArr[nDims+1,:] + InArr[nDims+2,:] .* 1im
    data = InArr[nDims+3,:] + InArr[nDims+4,:] .* 1im
    return qVecArr, data_λ, data
end

function readFortranΓ(dirName::String)
    files = readdir(dirName)
    ωₙ, freqBox, Γcharge0, Γspin0 = readFortran3FreqFile(dirName * "/" * files[1], sign=-1.0)
    Γcharge = Array{Complex{Float64}}(undef, length(files), size(Γspin0,1), size(Γspin0,2))
    Γspin   = Array{Complex{Float64}}(undef, length(files), size(Γspin0,1), size(Γspin0,2))
    Γcharge[1,:,:] = Γcharge0
    Γspin[1,:,:]   = Γspin0
    ω_min = ωₙ
    ω_max = ωₙ

    for (i,file) in enumerate(files[2:end])
        ωₙ, _, Γcharge_new, Γspin_new = readFortran3FreqFile(dirName * "/" * file, sign=-1.0)
        ω_min = if ωₙ < ω_min ωₙ else ω_min end
        ω_max = if ωₙ > ω_max ωₙ else ω_max end
        Γcharge[i+1,:,:] = Γcharge_new
        Γspin[i+1,:,:] = Γspin_new
    end
    freqBox = [ω_min ω_max; freqBox[1,1] freqBox[1,2]; freqBox[2,1] freqBox[2,2]]
    Γcharge = permutedims(Γcharge, [1,3,2])
    Γspin   = permutedims(Γspin, [1,3,2])
    return freqBox, Γcharge, Γspin
end

"""
    Returns χ_DMFT[ω, ν, ν']
"""
function readFortranχDMFT(dirName::String)
    files = readdir(dirName)
    _, _, χupup, χupdo = readFortran3FreqFile(dirName * "/" * files[1], sign=1.0, freqInteger=false)
    for file in files[2:end]
        _, _, χup_new, χupdo_new = readFortran3FreqFile(dirName * "/" * file, sign=1.0, freqInteger=false)
        χupup = cat(χupup, χup_new, dims=3)
        χupdo = cat(χupdo, χupdo_new, dims=3)
    end
    χupup = permutedims(χupup, [3,2,1])
    χupdo   = permutedims(χupdo, [3,2,1])
    χCharge = χupup .+ χupdo
    χSpin   = χupup .- χupdo
    return χCharge, χSpin
end

function readFortranχlDGA(dirName::String, nDims)
    files = readdir(dirName)
    qVecs, data_λ_i, data_i = readFortranqωFile(dirName * "/" * files[1], nDims, readq = true);
    data = Array{Complex{Float64}}(undef, length(files), 2, length(data_i))
    data[1,1,:] = data_λ_i
    data[1,2,:] = data_i
    for (i,file) in enumerate(files[2:end])
        _, data_λ_i, data_i = readFortranqωFile(dirName * "/" * file, nDims, readq = false);
        data[i+1,1,:] = data_λ_i
        data[i+1,2,:] = data_i
    end
    qVecs, data
end

function readFortranBubble(dirName::String, nBose, nFermi, nQ)
    bubble = Array{Complex{Float64}}(undef, nBose, nQ, nFermi)
    files_i = readdir(dirName)
    files = [dirName * "/" * f for f  in files_i]
    iBose = 1
    for file in files[1:end]
        InString = open(file, "r") do f
            readlines(f)
        end
        InArr = parse.(Float64,hcat(split.(InString)...)[4:end,:])
        bubble[iBose, :, :] = reshape(InArr[1,:] .+ InArr[2,:].*1im, (nQ, nFermi))
        iBose += 1
    end

    bubble = permutedims(bubble, [1,3,2])
    return bubble
end

function read_anderson_parameters(file)
    content = open(file) do f
        readlines(f)
    end

    in_epsk = false
    in_tpar = false
    ϵₖ = []
    Vₖ = []
    μ = 0
    for line in content
        if "Eps(k)" == strip(line)
            in_epsk = true
            continue
        elseif "tpar(k)" == strip(line)
            in_epsk = false
            in_tpar = true
            continue
        end

        if in_epsk
            push!(ϵₖ, parse(Float64, line))
        elseif in_tpar
            # skip last line, which is mu
            if length(Vₖ) < length(ϵₖ)
                push!(Vₖ, parse(Float64, line))
            else
                if occursin("#", line)
                    μ = parse(Float64, line[1:(findfirst("#", line))[1] - 1])
                else
                    μ = parse(Float64, line)
                end
            end
        end
    end
    return convert(Array{Float64,1}, ϵₖ), convert(Array{Float64,1}, Vₖ), μ
end

function readGImp(filename; only_positive=false)
    GFString = open(filename, "r") do f
        readlines(f)
    end

    tmp = parse.(Float64,hcat(split.(GFString)...)) # Construct a 2xN array of floats (re,im as 1st index)
    tmpG = tmp[2,:] .+ tmp[3,:].*1im
    tmpiνₙ = tmp[1,:] .* 1im
    if only_positive
        GImp = tmpG
        iνₙ  = tmpiνₙ
    else
        N = 2*size(tmpG,1)
        NH = size(tmpG,1)
        GImp = zeros(Complex{Float64}, N)
        iνₙ  = zeros(Complex{Float64}, N)
        GImp[1:(NH)] = reverse(conj.(tmpG[1:NH]))
        GImp[(NH+1):N] = tmpG[1:NH]
        iνₙ[1:(NH)] = conj.(reverse(tmpiνₙ[1:(NH)]))
        iνₙ[(NH+1):N] = tmpiνₙ[1:NH]
    end
    return iνₙ, GImp
end

function readEDAsymptotics(env, mP)
    χ_asympt = readdlm(env.inputDir * "/chi_asympt")
    χchAsympt = (χ_asympt[:,2] + χ_asympt[:,4]) / (2*mP.β*mP.β);
    χspAsympt = (χ_asympt[:,2] - χ_asympt[:,4]) / (2*mP.β*mP.β);
    #= _, χup, χdo = readFortranEDχ(env.inputDir * "/chi_dir", freqInteger = false) =#
    #= χchED = χup .+ χdo =#
    #= χspED = χup .- χdo =#
    save(env.asymptVars, "chi_ch_asympt", χchAsympt, "chi_sp_asympt", χspAsympt,
         compress=true, compatible=true)
end

function readFortranSymmGF(nFreq, filename; storedInverse, storeFull=false)
    GFString = open(filename, "r") do f
        readlines(f)
    end

    if size(GFString, 1)*(1 + 1*storeFull) < nFreq
        throw(BoundsError("nFermFreq in simulation parameters too large!"))
    end

    tmp = parse.(Float64,hcat(split.(GFString)...)[2:end,:]) # Construct a 2xN array of floats (re,im as 1st index)
    tmp = tmp[1,:] .+ tmp[2,:].*1im

    if storedInverse
        tmp = 1 ./ tmp
    end

    GF = Array{Complex{Float64}}(undef, nFreq)
    if storeFull
        NH = Int(nFreq/2)
        GF[1:(NH-1)] = reverse(conj.(tmp[2:NH]))
        GF[NH:nFreq] = tmp[1:NH]
    else
        GF = tmp[1:nFreq]
    end
    return GF
end

function readFortran3FreqFile(filename; sign = 1.0, freqInteger = true)
    InString = open(filename, "r") do f
        readlines(f)
    end
    NCols = 7
    NRows = length(InString)
    lineLen = length(InString[1])
    splitLength = floor(Int64,lineLen/NCols)
    if lineLen%NCols != 0
        @warn "   ---> Warning!! Could not find fixed column width!"
    end
    #InArr = sign .* parse.(Float64,hcat(split.(InString)...)[2:end,:])
    InArr = zeros(Float64, NRows, NCols-1)
    N = Int(sqrt(NRows))
    ErrorLine = []
    for i in 1:NRows
        row =  split_n(InString[i], splitLength, lineLen)[2:end]
        for j in 1:(NCols-1)
            InArr[i,j] = try
                sign * parse(Float64,row[j])
            catch e
                NaN
            end
        end
    end
    if freqInteger
        ωₙ = Int(parse(Float64, split(InString[1])[1]))
        tmpArr1 = filter(!isnan,InArr[:,1])
        tmpArr2 = filter(!isnan,InArr[:,2])
        middleFreqBox = Int64.([minimum(tmpArr1),maximum(tmpArr1)])
        innerFreqBox  = Int64.([minimum(tmpArr2),maximum(tmpArr2)])
        freqBox = permutedims([middleFreqBox innerFreqBox], [2,1])
    else
        #print("\rWarning, non integer frequencies in "*filename*" ignored!")
        ωₙ = 0
        freqBox = []
    end

    InCol1  = reshape(InArr[:,3] .+ InArr[:,4].*1im, (N, N))
    InCol2  = reshape(InArr[:,5] .+ InArr[:,6].*1im, (N, N))
    return ωₙ, freqBox, InCol1, InCol2
end

# ==================================================================================================== #
#                                                                                #
#                                            Text Output                                                #
#                                                                                #
# ==================================================================================================== #
function writeFortranΓ(dirName::String, fileName::String, simParams, inCol1, inCol2)
    simParams.n_iν+simParams.n_iω
    if !isdir(dirName)
        mkdir(dirName)
    end
    for (ωi,ωₙ) in enumerate(-simParams.n_iω:simParams.n_iω)
        filename = dirName * "/" * fileName * lpad(ωi-1,3,"0")
        open(filename, "w") do f
            for (νi,νₙ) in enumerate(-simParams.n_iν:(simParams.n_iν-1))
                for (ν2i,ν2ₙ) in enumerate(-simParams.n_iν:(simParams.n_iν-1))
                    @printf(f, "  %18.10f  %18.10f  %18.10f  %18.10f  %18.10f  %18.10f  %18.10f\n", ωₙ, νₙ, ν2ₙ,
                            real(inCol1[ωi, νi, ν2i]), imag(inCol1[ωi, νi, ν2i]),
                            real(inCol2[ωi, νi, ν2i]), imag(inCol2[ωi, νi, ν2i]))
                end
            end
        end
    end
end

function writeFortranΣ(dirName::String, Σ_ladder)
    res = zeros(size(Σ_ladder,1), 3)
    if !isdir(dirName)
        mkdir(dirName)
    end

    for ki in 1:size(Σ_ladder,2)
        fn = dirName * "/SELF_Q_" * lpad(ki,6,"0") * ".dat"
        open(fn, write=true) do f
            write(f, "header...\n")
            res[:,1] = (2 .*(0:size(Σ_ladder,1)-1) .+ 1) .* π ./ mP.β
            res[:,2] = real.(Σ_ladder[:,ki])
            res[:,3] = imag.(Σ_ladder[:,ki])
            writedlm(f,  rpad.(round.(res; digits=14), 22, " "), "\t")
        end
    end
end

function writeFortranχ(dirName::String, χ, χ_λ, qGrid, usable_ω)
    if !isdir(dirName)
        mkdir(dirName)
    end
    for ωi in 1:size(χ, 1)
        fn = dirName * "/chi" * lpad(ωi-1,3,"0")
        open(fn, write=true) do f
            for qi in 1:size(χ, 2)
                res = ωi in usable_ω ? χ_λ[ωi - first(usable_ω) + 1,1] : 0.0
                if length(qGrid[1]) == 3
                    @printf(f, "  %18.10f  %18.10f  %18.10f  %18.10f  %18.10f  %18.10f  %18.10f\n",
                        qGrid[qi][1], qGrid[qi][2], qGrid[qi][3],
                        real(χ[ωi, qi]), imag(χ[ωi, qi]),
                        real(res), imag(res))
                else
                    @printf(f, "  %18.10f  %18.10f  %18.10f  %18.10f  %18.10f  %18.10f\n",
                        qGrid[qi][1], qGrid[qi][2],
                        real(χ[ωi, qi]), imag(χ[ωi, qi]),
                        real(res), imag(res))
                end
            end
        end
    end
end

function writeFortranEnergies(E_Kin, E_Pot, β, dirName::String)
    if !isdir(dirName)
        mkdir(dirName)
    end
    @assert size(E_Kin) == size(E_Pot)
    νGrid = 0:(length(E_Kin)-1)
    iν_n = iν_array(β, νGrid)

    fn_DGA = dirName * "/energiesDGA.dat"
    open(fn_DGA, "w") do f
        for i in 1:length(E_Kin)
            @printf(f, "  %18.10f  %18.10f  %18.10f  %18.10f  %18.10f  %18.10f  %18.10f  %18.10f  %18.10f\n",
                    imag(iν_n[i]), 0.0, 0.0, real(E_Kin[i]), 0.0, 0.0, 0.0, real(E_Pot[i]), 0.0)
        end
    end
end

function writeFortranEnergies(E_Kin, E_Pot, E_Kin_ED, E_Pot_ED, β, dirName::String)
    if !isdir(dirName)
        mkdir(dirName)
    end
    @assert size(E_Kin) == size(E_Pot) == size(E_Kin_ED) == size(E_Pot_ED)
    νGrid = 0:(length(E_Kin)-1)
    iν_n = iν_array(β, νGrid)

    writeFortranEnergies(E_Kin, E_Pot, β, dirName)

    fn_DMFT = dirName * "/energiesDMFT.dat"
    open(fn_DMFT, "w") do f
        for i in 1:length(E_Kin)
            @printf(f, "  %18.10f  %18.10f  %18.10f  %18.10f  %18.10f  %18.10f  %18.10f  %18.10f  %18.10f\n",
                    imag(iν_n[i]), 0.0, 0.0, real(E_Kin_ED[i]), 0.0, 0.0, 0.0, real(E_Pot_ED[i]), 0.0)
        end
    end
end


text:
using SymPy
using Test

@testset "Latexify" begin
  @vars α
  @test_nowarn SymPy.Latexify.latexify(- 3/(8*α) + 1/(8*α^2))
  @test  occursin("\\cdot", SymPy.Latexify.latexify(- 3/(8*α) + 1/(8*α^2), cdot=true))
  @test !occursin("\\cdot", SymPy.Latexify.latexify(- 3/(8*α) + 1/(8*α^2), cdot=false))
end


text:
using QMTK
using Compat.Test

include("Base/Base.jl")
include("Lattice/Lattice.jl")
include("Hamiltonian/Hamiltonian.jl")
include("Samplers.jl")


text:
function verify_f2(f, p, q, pa, t, integrator, ::BAOABConstantCache)
    res = f(p, q, pa, t)
    res != p && throwex(integrator)
end
function verify_f2(f, res, p, q, pa, t, integrator, ::BAOABCache)
    f(res, p, q, pa, t)
    res != p && throwex(integrator)
end
function throwex(integrator)
  algn = typeof(integrator.alg)
  throw(ArgumentError("Algorithm $algn is not applicable if f2(p, q, t) != p"))
end

function initialize!(integrator, cache::BAOABConstantCache)
  @unpack t,dt,uprev,u,p,W = integrator
  du1 = integrator.uprev.x[1]
  u1 = integrator.uprev.x[2]

  verify_f2(integrator.f.f2, du1, u1, p, t, integrator, cache)
  cache.k .= integrator.f.f1(du1,u1,p,t)
end

function initialize!(integrator, cache::BAOABCache)
  @unpack t,dt,uprev,u,p,W = integrator
  du1 = integrator.uprev.x[1]
  u1 = integrator.uprev.x[2]

  verify_f2(integrator.f.f2, cache.k, du1, u1, p, t, integrator, cache)
  integrator.f.f1(cache.k,du1,u1,p,t)
end

@muladd function perform_step!(integrator,cache::BAOABConstantCache)
  @unpack t,dt,sqdt,uprev,u,p,W,f = integrator
  @unpack k, half, c1, c2 = cache
  du1 = uprev.x[1]
  u1 = uprev.x[2]

  # B
  du2 = du1 + half*dt*k

  # A
  u2 = u1 + half*dt*du2

  # O
  noise = integrator.g(u2,p,t+dt*half).*W.dW / sqdt
  du3 = c1*du2 + c2*noise

  # A
  u = u2 + half*dt*du3

  # B
  k .= f.f1(du3,u,p,t+dt)
  du = du3 + half*dt*k

  integrator.u = ArrayPartition((du, u))
end

@muladd function perform_step!(integrator,cache::BAOABCache)
  @unpack t,dt,sqdt,uprev,u,p,W,f = integrator
  @unpack utmp, dutmp, k, gtmp, noise, half, c1, c2 = cache
  du1 = uprev.x[1]
  u1 = uprev.x[2]

  # B
  @.. dutmp = du1 + half*dt*k

  # A
  @.. utmp = u1 + half*dt*dutmp

  # O
  integrator.g(gtmp,utmp,p,t+dt*half)
  @.. noise = gtmp*W.dW / sqdt
  @.. dutmp = c1*dutmp + c2*noise

  # A
  @.. u.x[2] = utmp + half*dt*dutmp

  # B
  f.f1(k,dutmp,u.x[2],p,t+dt)
  @.. u.x[1] = dutmp + half*dt*k
end

max_stars_repo_name: JuliaDiffEq/StochasticDiffEq


text:
module Grids

export RegularCartesianGrid, VerticallyStretchedCartesianGrid

using Oceananigans

include("grid_utils.jl")
include("regular_cartesian_grid.jl")
include("vertically_stretched_cartesian_grid.jl")

end


text:
module TestPlot
using GaussianDistributions
using Plots
using Kernel
using Test
@testset "plot" begin
    l = .3
    k = Kernel.Lengthscale(Kernel.EQ(), l)
    G = Gaussian(k)

    n = 8
    σ = 1e-2
    x = randn(n)
    y = @. sin(2π*x) + σ * randn()

    C = G | (x, y, σ^2)
    fig = plot(C)
    scatter!(x, y)
    gui()
    @test fig isa Plots.Plot
end
end


text:
__precompile__()

module EFIT

include("io.jl")
export readg, GEQDSKFile

include("utils.jl")
export triangularity, ellipticity, elongation, minor_radius, major_radius, aspect_ratio

const test_gfile = (@__DIR__)*"/../test/g000001.01000"

end # module


text:
using RecipesBase
import ColorTypes.HSV

@recipe function f(b::UniformRBFE)
    Nv = size(b.μ,1)
    dist = b.μ[end]-b.μ[1]
    v = range(b.μ[1]-0.1dist, b.μ[end]+0.1dist, length=200)
    seriestype := :line
    color_palette --> [HSV(i,1,0.7) for i in range(0, 255, length=Nv)]
    xguide --> "Scheduling signal \$v\$"
    title --> "Basis function expansion"
    label --> ""
    v, b(v)
end

@recipe function f(rbf::MultiUniformRBFE, style=:default)
    # color_palette --> [HSV(i,1,0.7) for i in range(0, 255, length=Nv)]
    xguide --> "Scheduling signal \$v_1\$"
    yguide --> "Scheduling signal \$v_2\$"
    title --> "Basis function expansion"
    label --> ""
    seriestype --> :surface

    c       = rbf.μ
    minb    = minimum(c,dims=2)
    maxb    = maximum(c,dims=2)
    dist    = maxb-minb
    Npoints = 50
    v = [range(mi, ma, length=Npoints) for (mi,ma) in zip(minb,maxb)]
    vg = meshgrid(v...)
    # v  = plotattributes[:seriestype] == :surface ? vg : v # not really sure what this is for
    if style == :default
        bg = map(vg...) do v1,v2
            maximum(rbf([v1,v2]))
        end
        (v..., bg)
    elseif style == :full
        colorbar --> false
        seriesalpha --> 0.5
        bg = map(vg...) do v1,v2
            rbf([v1,v2])
        end
        bg = cat(bg...,dims=3)
        bg = reshape(bg, size(bg,1), Npoints, Npoints)
        for i = 1:size(bg,1)
            @series (v..., bg[i,:,:])
        end
    end

end


text:
using RecipesBase

#=
!!!! IMPORTANT !!!!
Each code path of the @series macro must return the data for the series.
data = ... should be the last statement
=#
@recipe function abc_output_recipe(abco::ABCOutput;
      runs=nothing, params=nothing, population_colors=nothing)
      if params === nothing
            params = [i for i in 1:abco.n_params]
      end
      if runs === nothing
            runs = [i for i in 1:length(abco.population)]
      elseif isa(runs, Int)
            runs = [i for i in max(length(abco.population) - runs + 1, 1):length(abco.population)]
      end
      legend --> false
      layout := length(params) ^ 2
      for (i, par1) in enumerate(params)
            for (j, par2) in enumerate(params)
                  subplot := (i - 1) * length(params) + j
                  if i == j
                        @series begin
                              seriestype := :histogram
                              xlabel --> "Parameter $(par1)"
                              ylabel := ""
                              yaxis := false
                              bins := 50
                              @debug repr(population_colors)
                              # if population_colors !== nothing
                              #       idx = max_color_idx % length(population_colors)
                              #       if idx == 0
                              #             idx = length(population_colors)
                              #       end
                              #       seriescolor := population_colors[idx]
                              # end
                              if isa(abco.population, Vector)
                                    data = abco.population[end][:, par1]
                                    # max_color_idx = 1
                              else
                                    data = abco.population[:, par1]
                                    # max_color_idx = length(abco.population)
                              end # if Vector
                        end  # @series
                  elseif j < i
                        seriestype := :scatter
                        markerstrokecolor --> false
                        if isa(abco.population, Vector)
                              for r in runs
                                    @series begin
                                          xlabel --> "Parameter $(par2)"
                                          ylabel --> "Parameter $(par1)"
                                          if population_colors !== nothing
                                                idx = r % length(population_colors)
                                                if idx == 0
                                                      idx = length(population_colors)
                                                end
                                                seriescolor := population_colors[idx]
                                          end
                                          pop = abco.population[r]
                                          data = (pop[:, par2], pop[:, par1])
                                    end # @series
                              end # for r
                        else
                              @series begin
                                    data = (abco.population[:, par2], abco.population[:, par1])
                              end # @series
                        end # if Vector
                  else
                        @series begin
                              ylabel := ""
                              xlabel := ""
                              grid := false
                              xaxis := false
                              yaxis := false
                              data = []
                        end # @series
                  end # if/elseif
            end # for j, par2
      end # for i, par1
end # @recipe

# Plot recipe for mdoel selection output
@recipe function modelselection_plotrecipe(::Type{ModelSelectionOutput}, mso::ModelSelectionOutput)
    seriestype := :line
    xlabel --> "Population"
    ylabel --> "Number of accepted particles"
    labels --> [string("Model ", m) for m in 1:mso.M]
    data = [[mso.n_accepted[i][j] for i in 1:size(mso.n_accepted,1)] for j in 1:mso.M]
end


text:
## Plotting intro

## Step 1, add the plotting package.
# at command line, if you press "]"
# If you want to do it not at the command line, you can also use
using Pkg # this adds the Package manager to your Julia session
Pkg.add("Plots") # this adds the Plots package
## Let's create some data to plot.
"""
`brownian_motion`
=================

Create a simple brownian motion trail. A discrete-time Brownian motion
is the value that increases by a random, normally distributed variable at
each time step. So

    \$ x_{t+1} = x_{t} + N(0,1) \$

Returns: A Vector{Float64} with n entries that describes
    a Brownian motion over n time steps .
"""
function brownian_motion(n::Int)
    x = 0.0 # start our Brownian motion at x = 0.0
    B = zeros(0) # this is an empty vector of Float64's (or length 0)
    for t=1:n # for n time-steps
        x = x+randn() # randn gives you a normally distributed random value.
        push!(B,x) # add x to the "end" of the vector B.
    end
    return B
end
brownian_motion(5)

## Questions to explore on your own.
# - Modify the brownian_motion function to start at a different place
# - Modify the brownian_motion function to pre-allocate a length n array

## Other ways to implement Brownian motion
A = randn(5)
cumsum(A)

## This works
cumsum(randn(5)) # this does everything in one line
## TODO, figure out how to use the cumsum! function to do this in one step.
# cumsum!
# slighlty more annoying than I expected
## There is another way of defining a function in Julia
# I usually try and be explicit and define functions with the word function
# like the brownian_motion function.
# But you can also, just give it a name, like in math.
bmotion(n) = cumsum(randn(n))
# This syntax is equivalent to
# function bmotion(n)
#   return cumsum(randn(n))
# end
# but is much shorter.

# Sometimes, we use this for very simple functions.

##
B = brownian_motion(100)

##
using Plots # this adds the plotting functions to Julia's session
##
using Plots
plot(B) # this will plot a vector.
##
plot(brownian_motion(1000))

## Let's see multiple Plots
plot(brownian_motion(1000))
plot!(brownian_motion(1000))
##
plot!(title="two Brownian trails")
##
plot!(xlabel="Time") # alternative xlabel!("Time")
ylabel!("Value") # alternative plot!(ylabel="Value")
##
# Let's get rid of the legend
plot!(legend=false)
##
annotate!(200,15,"Max value")
##
savefig("brownian-trails.pdf") # save as a pdf
##
savefig("brownian-trails.png") # save as a png
# savefig("brownian-trails.eps") # save as a eps
# this doesn't work with GR :(
## We can do this all at once too.
A = brownian_motion(10000)
B = brownian_motion(9000)
plot(A,B) # won't give you what you expect ...
##
plot([A,B], lab=["Trail 1", "Trail 2"],  # this gives it a pair... and labels
    xlabel="Time", ylabel="Value", xscale=:log10)
title!("Two Brownian Trails")
savefig("brownian-trails-log.pdf")

## Plotting scattered data
# Don't use syntax like this without a good reason.
A, B = brownian_motion(10000), brownian_motion(10000) # all in one linef
##
scatter(A,B,alpha=0.5, markersize=1,markerstrokewidth=0,) # plot A as a function of B
plot!(title="2D Brownian motion")
## Playing around with options
x = brownian_motion(50)
y = brownian_motion(50)
##
scatter(x,y,marker=:x)
plot!(x,y,markersize=0) # this makes the line connecting them
##


text:
function kf_legacy(obs,lin1,lin2,cov1,cov2,T,initMean,initCov)
# Coded: Jan 25 2016 / T. Lienart
# ---------------------------------------------
    dx = size(lin1,1)
    dy = size(lin2,1)
    #
    kfMean = zeros(dx,T)
    kfCov  = zeros(dx,dx,T)
    # init
    kfMean[:,1]  = initMean
    kfCov[:,:,1] = initCov
    #
    kfMean_ = copy(kfMean)
    kfCov_  = copy(kfCov)
    #
    for t=2:T
        # intermediate (t|t-1)
        mu_ = lin1*kfMean[:,t-1]
        C_  = lin1*kfCov[:,:,t-1]*lin1'+cov1
        K_  = (C_*lin2') / (cov2+lin2*C_*lin2')
        # update (t|t)
        kfMean[:,t]  = mu_+K_*(obs[:,t]-lin2*mu_)
        kfCov[:,:,t] = (eye(dx)-K_*lin2)*C_
        # storage for smoothing
        kfMean_[:,t]  = mu_
        kfCov_[:,:,t] = C_
    end
    return(kfMean,kfCov,kfMean_,kfCov_)
end

function ks_legacy(obs,lin1,lin2,cov1,cov2,T,kfMean_,kfCov_)
# Coded: Jan 25 2016 / T. Lienart
# ---------------------------------------------
    dx = size(kfMean_,1)
    dy = size(obs,1)
    #
    ksMean = kfMean_*0.
    ksCov  = kfCov_*0.
    #
    A = lin1
    B = lin2
    C = cov1
    D = cov2
    # Pre-computations
    iCA = C\A
    iDB = D\B
    #
    Pbi_T = B'*iDB
    Pbi   = copy(Pbi_T)
    c     = iDB'*obs[:,T]
    #
    # Initialization
    ksCov[:,:,T] = inv((inv(kfCov_[:,:,T])) + Pbi_T)
    ksMean[:,T]  = ksCov[:,:,T]*(kfCov_[:,:,T]\kfMean_[:,T]+c)
    #
    # KS 2FS update
    for t = (T-1):-1:1
        K   = (eye(dx)+C*Pbi)\A
        Pbi = Pbi_T+iCA'*(A-K)
        c   = iDB'*obs[:,t]+K'*c
        #
        ksCov[:,:,t] = inv(inv(kfCov_[:,:,t])+Pbi)
        ksMean[:,t]  = ksCov[:,:,t]*(kfCov_[:,:,t]\kfMean_[:,t]+c)
    end
    return ksMean,ksCov
end


text:
## Progress printer
Ni = 1000

function print_progress(simulation)
    model = simulation.model
    i, t = model.clock.iteration, model.clock.time

    progress = 100 * (model.clock.time / end_time)

    umax = maximum(abs, model.velocities.u.data.parent)
    vmax = maximum(abs, model.velocities.v.data.parent)
    wmax = maximum(abs, model.velocities.w.data.parent)
    meanb = mean(view(model.tracers.b.data, 2:Nx-1, 2:Ny-1, 2:Nz-1))

    @printf("[%05.2f%%] i: %d, t: %.2e days, umax: (%6.3e, %6.3e, %6.3e) m/s, CFL: %6.4e, next Δt: %.1e s\n",
            progress, i, t / day, umax, vmax, wmax, cfl(model), Δt_wizard.Δt)
    println(" ")
    @printf("The mean value of b is %.2e", meanb)
    println(" ")
end

## Diagnostics
if write_slices
fields = Dict(
    "u" => model.velocities.u,
    "v" => model.velocities.v,
    "w" => model.velocities.w,
    "b" => model.tracers.b
)

surface_output_writer =
    NetCDFOutputWriter(model, fields, filepath = filename * "_surface.nc",
			           schedule=TimeInterval(slice_output_interval), field_slicer = FieldSlicer(k = Nz))

middepth_output_writer =
    NetCDFOutputWriter(model, fields, filepath = filename * "_middepth.nc",
                       schedule=TimeInterval(slice_output_interval), field_slicer = FieldSlicer(k = Int(floor(Nz/2))))

zonal_output_writer =
    NetCDFOutputWriter(model, fields, filepath = filename * "_zonal.nc",
                       schedule=TimeInterval(slice_output_interval), field_slicer = FieldSlicer(j = Int(floor(Ny/2))))

meridional_output_writer =
    NetCDFOutputWriter(model, fields, filepath = filename * "_meridional.nc",
                       schedule=TimeInterval(slice_output_interval), field_slicer = FieldSlicer(i = Int(floor(Nx/2))))
function debug_f(debug)
    if debug
        close(surface_output_writer)
        close(middepth_output_writer)
        close(zonal_output_writer)
        close(meridional_output_writer)
    end
    return nothing
end
end

## Zonal and Time Averages
if write_zonal
u, v, w = model.velocities
b = model.tracers.b

u_scratch = XFaceField(model.architecture, model.grid)
v_scratch = YFaceField(model.architecture, model.grid)
w_scratch = ZFaceField(model.architecture, model.grid)
b_scratch =  CellField(model.architecture, model.grid)

zonal_fields = Dict(
    :u => AveragedField( u, dims= 1),
    :v => AveragedField( v, dims= 1),
    :w => AveragedField( w, dims= 1),
    :b => AveragedField( b, dims= 1),
    :uu => AveragedField(u * u, dims= 1),
    :vv => AveragedField(v * v, dims= 1),
    :ww => AveragedField(w * w, dims= 1),
    :uv => AveragedField(u * v, dims= 1),
    :vw => AveragedField(w * v, dims= 1),
    :uw => AveragedField(w * u, dims= 1),
    :ub => AveragedField(u * b, dims= 1),
    :vb => AveragedField(v * b, dims= 1),
    :wb => AveragedField(w * b, dims= 1),
)
#
# AveragedTimeInterval(time_avg_window, window=zonal_output_interval, stride=1)
zonalschedule = AveragedTimeInterval(zonal_output_interval, window=time_avg_window, stride=5) # TimeInterval(checkpoint_interval)
zonal_statistics = JLD2OutputWriter(model, zonal_fields,
                                    schedule = zonalschedule,
                                    prefix = filename * "_zonal_averages", force = true)
end
## Checkpointer
checkpointer = Checkpointer(model, prefix = filename * "_checkpoint",
                            schedule = TimeInterval(checkpoint_interval),
                            force = true)

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
using PyPlot,LinearAlgebra,Random
"""
```
gen_ellipse(F1:Vector{Float64},F2=Vector{Float64},
	a::Float64,p::Int,npun::Int,ntot::Int,
	tmin::Float64,tmax::Float64)
```

Generate perturbed points in a ellipse given focus F1 and F2 major axis a. It is mandatory to define the number of points `ntot` the of exacts points `p` and the number of point to compose the ellipse `npun`.

`ntot-npun` is the number of random points
`npun-p` is the number of points of the ellipse which are normally-distributed.
`tmin` and `tmax` are parameter to define limits for the figure.

## Example

```julia-repl
julia> gen_ellipse([1.0,2.0],[3.0,5.0],5.0,0,100,1000,-10,10)
```

"""
function gen_ellipse(F1,F2,a,p,npun,ntot,tmin,tmax)
	c = norm(F2-F1,2)
	if a<=c
		error("a<=c, it isn't an elipse")
	end
	b = sqrt(a^2-c^2)
	# let us assume β is the rotation angle from (1,0)
	cosβ = (F2-F1)[1]/norm(F2-F1,2)
	sinβ = sqrt(1.0-cosβ)

	center = 0.5*(F1+F2)
	rng = MersenneTwister(1234)
	x = zeros(ntot)
	y = zeros(ntot)
	θ = shuffle(rng,[0:(2*π)/npun:2*π;])

	for i=1:p
		x[i] = center[1]*cosβ+b*cos(θ[i])*cosβ+center[2]*sinβ+a*sin(θ[i])*sinβ
		y[i] = -center[1]*sinβ -b*cos(θ[i])*sinβ+center[2]*cosβ+a*sin(θ[i])*cosβ
	end
	for i=p+1:npun
		anoise = a+0.1*randn()
		bnoise = b+0.1*randn()
		x[i] = center[1]*cosβ+bnoise*cos(θ[i])*cosβ+center[2]*sinβ+anoise*sin(θ[i])*sinβ
		y[i] = -center[1]*sinβ -bnoise*cos(θ[i])*sinβ+center[2]*cosβ+anoise*sin(θ[i])*cosβ
	end

	x[npun+1:end]= tmin.+(tmax-tmin)*rand(ntot-npun)
	y[npun+1:end]= tmin.+(tmax-tmin)*rand(ntot-npun)
	plot(x,y,".")

	return [x y]
end


text:
##file to include in ComplexModeler.jl which contains all files associated to Matpower

include("MatpowerGenerator.jl")
include("MatpowerLoad.jl")
include("MatpowerShunt.jl")
include("MatpowerLine_π.jl")
include("MatpowerLine_π_Smax.jl")
include("MatpowerLine_π_Imax.jl")
include("MatpowerVolt.jl")
include("read_matpower.jl")
include("read_matpower_Smax.jl")
include("read_matpower_simple_Imax.jl")
include("read_matpower_simple.jl")
include("read_matpower_graph.jl")


text:
output_files = ["../example_output/output.e"]

function test_on_output(output_file::String)
    output_file_name = abspath(output_file)
    test_name = rpad("Testing output file: $(basename(output_file_name))", 96)
    @testset "$test_name" begin
        exo_id = Exodus.open_exodus_database(output_file_name)
        init = Exodus.Initialization(exo_id)
        num_vars = Exodus.read_number_of_nodal_variables(exo_id)
        @test num_vars == 1
        variable_names = Exodus.read_nodal_variable_names(exo_id)
        @test variable_names == ["u"]
        num_time_steps = Exodus.read_number_of_time_steps(exo_id)
        @test num_time_steps == 2
        times = Exodus.read_times(exo_id)
        @test times == [0.0, 1.0]

        for (n, var_name) in enumerate(variable_names)
            for (n_t, time) in enumerate(times)
                values = Exodus.read_nodal_variable_values(exo_id, n_t, n, init.num_nodes)
            end
        end
    end
end

for output_file in output_files
    test_on_output(output_file)
end


text:
module Utils

export Resolution
export get_sinewave
export get_evolved
export get_evolved_slice
export commandlineparser
export commandlineparser_int
export check_time

using Burgers

function get_sinewave(x)
	# returns a sine wave with a period if 1 on [0, 1]
	return 0.35 .+ 0.25*sin.(2*pi*x)
end

struct Resolution{T, I}
	# contains information about the resolution of the 2D (x, t) grid
	# that will be used by the FunArr2D instances
	Ttot::T
	Nt::I
	Ncells::I
	Nghost::I
end

function check_time(res::Resolution{T}, t0::T) where {T}
	# makes sure that the given time is valid for the resolution
	return 0.0 <= t0 < res.Ttot
end

function get_evolved_slice(func::Function, res::Resolution{T, I}, t::T, sc::SlopeCalc)::FunArr{T, T} where {T, I<:Int}
	# evaluates 'func' on [0, 1] to initialize u(x, 0), and evolves in time,
	# then returns the function u(x) evaluated at a time t
	it = Int64(res.Nt*floor(t/res.Ttot))

	init_positions = FunArr{Float64, Float64}(func, res.Ncells, res.Nghost)
	u2fun          = evolve(init_positions, res.Nt, res.Ttot, sc, it)
	ufun_t         = get_timeslice(u2fun, t, res.Ncells, res.Nghost)

	return ufun_t
end

function get_evolved_slice(func::Function, res::Resolution{T, I}, it::I, sc::SlopeCalc)::FunArr{T, T} where {T, I<:Int}
	# evaluates 'func' on [0, 1] to initialize u(x, 0), and evolves in time,
	# then returns the function u(x) evaluated at the timeslice it
	init_positions = FunArr{Float64, Float64}(func, res.Ncells, res.Nghost)
	u2fun          = evolve(init_positions, res.Nt, res.Ttot, sc, it)
	uarr_t         = u2fun.u2dat[:,it]
	ufun_t         = FunArr{T, T}(uarr_t, res.Ncells, 0)

	return ufun_t
end

function get_evolved(func::Function, res::Resolution{T, I}, sc::SlopeCalc)::FunArr2D{T, T} where {T, I<:Int}
	# evaluates 'func' on [0, 1] to initialize u(x, 0), and evolves in time,
	# then returns the function u(x) evaluated at a time t
	init_positions = FunArr{Float64, Float64}(func, res.Ncells, res.Nghost)
	u2fun          = evolve(init_positions, res.Nt, res.Ttot, sc, res.Nt)

	return u2fun
end

function commandlineparser(check_int::Bool)::Tuple{Number, SlopeCalc, AbstractString}
	# parses the command line arguments
	if length(ARGS) != 3
		if check_int == true
			error("This program requires 3 arguments: (1) index for the evolved time, (2) the type of slope calculation to perform, and (3) the name of the file where the image will be saved.")
		else
			error("This program requires 3 arguments: (1) the evolved time, (2) the type of slope calculation to perform, and (3) the name of the file where the image will be saved.")
		end
	end

	t_str, sc_str, fname = ARGS

	# parse the slope calculation type
	sc =
	if sc_str == "minmod"
		minmod::SlopeCalc
	elseif sc_str == "linear"
		linear::SlopeCalc
	elseif sc_str == "constant"
		constant::SlopeCalc
	else
		error("the second argument must be one of: minmod linear constant")
	end

	# parse the time
	if check_int == true
		it = tryparse(Int64, t_str)
		if typeof(it) == Nothing
			error("the first argument must be an Int64 value.")
		end
		return it, sc, fname
	else
		t0 = tryparse(Float64, t_str)
		if typeof(t0) == Nothing
			error("the first argument must be a Float64 value.")
		end
		return t0, sc, fname
	end
end

end #module


text:
function overpass_polygon_buildings_query(geojson_polygons::Vector{Vector{Any}},
                                          metadata::Bool=false,
                                          download_format::Symbol=:osm
                                          )::String
    filters = ""
    for polygon in geojson_polygons
        polygon = map(x -> [x[2], x[1]], polygon) # switch lon-lat to lat-lon
        polygon_str = replace("$polygon", r"[\[,\]]" =>  "")
        filters *= """node["building"](poly:"$polygon_str");<;way["building"](poly:"$polygon_str");>;rel["building"](poly:"$polygon_str");>;"""
    end

    return overpass_query(filters, metadata, download_format)
end

"""
    osm_buildings_from_place_name(;place_name::String,
                                metadata::Bool=false,
                                download_format::Symbol=:osm
                                )::String

Downloads OpenStreetMap buildings using any place name string.

# Arguments
- `place_name::String`: Any place name string used as a search argument to the Nominatim API.
- `metadata::Bool=false`: Set true to return metadata.
- `download_format::Symbol=:osm`: Download format, either `:osm`, `:xml` or `json`.

# Return
- `String`: OpenStreetMap buildings data response string.
"""
function osm_buildings_from_place_name(;place_name::String,
                                       metadata::Bool=false,
                                       download_format::Symbol=:osm
                                       )::String
    geojson_polygons = polygon_from_place_name(place_name)
    query = overpass_polygon_buildings_query(geojson_polygons, metadata, download_format)
    return overpass_request(query)
end

"""
    osm_buildings_from_bbox(;minlat::AbstractFloat,
                          minlon::AbstractFloat,
                          maxlat::AbstractFloat,
                          maxlon::AbstractFloat,
                          metadata::Bool=false,
                          download_format::Symbol=:osm
                          )::String

Downloads OpenStreetMap buildings using bounding box coordinates.

# Arguments
- `minlat::AbstractFloat`: Bottom left bounding box latitude coordinate.
- `minlon::AbstractFloat`: Bottom left bounding box longitude coordinate.
- `maxlat::AbstractFloat`: Top right bounding box latitude coordinate.
- `maxlon::AbstractFloat`: Top right bounding box longitude coordinate.
- `metadata::Bool=false`: Set true to return metadata.
- `download_format::Symbol=:osm`: Download format, either `:osm`, `:xml` or `json`.

# Return
- `String`: OpenStreetMap buildings data response string.
"""
function osm_buildings_from_bbox(;minlat::Float64,
                                 minlon::Float64,
                                 maxlat::Float64,
                                 maxlon::Float64,
                                 metadata::Bool=false,
                                 download_format::Symbol=:osm
                                 )::String
    filters = """node["building"];<;way["building"];>;rel["building"];>;"""
    bbox = [minlat, minlon, maxlat, maxlon]
    query = overpass_query(filters, metadata, download_format, bbox)
    return overpass_request(query)
end

"""
    osm_buildings_from_point(;point::GeoLocation,
                           radius::Number,
                           metadata::Bool=false,
                           download_format::Symbol=:osm
                           )::String

Downloads OpenStreetMap buildings using bounding box coordinates calculated from a centroid point and radius (km).

# Arguments
- `point::GeoLocation`: Centroid point to draw the bounding box around.
- `radius::Number`: Distance (km) from centroid point to each bounding box corner.
- `metadata::Bool=false`: Set true to return metadata.
- `download_format::Symbol=:osm`: Download format, either `:osm`, `:xml` or `json`.

# Return
- `String`: OpenStreetMap buildings data response string.
"""
function osm_buildings_from_point(;point::GeoLocation,
                                  radius::Number,
                                  metadata::Bool=false,
                                  download_format::Symbol=:osm
                                  )::String
    bbox = bounding_box_from_point(point, radius)
    return osm_buildings_from_bbox(;bbox..., metadata=metadata, download_format=download_format)
end

function osm_buildings_downloader(download_method::Symbol)::Function
    if download_method == :place_name
        return osm_buildings_from_place_name
    elseif download_method == :bbox
        return osm_buildings_from_bbox
    elseif download_method == :point
        return osm_buildings_from_point
    else
        throw(ErrorException("OSM buildings downloader $download_method does not exist"))
    end
end

"""
    download_osm_buildings(download_method::Symbol;
                           metadata::Bool=false,
                           download_format::Symbol=:osm,
                           save_to_file_location::Union{String,Nothing}=nothing,
                           download_kwargs...
                           )::Union{XMLDocument,Dict{String,Any}}

Downloads OpenStreetMap buildings data by querying with a place name, bounding box, or centroid point.

# Arguments
- `download_method::Symbol`: Download method, choose from `:place_name`, `:bbox` or `:point`.
- `metadata::Bool=false`: Set true to return metadata.
- `download_format::Symbol=:osm`: Download format, either `:osm`, `:xml` or `json`.
- `save_to_file_location::Union{String,Nothing}=nothing`: Specify a file location to save downloaded data to disk.

# Required Download Kwargs

*`download_method=:place_name`*
- `place_name::String`: Any place name string used as a search argument to the Nominatim API.

*`download_method=:bbox`*
- `minlat::AbstractFloat`: Bottom left bounding box latitude coordinate.
- `minlon::AbstractFloat`: Bottom left bounding box longitude coordinate.
- `maxlat::AbstractFloat`: Top right bounding box latitude coordinate.
- `maxlon::AbstractFloat`: Top right bounding box longitude coordinate.

*`download_method=:point`*
- `point::GeoLocation`: Centroid point to draw the bounding box around.
- `radius::Number`: Distance (km) from centroid point to each bounding box corner.

# Return
- `Union{XMLDocument,Dict{String,Any}}`: OpenStreetMap buildings data parsed as either XML or Dictionary object depending on the download method.
"""
function download_osm_buildings(download_method::Symbol;
                                metadata::Bool=false,
                                download_format::Symbol=:osm,
                                save_to_file_location::Union{String,Nothing}=nothing,
                                download_kwargs...
                                )::Union{XMLDocument,Dict{String,Any}}
    downloader = osm_buildings_downloader(download_method)
    data = downloader(metadata=metadata, download_format=download_format; download_kwargs...)
    @info "Downloaded osm buildings data from $(["$k: $v" for (k, v) in download_kwargs]) in $download_format format"

    if !(save_to_file_location isa Nothing)
        save_to_file_location = validate_save_location(save_to_file_location, download_format)
        write(save_to_file_location, data)
        @info "Saved osm buildings data to disk: $save_to_file_location"
    end

    deserializer = string_deserializer(download_format)

    return deserializer(data)
end

is_building(tags::Dict)::Bool = haskey(tags, "building") ? true : false

function height(tags::Dict)::Number
    height = get(tags, "height", nothing)
    levels = get(tags, "building:levels", nothing) !== nothing ? tags["building:levels"] : get(tags, "level", nothing)

    if height !== nothing
        return height isa String ? max([remove_non_numeric(h) for h in split(height, r"[+^;,-]")]...) : height
    elseif levels !== nothing
        levels = levels isa String ? round(max([remove_non_numeric(l) for l in split(levels, r"[+^;,-]")]...)) : levels
        levels = levels == 0 ? rand(1:DEFAULT_MAX_BUILDING_LEVELS) : levels
    else
        levels = rand(1:DEFAULT_MAX_BUILDING_LEVELS)
    end

    return levels * DEFAULT_BUILDING_HEIGHT_PER_LEVEL
end

function parse_osm_buildings_dict(osm_buildings_dict::AbstractDict)::Dict{Integer,Building}
    T = DEFAULT_OSM_ID_TYPE

    nodes = Dict{T,Node{T}}()
    for node in osm_buildings_dict["node"]
        id = node["id"]
        nodes[id] = Node{T}(
            id,
            GeoLocation(node["lat"], node["lon"]),
            haskey(node, "tags") ? node["tags"] : nothing
        )
    end

    ways = Dict(way["id"] => way for way in osm_buildings_dict["way"]) # for lookup

    added_ways = Set{T}()
    buildings = Dict{T,Building{T}}()
    for relation in osm_buildings_dict["relation"]
        is_relation = true

        if haskey(relation, "tags") && is_building(relation["tags"])
            tags = relation["tags"]
            rel_id = relation["id"]
            members = relation["members"]

            polygons = Vector{Polygon{T}}()
            for member in members
                way_id = member["ref"]
                way = ways[way_id]

                haskey(way, "tags") && merge!(tags, way["tags"]) # could potentially overwrite some data
                push!(added_ways, way_id)

                is_outer = member["role"] == "outer" ? true : false
                nds = [nodes[n] for n in way["nodes"]]
                push!(polygons, Polygon(way_id, nds, is_outer))
            end

            tags["height"] = height(tags)
            sort!(polygons, by = x -> x.is_outer, rev=true) # sorting so outer polygon is always first
            buildings[rel_id] = Building{T}(rel_id, is_relation, polygons, tags)
        end
    end

    for (way_id, way) in ways
        is_relation = false
        is_outer = true

        if haskey(way, "tags") && is_building(way["tags"]) && !(way_id in added_ways)
            tags = way["tags"]
            tags["height"] = height(tags)
            nds = [nodes[n] for n in way["nodes"]]
            polygons = [Polygon(way_id, nds, is_outer)]
            buildings[way_id] = Building{T}(way_id, is_relation, polygons, tags)
        end
    end

    return buildings
end

"""
    buildings_from_object(buildings_xml_object::XMLDocument)::Dict{Integer,Building}

Creates `Building` objects from data downloaded with `download_osm_buildings`.

# Arguments
- `buildings_xml_object::XMLDocument`: Buildings data downloaded in `:xml` or `:osm` format.

# Return
- `Dict{Integer,Building}`: Mapping from building relation/way ids to `Building` objects.
"""
function buildings_from_object(buildings_xml_object::XMLDocument)::Dict{Integer,Building}
    dict_to_parse = osm_dict_from_xml(buildings_xml_object)
    return parse_osm_buildings_dict(dict_to_parse)
end

"""
    buildings_from_object(buildings_json_object::AbstractDict)::Dict{Integer,Building}

Creates `Building` objects from data downloaded with `download_osm_buildings`.

# Arguments
- `buildings_json_object::AbstractDict`: Buildings data downloaded in `:json` format.

# Return
- `Dict{Integer,Building}`: Mapping from building relation/way ids to `Building` objects.
"""
function buildings_from_object(buildings_json_object::AbstractDict)::Dict{Integer,Building}
    dict_to_parse = osm_dict_from_json(buildings_json_object)
    return parse_osm_buildings_dict(dict_to_parse)
end

"""
    buildings_from_file(file_path::String)::Dict{Integer,Building}

Creates `Building` objects from OpenStreetMap data file (either `:osm`, `:xml` or `:json` format).

# Arguments
- `file_path::String`: Path to OpenStreetMap data file.

# Return
- `Dict{Integer,Building}`: Mapping from building relation/way ids to `Building` objects.
"""
function buildings_from_file(file_path::String)::Dict{Integer,Building}
    !isfile(file_path) && throw(ArgumentError("File $file_path does not exist"))
    deserializer = file_deserializer(file_path)
    obj = deserializer(file_path)
    return buildings_from_object(obj)
end

"""
    buildings_from_download(download_method::Symbol;
                            metadata::Bool=false,
                            download_format::Symbol=:osm,
                            save_to_file_location::Union{String,Nothing}=nothing,
                            download_kwargs...
                            )::Dict{Integer,Building}

Downloads and Creates `Building` objects from OpenStreetMap APIs.

# Arguments
- `download_method::Symbol`: Download method, choose from `:place_name`, `:bbox` or `:point`.
- `metadata::Bool=false`: Set true to return metadata.
- `download_format::Symbol=:osm`: Download format, either `:osm`, `:xml` or `json`.
- `save_to_file_location::Union{String,Nothing}=nothing`: Specify a file location to save downloaded data to disk.

# Required Download Kwargs

*`download_method=:place_name`*
- `place_name::String`: Any place name string used as a search argument to the Nominatim API.

*`download_method=:bbox`*
- `minlat::AbstractFloat`: Bottom left bounding box latitude coordinate.
- `minlon::AbstractFloat`: Bottom left bounding box longitude coordinate.
- `maxlat::AbstractFloat`: Top right bounding box latitude coordinate.
- `maxlon::AbstractFloat`: Top right bounding box longitude coordinate.

*`download_method=:point`*
- `point::GeoLocation`: Centroid point to draw the bounding box around.
- `radius::Number`: Distance (km) from centroid point to each bounding box corner.

# Return
- `Dict{Integer,Building}`: Mapping from building relation/way ids to `Building` objects.
"""
function buildings_from_download(download_method::Symbol;
                                 metadata::Bool=false,
                                 download_format::Symbol=:osm,
                                 save_to_file_location::Union{String,Nothing}=nothing,
                                 download_kwargs...
                                 )::Dict{Integer,Building}
    obj = download_osm_buildings(download_method,
                                 metadata=metadata,
                                 download_format=download_format,
                                 save_to_file_location=save_to_file_location;
                                 download_kwargs...)
    return buildings_from_object(obj)
end


text:
@testset "umap tests" begin

    @testset "constructor" begin
        @testset "argument validation tests" begin
            data = rand(5, 10)
            @test_throws ArgumentError UMAP_([1. 1.]; n_neighbors=0) # n_neighbors error
            @test_throws ArgumentError UMAP_([1. 1.], 0; n_neighbors=1) # n_comps error
            @test_throws ArgumentError UMAP_([1. 1.], 2; n_neighbors=1) # n_comps error
            @test_throws ArgumentError UMAP_([1. 1.; 1. 1.; 1. 1.];
                    n_neighbors=1, min_dist = 0.) # min_dist error
        end
    end

    @testset "input type stability tests" begin
        data = rand(5, 100)
        umap_ = UMAP_(data; init=:random)
        @test umap_ isa UMAP_{Float64}
        @test size(umap_.graph) == (100, 100)
        @test size(umap_.embedding) == (2, 100)
        @test umap_.data === data

        data = rand(Float32, 5, 100)
        @test UMAP_(data; init=:random) isa UMAP_{Float32}
    end

    @testset "fuzzy_simpl_set" begin
        # knns = rand(1:50, 5, 50)
        # dists = rand(5, 50)
        knns = [2 3 2; 3 1 1]
        dists = [1.5 .5 .5; 2. 1.5 2.]
        k = 3
        umap_graph = fuzzy_simplicial_set(knns, dists, k, 3, 1, 1.)
        @test issymmetric(umap_graph)
        @test all(0. .<= umap_graph .<= 1.)
        @test size(umap_graph) == (3, 3)

        dists = convert.(Float32, dists)
        umap_graph = fuzzy_simplicial_set(knns, dists, k, 3, 1, 1.f0)
        @test issymmetric(umap_graph)
        @test eltype(umap_graph) == Float32

        umap_graph = fuzzy_simplicial_set(knns, dists, k, 200, 1, 1., false)
        @test all(0. .<= umap_graph .<= 1.)
        @test size(umap_graph) == (200, 3)
    end

    @testset "smooth_knn_dists" begin
        dists = [0., 1., 2., 3., 4., 5.]
        rho = 1
        k = 6
        local_connectivity = 1
        bandwidth = 1.
        niter = 64
        sigma = smooth_knn_dist(dists, rho, k, bandwidth, niter)
        psum(ds, r, s) = sum(exp.(-max.(ds .- r, 0.) ./ s))
        @test psum(dists, rho, sigma) - log2(k)*bandwidth < SMOOTH_K_TOLERANCE

        knn_dists = [0. 0. 0.;
                     1. 2. 3.;
                     2. 4. 5.;
                     3. 4. 5.;
                     4. 6. 6.;
                     5. 6. 10.]
        rhos, sigmas = smooth_knn_dists(knn_dists, k, local_connectivity)
        @test rhos == [1., 2., 3.]
        diffs = [psum(knn_dists[:,i], rhos[i], sigmas[i]) for i in 1:3] .- log2(6)
        @test all(diffs .< SMOOTH_K_TOLERANCE)

        knn_dists = [0. 0. 0.;
                     0. 1. 2.;
                     0. 2. 3.]
        rhos, sigmas = smooth_knn_dists(knn_dists, 2, 1)
        @test rhos == [0., 1., 2.]

        rhos, sigmas = smooth_knn_dists(knn_dists, 2, 1.5)
        @test rhos == [0., 1.5, 2.5]
    end

    @testset "compute_membership_strengths" begin
        knns = [1 2 3; 2 1 2]
        dists = [0. 0. 0.; 2. 2. 3.]
        rhos = [2., 1., 4.]
        sigmas = [1., 1., 1.]
        true_rows = [1, 2, 2, 1, 3, 2]
        true_cols = [1, 1, 2, 2, 3, 3]
        true_vals = [0., 1., 0., exp(-1.), 0., 1.]
        rows, cols, vals = compute_membership_strengths(knns, dists, rhos, sigmas)
        @test rows == true_rows
        @test cols == true_cols
        @test vals == true_vals
    end

    @testset "optimize_embedding" begin
        graph1 = sparse(Symmetric(sprand(6,6,0.4)))
        graph2 = sparse(sprand(5,3,0.4))
        graph3 = sparse(sprand(3,5,0.4))
        embedding = Vector{Float64}[[3, 1], [4, 5], [2, 3], [1, 7], [6, 3], [2, 6]]

        n_epochs = 1
        initial_alpha = 1.
        min_dist = 1.
        spread = 1.
        gamma = 1.
        neg_sample_rate = 5
        for graph in [graph1, graph2, graph3]
            ref_embedding = collect(eachcol(rand(2, size(graph, 1))))
            old_ref_embedding = deepcopy(ref_embedding)
            query_embedding = rand(2, size(graph, 2))
            query_embedding = [query_embedding[:, i] for i in 1:size(query_embedding, 2)]
            res_embedding = optimize_embedding(graph, query_embedding, ref_embedding, n_epochs, initial_alpha,
                                               min_dist, spread, gamma, neg_sample_rate, move_ref=false)
            @test res_embedding isa Array{Array{Float64, 1}, 1}
            @test length(res_embedding) == length(query_embedding)
            for i in 1:length(res_embedding)
                @test length(res_embedding[i]) == length(query_embedding[i])
            end
            @test isapprox(old_ref_embedding, ref_embedding, atol=1e-4)
        end
    end

    @testset "spectral_layout" begin
        A = sprand(10000, 10000, 0.001)
        B = dropzeros(A + A' - A .* A')
        layout = spectral_layout(B, 5)
        @test layout isa Array{Float64, 2}
        @inferred spectral_layout(B, 5)
        layout32 = spectral_layout(convert(SparseMatrixCSC{Float32}, B), 5)
        @test layout32 isa Array{Float32, 2}
        @inferred spectral_layout(convert(SparseMatrixCSC{Float32}, B), 5)
    end

    @testset "initialize_embedding" begin
        graph = [5 0 1 1;
                 2 4 1 1;
                 3 6 8 8] ./10
        ref_embedding = Float64[1 2 0;
                                0 2 -1]
        actual = [[9, 1], [8, 2], [3, -6], [3, -6]] ./10

        embedding = initialize_embedding(graph, ref_embedding)
        @test embedding isa AbstractVector{<:AbstractVector{Float64}}
        @test length(embedding) == length(actual)
        for i in 1:length(embedding)
            @test length(embedding[i]) == length(actual[i])
        end
        @test isapprox(embedding, actual, atol=1e-8)

        graph = Float16.(graph[:, [1,2]])
        graph[:, end] .= 0
        ref_embedding = Float16[1 2 0;
                                0 2 -1]
        actual = Vector{Float16}[[9, 1], [0, 0]] ./10

        embedding = initialize_embedding(graph, ref_embedding)
        @test embedding isa AbstractVector{<:AbstractVector{Float16}}
        @test length(embedding) == length(actual)
        for i in 1:length(embedding)
            @test length(embedding[i]) == length(actual[i])
        end
        @test isapprox(embedding, actual, atol=1e-2)
    end

    @testset "umap transform" begin
        @testset "argument validation tests" begin
            data = rand(5, 10)
            model = UMAP_(data, 2, n_neighbors=2, n_epochs=1)
            query = rand(5, 8)
            @test_throws ArgumentError transform(model, rand(6, 8); n_neighbors=3) # query size error
            @test_throws ArgumentError transform(model, query; n_neighbors=0) # n_neighbors error
            @test_throws ArgumentError transform(model, query; n_neighbors=15) # n_neighbors error
            @test_throws ArgumentError transform(model, query; n_neighbors=1, min_dist = 0.) # min_dist error

            model = UMAP_(model.graph, model.embedding, rand(6, 10), model.knns, model.dists)
            @test_throws ArgumentError transform(model, query; n_neighbors=3) # data size error
            model = UMAP_(model.graph, model.embedding, rand(5, 9), model.knns, model.dists)
            @test_throws ArgumentError transform(model, query; n_neighbors=3) # data size error
        end

        @testset "transform test" begin
            data = rand(5, 30)
            model = UMAP_(data, 2, n_neighbors=2, n_epochs=1)
            embedding = transform(model, rand(5, 10), n_epochs=5, n_neighbors=5)
            @test size(embedding) == (2, 10)
            @test typeof(embedding) == typeof(model.embedding)

            data = rand(Float32, 5, 30)
            model = UMAP_(data, 2, n_neighbors=2, n_epochs=1)
            embedding = @inferred transform(model, rand(5, 50), n_epochs=5, n_neighbors=5)
            @test size(embedding) == (2, 50)
            @test typeof(embedding) == typeof(model.embedding)
        end
    end
end


text:
# abstract monte carlo flavor definition
"""
Abstract definition of a Monte Carlo flavor.
"""
abstract type MonteCarloFlavor end

# MonteCarloFlavor interface: mandatory
"""
    init!(mc)

Initialize the Monte Carlo simulation. Has an alias function `reset!`.
"""
init!(mc::MonteCarloFlavor) = error("MonteCarloFlavor $(typeof(mc)) doesn't implement `init!`!")

"""
    run!(mc)

Run the Monte Carlo Simulation.
"""
run!(mc::MonteCarloFlavor) = error("MonteCarloFlavor $(typeof(mc)) doesn't implement `run!`!")

# general functions
"""
    reset!(mc::MonteCarloFlavor)

Resets the Monte Carlo simulation `mc`.
Previously set parameters will be retained.
"""
function reset!(mc::MonteCarloFlavor)
    th_meas = Dict{Symbol, AbstractMeasurement}([
        k => typeof(v)(mc, mc.model) for (k, v) in mc.thermalization_measurements
    ])

    meas = Dict{Symbol, AbstractMeasurement}([
        k => typeof(v)(mc, mc.model) for (k, v) in mc.measurements
    ])

    init!(mc, thermalization_measurements=th_meas, measurements=meas)
end


text:
triangular = [ div(n * (n + 1), 2) for n in 1:100000 ]
pentagonal = [ div(n * (3n - 1),  2) for n in 1:100000 ]
hexagonal = [ n * (2n - 1) for n in 1:100000 ]

final = intersect(triangular, pentagonal, hexagonal)

println(final)

max_stars_repo_name: vmchale/julia-problem-solving

max_issues_repo_name: vmchale/julia-problem-solving

max_forks_repo_name: vmchale/julia-problem-solving


text:
using ModelingToolkit, ModelingToolkitStandardLibrary, OrdinaryDiffEq
using ModelingToolkitStandardLibrary.Blocks
using ModelingToolkitStandardLibrary.Blocks: smooth_sin, smooth_cos, smooth_damped_sin, smooth_square, smooth_step, smooth_ramp, smooth_triangular, triangular, square

@parameters t

@testset "Constant" begin
    @named src = Constant(k=2)
    @named int = Integrator()
    @named iosys = ODESystem([
        connect(src.output, int.input),
        ],
        t,
        systems=[int, src],
    )
    sys = structural_simplify(iosys)

    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0))

    sol = solve(prob, Rodas4())
    @test sol.retcode == :Success
    @test sol[src.output.u][end] ≈ 2 atol=1e-3
end

@testset "Sine" begin
    sine(t, frequency, amplitude, phase, offset, start_time) = offset + ifelse(t < start_time, 0, amplitude* sin(2*pi*frequency*(t - start_time) + phase))

    frequency=1
    amplitude=2
    phase=0
    offset=1
    start_time=2
    δ = 1e-5
    @named int = Integrator()

    @named src = Sine(frequency=frequency, amplitude=amplitude, phase=phase, offset=offset, start_time=start_time)
    @named iosys = ODESystem([
        connect(src.output, int.input),
        ],
        t,
        systems=[int, src],
    )
    sys = structural_simplify(iosys)

    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0))

    sol = solve(prob, Rodas4())
    @test sol.retcode == :Success
    @test sol[src.output.u] ≈ sine.(sol.t, frequency, amplitude, phase, offset, start_time) atol=1e-3

    @named smooth_src = Sine(frequency=frequency,
        amplitude=amplitude,
        phase=phase,
        offset=offset,
        start_time=start_time,
        smooth=true)
    @named smooth_iosys = ODESystem([
        connect(smooth_src.output, int.input),
        ],
        t,
        systems=[int, smooth_src],
    )

    smooth_sys = structural_simplify(smooth_iosys)
    smooth_prob = ODEProblem(smooth_sys, Pair[int.x=>0.0], (0.0, 10.0))
    smooth_sol = solve(smooth_prob, Rodas4())

    @test sol.retcode == :Success
    @test smooth_sol[smooth_src.output.u] ≈ smooth_sin.(smooth_sol.t, δ, frequency, amplitude, phase, offset, start_time) atol=1e-3
end

@testset "Cosine" begin
    cosine(t, frequency, amplitude, phase, offset, start_time) =
        offset + ifelse(t < start_time, 0, amplitude* cos(2*pi*frequency*(t - start_time) + phase))

    frequency=1
    amplitude=2
    phase=0
    offset=1
    start_time=2
    δ = 1e-5
    @named int = Integrator()

    @named src = Cosine(frequency=frequency,
        amplitude=amplitude,
        phase=phase,
        offset=offset,
        start_time=start_time,
        smooth=false)
    @named iosys = ODESystem([
        connect(src.output, int.input),
        ],
        t,
        systems=[int, src],
    )

    sys = structural_simplify(iosys)
    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0))
    sol = solve(prob, Rodas4())

    @test sol.retcode == :Success
    @test sol[src.output.u] ≈ cosine.(sol.t, frequency, amplitude, phase, offset, start_time) atol=1e-3

    @named smooth_src = Cosine(frequency=frequency,
        amplitude=amplitude,
        phase=phase,
        offset=offset,
        start_time=start_time,
        smooth=true)
    @named smooth_iosys = ODESystem([
        connect(smooth_src.output, int.input),
        ],
        t,
        systems=[int, smooth_src],
    )

    smooth_sys = structural_simplify(smooth_iosys)
    smooth_prob = ODEProblem(smooth_sys, Pair[int.x=>0.0], (0.0, 10.0))
    smooth_sol = solve(smooth_prob, Rodas4())

    @test smooth_sol.retcode == :Success
    @test smooth_sol[smooth_src.output.u] ≈ smooth_cos.(smooth_sol.t, δ, frequency, amplitude, phase, offset, start_time) atol=1e-3
end

@testset "ContinuousClock" begin
    cont_clock(t, offset, start_time) = offset + ifelse(t < start_time, 0, t - start_time)

    offset, start_time = 1, 0

    @named src = ContinuousClock(offset=offset, start_time=start_time)
    @named int = Integrator()
    @named iosys = ODESystem([
        connect(src.output, int.input),
        ],
        t,
        systems=[int, src],
    )
    sys = structural_simplify(iosys)

    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0))

    sol = solve(prob, Rodas4())
    @test sol.retcode == :Success
    @test sol[src.output.u] ≈ cont_clock.(sol.t, offset, start_time) atol=1e-3
end

@testset "Ramp" begin
    ramp(t, offset, height, duration, start_time) = offset + ifelse(t < start_time, 0, ifelse(t < (start_time + duration), (t - start_time) * height / duration, height))

    offset, height, duration, start_time, δ = 1, 2, 2, 0, 1e-5
    @named int = Integrator()

    @named src = Ramp(offset=offset, height=height, duration=duration, start_time=start_time)
    @named iosys = ODESystem([
        connect(src.output, int.input),
        ],
        t,
        systems=[int, src],
    )

    sys = structural_simplify(iosys)
    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0))
    sol = solve(prob, Rodas4())

    @test sol.retcode == :Success
    @test sol[src.output.u] ≈ ramp.(sol.t, offset, height, duration, start_time) atol=1e-3

    start_time = 2
    @named smooth_src = Ramp(offset=offset, height=height, duration=duration, start_time=start_time, smooth=true)
    @named smooth_iosys = ODESystem([
        connect(smooth_src.output, int.input),
        ],
        t,
        systems=[int, smooth_src],
    )

    smooth_sys = structural_simplify(smooth_iosys)
    smooth_prob = ODEProblem(smooth_sys, Pair[int.x=>0.0], (0.0, 10.0))
    smooth_sol = solve(smooth_prob, Rodas4())

    @test smooth_sol.retcode == :Success
    @test smooth_sol[smooth_src.output.u] ≈ smooth_ramp.(smooth_sol.t, δ, height, duration, offset, start_time) atol=1e-3
end

@testset "Step" begin
    step(t, offset, height, start_time) = offset + ifelse(t < start_time, 0, height)

    offset, height, start_time, δ = 1, 2, 5, 1e-5
    @named int = Integrator()

    @named src = Step(offset=offset, height=height, start_time=start_time, smooth=false)
    @named iosys = ODESystem([
        connect(src.output, int.input),
        ],
        t,
        systems=[int, src],
    )

    sys = structural_simplify(iosys)
    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0))
    sol = solve(prob, Rodas4())

    @test sol.retcode == :Success
    @test sol[src.output.u] ≈ step.(sol.t, offset, height, start_time) atol=1e-2

    # test with duration
    duration = 1.2
    @named src = Step(offset=offset, height=height, start_time=start_time, duration=duration, smooth=false)
    @named iosys = ODESystem([
        connect(src.output, int.input),
        ],
        t,
        systems=[int, src],
    )

    sys = structural_simplify(iosys)
    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0))
    sol = solve(prob, Rodas4(), dtmax=0.1) # set dtmax to prevent the solver from overstepping the entire step disturbance

    @test sol.retcode == :Success
    @test sol[src.output.u] ≈ step.(sol.t, offset, height, start_time) - step.(sol.t, 0, height, start_time+duration) atol=1e-2

    @named smooth_src = Step(offset=offset, height=height, start_time=start_time, smooth=true)
    @named smooth_iosys = ODESystem([
        connect(smooth_src.output, int.input),
        ],
        t,
        systems=[int, smooth_src],
    )

    smooth_sys = structural_simplify(smooth_iosys)
    smooth_prob = ODEProblem(smooth_sys, Pair[int.x=>0.0], (0.0, 10.0))
    smooth_sol = solve(smooth_prob, Rodas4(), dtmax=0.1) # set dtmax to prevent the solver from overstepping the entire step disturbance)

    @test smooth_sol.retcode == :Success
    @test smooth_sol[smooth_src.output.u] ≈ smooth_step.(smooth_sol.t, δ, height, offset, start_time)

    # with duration
    @named smooth_src = Step(offset=offset, height=height, start_time=start_time, smooth=true, duration=duration)
    @named smooth_iosys = ODESystem([
        connect(smooth_src.output, int.input),
        ],
        t,
        systems=[int, smooth_src],
    )

    smooth_sys = structural_simplify(smooth_iosys)
    smooth_prob = ODEProblem(smooth_sys, Pair[int.x=>0.0], (0.0, 10.0))
    smooth_sol = solve(smooth_prob, Rodas4())

    @test smooth_sol.retcode == :Success
    @test smooth_sol[smooth_src.output.u] ≈ smooth_step.(smooth_sol.t, δ, height, offset, start_time) - smooth_step.(smooth_sol.t, δ, height, 0, start_time+ duration)
end

@testset "Square" begin
    frequency=1
    amplitude=2
    offset=1
    start_time=2.5
    δ = 1e-5
    @named int = Integrator()

    @named src = Square(frequency=frequency, amplitude=amplitude,
        offset=offset, start_time=start_time)
    @named iosys = ODESystem([
        connect(src.output, int.input),
        ],
        t,
        systems=[int, src],
    )

    sys = structural_simplify(iosys)
    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0))
    sol = solve(prob, Rodas4())

    @test sol.retcode == :Success
    @test sol[src.output.u] ≈ square.(sol.t, frequency, amplitude, offset, start_time) atol=1e-3

    @named smooth_src = Square(frequency=frequency, amplitude=amplitude,
    offset=offset, start_time=start_time, smooth=true)
    @named smooth_iosys = ODESystem([
        connect(smooth_src.output, int.input),
        ],
        t,
        systems=[int, smooth_src]
    )

    smooth_sys = structural_simplify(smooth_iosys)
    smooth_prob = ODEProblem(smooth_sys, Pair[int.x=>0.0], (0.0, 10.0))
    smooth_sol = solve(smooth_prob, Rodas4())

    @test smooth_sol.retcode == :Success
    @test smooth_sol[smooth_src.output.u] ≈ smooth_square.(smooth_sol.t, δ, frequency, amplitude, offset, start_time) atol=1e-3
end

@testset "Triangular" begin
    frequency=5
    amplitude=1
    offset = 2
    start_time = 1
    δ = 1e-5
    @named int = Integrator()

    @named src = Triangular(frequency=frequency, amplitude=amplitude,
        offset=offset, start_time=start_time)
    @named iosys = ODESystem([
        connect(src.output, int.input),
        ],
        t,
        systems=[int, src],
    )

    sys = structural_simplify(iosys)
    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 4.0))
    sol = solve(prob, Rodas4(), saveat=0.01)

    @test sol.retcode == :Success
    @test sol[src.output.u] ≈ triangular.(sol.t, frequency, amplitude, offset, start_time) atol=1e-3

    @named smooth_src = Triangular(frequency=frequency, amplitude=amplitude,
        offset=offset, start_time=start_time, smooth=true)
    @named smooth_iosys = ODESystem([
        connect(smooth_src.output, int.input),
        ],
        t,
        systems=[int, smooth_src],
    )

    smooth_sys = structural_simplify(smooth_iosys)
    smooth_prob = ODEProblem(smooth_sys, Pair[int.x=>0.0], (0.0, 4.0))
    smooth_sol = solve(smooth_prob, Rodas4(), saveat=0.01)

    @test smooth_sol.retcode == :Success
    @test smooth_sol[smooth_src.output.u] ≈ smooth_triangular.(smooth_sol.t, δ, frequency, amplitude, offset, start_time) atol=1e-3
end

@testset "ExpSine" begin
    exp_sine(t, amplitude, frequency, damping, phase, start_time) = offset + ifelse(t < start_time, 0, amplitude * exp(-damping * (t - start_time)) * sin(2*pi*frequency*(t - start_time) + phase))

    frequency, amplitude, damping = 3, 2, 0.10
    phase, offset, start_time, δ =  0, 0, 0, 1e-5
    @named src = ExpSine(frequency=frequency, amplitude=amplitude, damping=damping, phase=phase, offset=offset, start_time=start_time)
    @named int = Integrator()
    @named iosys = ODESystem([
        connect(src.output, int.input),
        ],
        t,
        systems=[int, src],
    )
    sys = structural_simplify(iosys)
    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0))
    sol = solve(prob, Rodas4())

    @test sol.retcode == :Success
    @test sol[src.output.u] ≈ exp_sine.(sol.t, amplitude, frequency, damping, phase, start_time) atol=1e-3

    offset, start_time = 1, 2
    @named smooth_src = ExpSine(frequency=frequency, amplitude=amplitude, damping=damping, phase=phase, offset=offset, start_time=start_time, smooth=true)
    @named smooth_iosys = ODESystem([
        connect(smooth_src.output, int.input),
        ],
        t,
        systems=[int, smooth_src],
    )
    smooth_sys = structural_simplify(smooth_iosys)
    smooth_prob = ODEProblem(smooth_sys, Pair[int.x=>0.0], (0.0, 10.0))
    smooth_sol = solve(smooth_prob, Rodas4())

    @test smooth_sol.retcode == :Success
    @test smooth_sol[smooth_src.output.u] ≈ smooth_damped_sin.(smooth_sol.t, δ, frequency, amplitude, damping, phase, offset, start_time) atol=1e-3
end