text:
"createHDF5: Creates or opens HDF5 file."
function createHDF5(sol::Solution, file::AbstractString; overwrite=true)
    if overwrite
        flag = "w"
        get_config(:verbosity) > 1 ? @info("Creating HDF5 file ", file) : nothing
        isfile(file) ? @warn("Overwriting existing HDF5 file.") : nothing
    else
        flag = "r+"
        get_config(:verbosity) > 1 ? @info("Opening HDF5 file ", file) : nothing
    end

    # create or open HDF5 file
    h5 = h5open(file, flag)

    return h5
end

"save_attributes: Saves common attributes of Solution to HDF5 file."
function save_attributes(sol::Solution)
    save_attributes(sol, hdf5(sol))
end

"save_attributes: Saves attributes of Deterministic Solutions to HDF5 file."
function save_attributes(sol::DeterministicSolution, h5::HDF5File)
    # save attributes
    attrs(h5)["ntime"] = ntime(sol)
    attrs(h5)["nsave"] = nsave(sol)
end

"save_attributes: Saves attributes of Stochastic Solutions to HDF5 file."
function save_attributes(sol::StochasticSolution, h5::HDF5File)
    attrs(h5)["ntime"] = ntime(sol)
    attrs(h5)["nsave"] = nsave(sol)
    attrs(h5)["conv"]  = string(conv(sol))
    attrs(h5)["nd"] = sol.nd
    attrs(h5)["nm"] = sol.nm
    attrs(h5)["ns"] = sol.ns
    attrs(h5)["K"]  = sol.K
end

"write_to_hdf5: Wrapper for saving Solution to HDF5 file."
function CommonFunctions.write_to_hdf5(solution::Solution)
    write_to_hdf5(solution, hdf5(solution), offset(solution))
end

"Creates HDF5 file, writes solution to file, and closes file."
function CommonFunctions.write_to_hdf5(solution::Solution, file::AbstractString)
    h5 = create_hdf5(solution, file)
    write_to_hdf5(solution, h5)
    close(h5)
end

"Creates HDF5 file and initialises datasets for solution object."
function create_hdf5!(solution::Solution, file::AbstractString; kwargs...)
    solution.h5 = create_hdf5(solution, file; kwargs...)
end

"Creates HDF5 file and initialises datasets for deterministic solution object."
function create_hdf5(solution::DeterministicSolution, file::AbstractString)
    # create HDF5 file and save attributes and common parameters
    h5 = createHDF5(solution, file)

    # save attributes
    save_attributes(solution, h5)

    # create dataset
    init_timeteps_in_hdf5(solution, h5)
    init_solution_in_hdf5(solution, h5)

    return h5
end

"Creates HDF5 file and initialises datasets for stochastic solution object."
function create_hdf5(solution::StochasticSolution, file::AbstractString; save_W=true)
    # create HDF5 file and save attributes and common parameters
    h5 = createHDF5(solution, file)

    # save attributes
    save_attributes(solution, h5)

    # create dataset
    init_timeteps_in_hdf5(solution, h5)
    init_solution_in_hdf5(solution, h5)

    if save_W
        init_increments_in_hdf5(solution, h5)
    end

    return h5
end

function init_timeteps_in_hdf5(solution::Solution{DT,TT}, h5::HDF5File) where {DT,TT}
    t = d_create(h5, "t", TT, ((solution.nt+1,), (-1,)), "chunk", (1,))
    t[1] = timesteps(solution)[0]
end

function init_solution_in_hdf5(solution::SolutionODE{DT,TT,2}, h5::HDF5File) where {DT,TT}
    q = d_create(h5, "q", DT, ((solution.nd, solution.nt+1), (solution.nd, -1)), "chunk", (solution.nd,1))
    q[:, 1] = solution.q[:, 0]
end

function init_solution_in_hdf5(solution::SolutionDAE{DT,TT,2}, h5::HDF5File) where {DT,TT}
    q = d_create(h5, "q", DT, ((solution.nd, solution.nt+1), (solution.nd, -1)), "chunk", (solution.nd,1))
    λ = d_create(h5, "λ", DT, ((solution.nm, solution.nt+1), (solution.nm, -1)), "chunk", (solution.nm,1))
    q[:, 1] = solution.q[:, 0]
    λ[:, 1] = solution.λ[:, 0]
end

function init_solution_in_hdf5(solution::SolutionSDE{DT,TT,2}, h5::HDF5File) where {DT,TT}
    q = d_create(h5, "q", DT, ((solution.nd, solution.nt+1), (solution.nd, -1)), "chunk", (solution.nd,1))
    q[:,1] = solution.q[:,0]
end

function init_solution_in_hdf5(solution::SolutionODE{DT,TT,3}, h5::HDF5File) where {DT,TT}
    q = d_create(h5, "q", DT, ((solution.nd, solution.nt+1, solution.ni),(solution.nd, -1, solution.ni)), "chunk", (solution.nd,1,1))
    q[:, 1, :] = solution.q[:, 0, :]
end

function init_solution_in_hdf5(solution::SolutionDAE{DT,TT,3}, h5::HDF5File) where {DT,TT}
    q = d_create(h5, "q", DT, ((solution.nd, solution.nt+1, solution.ni),(solution.nd, -1, solution.ni)), "chunk", (solution.nd,1,1))
    λ = d_create(h5, "λ", DT, ((solution.nm, solution.nt+1, solution.ni),(solution.nm, -1, solution.ni)), "chunk", (solution.nm,1,1))
    q[:, 1, :] = solution.q[:, 0, :]
    λ[:, 1, :] = solution.λ[:, 0, :]
end

function init_solution_in_hdf5(solution::SolutionSDE{DT,TT,3}, h5::HDF5File) where {DT,TT}
    q = d_create(h5, "q", DT, ((solution.nd, solution.nt+1, solution.ns),(solution.nd, -1, solution.ns)), "chunk", (solution.nd,1,1))
    q[:,1,:] = solution.q[:,0,:]
end

function init_solution_in_hdf5(solution::SolutionPODE{DT,TT,2}, h5::HDF5File) where {DT,TT}
    q = d_create(h5, "q", DT, ((solution.nd, solution.nt+1), (solution.nd, -1)), "chunk", (solution.nd,1))
    p = d_create(h5, "p", DT, ((solution.nd, solution.nt+1), (solution.nd, -1)), "chunk", (solution.nd,1))
    q[:, 1] = solution.q[:, 0]
    p[:, 1] = solution.p[:, 0]
end

function init_solution_in_hdf5(solution::SolutionPDAE{DT,TT,2}, h5::HDF5File) where {DT,TT}
    q = d_create(h5, "q", DT, ((solution.nd, solution.nt+1), (solution.nd, -1)), "chunk", (solution.nd,1))
    p = d_create(h5, "p", DT, ((solution.nd, solution.nt+1), (solution.nd, -1)), "chunk", (solution.nd,1))
    λ = d_create(h5, "λ", DT, ((solution.nm, solution.nt+1), (solution.nm, -1)), "chunk", (solution.nm,1))
    q[:, 1] = solution.q[:, 0]
    p[:, 1] = solution.p[:, 0]
    λ[:, 1] = solution.λ[:, 0]
end

function init_solution_in_hdf5(solution::SolutionPSDE{DT,TT,2}, h5::HDF5File) where {DT,TT}
    q = d_create(h5, "q", DT, ((solution.nd, solution.nt+1), (solution.nd, -1)), "chunk", (solution.nd,1))
    p = d_create(h5, "p", DT, ((solution.nd, solution.nt+1), (solution.nd, -1)), "chunk", (solution.nd,1))
    q[:,1] = solution.q[:,0]
    p[:, 1] = solution.p[:, 0]
end

function init_solution_in_hdf5(solution::SolutionPODE{DT,TT,3}, h5::HDF5File) where {DT,TT}
    q = d_create(h5, "q", DT, ((solution.nd, solution.nt+1, solution.ni),(solution.nd, -1, solution.ni)), "chunk", (solution.nd,1,1))
    p = d_create(h5, "p", DT, ((solution.nd, solution.nt+1, solution.ni),(solution.nd, -1, solution.ni)), "chunk", (solution.nd,1,1))
    q[:, 1, :] = solution.q[:, 0, :]
    p[:, 1, :] = solution.p[:, 0, :]
end

function init_solution_in_hdf5(solution::SolutionPDAE{DT,TT,3}, h5::HDF5File) where {DT,TT}
    q = d_create(h5, "q", DT, ((solution.nd, solution.nt+1, solution.ni),(solution.nd, -1, solution.ni)), "chunk", (solution.nd,1,1))
    p = d_create(h5, "p", DT, ((solution.nd, solution.nt+1, solution.ni),(solution.nd, -1, solution.ni)), "chunk", (solution.nd,1,1))
    λ = d_create(h5, "λ", DT, ((solution.nm, solution.nt+1, solution.ni),(solution.nm, -1, solution.ni)), "chunk", (solution.nm,1,1))
    q[:, 1, :] = solution.q[:, 0, :]
    p[:, 1, :] = solution.p[:, 0, :]
    λ[:, 1, :] = solution.λ[:, 0, :]
end

function init_solution_in_hdf5(solution::SolutionPSDE{DT,TT,3}, h5::HDF5File) where {DT,TT}
    q = d_create(h5, "q", DT, ((solution.nd, solution.nt+1, solution.ns),(solution.nd, -1, solution.ns)), "chunk", (solution.nd,1,1))
    p = d_create(h5, "p", DT, ((solution.nd, solution.nt+1, solution.ns),(solution.nd, -1, solution.ns)), "chunk", (solution.nd,1,1))
    q[:,1,:] = solution.q[:,0,:]
    p[:, 1, :] = solution.p[:, 0, :]
end

function init_increments_in_hdf5(solution::StochasticSolution{DT,TT,NQ,2}, h5::HDF5File) where {DT,TT,NQ}
    dW = d_create(h5, "ΔW", DT, ((solution.nm, solution.ntime),(solution.nm, -1)), "chunk", (solution.nm,1))
    dZ = d_create(h5, "ΔZ", DT, ((solution.nm, solution.ntime),(solution.nm, -1)), "chunk", (solution.nm,1))
end

function init_increments_in_hdf5(solution::StochasticSolution{DT,TT,NQ,3}, h5::HDF5File) where {DT,TT,NQ}
    dW = d_create(h5, "ΔW", DT, ((solution.nm, solution.ntime, solution.ns),(solution.nm, -1, solution.ns)), "chunk", (solution.nm,1,1))
    dZ = d_create(h5, "ΔZ", DT, ((solution.nm, solution.ntime, solution.ns),(solution.nm, -1, solution.ns)), "chunk", (solution.nm,1,1))
end

function copy_timeteps_to_hdf5(sol::Solution, h5::HDF5File, j1, j2, n1, n2)
    if size(h5["t"],1) < j2
        set_dims!(h5["t"], (j2,))
    end
    h5["t"][j1:j2] = timesteps(sol)[n1:n2]
end

function copy_solution_to_hdf5(solution::Union{SolutionODE{DT,TT,2},SolutionDAE{DT,TT,2},SolutionSDE{DT,TT,2}}, h5::HDF5File, j1, j2, n1, n2) where {DT,TT}
    if size(h5["q"],2) < j2
        set_dims!(h5["q"], (size(h5["q"],1), j2))
    end
    h5["q"][:, j1:j2] = solution.q[:, n1:n2]
end

function copy_solution_to_hdf5(solution::Union{SolutionODE{DT,TT,3},SolutionDAE{DT,TT,3},SolutionSDE{DT,TT,3}}, h5::HDF5File, j1, j2, n1, n2) where {DT,TT}
    if size(h5["q"],2) < j2
        set_dims!(h5["q"], (size(h5["q"],1), j2, size(h5["q"],3)))
    end
    h5["q"][:, j1:j2, :] = solution.q[:, n1:n2, :]
end

function copy_solution_to_hdf5(solution::Union{SolutionPODE{DT,TT,2},SolutionPDAE{DT,TT,2},SolutionPSDE{DT,TT,2}}, h5::HDF5File, j1, j2, n1, n2) where {DT,TT}
    if size(h5["q"],2) < j2
        set_dims!(h5["q"], (size(h5["q"],1), j2))
    end
    if size(h5["p"],2) < j2
        set_dims!(h5["p"], (size(h5["p"],1), j2))
    end
    h5["q"][:, j1:j2] = solution.q[:, n1:n2]
    h5["p"][:, j1:j2] = solution.p[:, n1:n2]
end

function copy_solution_to_hdf5(solution::Union{SolutionPODE{DT,TT,3},SolutionPDAE{DT,TT,3},SolutionPSDE{DT,TT,3}}, h5::HDF5File, j1, j2, n1, n2) where {DT,TT}
    if size(h5["q"],2) < j2
        set_dims!(h5["q"], (size(h5["q"],1), j2, size(h5["q"],3)))
    end
    if size(h5["p"],2) < j2
        set_dims!(h5["p"], (size(h5["p"],1), j2, size(h5["p"],3)))
    end
    h5["q"][:, j1:j2, :] = solution.q[:, n1:n2, :]
    h5["p"][:, j1:j2, :] = solution.p[:, n1:n2, :]
end

function copy_multiplier_to_hdf5(solution::Union{SolutionDAE{DT,TT,2},SolutionPDAE{DT,TT,2}}, h5::HDF5File, j1, j2, n1, n2) where {DT,TT}
    if size(h5["λ"],2) < j2
        set_dims!(h5["λ"], (size(h5["λ"],1), j2))
    end
    h5["λ"][:, j1:j2] = solution.λ[:, n1:n2]
end

function copy_multiplier_to_hdf5(solution::Union{SolutionDAE{DT,TT,3},SolutionPDAE{DT,TT,3}}, h5::HDF5File, j1, j2, n1, n2) where {DT,TT}
    if size(h5["λ"],2) < j2
        set_dims!(h5["λ"], (size(h5["λ"],1), j2, size(h5["λ"],3)))
    end
    h5["λ"][:, j1:j2, :] = solution.λ[:, n1:n2, :]
end

function copy_increments_to_hdf5(solution::StochasticSolution{DT,TT,NQ,2}, h5::HDF5File, j1, j2, n1, n2) where {DT,TT,NQ}
    if size(h5["ΔW"],2) < j2
        set_dims!(h5["ΔW"], (size(h5["ΔW"],1), j2))
    end
    if size(h5["ΔZ"],2) < j2
        set_dims!(h5["ΔZ"], (size(h5["ΔZ"],1), j2))
    end
    h5["ΔW"][:,j1:j2] = solution.W.ΔW[:,n1:n2]
    h5["ΔZ"][:,j1:j2] = solution.W.ΔZ[:,n1:n2]
end

function copy_increments_to_hdf5(solution::StochasticSolution{DT,TT,NQ,3}, h5::HDF5File, j1, j2, n1, n2) where {DT,TT,NQ}
    if size(h5["ΔW"],2) < j2
        set_dims!(h5["ΔW"], (size(h5["ΔW"],1), j2, size(h5["ΔW"],3)))
    end
    if size(h5["ΔZ"],2) < j2
        set_dims!(h5["ΔZ"], (size(h5["ΔZ"],1), j2, size(h5["ΔZ"],3)))
    end
    h5["ΔW"][:,j1:j2,:] = solution.W.ΔW[:,n1:n2,:]
    h5["ΔZ"][:,j1:j2,:] = solution.W.ΔZ[:,n1:n2,:]
end

"Append solution to HDF5 file."
function CommonFunctions.write_to_hdf5(solution::DeterministicSolution, h5::HDF5File=hdf5(solution), offset=offset(solution))
    # set convenience variables and compute ranges
    j1 = offset+2
    j2 = offset+1+solution.nt

    # copy data from solution to HDF5 dataset
    copy_timeteps_to_hdf5(solution, h5, j1, j2, 1, solution.nt)
    copy_solution_to_hdf5(solution, h5, j1, j2, 1, solution.nt)

    return nothing
end

"Append solution to HDF5 file."
function CommonFunctions.write_to_hdf5(solution::Union{SolutionDAE,SolutionPDAE}, h5::HDF5File=hdf5(solution), offset=offset(solution))
    # set convenience variables and compute ranges
    j1 = offset+2
    j2 = offset+1+solution.nt

    # copy data from solution to HDF5 dataset
    copy_timeteps_to_hdf5(solution, h5, j1, j2, 1, solution.nt)
    copy_solution_to_hdf5(solution, h5, j1, j2, 1, solution.nt)
    copy_multiplier_to_hdf5(solution, h5, j1, j2, 1, solution.nt)

    return nothing
end
"""
Append solution to HDF5 file.
  soffset - start writing the solution q at the position soffset+2
  woffset - start writing the increments ΔW, ΔZ at the position woffset+1
"""
function CommonFunctions.write_to_hdf5(solution::StochasticSolution, h5::HDF5File=hdf5(solution), soffset=offset(solution), woffset=ioffset(solution))
    # set convenience variables and compute ranges
    js1 = soffset+2
    js2 = soffset+1+solution.nt
    jw1 = woffset+1
    jw2 = woffset+solution.nwrite

    # copy data from solution to HDF5 dataset
    copy_timeteps_to_hdf5(solution, h5, js1, js2, 1, solution.nt)
    copy_solution_to_hdf5(solution, h5, js1, js2, 1, solution.nt)

    if exists(h5, "ΔW") && exists(h5, "ΔZ")
        # copy the Wiener process increments from solution to HDF5 dataset
        copy_increments_to_hdf5(solution, h5, jw1, jw2, 1, solution.nwrite)
    end

    return nothing
end


text:
using DeterminantalPointProcesses
using Test
using Random, LinearAlgebra, IterTools
using KernelFunctions
import Combinatorics: combinations

rng = MersenneTwister(42)

n = 5
k = 2
A = rand(rng, n, n)
L = Symmetric(A * A')
kernel = SqExponentialKernel()
X = [rand(2) for _ in 1:3]
dpp = DPP(L)

@testset "Ensure correct pdf and logpdf computation" begin
    @testset "Kernel API" begin
        K = kernelmatrix(kernel, X)
        @test DPP(K).L ≈ DPP(kernel, X).L
        @test DPP(K).L ≈ DPP(kernel, reduce(hcat, X); obsdim=2).L
    end
    @testset "For DPP" begin
        # compute the true distribution
        true_pdf = Float64[]
        true_logpdf = Float64[]
        for z in subsets(1:n)
            push!(true_pdf, pdf(dpp, z))
            push!(true_logpdf, logpdf(dpp, z))
        end
        @test sum(true_pdf) ≈ 1.0
        @test all(true_pdf .<= 1.0)
        @test all(true_logpdf .<= 0.0)
    end

    @testset "For k-DPP" begin
        # compute the true distribution
        true_pdf = Float64[]
        true_logpdf = Float64[]
        for z in combinations(1:n, k)
            push!(true_pdf, pdf(dpp(k), z))
            push!(true_logpdf, logpdf(dpp(k), z))
        end
        @test sum(true_pdf) ≈ 1.0
        @test all(true_pdf .<= 1.0)
        @test all(true_logpdf .<= 0.0)
    end
end

@testset "Ensure correct sampling from DPP" begin
    # compute the true distribution
    all_subsets = []
    true_pdf = Float64[]
    for (i, z) in enumerate(subsets(1:n))
        push!(true_pdf, pdf(dpp, z))
        push!(all_subsets, (z, i))
    end
    global subset_to_idx = Dict(all_subsets)

    @testset "Exact sampling" begin
        nb_samples = 1000
        samples = rand(dpp, nb_samples)

        # compute the empirical distribution
        empirical_pdf = zeros(Float64, length(true_pdf))
        for z in samples
            empirical_pdf[subset_to_idx[z]] += 1
        end
        empirical_pdf ./= nb_samples

        # ensure that the empirical pdf is close to the true pdf
        total_variation = maximum(abs.(true_pdf - empirical_pdf))
        @test total_variation ≈ 0.0 atol = 1e-1
    end

    @testset "MCMC sampling" begin
        nb_samples = 1000
        samples, state = randmcmc(dpp, nb_samples; return_final_state=true)

        # ensure that L_z_inv makes sense (i.e., noise did not accumulate)
        z, L_z_inv = state
        @test sum(L_z_inv[z, z] * dpp.L[z, z] - I) ≈ 0 atol=1e-10

        # compute the empirical distribution
        empirical_pdf = zeros(Float64, length(true_pdf))
        for z in samples
            empirical_pdf[subset_to_idx[z]] += 1
        end
        empirical_pdf ./= nb_samples

        # ensure that the empirical pdf is close to the true pdf
        total_variation = maximum(abs.(true_pdf - empirical_pdf))
        @test total_variation ≈ 0.0 atol = 1e-1
    end
end

@testset "Ensure correct sampling from k-DPP" begin
    # compute the true distribution
    all_k_subsets = []
    true_pdf = Float64[]
    for (i, z) in enumerate(combinations(1:n, k))
        push!(true_pdf, pdf(dpp, z))
        push!(all_k_subsets, (z, i))
    end
    true_pdf ./= sum(true_pdf)
    k_subset_to_idx = Dict(all_k_subsets)

    @testset "Exact sampling" begin
        nb_samples = 10000
        samples = rand(dpp(k), nb_samples)

        # ensure that samples are of proper cardinality
        @test all(map(length, samples) .== k)

        # compute the empirical distribution
        empirical_pdf = zeros(Float64, length(true_pdf))
        for z in samples
            empirical_pdf[k_subset_to_idx[z]] += 1
        end
        empirical_pdf ./= nb_samples

        # ensure that the empirical pdf is close to the true pdf
        total_variation = maximum(abs.(true_pdf - empirical_pdf))
        @test total_variation ≈ 0.0 atol = 1e-1
    end

    @testset "MCMC sampling" begin
        nb_samples = 1000
        samples, state = randmcmc(dpp(k), nb_samples; return_final_state=true)

        # ensure that L_z_inv makes sense (i.e., noise did not accumulate)
        z, L_z_inv = state
        @test sum(L_z_inv[z, z] * dpp.L[z, z] - I) ≈ 0 atol = 1e-1

        # ensure that samples are of proper cardinality
        @test all(map(length, samples) .== k)

        # compute the empirical distribution
        empirical_pdf = zeros(Float64, length(true_pdf))
        for z in samples
            empirical_pdf[k_subset_to_idx[z]] += 1
        end
        empirical_pdf ./= nb_samples

        # ensure that the empirical pdf is close to the true pdf
        total_variation = maximum(abs.(true_pdf - empirical_pdf))
        @test_broken total_variation ≈ 0.0
        atol = 1e-1
    end
end

ext: jl


text:
module TestEnvironment
    import Benchmarks
    using Base.Test
    using Compat

    e = Benchmarks.Environment()

    io = IOBuffer()
    show(io, e)

    path = tempname()
    writecsv(path, e)
    bytes = readstring(path)
    rm(path)
end


text:
type FdG2SwaptionEngine{Y <: TermStructure} <: PricingEngine
  model::G2
  tGrid::Int
  xGrid::Int
  yGrid::Int
  dampingSteps::Int
  invEps::Float64
  schemeDesc::FdmSchemeDesc
  ts::Y
end

FdG2SwaptionEngine(model::G2, tGrid::Int = 100, xGrid::Int = 50, yGrid::Int = 50, dampingSteps::Int = 0, invEps::Float64 = 1e-5,
                  schemeDesc::FdmSchemeDesc = FdmSchemeDesc(Hundsdorfer())) =
                  FdG2SwaptionEngine(model, tGrid, xGrid, yGrid, dampingSteps, invEps, schemeDesc, model.ts)

# methods #
function _calculate!(pe::FdG2SwaptionEngine, swaption::Swaption)
  # 1. Term structure
  ts = pe.ts

  # 2. Mesher
  dc = ts.dc
  refDate = reference_date(ts)
  maturity = year_fraction(dc, refDate, swaption.exercise.dates[end])

  process1 = OrnsteinUhlenbeckProcess(get_a(pe.model), get_sigma(pe.model))
  process2 = OrnsteinUhlenbeckProcess(get_b(pe.model), get_eta(pe.model))

  xMesher = FdmSimpleProcess1dMesher(pe.xGrid, process1, maturity, 1, pe.invEps)
  yMesher = FdmSimpleProcess1dMesher(pe.yGrid, process2, maturity, 1, pe.invEps)
  mesher = FdmMesherComposite(xMesher, yMesher)

  # 3. Inner Value calculator
  exerciseDates = swaption.exercise.dates
  t2d = Dict{Float64, Date}()
  for i = 1:length(exerciseDates)
    t = year_fraction(dc, refDate, exerciseDates[i])
    t2d[t] = exerciseDates[i]
  end

  disTs = pe.model.ts
  fwdTs = swaption.swap.iborIndex.ts

  fwdModel = G2(fwdTs, get_a(pe.model), get_sigma(pe.model), get_b(pe.model), get_eta(pe.model), get_rho(pe.model))
  calculator = FdmAffineModelSwapInnerValue(pe.model, fwdModel, swaption.swap, t2d, mesher, 1)

  # 4. Step Conditions
  conditions = vanilla_FdmStepConditionComposite(DividendSchedule(), swaption.exercise, mesher, calculator, refDate, dc)

  # 5. Boundary conditions
  boundaries = FdmBoundaryConditionSet()

  # 6. Solver
  solverDesc = FdmSolverDesc(mesher, boundaries, conditions, calculator, maturity, pe.tGrid, pe.dampingSteps)
  solver = FdmG2Solver(pe.model, solverDesc, pe.schemeDesc)
  swaption.results.value = value_at(solver, 0.0, 0.0)
end


text:
include("./model/modelMain.jl")
include("./heuristics/heuristicsMain.jl")

include("./relaxed/relaxedMain.jl")
include("./restricted/restrictedMain.jl")
include("./branchAndBound/branchAndBoundMain.jl")

"""
    checkForSolution(dd)

A debug method that checks a dd for a hard coded solution in the form of Vector{Int}

# Arguments
- `dd::Vector{Vector{RelaxedSequencingNode}}`: the dd to check
"""
function checkForSolution(dd::Vector{Vector{RelaxedSequencingNode}})
    solution = [1, 5, 2, 10, 3, 6, 4, 7, 11, 8, 9, 12, 13]
    solLength = length(dd[1][1])
    startingNode = nothing

    for node in dd[1]
        if node.state == solution[solLength]
            secondCheck = true
            for e in 1:solLength
                if !(solution[e] in node.someDown)
                    secondCheck = false
                    break
                end
            end
            if secondCheck
                startingNode = node
                break
            end
        end
    end
    if isnothing(startingNode)
        return false
    end

    for i in solLength+1:length(solution)
        nextNum = solution[i]

        for arc in startingNode.arcsOut
            if arc.assignment==nextNum
                startingNode = arc.destination
                break
            end
        end

        if length(startingNode) != i
            println("_______________check for solution_____________")
            println(startingNode)
            println("_______________checked for solution_____________")
            return false
        end
    end

    return true
end


text:
"""
Annealed Gibbs sampling, to remove small amplitude events.
"""
function annealed_gibbs!(
        model::Union{SeqModel,DistributedSeqModel},
        spikes::Vector{Spike},
        initial_assignments::Vector{Int64},
        num_anneals::Int64,
        samples_per_anneal::Int64,
        max_temperature::Float64,
        extra_split_merge_moves::Int64,
        split_merge_window::Float64,
        save_every::Int64;
        verbose::Bool=false
    )

    # Initialize storage.
    assignment_hist = zeros(Int64, length(spikes), 0)
    log_p_hist = Float64[]
    latent_event_hist = Vector{EventSummaryInfo}[]
    globals_hist = SeqGlobals[]

    # Return early if no samples.
    if num_anneals == 0
        return (
            initial_assignments,
            assignment_hist,
            log_p_hist,
            latent_event_hist,
            globals_hist
        )
    end

    # Final amplitude for anneal.
    target_mean = mean(priors(model).seq_event_amplitude)
    target_var = var(priors(model).seq_event_amplitude)

    # Begin annealing.
    temperatures = exp10.(range(log10(max_temperature), 0, length=num_anneals))
    assignments = initial_assignments

    for temp in temperatures

        # Print progress.
        verbose && println("TEMP:  ", temp)
        flush(stdout)

        # Anneal prior on sequence amplitude.
        prior = priors(model)
        set_priors!(
            model,
            SeqPriors(
                prior.seq_event_rate,
                prior.seq_type_proportions,
                specify_gamma(target_mean, target_var * temp), # Anneal prior.
                prior.neuron_response_proportions,
                prior.neuron_response_profile,
                prior.bkgd_amplitude,
                prior.bkgd_proportions,
                prior.warp_values,
                prior.warp_log_proportions
            )
        )

        # Recompute probability of introducing a new cluster.
        #  ==> TODO: maybe set_priors! should do this automatically?
        max_time = get_max_time(model)
        prior = priors(model)
        α = prior.seq_event_amplitude.α
        β = prior.seq_event_amplitude.β
        λ = prior.seq_event_rate
        set_new_cluster_log_prob!(
            model,
            log(α) + log(λ) + log(max_time) + α * (log(β) - log(1 + β))
        )
        # TODO fix this, can't use primary model in not distributed code
        set_bkgd_log_prob!(
            model,
            (log(get_bkgd_amplitude(model))
            + log(max_time)
            + log(1 + β))
        )

        # Draw gibbs samples.
        (
            assignments,
            _assgn,
            _logp,
            _latents,
            _globals
        ) =
        gibbs_sample!(
            model,
            spikes,
            assignments,
            samples_per_anneal,
            extra_split_merge_moves,
            split_merge_window,
            save_every,
            verbose=verbose
        )

        # Save samples.
        assignment_hist = cat(assignment_hist, _assgn, dims=2)
        append!(log_p_hist, _logp)
        append!(latent_event_hist, _latents)
        append!(globals_hist, _globals)

    end

    return (
        assignments,
        assignment_hist,
        log_p_hist,
        latent_event_hist,
        globals_hist
    )

end


text:
using BinaryBuilder, Pkg

name = "SSGraphBLAS"
version = v"4.0.3"

# Collection of sources required to build SuiteSparse:GraphBLAS
sources = [
    ArchiveSource("https://github.com/DrTimothyAldenDavis/GraphBLAS/archive/v4.0.3.tar.gz",
        "43519783625f1a0a631158603850cfcf0d9681646dbb5c64ae2eaf27e1444b90")
]

# Bash recipe for building across all platforms
script = raw"""
# Compile GraphBLAS
cd $WORKSPACE/srcdir/GraphBLAS-*/build
cmake .. -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_BUILD_TYPE=Release -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DUSER_NONE=1

make -j${nproc} install

if [[ ! -f "${libdir}/libgraphblas.${dlext}" ]]; then
    # For mysterious reasons, the shared library is not installed
    # when building for Windows
    mkdir -p "${libdir}"
    cp "libgraphblas.${dlext}" "${libdir}"
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libgraphblas", :libgraphblas),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"6")

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
function secret_handshake(code)
    handshake = []
    code & 1 != 0 && (handshake = ["wink"])
    code & 2 != 0 && (push!(handshake, "double blink"))
    code & 4 != 0 && (push!(handshake, "close your eyes"))
    code & 8 != 0 && (push!(handshake, "jump"))
    code & 16 != 0 && (reverse!(handshake))
    return handshake

end

max_stars_repo_name: Rigante-pl/Exercism-Julia

max_issues_repo_name: Rigante-pl/Exercism-Julia

max_forks_repo_name: Rigante-pl/Exercism-Julia


text:
"""
    struct QCurlGradMonomialBasis{...} <: Field

This type implements a multivariate vector-valued polynomial basis
spanning the space needed for Raviart-Thomas reference elements on n-cubes.
The type parameters and fields of this `struct` are not public.
This type fully implements the [`Field`](@ref) interface, with up to first order
derivatives.
"""
struct QCurlGradMonomialBasis{D,T} <: Field
  qgrad::QGradMonomialBasis{D,T}
  function QCurlGradMonomialBasis(::Type{T},order::Int,terms::CartesianIndices{D},perms::Matrix{Int}) where {D,T}
    qgrad = QGradMonomialBasis(T,order,terms,perms)
    new{D,T}(qgrad)
  end
end

"""
    QCurlGradMonomialBasis{D}(::Type{T},order::Int) where {D,T}

Returns a `QCurlGradMonomialBasis` object. `D` is the dimension
of the coordinate space and `T` is the type of the components in the vector-value.
"""
function QCurlGradMonomialBasis{D}(::Type{T},order::Int) where {D,T}
  @assert T<:Real "T needs to be <:Real since represents the type of the components of the vector value"
  _t = tfill(order,Val{D-1}())
  t = (order+1,_t...)
  terms = CartesianIndices(t)
  perms = _prepare_perms(D)
  QCurlGradMonomialBasis(T,order,terms,perms)
end

function field_cache(f::QCurlGradMonomialBasis,x)
  field_cache(f.qgrad,x)
end

@inline function evaluate_field!(cache,f::QCurlGradMonomialBasis,x)
  evaluate_field!(cache,f.qgrad,x)
end

function gradient_cache(f::QCurlGradMonomialBasis,x)
  gradient_cache(f.qgrad,x)
end

@inline function evaluate_gradient!(cache,f::QCurlGradMonomialBasis,x)
  evaluate_gradient!(cache,f.qgrad,x)
end


text:
mutable struct SingleOCFUpdater <: Updater
    pomdp::SingleOCFPOMDP
end

const SingleOCFBelief = SparseCat{Vector{SingleOCFState},Vector{Float64}}

function POMDPs.update(pomdp::SingleOCFPOMDP, up::SingleOCFUpdater, bold::Dict{Int64, SingleOCFBelief}, a::SingleOCFAction, o::Dict{Int64, SingleOCFObs})
    bnew = Dict{Int64, SingleOCFBelief}()

    for oid in keys(o)

        if haskey(bold, oid) && oid != PEDESTRIAN_OFF_KEY  # old measurment
           bnew[oid] = update(pomdp, up, bold[oid], a, o[oid])

        elseif oid == PEDESTRIAN_OFF_KEY  # absent state
            bnew[PEDESTRIAN_OFF_KEY] = initBeliefAbsentPedestrian(pomdp, o[oid].ego_y, o[oid].ego_v)
           # bnew[PEDESTRIAN_OFF_KEY] = update(pomdp, up, bold[PEDESTRIAN_OFF_KEY], a, o[oid])
           # bnew[PEDESTRIAN_OFF_KEY] = initBeliefAbsentPedestrianBorder(pomdp, o[oid].ego_y, o[oid].ego_v)
           b_states = []
           b_prob = []
           for (s, prob) in weighted_iterator(bnew[PEDESTRIAN_OFF_KEY] )
               if ( prob > 1e-4)
                   push!(b_states, s)
                   push!(b_prob, prob)
               end
           end
           normalize!(b_prob, 1)

           bnew[PEDESTRIAN_OFF_KEY] =  SingleOCFBelief(b_states, b_prob)
           #println(bnew[PEDESTRIAN_OFF_KEY])
        else # ped appeared
            bnew[oid] = update(pomdp, up, bold[PEDESTRIAN_OFF_KEY], a, o[oid])
        end
    end

    return bnew
end

function POMDPs.update(pomdp::SingleOCFPOMDP, up::SingleOCFUpdater, b::SingleOCFBelief, a::SingleOCFAction, o::SingleOCFObs)

    states_p = SingleOCFState[]
    sizehint!(states_p, 1000);

    bp = Float64[]
    sizehint!(bp, 1000);

    (ego_y_state_space, ego_v_state_space) = getEgoDataInStateSpace(pomdp, o.ego_y, o.ego_v)

    # object is outside the defined state space
    if ( is_observation_absent(pomdp, o) )
        o = get_state_absent(pomdp,o)
        #println("Pedestrian is absent: ", o)
    else
       # println("Pedestrian is visible: ", o)
    end

    bp_sum = 0.0   # to normalize the distribution
    for sp_ped in pomdp.state_space_ped

        if ( sp_ped != pomdp.state_space_ped[end])  # not absent state
            sp = SingleOCFState(ego_y_state_space, ego_v_state_space, sp_ped.ped_s, sp_ped.ped_T, sp_ped.ped_theta ,sp_ped.ped_v)
        else # absent state
            sp = get_state_absent(pomdp, ego_y_state_space, ego_v_state_space)
        end
        po = observation_weight(pomdp, sp, o)

        if po < 0.000000000001
            continue
        end

        b_sum = 0.0
        for (s, prob) in weighted_iterator(b)
            td = transition(pomdp, s, a)
            pp = pdf(td, sp)
            b_sum += pp * prob
        end

        if b_sum > 0. # && b_sum > 1e-14
            push!(states_p, sp)
            push!(bp, po * b_sum)
            bp_sum += po * b_sum
        end

    end

    if bp_sum == 0.0

        # pedestrian appeared inside the state space, initialize belief with the observation
        if ( !is_observation_absent(pomdp,o) && length(b) == 1)
            for (s,p) in weighted_iterator(b)
                if is_state_absent(pomdp,s)
                    return initBeliefPedestrian(pomdp, o)
                end
            end
        end
# TODO: not perfect solution
 return initBeliefPedestrian(pomdp, o)

 error("""
              Failed discrete belief update: new probabilities sum to zero.
            #  b = $b
              a = $a
              o = $o
              Failed discrete belief update: new probabilities sum to zero.
              """)

    else
        bp ./= bp_sum
    end

   # println("b-length: ", length(states_p))
  #  println(states_p)
    return SingleOCFBelief(states_p, bp)

#=
    b_states = []
    b_prob = []
    for (s, prob) in weighted_iterator(SingleOCFBelief(states_p, bp) )
        if ( prob > 1e-4)
            push!(b_states, s)
            push!(b_prob, prob)
        end
    end
    return SingleOCFBelief(b_states, b_prob)
 =#

end

function POMDPs.initialstate_distribution(pomdp::SingleOCFPOMDP)

    bp = ones(n_states(pomdp)) / n_states(pomdp)

    states = SingleOCFState[]
    sizehint!(states, n_states(pomdp))

    for state in pomdp.state_space
        push!(states,state)
    end

    return SingleOCFBelief(states, bp)
end

function initBeliefAbsentPedestrianBorder(pomdp::SingleOCFPOMDP, ego_y::Float64, ego_v::Float64)

    states = SingleOCFState[]
    sizehint!(states, 500);
    probs = Float64[]
    sizehint!(probs, 500);
    (ego_y_state_space,ego_v_state_space) = getEgoDataInStateSpace(pomdp, ego_y, ego_v)

    s_min = pomdp.S_MIN
    if ( pomdp.env.params.obstacles_visible )
        for i = 1:length(pomdp.env.params.obstacles)
            (obst_s, obst_T, right_side) = getObstructionCorner(pomdp, pomdp.env.params.obstacles[i])
            if ( right_side )
                s_min = clamp(obst_s, pomdp.S_MIN, pomdp.S_MAX)
            end
        end
    end

    # add states on the right side
    for ped_theta in pomdp.PED_THETA_RANGE
        for ped_v in pomdp.PED_V_RANGE
            for ped_s in pomdp.S_RANGE
                if ( ped_s > s_min)
                    push!(states,SingleOCFState(ego_y_state_space, ego_v_state_space, ped_s, pomdp.T_MIN, ped_theta, ped_v))
                end
            end
        end
    end

    # add absent state
    absent_state = get_state_absent(pomdp, ego_y_state_space, pomdp.ego_vehicle.state.v)
    push!(states, absent_state)

    probs = ones(length(states))
    probs[1:end - 1] .= pomdp.PROBABILITY_PEDESTRIAN_BIRTH / length(states)
    probs[end] = 1.0 - pomdp.PROBABILITY_PEDESTRIAN_BIRTH

    return SingleOCFBelief(states,probs)

end

function initBeliefAbsentPedestrian(pomdp::SingleOCFPOMDP, ego_y::Float64, ego_v::Float64)

    states = SingleOCFState[]
    sizehint!(states, 1000);
    probs = Float64[]
    sizehint!(probs, 1000);

    (ego_y_state_space, ego_v_state_space) = getEgoDataInStateSpace(pomdp, ego_y, ego_v)

    # add occluded states
    if ( pomdp.env.params.obstacles_visible )
        for i = 1:length(pomdp.env.params.obstacles)
            ego_pos = VecE2(pomdp.ego_vehicle.state.posG.x, pomdp.ego_vehicle.state.posG.y)
            #println("ego_pos: ", ego_pos)
            #println("pomdp.env.params.obstacles[i]: ", pomdp.env.params.obstacles[i])

            (obst_s, obst_T, right_side) = getObstructionCorner(pomdp, pomdp.env.params.obstacles[i])
            #println("obst_s: ", obst_s, " obst_T: ", obst_T , " right_side: ", right_side)
            if ( right_side )
                occluded_positions = calulateHiddenPositionsRightSide(pomdp, obst_s, obst_T)

                for ped_theta in pomdp.PED_THETA_RANGE
                    for ped_v = 0.:1.:2.# in pomdp.PED_V_RANGE
                        for (hidden_s, hidden_T) in occluded_positions
                            push!(states,SingleOCFState(ego_y_state_space, ego_v_state_space, hidden_s, hidden_T, ped_theta, ped_v))
                        end
                    end
                end
            end
        end
    end

    # add absent state
    absent_state = get_state_absent(pomdp, ego_y_state_space, ego_v_state_space)
    push!(states,absent_state)

    probs = ones(length(states))
    probs[1:end - 1] .= pomdp.PROBABILITY_PEDESTRIAN_BIRTH / length(states)
    probs[end] = 1.0 - pomdp.PROBABILITY_PEDESTRIAN_BIRTH

    return SingleOCFBelief(states,probs)
end

function initBeliefPedestrian(pomdp::SingleOCFPOMDP, o::SingleOCFObs)

    states = SingleOCFState[]
    sizehint!(states, 2500);
    probs = Float64[]
    sizehint!(probs, 2500);

    (ego_y_state_space, ego_v_state_space) = getEgoDataInStateSpace(pomdp, o.ego_y, o.ego_v)

    for ped_theta in pomdp.PED_THETA_RANGE
        for ped_v in pomdp.PED_V_RANGE
            for ds = -1.0:1.0:1.0
                for dT = -1.0:1.0:1.0
                    obs = SingleOCFState(ego_y_state_space, ego_v_state_space, o.ped_s+ds, o.ped_T+dT, ped_theta, ped_v)
                    obs_int = pomdp.state_space[stateindex(pomdp, obs)]
                    push!(states, obs_int)
                end
            end
        end
    end
    probs = ones(length(states)) / length(states)

    return SingleOCFBelief(states,probs)
end


text:
using CxxWrap
using CxxWrap.CxxWrapCore: SafeCFunction
using BenchmarkTools

if !isdefined(Main, :M)
    include("ibapi-mod.jl")
end
using .M

function main2()
    # jc = M.MyClient(
    #     @safe_cfunction(fuckCb, Nothing, (Int32,))
    # )
    # c = M.TestCppClient(jc)

    # c = M.TestCppClient(
    #     @safe_cfunction(fuckCb, Nothing, (Int32,))
    # )
    #
    # for i=1:3
    #     M.getFuck(c)
    #     println(fuckVal)
    # end

    # @btime M.getFuck($c)
end

function main()
    c = M.Client()
    M.getFuck(c)
end


text:
#=
Macklemore Problem:

We are going to the thrift shop and need 99 cents to buy something. What is the
lightest way to carry 99 cents using US currency?

This is a Knapsack problem.

We specify that you need at least 99 cents - does the answer change if you need exact change?
=#

using JuMP
using Cbc # Must support integer programming.

m = Model(solver=CbcSolver())

# Variables represent how many of each coin we want to carry
@defVar(m, pennies >= 0, Int)
@defVar(m, nickels >= 0, Int)
@defVar(m, dimes >= 0, Int)
@defVar(m, quarters >= 0, Int)

# We need at least 99 cents
@addConstraint(m, 1 * pennies + 5 * nickels + 10 * dimes + 25 * quarters >= 99)

# Minimize mass (Grams)
# (source: US Mint)
@setObjective(m, Min, 2.5 * pennies + 5 * nickels + 2.268 * dimes + 5.670 * quarters)

# Solve
status = solve(m)

println("Minimum mass: ", getObjectiveValue(m), " grams")
println("using:")
println(round(getValue(pennies)), " pennies") # "round" to cast as integer
println(round(getValue(nickels)), " nickels")
println(round(getValue(dimes)), " dimes")
println(round(getValue(quarters)), " quarters")

max_stars_repo_name: JuliaQuantum/Resources

max_issues_repo_name: JuliaQuantum/Resources

max_forks_repo_name: JuliaQuantum/Resources


text:
using Graphs

# Imports
include("entities/GraphsDFG.jl")
include("services/GraphsDFG.jl")

# Exports
export GraphsDFG
export exists
export getLabelDict, getDescription, setDescription, getAddHistory, getSolverParams, setSolverParams

export getAddHistory, getDescription, getLabelDict
export addVariable!, addFactor!
export ls, lsf, getVariables, getFactors, getVariableIds, getFactorIds
export getVariable, getFactor
export updateVariable!, updateFactor!, updateVariableSolverData!
export deleteVariable!, deleteFactor!
export getAdjacencyMatrix
export getNeighbors
export getSubgraphAroundNode
export getSubgraph
export isFullyConnected, hasOrphans
export toDot, toDotFile


text:
"""
# module Reciprocal

# Examples

```jldoctest
julia>
```
"""
module Reciprocal

using StaticArrays: SVector

export KPoint

struct KPoint
    coordinates::SVector{3, Float64}
    weight::Float64
end

end


text:
export PhasePartition
# Thermodynamic states
export ThermodynamicState, PhaseEquil,
  PhaseNonEquil, LiquidIcePotTempSHumEquil

"""
    PhasePartition

Represents the mass fractions of the moist air mixture.

# Constructors

    PhasePartition(q_tot::Real[, q_liq::Real[, q_ice::Real]])
    PhasePartition(ts::ThermodynamicState)

See also [`PhasePartition_equil`](@ref)

# Fields

$(DocStringExtensions.FIELDS)
"""
struct PhasePartition{DT<:Real}
  "total specific humidity"
  tot::DT
  "liquid water specific humidity (default: `0`)"
  liq::DT
  "ice specific humidity (default: `0`)"
  ice::DT
end

PhasePartition(q_tot::DT,q_liq::DT) where {DT<:Real} =
  PhasePartition(q_tot, q_liq, zero(DT))
PhasePartition(q_tot::DT) where {DT<:Real} =
  PhasePartition(q_tot, zero(DT), zero(DT))

"""
    ThermodynamicState{DT}

A thermodynamic state, which can be initialized for
various thermodynamic formulations (via its sub-types).
All `ThermodynamicState`'s have access to functions to
compute all other thermodynamic properties.
"""
abstract type ThermodynamicState{DT} end

"""
    PhaseEquil{DT} <: ThermodynamicState

A thermodynamic state assuming thermodynamic equilibrium (therefore, saturation adjustment
is needed).

# Constructors

    PhaseEquil(e_int, q_tot, ρ)

# Fields

$(DocStringExtensions.FIELDS)
"""
struct PhaseEquil{DT} <: ThermodynamicState{DT}
  "internal energy"
  e_int::DT
  "total specific humidity"
  q_tot::DT
  "density of air (potentially moist)"
  density::DT
  "temperature: computed via [`saturation_adjustment`](@ref)"
  T::DT
end
PhaseEquil(e_int, q_tot, ρ) =
  PhaseEquil(e_int, q_tot, ρ, saturation_adjustment(e_int, ρ, q_tot))

"""
    LiquidIcePotTempSHumEquil(θ_liq_ice, q_tot, ρ, p)

Constructs a [`PhaseEquil`](@ref) thermodynamic state from liquid-ice potential temperature.

 - `θ_liq_ice` - liquid-ice potential temperature
 - `q_tot` - total specific humidity
 - `ρ` - density
 - `p` - pressure
"""
function LiquidIcePotTempSHumEquil(θ_liq_ice, q_tot, ρ, p)
    T = saturation_adjustment_q_tot_θ_liq_ice(θ_liq_ice, q_tot, ρ, p)
    q = PhasePartition(T, ρ, q_tot)
    e_int = internal_energy(T, q)
    PhaseEquil(e_int, q_tot, ρ, T)
end

"""
    PhaseNonEquil{DT} <: ThermodynamicState

A thermodynamic state assuming thermodynamic non-equilibrium (therefore, temperature can
be computed directly).

# Constructors

    PhaseNonEquil(e_int, q::PhasePartition, ρ)

# Fields

$(DocStringExtensions.FIELDS)

"""
struct PhaseNonEquil{DT} <: ThermodynamicState{DT}
  "internal energy"
  e_int::DT
  "phase partition"
  q::PhasePartition{DT}
  "density of air (potentially moist)"
  density::DT
end

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
using ImageBinarization
using Test, TestImages
using Statistics
using ImageCore, ImageTransformations
using OffsetArrays
using ReferenceTests

include("testutils.jl")

@testset "ImageBinarization.jl" begin
    include("util.jl")

    include("algorithms/single_histogram_threshold.jl")
    include("algorithms/adaptive_threshold.jl")
    include("algorithms/niblack.jl")
    include("algorithms/polysegment.jl")
    include("algorithms/sauvola.jl")
end

nothing


text:
# Load the Interpolator in the FEM module
FEM.quadquadinterpmod()
using FEM.QuadQuadInterpMod

module QuadQuadMod

using FEM.LinTrigInterpMod

import FEM: createinterp, creategps, createstorage, get_field,
            Bmatrix, doftypes, get_ndofs, get_geotype, get_ref_area, dNdxmatrix

import FEM: AbstractMaterialStatus, AbstractElemStorage, AbstractFElement, FENode2
import FEM: Vertex6, Point2, GaussPoint2, Du, Dv, GeoTrig

export QuadQuad

type QuadQuadStorage <: ElemStorage
    B::Matrix{Float64}
    DeBe::Matrix{Float64}
    Ke::Matrix{Float64}
    ɛ::Vector{Float64}
    f_int::Vector{Float64}
    u_field::Vector{Float64}
end

QuadQuadStorage() = QuadQuadStorage(zeros(4, 16), zeros(4,16),
                                  zeros(16,16), zeros(4), zeros(16), zeros(16))

type QuadQuad{T <: AbstractMaterialStatus} <: AbstractFElement{T}
    vertices::Vertex6
    gps::Vector{GaussPoint2}
    n::Int
    interp::QuadQuadInterp
    storage::QuadQuadStorage
    matstats::Vector{T}
    temp_matstats::Vector{T}
end
gausspoints(elem::QuadQuad) = elem.gps

# Constructor
function QuadQuad{T <: AbstractMaterialStatus}(vertices::Vertex6, n, interp::QuadQuadInterp,
                 storage::QuadQuadStorage, gps::Vector{GaussPoint2}, matstat::T)
    matstats = T[]
    temp_matstats = T[]
    for i in 1:length(gps)
        push!(matstats, copy(matstat))
        push!(temp_matstats, copy(matstat))
    end
    QuadQuad(vertices, gps, n, interp, storage, matstats, temp_matstats)
end

get_ref_area(::QuadQuad) = 4.0
get_geoelem(ele::QuadQuad) = GeoQTrig(ele.n, ele.vertices)
get_geotype(::QuadQuad) = GeoQTrig

createstorage(::Type{QuadQuad}) = QuadQuadStorage()
createinterp(::Type{QuadQuad}) = QuadQuadInterp()

function creategps(::Type{QuadQuad})
end

@inline function get_ndofs(::QuadQuad)
    return 12
end

# Avoids allocation (remove this?)
const QUADTRIG_DOFTYPES = [Du, Dv]
function doftypes(::QuadQuad, ::Int)
    return return QUADTRIG_DOFTYPES
end

function Bmatrix(elem::QuadQuad, gp::GaussPoint2, nodes::Vector{FENode2})
    dNdx = dNdxmatrix(elem.interp, gp.local_coords, elem.vertices, nodes)
    B = elem.storage.B
    for i in 1:6
        B[1, 2*i - 1] = dNdx[i, 1]
        B[2, 2*i]     = dNdx[i, 2]

        B[4, 2*i - 1] = dNdx[i, 2]
        B[4, 2*i]     = dNdx[i, 1]
    end
    return B
end

# Get the stress/strain in gausspoint i
get_field(elem::QuadQuad, ::Type{Stress}, i::Int) = elem.matstats[i].stress
get_field(elem::QuadQuad, ::Type{Strain}, i::Int) = elem.matstats[i].strain


text:
# Quadrotor Simulation Environment - Gabriel Barsi Haberfeld
# University of Illinois at Urbana-Champaign, 2019
# Heavily based on original MATLAB version by William Selby:
# https://www.wilselby.com/research/arducopter/simulation-environment/

# This script implements a basic PID position control and associated plots
# and animation. The code follows sequentially and everyrhing lives inside the
# parameter struct QuadParam. It initializes with default values.

using Parameters
using BenchmarkTools
using Plots

# If the struct can be non-mutable, restart julia to redefine
@with_kw mutable struct QuadParam @deftype Float64

        init = 0;     # used in initilization
        Ts = .01;     # Sampling time (100 Hz)
        sim_time = 10; # Simulation time (seconds)
        counter = 1;   # the counter that holds the time value

        # Plotting Variables
        t_plot::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}} = 0:Ts:sim_time-Ts;       # the time values

        # Environmental Parametes
        g = 9.81;     # Gravity (m/s^2)

        # Quadrotor Physical Parameters
        m = 1.4;      # Quadrotor mass (kg)
        l = .56;     # Distance from the center of mass to the each motor (m)
        t = .02;   #Thickness of the quadrotor's arms for drawing purposes (m)
        rot_rad = .1;   #Radius of the propellor (m)
        Kd = 1.3858e-6;    # Drag torque coeffecient (kg-m^2)

        Kdx = 0.16481;    # Translational drag force coeffecient (kg/s)
        Kdy = 0.31892;    # Translational drag force coeffecient (kg/s)
        Kdz = 1.1E-6;    # Translational drag force coeffecient (kg/s)

        Jx = .05;     # Moment of inertia about X axis (kg-m^2)
        Jy = .05;     # Moment of inertia about Y axis (kg-m^2)
        Jz = .24;    # Moment of inertia about Z axis (kg-m^2)

        # Quadrotor Sensor Paramaters
        GPS_freq = (1/Ts)/1;
        X_error = .01;  #+/- m
        Y_error = .01;  #+/- m
        Z_error = .01;  #+/- m
        phi_error = .01;
        theta_error = .01;
        psi_error = .01;

        x_acc_bias = 0.16594;  # m/s^2
        x_acc_sd = 0.0093907;
        y_acc_bias = 0.31691;  # m/s^2
        y_acc_sd = 0.011045;
        z_acc_bias = -8.6759;  # m/s^2
        z_acc_sd = 0.016189;

        x_gyro_bias = 0.00053417;  # rad/s
        x_gyro_sd = 0.00066675;
        y_gyro_bias = -0.0011035;  # rad/s
        y_gyro_sd = 0.00053642;
        z_gyro_bias = 0.00020838;  # rad/s
        z_gyro_sd = 0.0004403;

        ground_truth = 1;  # Use perfect sensor measurements
        sensor_unfiltered = 0; # Use sensor errors, no filter
        sensor_kf = 0;     # Use sensor error, Kalman Filter

        # Motor Parameters
        KT = 1.3328e-5;    # Thrust force coeffecient (kg-m)
        Jp = 0.044;     # Moment of Intertia of the rotor (kg-m^2)
        max_motor_speed = 925; # motors upper limit (rad/s)
        min_motor_speed = 0; #-1*((400)^2); # motors lower limit (can't spin in reverse)

        Obar = 0;     # sum of motor speeds (O1-O2+O3-O4, N-m)
        O1 = 0;       # Front motor speed (raidans/s)
        O2 = 0;       # Right motor speed (raidans/s)
        O3 = 0;       # Rear motor speed (raidans/s)
        O4 = 0;       # Left motor speed (raidans/s)

        #Translational Positions
        X = 0;        # Initial position in X direction GF (m)
        Y = 0;        # Initial position in Y direction GF (m)
        Z = 0;        # Initial position in Z direction GF (m)
        X_BF = 0;     # Initial position in X direction BF (m)
        Y_BF = 0;     # Initial position in Y direction BF (m)
        Z_BF = 0;     # Initial position in the Z direction BF (m)

        #Translational Velocities
        X_dot = 0;    # Initial velocity in X direction GF (m/s)
        Y_dot = 0;    # Initial velocity in Y direction GF (m/s)
        Z_dot = 0;    # Initial velocity in Z direction GF (m/s)
        X_dot_BF = 0;    # Initial velocity in X direction BF (m/s)
        Y_dot_BF = 0;    # Initial velocity in Y direction BF (m/s)
        Z_dot_BF = 0;    # Initial velocity in Y direction BF (m/s)

        #Angular Positions
        phi = 0;      # Initial phi value (rotation about X GF, roll,  radians)
        theta = 0;    # Initial theta value (rotation about Y GF, pitch, radians)
        psi = 0;      # Initial psi value (rotation about Z GF, yaw, radians)

        #Angular Velocities
        p = 0;        # Initial p value (angular rate rotation about X BF, radians/s)
        q = 0;        # Initial q value (angular rate rotation about Y BF, radians/s)
        r = 0;        # Initial r value (angular rate rotation about Z BF, radians/s)

        # Desired variables
        X_des_GF = 5;         # desired value of X in Global frame
        Y_des_GF = 5;         # desired value of Y in Global frame
        Z_des_GF = -0.5;      # desired value of Z in Global frame
        X_des = 0;            # desired value of X in Body frame
        Y_des = 0;            # desired value of Y in Body frame
        Z_des = 0;            # desired value of Z in Body frame

        phi_des = 0;          # desired value of phi (radians)
        theta_des = 0;        # desired value of theta (radians)
        psi_des = 0;          # desired value of psi (radians)

        p_des = 0;
        q_des = 0;
        r_des = 0;

        # Measured variables
        X_meas = 0;
        Y_meas = 0;
        Z_meas = 0;
        phi_meas = 0;
        theta_meas = 0;
        psi_meas = 0;

        # Disturbance Variables
        Z_dis = 0;            # Disturbance in Z direction
        X_dis = 0;            # Disturbance in X direction
        Y_dis = 0;            # Ddisturbance in Y direction
        phi_dis = 0;            # Disturbance in Yaw direction
        theta_dis = 0;            # Disturbance in Pitch direction
        psi_dis = 0;            # Disturbance in Roll direction

        # Control Inputs
        U1 = 0;       # Total thrust (N)
        U2 = 0;       # Torque about X axis BF (N-m)
        U3 = 0;       # Torque about Y axis BF (N-m)
        U4 = 0;       # Torque about Z axis BF (N-m)

        # Control Limits (update values)
        U1_max = 43.5;   # KT*4*max_motor_speed^2
        U1_min = 0;      #
        U2_max = 6.25;  # KT*l*max_motor_speed^2
        U2_min = -6.25; # KT*l*max_motor_speed^2
        U3_max = 6.25;  # KT*l*max_motor_speed^2
        U3_min = -6.25; # KT*l*max_motor_speed^2
        U4_max = 2.25;  # Kd*2*max_motor_speed^2
        U4_min = -2.25; # Kd*2*max_motor_speed^2

        # PID errors

        x_error_sum = 0;
        y_error_sum = 0
        z_error_sum = 0;
        phi_error_sum = 0;
        theta_error_sum = 0;
        psi_error_sum = 0;
        p_error_sum = 0;
        q_error_sum = 0;
        r_error_sum = 0;

        # PID parameters
        X_KP = .35;          # KP value in X position control
        X_KI = .25;            # KI value in X position control
        X_KD = -.35;         # KD value in X position control
        X_KI_lim = .25;         # Error to start calculating integral term

        Y_KP = .35;          # KP value in Y position control
        Y_KI = .25;            # KI value in Y position control
        Y_KD = -.35;         # KD value in Y position control
        Y_KI_lim = .25;         # Error to start calculating integral term

        Z_KP = 10/1.7;    # KP value in altitude control
        Z_KI = 0*3;    # KI value in altitude control
        Z_KD = -100/1.980;  # KD value in altitude control
        Z_KI_lim = .25;         # Error to start calculating integral term

        phi_KP = 4.5;      # KP value in roll control 2
        phi_KI = 0;       # KI value in roll control   1
        phi_KD = 0;     # KD value in roll control  -.5
        phi_max = pi/8;   # Maximum roll angle commanded
        phi_KI_lim = 2*(2*pi/360);  # Error to start calculating integral

        theta_KP = 4.5;    # KP value in pitch control 2
        theta_KI = 0;     # KI value in pitch control 1
        theta_KD = 0;   # KD value in pitch control -.5
        theta_max = pi/8; # Maximum pitch angle commanded
        theta_KI_lim = 2*(2*pi/360);  # Error to start calculating integral

        psi_KP = 10;     # KP value in yaw control
        psi_KI = 0;     # KI value in yaw control .75
        psi_KD = 0;     # KD value in yaw control -.5
        psi_KI_lim = 8*(2*pi/360);  # Error to start calculating integral

        p_KP = 2.7;    # KP value in pitch control 2
        p_KI = 1;     # KI value in pitch control
        p_KD = -.01;   # KD value in pitch control -.5
        p_max = 500*(2*pi/360); # Maximum pitch angle commanded
        p_KI_lim = 10*(2*pi/360);  # Error to start calculating integral

        q_KP = 2.7;    # KP value in pitch control
        q_KI = 1;     # KI value in pitch control
        q_KD = -.01;   # KD value in pitch control -.5
        q_max = 500*(2*pi/360); # Maximum pitch angle commanded
        q_KI_lim = 10*(2*pi/360);  # Error to start calculating integral

        r_KP = 2.7;    # KP value in pitch control
        r_KI = 1;     # KI value in pitch control
        r_KD = -.01;   # KD value in pitch control
        r_max = 50*(2*pi/360); # Maximum pitch angle commanded
        r_KI_lim = 10*(2*pi/360);  # Error to start calculating integral

        #dynamics
        X_ddot = (-cos(phi)*sin(theta)*cos(psi)+sin(phi)*sin(psi)*U1-Kdx*X_dot)/m;
        Y_ddot = (-cos(phi)*sin(psi)*sin(theta)-cos(psi)*sin(phi)*U1-Kdy*Y_dot)/m;
        Z_ddot = (-cos(phi)*cos(theta)*U1-Kdz*Z_dot)/m+g;

        p_dot = (q*r*(Jy - Jz) - Jp*p*Obar + l*U2)/Jx;
        q_dot = (p*r*(Jz - Jx) + Jp*q*Obar + l*U3)/Jy;
        r_dot = (p*q*(Jx - Jy) + U4)/Jz;

        phi_dot   = p + sin(phi)*tan(theta)*q + cos(phi)*tan(theta)*r;
        theta_dot = cos(phi)*q - sin(phi)*r;
        psi_dot   = sin(phi)/cos(theta)*q + cos(phi)/cos(theta)*r;
end

# Nonlinear Dynamics
function QuadDynamicsUpdate(Quad)

        ## Update Accelerations

        Quad.X_ddot = (-(cos(Quad.phi)*sin(Quad.theta)*cos(Quad.psi)+sin(Quad.phi)*sin(Quad.psi))*Quad.U1-Quad.Kdx*Quad.X_dot)/Quad.m;
        Quad.Y_ddot = (-(cos(Quad.phi)*sin(Quad.psi)*sin(Quad.theta)-cos(Quad.psi)*sin(Quad.phi))*Quad.U1-Quad.Kdy*Quad.Y_dot)/Quad.m;
        Quad.Z_ddot = (-(cos(Quad.phi)*cos(Quad.theta))*Quad.U1-Quad.Kdz*Quad.Z_dot)/Quad.m+Quad.g;

        Quad.p_dot = (Quad.q*Quad.r*(Quad.Jy - Quad.Jz) - Quad.Jp*Quad.p*Quad.Obar + Quad.l*Quad.U2)/Quad.Jx;
        Quad.q_dot = (Quad.p*Quad.r*(Quad.Jz - Quad.Jx) + Quad.Jp*Quad.q*Quad.Obar + Quad.l*Quad.U3)/Quad.Jy;
        Quad.r_dot = (Quad.p*Quad.q*(Quad.Jx - Quad.Jy) + Quad.U4)/Quad.Jz;

        Quad.phi_dot   = Quad.p + sin(Quad.phi)*tan(Quad.theta)*Quad.q + cos(Quad.phi)*tan(Quad.theta)*Quad.r;
        Quad.theta_dot = cos(Quad.phi)*Quad.q - sin(Quad.phi)*Quad.r;
        Quad.psi_dot   = sin(Quad.phi)/cos(Quad.theta)*Quad.q + cos(Quad.phi)/cos(Quad.theta)*Quad.r;

        ## Disturbance model

        Quad.X_ddot = Quad.X_ddot + Quad.X_dis/Quad.m;
        Quad.Y_ddot = Quad.Y_ddot + Quad.Y_dis/Quad.m;
        Quad.Z_ddot = Quad.Z_ddot + Quad.Z_dis/Quad.m;
        Quad.phi_dot = Quad.phi_dot + Quad.phi_dis/Quad.Jx*Quad.Ts;
        Quad.theta_dot = Quad.theta_dot + Quad.theta_dis/Quad.Jy*Quad.Ts;
        Quad.psi_dot = Quad.psi_dot + Quad.psi_dis/Quad.Jz*Quad.Ts;

        ## Update Velocities and Positions

        # Calculating the Z velocity & position
        Quad.Z_dot = Quad.Z_ddot*Quad.Ts + Quad.Z_dot;
        Quad.Z = Quad.Z_dot*Quad.Ts + Quad.Z;

        # Calculating the X velocity & position
        Quad.X_dot = Quad.X_ddot*Quad.Ts + Quad.X_dot;
        Quad.X = Quad.X_dot*Quad.Ts + Quad.X;

        # Calculating the Y velocity & position
        Quad.Y_dot = Quad.Y_ddot*Quad.Ts + Quad.Y_dot;
        Quad.Y = Quad.Y_dot*Quad.Ts + Quad.Y;

        # Calculating p,q,r
        Quad.p = Quad.p_dot*Quad.Ts+Quad.p;
        Quad.q = Quad.q_dot*Quad.Ts+Quad.q;
        Quad.r = Quad.r_dot*Quad.Ts+Quad.r;

        # Calculating angular velocity and position
        Quad.phi = Quad.phi_dot*Quad.Ts + Quad.phi;
        Quad.theta = Quad.theta_dot*Quad.Ts+Quad.theta;
        Quad.psi = Quad.psi_dot*Quad.Ts+Quad.psi;

        return Quad
end

function QuadSensorMeasurements(Quad)

        if Quad.ground_truth == 1
            noise = 0;
        else
            noise = 1;
        end

        ## VICON
        Quad.X_meas = Quad.X + noise*(randn()*Quad.X_error);
        Quad.Y_meas = Quad.Y + noise*(randn()*Quad.Y_error);
        Quad.Z_meas = Quad.Z + noise*(randn()*Quad.Z_error);
        Quad.phi_meas = Quad.phi + noise*(randn()*Quad.phi_error);
        Quad.theta_meas = Quad.theta + noise*(randn()*Quad.theta_error);
        Quad.psi_meas = Quad.psi + noise*(randn()*Quad.psi_error);

        ## IMU Measurements

        Quad.X_ddot = Quad.X_ddot + noise*(Quad.x_acc_bias + Quad.x_acc_sd*randn());
        Quad.Y_ddot = Quad.Y_ddot + noise*(Quad.y_acc_bias + Quad.y_acc_sd*randn());
        Quad.Z_ddot = Quad.Z_ddot + noise*(Quad.z_acc_bias + Quad.z_acc_sd*randn());

        Quad.p = Quad.p + noise*(Quad.x_gyro_bias + Quad.x_gyro_sd*randn());
        Quad.q = Quad.q + noise*(Quad.y_gyro_bias + Quad.y_gyro_sd*randn());
        Quad.r = Quad.r + noise*(Quad.z_gyro_bias + Quad.z_gyro_sd*randn());

        return Quad
end

function QuadPositionPID(Quad)

        x = Quad.X_meas;
        y = Quad.Y_meas;
        z = Quad.Z_meas;
        phi = Quad.phi_meas;
        theta = Quad.theta_meas;
        psi = Quad.psi_meas;

        # Rotate Desired Position from GF to BF (Z axis rotation only maybe)
        pts = [Quad.X_des_GF Quad.Y_des_GF Quad.Z_des_GF]'
        pts = rotateGFtoBF(pts,1*phi,1*theta,psi);
        Quad.X_des = pts[1]
        Quad.Y_des = pts[2]
        Quad.Z_des = pts[3]

        # Rotate Current Position from GF to BF
        pts = [x y z]'
        pts = rotateGFtoBF(pts,phi,theta,psi);
        Quad.X_BF = pts[1]
        Quad.Y_BF = pts[2]
        Quad.Z_BF = pts[3]

        # Rotate Current Velocity from GF to BF
        pts = [Quad.X_dot Quad.Y_dot Quad.Z_dot]'
        pts = rotateGFtoBF(pts,phi,theta,psi);
        Quad.X_dot_BF = pts[1]
        Quad.Y_dot_BF = pts[2]
        Quad.Z_dot_BF = pts[3]

        # X Position PID controller
        x_error = Quad.X_des - Quad.X_BF;
        if(abs(x_error) < Quad.X_KI_lim) #only use KI for fine tracking
            Quad.x_error_sum = Quad.x_error_sum + x_error;
        end
        cp = Quad.X_KP*x_error;    #Proportional term
        ci = Quad.X_KI*Quad.Ts*Quad.x_error_sum;
        ci = min(Quad.theta_max, max(-Quad.theta_max, ci));    #Saturate ci
        cd = Quad.X_KD*Quad.X_dot_BF;                    #Derivative term
        Quad.theta_des =  - (cp + ci + cd);   #Theta and X inversely related
        Quad.theta_des = min(Quad.theta_max, max(-Quad.theta_max, Quad.theta_des));

        # Y Position PID controller
        y_error = Quad.Y_des - Quad.Y_BF;
        if(abs(y_error) < Quad.Y_KI_lim)
            Quad.y_error_sum = Quad.y_error_sum + y_error;
        end
        cp = Quad.Y_KP*y_error;    #Proportional term
        ci = Quad.Y_KI*Quad.Ts*Quad.y_error_sum;
        ci = min(Quad.phi_max, max(-Quad.phi_max, ci));    #Saturate ci
        cd = Quad.Y_KD*Quad.Y_dot_BF;                      #Derivative term
        Quad.phi_des = cp + ci + cd;
        Quad.phi_des = min(Quad.phi_max, max(-Quad.phi_max, Quad.phi_des));

        ## Z Position PID Controller/Altitude Controller
        z_error = Quad.Z_des_GF-Quad.Z_BF;
        if(abs(z_error) < Quad.Z_KI_lim)
            Quad.z_error_sum = Quad.z_error_sum + z_error;
        end
        cp = Quad.Z_KP*z_error;         #Proportional term
        ci = Quad.Z_KI*Quad.Ts*Quad.z_error_sum; #Integral term
        ci = min(Quad.U1_max, max(Quad.U1_min, ci));    #Saturate ci
        cd = Quad.Z_KD*Quad.Z_dot;                  #Derivative term
        Quad.U1 = -(cp + ci + cd)/(cos(theta)*cos(phi)) + (Quad.m * Quad.g)/(cos(theta)*cos(phi));   #Negative since Thurst and Z inversely related
        Quad.U1 = min(Quad.U1_max, max(Quad.U1_min, Quad.U1));

        return Quad
end

function QuadAttitudePID(Quad)

        phi = Quad.phi_meas;
        theta = Quad.theta_meas;
        psi = Quad.psi_meas;

        # Roll PID Controller
        phi_error = Quad.phi_des - phi;
        if(abs(phi_error) < Quad.phi_KI_lim)
            Quad.phi_error_sum = Quad.phi_error_sum + phi_error;
        end
        cp = Quad.phi_KP*phi_error;
        ci = Quad.phi_KI*Quad.Ts*Quad.phi_error_sum;
        ci = min(Quad.p_max, max(-Quad.p_max, ci));
        cd = Quad.phi_KD*Quad.p;
        Quad.p_des = cp + ci + cd;
        Quad.p_des = min(Quad.p_max, max(-Quad.p_max, Quad.p_des));

        # Pitch PID Controller
        theta_error = Quad.theta_des - theta;
        if(abs(theta_error) < Quad.theta_KI_lim)
            Quad.theta_error_sum = Quad.theta_error_sum + theta_error;
        end
        cp = Quad.theta_KP*theta_error;
        ci = Quad.theta_KI*Quad.Ts*Quad.theta_error_sum;
        ci = min(Quad.q_max, max(-Quad.q_max, ci));
        cd = Quad.theta_KD*Quad.q;
        Quad.q_des = cp + ci + cd;
        Quad.q_des = min(Quad.q_max, max(-Quad.q_max, Quad.q_des));

        # Yaw PID Controller
        psi_error = Quad.psi_des - psi;
        if(abs(psi_error) < Quad.psi_KI_lim)
            Quad.psi_error_sum = Quad.psi_error_sum + psi_error;
        end
        cp = Quad.psi_KP*psi_error;
        ci = Quad.psi_KI*Quad.Ts*Quad.psi_error_sum;
        ci = min(Quad.r_max, max(-Quad.r_max, ci));
        cd = Quad.psi_KD*Quad.r;
        Quad.r_des = cp + ci + cd;
        Quad.r_des = min(Quad.r_max, max(-Quad.r_max, Quad.r_des));

        return Quad
end

function QuadRatePID(Quad)

        p = Quad.p;
        q = Quad.q;
        r = Quad.r;

        ## Angular Rate Controller

        # Roll PID Controller
        p_error = Quad.p_des - p;
        if(abs(p_error) < Quad.p_KI_lim)
            Quad.p_error_sum = Quad.p_error_sum + p_error;
        end
        cp = Quad.p_KP*p_error;
        ci = Quad.p_KI*Quad.Ts*Quad.p_error_sum;
        ci = min(Quad.U2_max, max(Quad.U2_min, ci));
        cd = Quad.p_KD*Quad.p_dot;
        Quad.U2 = cp + ci + cd;
        Quad.U2 = min(Quad.U2_max, max(Quad.U2_min, Quad.U2));

        # Pitch PID Controller
        q_error = Quad.q_des - q;
        if(abs(q_error) < Quad.q_KI_lim)
            Quad.q_error_sum = Quad.q_error_sum + q_error;
        end
        cp = Quad.q_KP*q_error;
        ci = Quad.q_KI*Quad.Ts*Quad.q_error_sum;
        ci = min(Quad.U3_max, max(Quad.U3_min, ci));
        cd = Quad.q_KD*Quad.q_dot;
        Quad.U3 = cp + ci + cd;
        Quad.U3 = min(Quad.U3_max, max(Quad.U3_min, Quad.U3));

        # Yaw PID Controller
        r_error = Quad.r_des - r;
        if(abs(r_error) < Quad.r_KI_lim)
            Quad.r_error_sum = Quad.r_error_sum + r_error;
        end
        cp = Quad.r_KP*r_error;
        ci = Quad.r_KI*Quad.Ts*Quad.r_error_sum;
        ci = min(Quad.U4_max, max(Quad.U4_min, ci));
        cd = Quad.r_KD*Quad.r_dot;
        Quad.U4 = cp + ci + cd;
        Quad.U4 = min(Quad.U4_max, max(Quad.U4_min, Quad.U4));

        return Quad
end

function rotateGFtoBF(pts,phi,theta,psi)

        # define rotation matrix
          R_roll = [
                  1     0           0;
                  0     cos(phi)    -sin(phi);
                  0     sin(phi)    cos(phi)];
          R_pitch = [
                  cos(theta)    0       sin(theta);
                  0             1       0;
                  -sin(theta)   0       cos(theta)];
          R_yaw = [
                  cos(psi)      -sin(psi)       0;
                  sin(psi)      cos(psi)        0;
                  0             0               1];

          R = R_yaw*R_pitch*R_roll;
          pts = pts'*R;

        return pts
end

function QuadMotorSpeeds(Quad)

        # Calculate motor speeds (rad/s)^2
        w1 = Quad.U1/(4*Quad.KT) + Quad.U3/(2*Quad.KT*Quad.l) + Quad.U4/(4*Quad.Kd);
        w2 = Quad.U1/(4*Quad.KT) - Quad.U2/(2*Quad.KT*Quad.l) - Quad.U4/(4*Quad.Kd);
        w3 = Quad.U1/(4*Quad.KT) - Quad.U3/(2*Quad.KT*Quad.l) + Quad.U4/(4*Quad.Kd);
        w4 = Quad.U1/(4*Quad.KT) + Quad.U2/(2*Quad.KT*Quad.l) - Quad.U4/(4*Quad.Kd);

        # Apply realistic motor speed limits
        if w1 > Quad.max_motor_speed^2
            w1 = Quad.max_motor_speed^2;
        end
        if w1 < Quad.min_motor_speed^2
            w1 = Quad.min_motor_speed^2;
        end

        if w2 > Quad.max_motor_speed^2
            w2 = Quad.max_motor_speed^2;
        end
        if w2 < Quad.min_motor_speed^2
            w2 = Quad.min_motor_speed^2;
        end

        if w3 > Quad.max_motor_speed^2
            w3 = Quad.max_motor_speed^2;
        end
        if w3 < Quad.min_motor_speed^2
            w3 = Quad.min_motor_speed^2;
        end

        if w4 > Quad.max_motor_speed^2
            w4 = Quad.max_motor_speed^2;
        end
        if w4 < Quad.min_motor_speed^2
            w4 = Quad.min_motor_speed^2;
        end

        Quad.O1 = sqrt(w1);    # Front M
        Quad.O2 = sqrt(w2);    # Right M
        Quad.O3 = sqrt(w3);    # Rear M
        Quad.O4 = sqrt(w4);    # Left M

        ## Re-compute traditional control inputs

        Quad.U1 = Quad.KT*(Quad.O1^2 + Quad.O2^2 + Quad.O3^2 + Quad.O4^2);
        Quad.U2 = Quad.KT*Quad.l*(Quad.O4^2 - Quad.O2^2);
        Quad.U3 = Quad.KT*Quad.l*(Quad.O1^2 - Quad.O3^2);
        Quad.U4 = Quad.Kd*(Quad.O1^2 + Quad.O3^2 - Quad.O2^2 - Quad.O4^2);
        Quad.Obar = (Quad.O1 - Quad.O2 + Quad.O3 - Quad.O4);

        return Quad
end

function QuadTorqueDisturbances(Quad,t)

        Tx = 1.0*randn() + 0.5*cos(t*5)
        Ty = 1.0*randn() + 0.5*sin(t*3)
        Tz = 1.0*randn() + 0.5*sin(t*1)

        Quad.U2 += Tx
        Quad.U3 += Ty
        Quad.U4 += Tz

        return Quad
end

function QuadSim(Xd,Yd,Zd,ground_truth)

        Quad = QuadParam()
        Quad.X_des_GF = Xd
        Quad.Y_des_GF = Yd
        Quad.Z_des_GF = Zd
        Quad.ground_truth = ground_truth
        Quad.sim_time = 10
        Quad.t_plot::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}} = 0:Quad.Ts:Quad.sim_time-Quad.Ts;
        l = size(Quad.t_plot,1) + 1
        U1 = zeros(l)
        U2 = zeros(l)
        U3 = zeros(l)
        U4 = zeros(l)
        X = zeros(l)
        Y = zeros(l)
        Z = zeros(l)
        phi = zeros(l)
        theta = zeros(l)
        psi = zeros(l)

        t_plot = 0:Quad.Ts:Quad.sim_time
        c = 1

        for t = 0:Quad.Ts:Quad.sim_time
                # #Spiral Trajectory
                # Quad.X_des_GF = 4*sin(3*t/pi)
                # Quad.Y_des_GF = 4*cos(3*t/pi)
                # Quad.Z_des_GF = -2*sin(2*t/pi)-1
                # Figure Eight
                Quad.X_des_GF = 4*sin(2*t/pi)
                Quad.Y_des_GF = 4*sin(2*t/pi)*cos(2*t/pi)
                Quad.Z_des_GF = Zd
                QuadSensorMeasurements(Quad)
                QuadPositionPID(Quad)
                QuadAttitudePID(Quad)
                QuadRatePID(Quad)
                QuadTorqueDisturbances(Quad,t)
                #QuadMotorSpeeds(Quad) # this is broken
                QuadDynamicsUpdate(Quad)
                U1[c] = Quad.U1
                U2[c] = Quad.U2
                U3[c] = Quad.U3
                U4[c] = Quad.U4
                X[c] = Quad.X
                Y[c] = Quad.Y
                Z[c] = -Quad.Z # flip for plotting (NED)
                phi[c] = Quad.phi
                theta[c] = Quad.theta
                psi[c] = Quad.psi
                c = c + 1
        end

        return Quad, t_plot, U1, U2, U3, U4, X, Y, Z, phi, theta, psi
end

function QuadPlot(j, plotlim, qb)
        for i = 1:size(qb,1)
                qb[i,:] = rotateGFtoBF(qb[i,:], phi[j], theta[j], psi[j])
        end
        qbx = qb[:,1] .+ X[j]
        qby = qb[:,2] .+ Y[j]
        qbz = qb[:,3] .+ Z[j]
        plot(qbx,qby,qbz,lw=2 )
        plot!(X[1:j], Y[1:j], Z[1:j],lw=2)
        plot!(X[1:j], Y[1:j], 0 .*Z[1:j],lw=2, color=:black)
        scatter!(qbx[1:4],qby[1:4],qbz[1:4], markersize=10, markercolor=:blue, markeralpha=0.5)
        scatter!(qbx[3:4],qby[3:4],qbz[3:4], markersize=10, markercolor=:red, markeralpha=0.5)
        plot!(aspect_ratio=:equal, xlim=(-plotlim, plotlim), ylim=(-plotlim, plotlim), zlim=(-0.1*plotlim, plotlim), legend=:false)
        plot!(xlabel = "x", ylabel = "y")
        plot!(size = (1064, 768))
end

Quad, t_plot, U1, U2, U3, U4, X, Y, Z, phi, theta, psi  = QuadSim(-4,-3,-1,0)
@show Quad.X Quad.Y Quad.Z

l = Quad.l

frames = 30
animtime = Quad.sim_time
numframes = animtime*frames
skip = Int64(div(size(X,1),numframes))
plotlim = 5

if true

        anim = @animate for j = 1:skip:size(X,1)
                qb = [  l l 0 ;
                        -l l 0 ;
                        -l -l 0 ;
                        l -l 0 ;
                        l l 0 ]
                QuadPlot(j, plotlim, qb)
                progress = div((j+2)*100,size(X,1))
                @show progress
        end
        gif(anim, "quad.gif", fps = frames)

end
p1 = plot(t_plot,U1,label="F")
p1 = plot!(t_plot,U2,label="Tx")
p1 = plot!(t_plot,U3,label="Ty")
p1 = plot!(t_plot,U4,label="Tz")

p2 = plot(t_plot,X,label="X")
p2 = plot!(t_plot,Y,label="Y")
p2 = plot!(t_plot,Z,label="Z")

p3 = plot(t_plot, phi, label = "roll")
p3 = plot!(t_plot, theta, label = "pitch")
p3 = plot!(t_plot, psi, label = "yaw")

plot(p1, p2, p3, layout=(3,1))


text:
path = @__DIR__
ml = CellModel(path * "/../models/noble_1962/Noble_1962.cellml")

prob = ODEProblem(ml, (0,10000.0))
sol1 = solve(prob, Euler(), dt=0.01, saveat=1.0)
sol2 = solve(prob, TRBDF2(), dtmax=0.5, saveat=1.0)
V1 = map(x -> x[2], sol1.u)
V2 = map(x -> x[2], sol2.u)
err = sum(abs.(V1 .- V2)) / length(V1)
@test err < 1.0


text:
# Test reading in OSM data
module TestReadData

using OpenStreetMap
using Base.Test

MAP_FILENAME = "tech_square.osm"

if !isfile(MAP_FILENAME)
    url = "https://dl.dropboxusercontent.com/u/8297575/$MAP_FILENAME"
    download(url, MAP_FILENAME)
end

nodes, hwys, builds, feats = getOSMData(MAP_FILENAME)

@test length(nodes) == 1410
@test length(hwys) == 55
@test length(builds) == 32
@test length(feats) == 4

map = parseMapXML(MAP_FILENAME)
bounds = getBounds(map)

@test_approx_eq bounds.min_y 42.3626000
@test_approx_eq bounds.max_y 42.3659000
@test_approx_eq bounds.min_x -71.0939000
@test_approx_eq bounds.max_x -71.0891000

end # module TestReadData


text:
#####
##### Initial value problem
#####

struct InitialValueProblem{FT, IC, BC} <: AbstractSimpleBoxProblem
    Lˣ::FT
    Lʸ::FT
    H::FT
    initial_conditions::IC
    boundary_conditions::BC

    """
        InitialValueProblem{FT}(; dimensions, initial_conditions=InitialConditions(),
                                boundary_conditions = (OceanBC(Impenetrable(FreeSlip()), Insulating()),
                                                       OceanBC(Penetrable(FreeSlip()), Insulating())))

    Returns an `InitialValueProblem` with `dimensions = (Lˣ, Lʸ, H)`, `initial_conditions`,
    and `boundary_conditions`.

    The default `initial_conditions` are resting with no temperature perturbation;
    the default list of `boundary_conditions` provide implementations for an impenetrable, insulating
    boundary, and a penetrable, insulating boundary.
    """
    function InitialValueProblem{FT}(;
        dimensions,
        initial_conditions = InitialConditions(),
        boundary_conditions = (
            OceanBC(Impenetrable(FreeSlip()), Insulating()),
            OceanBC(Penetrable(FreeSlip()), Insulating()),
        ),
    ) where {FT}

        return new{FT, typeof(initial_conditions), typeof(boundary_conditions)}(
            FT.(dimensions)...,
            initial_conditions,
            boundary_conditions,
        )
    end

end

#####
##### Initial conditions
#####

resting(x, y, z) = 0

struct InitialConditions{U, V, T, E}
    u::U
    v::V
    θ::T
    η::E
end

"""
    InitialConditions(; u=resting, v=resting, θ=resting, η=resting)

Stores initial conditions for each prognostic variable provided as functions of `x, y, z`.

Example
=======

# A Gaussian surface perturbation
a = 0.1 # m, amplitude
L = 1e5 # m, horizontal scale of the perturbation

ηᵢ(x, y, z) = a * exp(-(x^2 + y^2) / 2L^2)

ics = InitialConditions(η=ηᵢ)
"""
InitialConditions(; u = resting, v = resting, θ = resting, η = resting) =
    InitialConditions(u, v, θ, η)

"""
    ocean_init_state!(::HydrostaticBoussinesqModel, ic::InitialCondition, state, aux, local_geometry, time)

Initialize the state variables `u = (u, v)` (a vector), `θ`, and `η`. Mutates `state`.

This function is called by `init_state_prognostic!(::HydrostaticBoussinesqModel, ...)`.
"""
function ocean_init_state!(
    ::HydrostaticBoussinesqModel,
    ivp::InitialValueProblem,
    state,
    aux,
    local_geometry,
    time,
)

    ics = ivp.initial_conditions
    x, y, z = local_geometry.coord

    state.u = @SVector [ics.u(x, y, z), ics.v(x, y, z)]
    state.θ = ics.θ(x, y, z)
    state.η = ics.η(x, y, z)

    return nothing
end


text:
# Copyright 2015 Martin Holters
# See accompanying license file.

using ACME

function sallenkey(::Type{Circuit})
    r1 = resistor(10e3)
    r2 = resistor(10e3)
    c1 = capacitor(10e-9)
    c2 = capacitor(10e-9)
    j_in = voltagesource()
    u1 = opamp()
    #u1 = opamp(Val{:macak}, 1000, -4, 4)
    j_out = voltageprobe()

    circ = Circuit()
    connect!(circ, j_in[:-], :gnd)
    connect!(circ, j_in[:+], r1[1])
    connect!(circ, r1[2], r2[1], c1[1])
    connect!(circ, r2[2], c2[1], u1["in+"])
    connect!(circ, c2[2], u1["out-"], j_out[:-], :gnd)
    connect!(circ, u1["in-"], u1["out+"], j_out[:+])

    return circ
end

sallenkey{T<:DiscreteModel}(::Type{T}=DiscreteModel; fs=44100) =
    T(sallenkey(Circuit), 1/fs)


text:
@testset "hartree_fock.extract_tij_Uijlk" begin
    ofermion = pyimport("openfermion")
    FermionOperator = ofermion.ops.operators.FermionOperator
    ham = FermionOperator("1^ 0^ 1 0") + FermionOperator("0^ 0") + FermionOperator("1^ 1", -1.0)
    tij, Uijlk = extract_tij_Uijlk(ham)

    tij_ref = [1.0 0.0; 0.0 -1.0]
    Uijlk_ref = zeros(Float64, 2, 2, 2, 2)
    Uijlk_ref[2, 1, 2, 1] = 2.0
    @test all(tij .== tij_ref)
    @test all(Uijlk .== Uijlk_ref)
end

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
struct UDFRAcomp
  passengerid::String
  flightid::String
  originairport::String
  dstairport::String
  departuretime::DateTime
  arrivaltime::DateTime
  totalflighttime::Int16
end

function UDFAirportName{T<:AbstractString}(value::T)
  value = uppercase(value)
  length(value) < 3 && throw(UDFException("Airport Name must be longer than 3 characters: $value"))
  length(value) > 20 && throw(UDFException("Airport Name must be shorter than 20 symbols: $value"))
  value
end

function UDFAirportCode{T<:AbstractString}(value::T)
  value = uppercase(value)
  m = match(r"^[A-Z]{3}$", value)
  typeof(m) === Void && throw(UDFException("Airport id must be comprised of three capital letters. Trying to initialise with $value instead."))
  value
end

function UDFFlightId{T<:AbstractString}(value::T)
  value = uppercase(value)
  m = match(r"^[A-Z]{3}[0-9]{4}[A-Z]$", value)
  typeof(m) === Void && throw(UDFException("Invalid flight id: $value"))
  value
end

function UDFPassengerId{T<:AbstractString}(value::T)
  value = uppercase(value)
  m = match(r"^[A-Z]{3}[0-9]{4}[A-Z]{2}[0-9]$", value)
  typeof(m) === Void && throw(UDFException("Invalid passenger id: $value"))
  value
end

UDFDepartureTime{T<:Real}(value::T) = Dates.unix2datetime(value)
function UDFDepartureTime{T<:AbstractString}(value::T)
  m = match(r"^[0-9]+$", value)
  typeof(m) === Void && throw(UDFException("Invalid DepartureTime: $value"))
  UDFDepartureTime(parse(value))
end

function UDFTotalFlightTime{T<:AbstractString}(value::T)
  m = match(r"^[0-9]+$", value)
  typeof(m) === Void && throw(UDFException("Invalid flight time, cannot parse to integer: $value"))
  x = parse(Int16, value)
  x < 0 && throw(UDFException("total flight time must be a positive integer"))
  x > 1800 && warn("Total flight time $x minutes is more than 30 hours")
  x
end


text:
using ModelingToolkit, MethodOfLines, LinearAlgebra, OrdinaryDiffEq
using ModelingToolkit: operation, istree, arguments
using DomainSets
using NonlinearSolve

# # Define some variables
@testset "Heat Equation 1D 2 variables" begin
    @parameters t x
    @variables u(..) v(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dxx = Differential(x)^2
    eqs = [Dt(u(t, x)) ~ Dxx(u(t, x)),
        Dt(v(t, x)) ~ Dxx(v(t, x))]
    bcs = [u(0, x) ~ -x * (x - 1) * sin(x),
        v(0, x) ~ -x * (x - 1) * sin(x),
        u(t, 0) ~ 0.0, u(t, 1) ~ 0.0,
        v(t, 0) ~ 0.0, v(t, 1) ~ 0.0]

    domains = [t ∈ Interval(0.0, 1.0),
        x ∈ Interval(0.0, 1.0)]

    @named pdesys = PDESystem(eqs, bcs, domains, [t, x], [u(t, x), v(t, x)])
    discretization = MOLFiniteDifference([x => 0.1], t; grid_align=edge_align)
    prob = discretize(pdesys, discretization) # This gives an ODEProblem since it's time-dependent
    sol = solve(prob, Tsit5())
end

@testset "Heat Equation 2D 1 variable" begin
    @parameters t x y
    @variables u(..)
    Dxx = Differential(x)^2
    Dyy = Differential(y)^2
    Dt = Differential(t)
    t_min = 0.0
    t_max = 2.0
    x_min = 0.0
    x_max = 2.0
    y_min = 0.0
    y_max = 2.0

    # 3D PDE
    eq = Dt(u(t, x, y)) ~ Dxx(u(t, x, y)) + Dyy(u(t, x, y))

    analytic_sol_func(t, x, y) = exp(x + y) * cos(x + y + 4t)
    # Initial and boundary conditions
    bcs = [u(t_min, x, y) ~ analytic_sol_func(t_min, x, y),
        u(t, x_min, y) ~ analytic_sol_func(t, x_min, y),
        u(t, x_max, y) ~ analytic_sol_func(t, x_max, y),
        u(t, x, y_min) ~ analytic_sol_func(t, x, y_min),
        u(t, x, y_max) ~ analytic_sol_func(t, x, y_max)]

    # Space and time domains
    domains = [t ∈ Interval(t_min, t_max),
        x ∈ Interval(x_min, x_max),
        y ∈ Interval(y_min, y_max)]
    @named pdesys = PDESystem([eq], bcs, domains, [t, x, y], [u(t, x, y)])
    # Method of lines discretization
    dx = 0.1
    dy = 0.2
    discretization = MOLFiniteDifference([x => dx, y => dy], t)
    prob = ModelingToolkit.discretize(pdesys, discretization)
    sol = solve(prob, Tsit5())
end

# Diffusion in a sphere
@testset "Spherical Diffusion" begin

    @parameters t r
    @variables u(..)
    Dt = Differential(t)
    Dr = Differential(r)
    Drr = Dr^2
    eq = Dt(u(t, r)) ~ 1 / r^2 * Dr(r^2 * Dr(u(t, r)))
    bcs = [u(0, r) ~ -r * (r - 1) * sin(r),
        Dr(u(t, 0)) ~ 0.0, u(t, 1) ~ sin(1)]

    domains = [t ∈ Interval(0.0, 1.0),
        r ∈ Interval(0.0, 1.0)]

    @named pdesys = PDESystem(eq, bcs, domains, [t, r], [u(t, r)])
    discretization = MOLFiniteDifference([r => 0.1], t)
    prob = discretize(pdesys, discretization) # This gives an ODEProblem since it's time-dependent
    sol = solve(prob, Tsit5())
end

@testset "RHS = 0" begin
    @parameters x, t
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dx2 = Differential(x)^2
    Dx3 = Differential(x)^3
    Dx4 = Differential(x)^4

    α = 1
    β = 4
    γ = 1
    eq = Dt(u(x, t)) + u(x, t) * Dx(u(x, t)) + α * Dx2(u(x, t)) + β * Dx3(u(x, t)) + γ * Dx4(u(x, t)) ~ 0

    du(x, t; z=-x / 2 + t) = 15 / 2 * (tanh(z) + 1) * (3 * tanh(z) - 1) * sech(z)^2

    bcs = [u(x, 0) ~ x^2,
        Dx(u(-10, t)) ~ du(-10, t),
        Dx(u(10, t)) ~ du(10, t)]

    # Space and time domains
    domains = [x ∈ Interval(-10.0, 10.0),
        t ∈ Interval(0.0, 0.5)]
    # Discretization
    dx = 0.4
    dt = 0.2

    discretization = MOLFiniteDifference([x => dx], t; approx_order=4)
    @named pdesys = PDESystem(eq, bcs, domains, [x, t], [u(x, t)])
    prob = discretize(pdesys, discretization)
end

@testset "Wave Equation" begin
    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    Dt = Differential(t)             # Required in ICs and equation
    Dtt = Differential(t)^2           # required in equation
    Dxxxx = Differential(x)^4       # required in equation
    Dxx = Differential(x)^2         # required in BCs
    # some parameters
    EI = 291.6667
    m = 1.3850
    c = 1.0
    p = 0.01
    L = 2.0

    # 1D PDE and boundaru conditions
    eq = m * Dtt(u(t, x)) + c * Dt(u(t, x)) + EI * Dxxxx(u(t, x)) ~ 0

    ic_bc = [u(0, x) ~ (p * x * (x^3 + L^3 - 2 * L * x^2) / (24 * EI)), #for all 0 < u < L
        Dt(u(0, x)) ~ 0.0,        # for all 0 < u < L
        u(t, 0) ~ 0.0,            # for all t > 0,, displacement zero at u=0
        u(t, 2) ~ 0.0,            # for all t > 0,, displacement zero at u=L
        Dxx(u(t, 0)) ~ 0.0,       # for all t > 0,, curvature zero at u=0
        Dxx(u(t, 2)) ~ 0.0]       # for all t > 0,, curvature zero at u=L

    # Space and time domains
    domains = [t ∈ Interval(0.0, 2.0),
        x ∈ Interval(0.0, L)]

    dt = 0.1   # dt related to saving the data.. not actual dt
    # PDE sustem
    @named pdesys = PDESystem(eq, ic_bc, domains, [t, x], [u(t, x)])

    # Method of lines discretization
    dx = 0.1
    order = 2
    discretization = MOLFiniteDifference([x => dx, t => dt], approx_order=order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys, discretization)

    # Solve the ODE problem
    sol = NonlinearSolve.solve(prob, NewtonRaphson())
end

@testset "Rearranged Robin" begin
    @parameters x t

    @variables c(..)

    ∂t = Differential(t)
    ∂x = Differential(x)
    ∂²x = Differential(x)^2

    D₀ = 1.0
    R = 0.5
    cₑ = 2.5
    ℓ = 2.0
    α = 1 / 3

    diff_eq = ∂t(c(x, t)) ~ ∂x((D₀ + α * c(x, t)) * ∂x(c(x, t)))

    bcs = [c(x, 0) ~ 0.0,         # initial condition
        R * (D₀ + α * c(0, t)) * ∂x(c(0, t)) ~ c(0, t) - cₑ, # Robin
        ∂x(c(ℓ, t)) ~ 0.0]   # no flux

    domains = [t ∈ Interval(0.0, 10.0),
        x ∈ Interval(0.0, ℓ)]

    @named pdesys = PDESystem(diff_eq, bcs, domains, [t, x], [c(x, t)])

    # Method of lines discretization
    dx = 0.1
    order = 2
    discretization = MOLFiniteDifference([x => dx], t, approx_order=order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys, discretization)

    sol = solve(prob, Rodas4())
end

@testset "Array u" begin
    # Dependencies
    N = 6 # number of dependent variables

    # Variables, parameters, and derivatives
    @parameters x
    @variables u[1:N](..)
    Dx = Differential(x)
    Dxx = Differential(x)^2

    # Domain edges
    x_min = 0.0
    x_max = 1.0

    # Discretization parameters
    dx = 0.1
    order = 2

    #u = collect(u)

    # Equations
    eqs = Vector{ModelingToolkit.Equation}(undef, N)
    for i = 1:N
        eqs[i] = Dxx(u[i](x)) ~ u[i](x)
    end

    # Initial and boundary conditions
    bcs = Vector{ModelingToolkit.Equation}(undef, 2 * N)
    for i = 1:N
        bcs[i] = Dx(u[i](x_min)) ~ 0.0
    end

    for i = 1:N
        bcs[i+N] = u[i](x_max) ~ rand()
    end

    # Space and time domains
    domains = [x ∈ Interval(x_min, x_max)]

    # PDE system
    @named pdesys = PDESystem(eqs, bcs, domains, [x], collect([u[i](x) for i = 1:N]))

    # Method of lines discretization
    discretization = MOLFiniteDifference([x => dx], nothing, approx_order=order)
    prob = ModelingToolkit.discretize(pdesys, discretization)

    # # Solution of the ODE system
    sol = NonlinearSolve.solve(prob, NewtonRaphson())
end

@testset "Non-Uniform Grid" begin
    @parameters x y t
    @variables u(..) v(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dy = Differential(y)
    Dxx = Differential(x)^2
    Dyy = Differential(y)^2

    ∇²(u) = Dxx(u) + Dyy(u)

    brusselator_f(x, y, t) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.0

    x_min = y_min = t_min = 0.0
    x_max = y_max = 1.0
    t_max = 11.5

    α = 10.0

    u0(x, y, t) = 22(y * (1 - y))^(3 / 2)
    v0(x, y, t) = 27(x * (1 - x))^(3 / 2)

    eq = [Dt(u(x, y, t)) ~ 1.0 + Dy(v(x, y, t)) * u(x, y, t)^2 - 4.4 * u(x, y, t) + α * ∇²(u(x, y, t)) + brusselator_f(x, y, t),
        Dt(v(x, y, t)) ~ -3.4 * Dx(u(x, y, t)) - v(x, y, t) * u(x, y, t)^2 + α * ∇²(v(x, y, t))]

    domains = [x ∈ Interval(x_min, x_max),
        y ∈ Interval(y_min, y_max),
        t ∈ Interval(t_min, t_max)]

    bcs = [u(x, y, 0) ~ u0(x, y, 0),
        u(0, y, t) ~ 0.0,
        u(1, y, t) ~ 0.0,
        u(x, 0, t) ~ 0.0,
        u(x, 1, t) ~ 0.0,
        v(x, y, 0) ~ v0(x, y, 0),
        v(0, y, t) ~ 0.0,
        v(1, y, t) ~ 0.0,
        v(x, 0, t) ~ 0.0,
        v(x, 1, t) ~ 0.0]

    @named pdesys = PDESystem(eq, bcs, domains, [x, y, t], [u(x, y, t), v(x, y, t)])

    # Method of lines discretization
    N = 5

    dx = [0.0, cumsum(rand(0.05:0.01:0.2, 5))...]
    if dx[end] < 1.0
        push!(dx, 1.0)
    end

    dy = [0.0, cumsum(rand(0.05:0.01:0.2, 5))...]
    if dy[end] < 1.0
        push!(dy, 1.0)
    end

    order = 2

    discretization = MOLFiniteDifference([x => dx, y => dy], t, approx_order=order)

    prob = discretize(pdesys, discretization)

    sol = solve(prob, TRBDF2())

end

@testset "2D variable connected to 1D variable at boundary #33" begin
    @parameters t x r
    @variables u(..) v(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dxx = Differential(x)^2
    Dr = Differential(r)
    Drr = Differential(r)^2

    s = u(t, x) + v(t, x, 1)

    eqs = [Dt(u(t, x)) ~ Dxx(u(t, x)) + s,
        Dt(v(t, x, r)) ~ Drr(v(t, x, r))]
    bcs = [u(0, x) ~ 1,
        v(0, x, r) ~ 1,
        Dx(u(t, 0)) ~ 0.0, Dx(u(t, 1)) ~ 0.0,
        Dr(v(t, x, 0)) ~ 0.0, Dr(v(t, x, 1)) ~ s]

    domains = [t ∈ Interval(0.0, 1.0),
        x ∈ Interval(0.0, 1.0),
        r ∈ Interval(0.0, 1.0)]

    @named pdesys = PDESystem(eqs, bcs, domains, [t, x, r], [u(t, x), v(t, x, r)])

    # Method of lines discretization
    dx = 0.1
    dr = 0.1
    order = 2
    discretization = MOLFiniteDifference([x => dx, r => dr], t, approx_order=order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys, discretization)

    sol = solve(prob, Tsit5())
end

@testset "Testing discretization of varied systems" begin
    @parameters x t

    @variables c(..)

    ∂t = Differential(t)
    ∂x = Differential(x)
    ∂²x = Differential(x)^2

    D₀ = 1.5
    α = 0.15
    χ = 1.2
    R = 0.1
    cₑ = 2.0
    ℓ = 1.0
    Δx = 0.1

    bcs = [
        # initial condition
        c(x, 0) ~ 0.0,
        # Robin BC
        ∂x(c(0.0, t)) / (1 + exp(α * (c(0.0, t) - χ))) * R * D₀ + cₑ - c(0.0, t) ~ 0.0,
        # no flux BC
        ∂x(c(ℓ, t)) ~ 0.0]

    # define space-time plane
    domains = [x ∈ Interval(0.0, ℓ), t ∈ Interval(0.0, 5.0)]

    @testset "Test 01: ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))" begin
        D = D₀ / (1.0 + exp(α * (c(x, t) - χ)))
        diff_eq = ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))
        @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)])
        discretization = MOLFiniteDifference([x => Δx], t)
    end

    @testset "Test 02: ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))" begin
        D = 1.0 / (1.0 + exp(α * (c(x, t) - χ)))
        diff_eq = ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))
        @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)])
        discretization = MOLFiniteDifference([x => Δx], t)
    end

    @testset "Test 03: ∂t(c(x, t)) ~ ∂x(1.0 / (1.0/D₀ + exp(α * (c(x, t) - χ))/D₀) * ∂x(c(x, t)))" begin
        diff_eq = ∂t(c(x, t)) ~ ∂x(1.0 / (1.0 / D₀ + exp(α * (c(x, t) - χ)) / D₀) * ∂x(c(x, t)))
        @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)])
        discretization = MOLFiniteDifference([x => Δx], t)
    end

    @testset "Test 04: ∂t(c(x, t)) ~ ∂x(D₀ / (1.0 + exp(α * (c(x, t) - χ))) * ∂x(c(x, t)))" begin
        diff_eq = ∂t(c(x, t)) ~ ∂x(D₀ / (1.0 + exp(α * (c(x, t) - χ))) * ∂x(c(x, t)))
        @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)])
        discretization = MOLFiniteDifference([x => Δx], t)
    end

    @testset "Test 05: ∂t(c(x, t)) ~ ∂x(1/x * ∂x(c(x, t)))" begin
        diff_eq = ∂t(c(x, t)) ~ ∂x(1 / x * ∂x(c(x, t)))
        @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)])
        discretization = MOLFiniteDifference([x => Δx], t)
    end

    @testset "Test 06: ∂t(c(x, t)) ~ ∂x(x*∂x(c(x, t)))/c(x,t)" begin
        diff_eq = ∂t(c(x, t)) ~ ∂x(x * ∂x(c(x, t))) / c(x, t)
        @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)])
        discretization = MOLFiniteDifference([x => Δx], t)
    end

    @testset "Test 07: ∂t(c(x, t)) ~ ∂x(1/(1+c(x,t)) ∂x(c(x, t)))" begin
        diff_eq = ∂t(c(x, t)) ~ ∂x(1 / (1 + c(x, t)) * ∂x(c(x, t)))
        @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)])
        discretization = MOLFiniteDifference([x => Δx], t)
    end

    @testset "Test 08: ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x,t) * ∂x(c(x, t)))" begin
        diff_eq = ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x, t) * ∂x(c(x, t)))
        @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)])
        discretization = MOLFiniteDifference([x => Δx], t)
    end

    @testset "Test 09: ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x,t) * ∂x(c(x, t)))/(1+c(x,t))" begin
        diff_eq = c(x, t) * ∂x(c(x, t) * ∂x(c(x, t))) / (1 + c(x, t))
        @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)])
        discretization = MOLFiniteDifference([x => Δx], t)
    end

    @testset "Test 10: ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x,t) * ∂x(c(x, t)))/(1+c(x,t))" begin
        diff_eq = c(x, t) * ∂x(c(x, t)^(-1) * ∂x(c(x, t)))
        @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)])
        discretization = MOLFiniteDifference([x => Δx], t)
    end

    @testset "Test 11: ∂t(c(x, t)) ~ ∂x(1/(1+c(x,t)^2) ∂x(c(x, t)))" begin
        diff_eq = ∂t(c(x, t)) ~ ∂x(1 / (1 + c(x, t)^2) * ∂x(c(x, t)))
        @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)])
        discretization = MOLFiniteDifference([x => Δx], t)
    end

end