text:
using Scratch, Test, Dates, Pkg
include("utils.jl")

# Set to true for verbose Pkg output
const verbose = false
global const pkgio = verbose ? stderr : (VERSION < v"1.6.0-DEV.254" ? mktemp()[2] : devnull)

@testset "Scratch Space Basics" begin
    # Run everything in a separate depot, so that we can test GC'ing and whatnot
    temp_pkg_dir() do project_dir
        # Create a global scratch space, ensure it exists and is writable
        dir = get_scratch!("test")
        @test isdir(dir)
        @test startswith(dir, scratch_dir())
        touch(joinpath(dir, "foo"))
        @test readdir(dir) == ["foo"]

        # Test that this created a `scratch_usage.toml` file, and that accessing it
        # again does not increase the size of the scratch_usage.toml file, since we
        # only mark usage once every so often per julia session.
        usage_path = Scratch.usage_toml()
        @test isfile(usage_path)
        size = filesize(usage_path)
        dir = get_scratch!("test")
        @test size == filesize(usage_path)

        # But accessing it from a new Julia instance WILL increase its size:
        code = "import Scratch; Scratch.get_scratch!(\"test\")"
        run(setenv(
            `$(Base.julia_cmd()) --project=$(dirname(@__DIR__)) -e $code`,
            "JULIA_DEPOT_PATH" => first(Base.DEPOT_PATH),
        ))
        @test size < filesize(usage_path)

        # Delete the scratch space, ensure it's gone.
        delete_scratch!("test")
        @test !isdir(dir)

        # Key verification
        @test isdir(get_scratch!("abcABC123._-"))
        @test_throws ArgumentError get_scratch!("")
        @test_throws ArgumentError get_scratch!("hello/world")
        @test_throws ArgumentError get_scratch!("hello\\world")
    end
end

@testset "Scratch Space Namespacing" begin
    temp_pkg_dir() do project_dir
        # ScratchUsage UUID
        su_uuid = "93485645-17f1-6f3b-45bc-419db53815ea"
        # The UUID that gets used when no good UUIDs are available
        global_uuid = string(Base.UUID(UInt128(0)))

        # Touch the spaces of a ScratchUsage v1.0.0
        install_test_ScratchUsage(project_dir, v"1.0.0")

        # Ensure that the files were created for v1.0.0
        @test isfile(scratch_dir(su_uuid, "1.0.0", "ScratchUsage-1.0.0"))
        @test length(readdir(scratch_dir(su_uuid, "1.0.0"))) == 1
        @test isfile(scratch_dir(su_uuid, "1", "ScratchUsage-1.0.0"))
        @test length(readdir(scratch_dir(su_uuid, "1"))) == 1
        @test isfile(scratch_dir(global_uuid, "GlobalSpace", "ScratchUsage-1.0.0"))
        @test length(readdir(scratch_dir(global_uuid, "GlobalSpace"))) == 1

        # Next, do the same but for more versions
        install_test_ScratchUsage(project_dir, v"1.1.0")
        install_test_ScratchUsage(project_dir, v"2.0.0")

        # Check the spaces were shared when they should have been, and not when they shouldn't
        @test isfile(scratch_dir(su_uuid, "1.0.0", "ScratchUsage-1.0.0"))
        @test length(readdir(scratch_dir(su_uuid, "1.0.0"))) == 1
        @test isfile(scratch_dir(su_uuid, "1.1.0", "ScratchUsage-1.1.0"))
        @test length(readdir(scratch_dir(su_uuid, "1.1.0"))) == 1
        @test isfile(scratch_dir(su_uuid, "2.0.0", "ScratchUsage-2.0.0"))
        @test length(readdir(scratch_dir(su_uuid, "2.0.0"))) == 1
        @test isfile(scratch_dir(su_uuid, "1", "ScratchUsage-1.0.0"))
        @test isfile(scratch_dir(su_uuid, "1", "ScratchUsage-1.1.0"))
        @test length(readdir(scratch_dir(su_uuid, "1"))) == 2
        @test isfile(scratch_dir(su_uuid, "2", "ScratchUsage-2.0.0"))
        @test length(readdir(scratch_dir(su_uuid, "2"))) == 1
        @test isfile(scratch_dir(global_uuid, "GlobalSpace", "ScratchUsage-1.0.0"))
        @test isfile(scratch_dir(global_uuid, "GlobalSpace", "ScratchUsage-1.1.0"))
        @test isfile(scratch_dir(global_uuid, "GlobalSpace", "ScratchUsage-2.0.0"))
        @test length(readdir(scratch_dir(global_uuid, "GlobalSpace"))) == 3

        clear_scratchspaces!(Base.UUID(su_uuid))
        @test !isdir(scratch_dir(su_uuid))

        # UUID lookup for active project when running in Main
        ## Project.toml without UUID
        project = temp_project_file()
        with_active_project(project) do
            @test (@__MODULE__) == Main
            ## Project.toml without UUID
            path = @get_scratch!("project-no-uuid")
            @test isdir(path)
            @test path == scratch_dir(global_uuid, "project-no-uuid")
            @test path === get_scratch!(@__MODULE__, "project-no-uuid")
        end
        ## Project.toml with UUID
        project_uuid = Base.UUID("69386cca-e009-4a96-a0ae-829213699cfc")
        project = temp_project_file(project_uuid)
        with_active_project(project) do
            path = @get_scratch!("project-uuid")
            @test isdir(path)
            @test path == scratch_dir(string(project_uuid), "project-uuid")
            @test path === get_scratch!(@__MODULE__, "project-uuid")
        end # do

        # Cross-package scratch usage: Test that the scratch space is namespaced
        # to the other package, but tracked to me when it comes to lifecycling
        # the other package uuid, but the scratchspace is
        other_uuid = Base.UUID("6dc9890c-246d-42f7-b07c-3ce39ca50d56")
        ## Activate test-project from above again
        with_active_project(project) do
            path = get_scratch!(other_uuid, "hello-there", project_uuid)
            ## Test that the path is namespaced to other_uuid, but lifecycled with me
            @test path == scratch_dir(string(other_uuid), "hello-there")
            usage_path = usage_path = Scratch.usage_toml()
            usage = Pkg.TOML.parsefile(usage_path)
            @test project ∈ usage[path][1]["parent_projects"]
        end

        # Internal scratch access tracking
        empty!(Scratch.scratch_access_timers)
        usage_path = usage_path = Scratch.usage_toml()
        project2 = temp_project_file(other_uuid)
        with_active_project(project) do
            ## This should track project_uuid as the user
            path = get_scratch!(project_uuid, "general-kenobi")
            @test project_uuid ∈ first.(keys(Scratch.scratch_access_timers))
            usage = Pkg.TOML.parsefile(usage_path)
            @test any(project ∈ record["parent_projects"] for record in usage[path])
        end
        with_active_project(project2) do
            ## Reaching for the same space as another owner should (i) track this
            ## usage and (ii) write to scratch_usage.toml since this belongs to
            ## another project file
            path = get_scratch!(project_uuid, "general-kenobi", other_uuid)
            @test other_uuid ∈ first.(keys(Scratch.scratch_access_timers))
            usage = Pkg.TOML.parsefile(usage_path)
            @test any(project2 ∈ record["parent_projects"] for record in usage[path])
        end
        ## Deleting the path should remove both UUIDs from the timers
        delete_scratch!(project_uuid, "general-kenobi")
        @test project_uuid ∉ first.(keys(Scratch.scratch_access_timers))
        @test other_uuid ∉ first.(keys(Scratch.scratch_access_timers))
    end
end

# Run GC tests only on Julia >1.6
@testset "Scratch Space Lifecycling" begin; if VERSION >= v"1.6.0-DEV.676"
    temp_pkg_dir() do project_dir
        # First, install ScratchUsage
        su_uuid = "93485645-17f1-6f3b-45bc-419db53815ea"
        global_uuid = string(Base.UUID(UInt128(0)))
        install_test_ScratchUsage(project_dir, v"1.0.0")

        # Ensure that a few files were created
        @test isfile(scratch_dir(su_uuid, "1.0.0", "ScratchUsage-1.0.0"))
        @test length(readdir(scratch_dir(su_uuid, "1.0.0"))) == 1
        @test isfile(scratch_dir(global_uuid, "GlobalSpace", "ScratchUsage-1.0.0"))
        @test length(readdir(scratch_dir(global_uuid, "GlobalSpace"))) == 1

        # Test that a gc() doesn't remove anything, and that there is no orphanage
        Pkg.gc(; io=pkgio)
        orphaned_path = joinpath(first(Base.DEPOT_PATH), "logs", "orphaned.toml")
        @test isfile(scratch_dir(su_uuid, "1.0.0", "ScratchUsage-1.0.0"))
        @test isfile(scratch_dir(global_uuid, "GlobalSpace", "ScratchUsage-1.0.0"))
        @test !isfile(orphaned_path)

        # Remove ScrachUsage, which causes the package (but not the scratch dirs)
        # to move to the orphanage
        Pkg.rm("ScratchUsage"; io=pkgio)
        rm(joinpath(project_dir, "ScratchUsage"); force=true, recursive=true)
        Pkg.gc(; io=pkgio)

        @test isfile(scratch_dir(su_uuid, "1.0.0", "ScratchUsage-1.0.0"))
        @test isfile(scratch_dir(global_uuid, "GlobalSpace", "ScratchUsage-1.0.0"))
        @test isfile(orphaned_path)
        orphanage = Pkg.TOML.parse(String(read(orphaned_path)))
        @test haskey(orphanage, scratch_dir(su_uuid, "1.0.0"))
        @test haskey(orphanage, scratch_dir(su_uuid, "1"))
        @test !haskey(orphanage, scratch_dir(global_uuid, "GlobalSpace"))

        # Run a GC, forcing collection to ensure that everything in the SpaceUsage
        # namespace gets removed (but still appears in the orphanage)
        sleep(0.2)
        Pkg.gc(;collect_delay=Millisecond(100), io=pkgio)
        @test !isdir(scratch_dir(su_uuid))
        @test isdir(scratch_dir(global_uuid, "GlobalSpace"))
        orphanage = Pkg.TOML.parse(String(read(orphaned_path)))
        @test haskey(orphanage, scratch_dir(su_uuid, "1.0.0"))
        @test haskey(orphanage, scratch_dir(su_uuid, "1"))
        @test !haskey(orphanage, scratch_dir(global_uuid, "GlobalSpace"))
    end
end end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/09/2014

module SimplePilotResponseImpl

export
    initialize,
    update,

    updatePilotResponse,

    SimplePilotResponse,
    SimplePRResolutionAdvisory,
    SimplePRCommand

using AbstractPilotResponseImpl
using AbstractPilotResponseInterfaces
using CommonInterfaces

import CommonInterfaces.initialize
import CommonInterfaces.update
import AbstractPilotResponseInterfaces.updatePilotResponse

type SimplePRResolutionAdvisory
  h_d::Float64
end

type SimplePRCommand

  t::Float64
  v_d::Float64
  h_d::Float64
  psi_d::Float64

  logProb::Float64
end

type SimplePilotResponse <: AbstractPilotResponse

  b_CAS_activated::Bool
  RA::SimplePRResolutionAdvisory

  function SimplePilotResponse()

    obj = new()

    obj.b_CAS_activated = false

    return obj
  end
end

function updatePilotResponse(pr::SimplePilotResponse, command::SimplePRCommand, RA::Union{SimplePRResolutionAdvisory,Void})

  t, v_d, h_d, psi_d = command.t, command.v_d, command.h_d, command.psi_d

  if RA != nothing && pr.b_CAS_activated == false
    pr.RA = RA
    pr.b_CAS_activated = true
  end

  if pr.b_CAS_activated
    h_d = pr.RA.h_d
  end

  return SimplePRCommand(t, v_d, h_d, psi_d, 0.0)
end

update(pr::SimplePilotResponse, command, RA) = update(pr,convert(SimplePRCommand, command), convert(SimplePRResolutionAdvisory,RA))

update(pr::SimplePilotResponse, command::SimplePRCommand, RA::Union{SimplePRResolutionAdvisory,Void}) = updatePilotResponse(pr, command, RA)

function initialize(pr::SimplePilotResponse)

  pr.b_CAS_activated = false
end

end

ext: jl


text:
# search for counter examples to bapat sunder to show that it is hard
"""
    search_until_user_stop(search_function)

Runs `search_function` until user-stop (Ctrl+C).
"""
function search_until_user_stop(search_function)

    n_trials = 1

    try

        while true
            search_function()
            n_trials += 1
        end

    catch err

        if isa(err, InterruptException)
            print("no counter example found after ")
        end
    finally
        println("$n_trials trials")
    end
end

"""
    random_search_counter_example_bapat_sunder(;m,n,r, physical_H = true)

Brute-force search of counter-examples of rank `r`.
"""
function random_search_counter_example_bapat_sunder(;m,n,r, physical_H = true)

    S = rand_gram_matrix_rank(n,r)

    if physical_H
        U = rand_haar(m)

        input_state = [i <= n ? 1 : 0 for i in 1:m]
        partition = [i <= r ? 1 : 0 for i in 1:m]

        H = H_matrix(U, input_state, partition)
    else
        H = rand_gram_matrix_rank(n)
    end

    if violates_bapat_sunder(H,S)
        println("counter example found : ")
        println(H)
        println(S)
    end

end


text:
module SpecialMatricesModule
export computeMatrixKinetic, computeMatrixNuclearAttraction, computeMatrixOverlap
using ..IntegralsModule
using ..BasisModule
using ..Geometry

function computeMatrixKinetic(basis::GaussianBasis)
  return [KineticIntegral(cgb1,cgb2) for cgb1 in basis.contractedBFs, cgb2 in basis.contractedBFs]
end

function computeMatrixOverlap(basis::GaussianBasis)
  return [Overlap(cgb1,cgb2) for cgb1 in basis.contractedBFs, cgb2 in basis.contractedBFs]
end

function computeMatrixNuclearAttraction(basis::GaussianBasis,geo::Geometry)
  return [sum([NuclearAttractionIntegral(cgb1,cgb2,atom) for atom in geo.atoms]) for cgb1 in basis.contractedBFs, cgb2 in basis.contractedBFs]
end

end

max_stars_repo_name: tkelman/QuantumLab

max_issues_repo_name: tkelman/QuantumLab

max_forks_repo_name: tkelman/QuantumLab


text:
mutable struct Tikhonov{T} <: AbstractRegularization{T}
    α::T
end
Tikhonov() =  Tikhonov

get_regularization_matrix(regop::Tikhonov{T},N,ds) where {T} =   Matrix(I(N))

mutable struct SecondDerivative{T} <: AbstractRegularization{T}
    α::T
end
SecondDerivative()=SecondDerivative

function get_regularization_matrix(regop::SecondDerivative{T},N,ds) where {T}
    α=regop.α
    L=zeros( N+2,N)
    for j in 3:N
        L[j,j]= 1/(ds[j]*(ds[j]+ds[j-1]))#
        L[j,j-1]= -2*1/(ds[j-1]*(ds[j]))#
        L[j,j-2]= 1/(ds[j-1]*(ds[j]+ds[j-1]))#
    end

    # if regop.lower_bc
         L[1,1]=1/ds[1]
         L[2,1]=-2*1/(ds[1]*(ds[2]))
         L[2,2]=1/ds[2]*(ds[1]+ds[2])
    # end
    # if regop.upper_bc
         L[N+1,N-1]=1/ds[N-1]*(ds[N]+ds[N-1])
         L[N+1,N]=-2*1/((ds[N]*ds[N-1]))
         L[N+2,N]=1/ds[N]
    # end
    return L
    #L[:,end] .=0  #Do not regularize on the y offset
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -2.9922365365158297 <= q <= 1.3959867501971397)

                     add_NL_constraint(m, :(swish(-0.7172866311917714 + 0.8530094465195428*swish(0.49975691401319544 + 0.817850285997455*swish(-0.9044530867886094 + -0.3504940249833268*swish(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + -0.925387475816613*swish(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))) + -0.15011340339301915*swish(-0.49866226232151334 + 0.31623027699891626*swish(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + 0.1823838704363938*swish(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4])))) + -0.2302663292873568*swish(0.6056562038214484 + 0.017530899400252764*swish(-0.9044530867886094 + -0.3504940249833268*swish(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + -0.925387475816613*swish(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))) + 0.9141090581590978*swish(-0.49866226232151334 + 0.31623027699891626*swish(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + 0.1823838704363938*swish(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))) + swish(0.11802092343800519 + -0.23198020992542823*swish(0.49975691401319544 + 0.817850285997455*swish(-0.9044530867886094 + -0.3504940249833268*swish(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + -0.925387475816613*swish(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))) + -0.15011340339301915*swish(-0.49866226232151334 + 0.31623027699891626*swish(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + 0.1823838704363938*swish(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4])))) + 0.3694431087223533*swish(0.6056562038214484 + 0.017530899400252764*swish(-0.9044530867886094 + -0.3504940249833268*swish(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + -0.925387475816613*swish(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))) + 0.9141090581590978*swish(-0.49866226232151334 + 0.31623027699891626*swish(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + 0.1823838704363938*swish(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
export treedepth

"""
    treedepth(jlt::AbstractJLBoostTree)

Return the depth of the `jlt` tree
"""
function treedepth(jlt::AbstractJLBoostTree)
    if has_children(jlt)
        return 1 + maximum(treedepth, children(jlt))
    else
        return 0
    end
end


text:
using BandwidthBenchmark
using Test

@testset "BandwidthBenchmark.jl" begin
    # Write your tests here.
end

max_stars_repo_name: carstenbauer/BandwidthBenchmark

max_issues_repo_name: carstenbauer/BandwidthBenchmark

max_forks_repo_name: carstenbauer/BandwidthBenchmark


text:
#---------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Model of a stoichiometric reactor
#*----------------------------------------------------------------------
#*
#*   Description:
#*       Modeling of a reactor based on a stoichiometric approach.
#*
#*   Assumptions:
#*		* single- and two-phases involved
#*		* steady-state
#*
#*	Specify:
#*		* inlet stream
#*		* extent of reactions or
#*		* conversion of a key component
#*
#*----------------------------------------------------------------------
#* Author: Rodolfo Rodrigues
#* $Id$
#*--------------------------------------------------------------------
#---------------------------------------------------------------------
#*	only vapour-phase
#*--------------------------------------------------------------------
type stoic_vap
	stoic_vap()=begin
		new(
			tank_vap(),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of reactions",
				:Default=>1
			)),
			fill(DanaReal(Dict{Symbol,Any}(
				:Brief=>"Stoichiometric matrix",
				:Symbol=>"\\nu"
			)),(NComp,NReac)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Outlet stream",
				:PosX=>1,
				:PosY=>1,
				:Symbol=>"_{out}"
			)),
			fill(reaction_mol(Dict{Symbol,Any}(
				:Brief=>"Overall component rate of reaction"
			)),(NComp)),
			fill(DanaReal(Dict{Symbol,Any}(
				:Brief=>"Fractional conversion of component",
				:Symbol=>"X",
				:Default=>0
			)),(NComp)),
			[
				:(Outlet.F*Outlet.z = _base_1._base_1.Outletm.F*_base_1._base_1.Outletm.z + rate*_base_1._base_1.Tank.V),
				:(Outlet.P = _base_1._base_1.Outletm.P),
				:(Outlet.F*Outlet.h = _base_1._base_1.Outletm.F*_base_1._base_1.Outletm.h),
				:(Outlet.F = _base_1._base_1.Outletm.F),
				:(Outlet.F*Outlet.z([1:_base_1._base_1.NComp]) = _base_1._base_1.Outletm.F*_base_1._base_1.Outletm.z([1:_base_1._base_1.NComp])*(1 - conv([1:_base_1._base_1.NComp]))),
				:(conv([1:_base_1._base_1.NComp]) = 1),
				:(conv([1:_base_1._base_1.NComp]) = 0),
			],
			[
				"Outlet stream","Mechanical equilibrium","Energy balance","Steady-state","Molar conversion","Molar conversion","Molar conversion",
			],
			[:NReac,:stoic,],
			[:Outlet,:rate,:conv,]
		)
	end
	_base_1::tank_vap
	NReac::DanaInteger
	stoic::Array{DanaReal}
	Outlet::vapour_stream
	rate::Array{reaction_mol}
	conv::Array{DanaReal}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export stoic_vap
function setEquationFlow(in::stoic_vap)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	if (Outletm.z(i) > 1e-16)
		addEquation(5)
	else
		if (Outlet.z(i) > 0)
			addEquation(6)
			# ?

		else
			addEquation(7)
			# ?

		end
	end
end
function atributes(in::stoic_vap,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Brief]="Basic model for a vapour-phase stoichiometric CSTR"
	fields[:Info]="
== Assumptions ==
* only vapour-phase
* steady-state
"
	drive!(fields,_)
	return fields
end
stoic_vap(_::Dict{Symbol,Any})=begin
	newModel=stoic_vap()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using Test
using MultiPoly

include("mpolycalculus.jl")
include("mpoly.jl")
include("mpolyarithmetic.jl")
include("mpolyprinting.jl")


text:
"""
    function clean(a)

Cleans the complex print of an expression to make it more clear to see.

Code reused from comments on [Julia](https://discourse.julialang.org/t/code-generation-unnecessary-comment-lines-when-using-quote/398/2)
"""
function clean(ex::Expr)
    filter!(ex.args) do e
        isa(e, LineNumberNode) && return false
        if isa(e, Expr)
            (e::Expr).head === :line && return false
            clean(e::Expr)
        end
        return true
    end
    return ex
end

function clean2(ex::Expr)
    filter!(e->typeof(e)!=LineNumberNode,ex.args)
    for i in ex.args
        if typeof(i) == Expr
            clean2(i)
        end
    end
    return ex
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
"""
    euler_totient(n::Int64)

Compute the positive integers up to a given integer n that are relatively prime to n; also
known as Euler's totient. For more information see: [https://en.wikipedia.org/wiki/Euler%27s_totient_function](https://en.wikipedia.org/wiki/Euler%27s_totient_function)

# Arguments
- `n::Int64`: Number of elements

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> ClassicAlgorithmsCollections.euler_totient(20)
[1, 1, 2, 2, 4, 2, 6, 4, 6, 4, 10, 4, 12, 6, 8, 8, 16, 6, 18, 8]
```
"""
function euler_totient(n::Int64)

    result = collect(Int64, 1:n)

    for i in 1:n
        m = i
        p = 2
        while m >= p^2
            if (m % p == 0)
                while (m % p == 0)
                    m = convert(Int64, m / p)
                end
                result[i] -= convert(Int64, result[i] / p)
            end
            p += 1
        end
        if m > 1
            result[i] -= convert(Int64, result[i] / m)
        end
    end
    return result
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
"""
Bridges `CP.SymmetricAllDifferent` to `CP.AllDifferent` and `CP.Inverse`.
"""
struct SymmetricAllDifferent2AllDifferentInverseBridge{T} <: MOIBC.AbstractBridge
    con_all_diff::MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.AllDifferent}
    con_inverse::MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Inverse}
end

function MOIBC.bridge_constraint(
    ::Type{SymmetricAllDifferent2AllDifferentInverseBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.SymmetricAllDifferent,
) where {T}
    return MOIBC.bridge_constraint(
        SymmetricAllDifferent2AllDifferentInverseBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{SymmetricAllDifferent2AllDifferentInverseBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.SymmetricAllDifferent,
) where {T}
    f_scalars = MOIU.scalarize(f)
    dim = MOI.output_dimension(f)

    con_all_diff = MOI.add_constraint(
        model,
        f,
        CP.AllDifferent(dim)
    )

    con_inverse = MOI.add_constraint(
        model,
        MOIU.vectorize([f_scalars..., f_scalars...]),
        CP.Inverse(dim)
    )

    return SymmetricAllDifferent2AllDifferentInverseBridge(con_all_diff, con_inverse)
end

function MOI.supports_constraint(
    ::Type{SymmetricAllDifferent2AllDifferentInverseBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.SymmetricAllDifferent},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{SymmetricAllDifferent2AllDifferentInverseBridge{T}}) where {T}
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(::Type{SymmetricAllDifferent2AllDifferentInverseBridge{T}}) where {T}
    return [
        (MOI.VectorAffineFunction{T}, CP.AllDifferent),
        (MOI.VectorAffineFunction{T}, CP.Inverse),
    ]
end

function MOI.get(
    ::SymmetricAllDifferent2AllDifferentInverseBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, CP.AllDifferent,
    },
) where {T}
    return 1
end

function MOI.get(
    ::SymmetricAllDifferent2AllDifferentInverseBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, CP.Inverse,
    },
) where {T}
    return 1
end

function MOI.get(
    b::SymmetricAllDifferent2AllDifferentInverseBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, CP.AllDifferent,
    },
) where {T}
    return [b.con_all_diff]
end

function MOI.get(
    b::SymmetricAllDifferent2AllDifferentInverseBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, CP.Inverse,
    },
) where {T}
    return [b.con_inverse]
end


text:
module OutlierDetection
    using OutlierDetectionInterface

    import MLJBase
    const MLJ = MLJBase

    # re-export from OutlierDetectionInterface
    export CLASS_NORMAL, CLASS_OUTLIER, DEFAULT_THRESHOLD

    export ScoreTransformer,
           ClassTransformer,
           scale_minmax,
           scale_unify,
           combine_mean,
           combine_median,
           combine_max,
           classify_percentile,
           to_categorical,
           to_univariate_finite,
           ProbabilisticDetector,
           DeterministicDetector

    # utilities
    include("normalization.jl")
    include("classification.jl")
    include("combination.jl")

    # extension
    include("mlj_helpers.jl")
    include("mlj_transformers.jl")
    include("mlj_wrappers.jl")
end


text:
export AiryU0

struct AiryU0 <: ContinuousMatrixDistribution
    beta::Real
end

function rand(d::AiryU0, trange::AbstractRange{T}) where T <: Real
    t_grid = collect(trange)
    dt = step(trange)

    B = rand_brownian(t_grid, dt)

    return make_U0(t_grid, dt, B, d.beta)
end

##

function rand_brownian(t_grid, dt)
    grid_len = length(t_grid)

    B_diff = randn(grid_len) * sqrt(abs(dt))
    cumsum(B_diff)
end

function make_U0(t_grid, dt, B, beta)
    U0_row = t_grid.^2/2 + B * sqrt(4/beta)
    U0_row .- U0_row'
end


text:
using QuerySQLite

using DataValues: DataValue
using Dates: Date, DateTime, format, Time
using Documenter: doctest
using Query
using QueryTables
using SQLite: DB, drop!, execute, Stmt
using Statistics: mean
using Test

doctest(QuerySQLite)

@testset "SQLite tutorial" begin

    filename = joinpath(@__DIR__, "Chinook_Sqlite.sqlite")
    database = Database(filename)
    database2 = Database(DB(filename))

    @test (database.Track |>
     @filter(_.AlbumId == 1) |>
     collect |>
     first).Name == "For Those About To Rock (We Salute You)"

    @test database.Track |>
    @map({_.TrackId, _.Name, _.Composer, _.UnitPrice}) |>
    collect |>
    first |>
    propertynames == (:TrackId, :Name, :Composer, :UnitPrice)

    @test (database2.Customer |>
    @map({_.City, _.Country}) |>
    @orderby(_.Country) |>
    DataTable).Country[1] == "Argentina"

    @test database.Customer |>
    @map({_.City}) |>
    @unique() |>
    collect |>
    length == 53

    @test database.Track |>
    @map({_.TrackId, _.Name}) |>
    @take(10) |>
    collect |>
    length == 10

    @test first((database.Track |>
    @map({_.TrackId, _.Name}) |>
    @take(10) |>
    @drop(10) |>
    DataTable).Name)  == "C.O.D."

    @test first((
    @from i in database.Track begin
        @orderby descending(i.Bytes), i.Name
        @select {i.TrackId, i.Name, i.Bytes}
        @collect DataTable
    end
    ).Bytes) == 1059546140

    @test first((
    @from i in database.Track begin
        @orderby i.Bytes, descending(i.Name)
        @select {i.TrackId, i.Name, i.Bytes}
        @collect DataTable
    end
    ).Bytes) == 38747

    @test database.Artist |>
    @join(database.Album, _.ArtistId, _.ArtistId, {_.ArtistId, __.AlbumId}) |>
    DataTable |>
    length == 347

    @test (database.Track |>
    @map({_.Name, _.Milliseconds, _.Bytes, _.AlbumId}) |>
    @filter(_.AlbumId == 1) |>
    collect |>
    first).Name == "For Those About To Rock (We Salute You)"

    group_by_row =
    database.Track |>
    @groupby(_.AlbumId) |>
    @map({
        AlbumId = key(_),
        length = length(_.AlbumId),
        sum = sum(_.Milliseconds),
        min = minimum(_.Milliseconds),
        max = maximum(_.Milliseconds),
        mean = mean(_.Milliseconds)
    }) |>
    collect |>
    first

    @test group_by_row.AlbumId == 1
    @test group_by_row.length == 10
    @test group_by_row.sum == 2400415
    @test group_by_row.min == 199836
    @test group_by_row.mean == 240041.5

    @test collect(
    database.Track |>
    @select(:AlbumId, :Composer) |>
    @mutate(bar = _.AlbumId * 2)
)[1].AlbumId == 1

end

@testset "Systematic tests" begin

    filename = joinpath(@__DIR__, "tmp", "test.sqlite")
    isfile(filename) && rm(filename)
    cp(joinpath(@__DIR__, "test.sqlite"), filename)

    connection = DB(filename)
    execute(Stmt(connection, """DROP TABLE IF EXISTS test"""))
    execute(Stmt(connection, """
    CREATE TABLE test (
        zero Int,
        one Int,
        negative_one Int,
        ab Text,
        null_column Int,
        A_code Int,
        a Text,
        b Text,
        a_space Text,
        a_underscore Text,
        a_wild Text,
        one_point_one_one Real,
        datetime_text Text,
        format Text
    )"""))
    execute(Stmt(connection, """
    INSERT INTO test VALUES(0, 1, -1, "ab", NULL, 65, "a", "b", " a ", "_a_", "a%", 1.11, "2019-12-08T11:09:00", "%Y-%m-%d %H:%M:%S")
"""))
    small = Database(connection)

    @test (small.test |>
    @groupby(_.zero) |>
    @map({
        join_test = join(_.ab)
    }) |>
    collect |>
    first).join_test == "ab"

    result =
    small.test |>
    @map({
        equals_test = _.zero == _.one,
        not_equals_test = _.zero != _.one,
        not_test = !(_.zero),
        and_test = _.zero & _.one,
        or_test = _.zero | _.one,
        times_test = _.zero * _.one,
        divide_test = _.zero / _.one,
        plus_test = _.zero + _.one,
        minus_test = _.one - _.zero,
        mod_test = _.zero % _.one,
        abs_test = abs(_.negative_one),
        in_test = _.zero in (0, 1),
        coalesce_test = coalesce(_.null_column, _.zero),
        if_else_test_1 = if_else(_.one, 1, 0),
        if_else_test_2 = if_else(1, _.one, 0),
        if_else_test_3 = if_else(0, 1, _.zero),
        if_else_test_4 = if_else(0, _.one, _.zero),
        if_else_test_5 = if_else(_.zero, 1, _.zero),
        if_else_test_6 = if_else(_.one, _.one, 0),
        if_else_test_7 = if_else(_.one, _.one, _.zero),
        ismissing_test = ismissing(_.null_column),
        lowercase_test = lowercase(_.ab),
        max_test = max(_.one, 0),
        min_test = min(_.zero, 1),
        occursin_test = occursin("a%", _.ab),
        occursin_test_2 = occursin(_.a_wild, "ab"),
        occursin_test_3 = occursin(_.a_wild, _.ab),
        uppercase_test = uppercase(_.ab),
        char_test = char(_.A_code),
        instr_test_1 = instr(_.ab, "b"),
        instr_test_2 = instr("ab", _.b),
        instr_test_3 = instr(_.ab, _.b),
        hex_test = hex(_.ab),
        strip_test = strip(_.a_space),
        strip_test_2 = strip(_.a_underscore, '_'),
        repr_test = repr(_.ab),
        replace_test = replace(_.ab, "b" => "a"),
        round_test_1 = round(_.one_point_one_one, digits = 1),
        round_test_2 = round(_.one_point_one_one),
        SubString_test_1 = SubString(_.ab, 2, 2),
        SubString_test_2 = SubString(_.ab, 2),
        random_test = rand(BySQL(_), Int),
        # TODO: fix
        # randstring_test = randstring(BySQL(_), 4),
        # randstring_test2 = randstring(_.one),
        format_test = format(_.datetime_text, "%Y-%m-%d %H:%M:%S"),
        format_test_2 = format("2019-12-08T11:09:00", _.format),
        format_test_3 = format(_.datetime_text, _.format),
        type_of_test = type_of(_.zero),
        convert_test = convert(Int, _.b),
        string_test = string(_.a, _.b)
    }) |>
    collect |>
    first

    @test result.equals_test == 0
    @test result.not_equals_test == 1
    @test result.not_test == 1
    @test result.and_test == 0
    @test result.or_test == 1
    @test result.times_test == 0
    @test result.divide_test == 0
    @test result.plus_test == 1
    @test result.minus_test == 1
    @test result.mod_test == 0
    @test result.abs_test == 1
    @test result.in_test == 1
    @test result.coalesce_test == 0
    @test result.if_else_test_1 == 1
    @test result.if_else_test_2 == 1
    @test result.if_else_test_3 == 0
    @test result.if_else_test_4 == 0
    @test result.if_else_test_5 == 0
    @test result.if_else_test_6 == 1
    @test result.if_else_test_7 == 1
    @test result.ismissing_test == 1
    @test result.max_test == 1
    @test result.min_test == 0
    @test result.occursin_test == 1
    @test result.occursin_test_2 == 1
    @test result.occursin_test_3 == 1
    @test result.uppercase_test == "AB"
    @test result.lowercase_test == "ab"
    @test result.char_test == "A"
    @test result.instr_test_1 == 2
    @test result.instr_test_2 == 2
    @test result.instr_test_3 == 2
    @test result.hex_test == "6162"
    @test result.strip_test == "a"
    @test result.strip_test_2 == "a"
    @test result.repr_test == "\'ab\'"
    @test result.replace_test == "aa"
    @test result.round_test_1 == 1.1
    @test result.round_test_2 == 1.0
    @test result.SubString_test_1 == "b"
    @test result.SubString_test_2 == "b"
    @test result.random_test isa DataValue{Int}
    @test_broken length(result.randomstring_test) == 4
    @test_broken length(result.randomstring_test2) == 1
    @test result.format_test == "2019-12-08 11:09:00"
    @test result.format_test_2 == "2019-12-08 11:09:00"
    @test result.format_test_3 == "2019-12-08 11:09:00"
    @test result.type_of_test == "integer"
    @test result.string_test == "ab"

    drop!(connection, "test")

    @test_throws ArgumentError Database("file.not_sqlite")

end


text:
module TestDoctest

import PublicAPI
using Documenter: doctest
using Test

function test()
    VERSION ≥ v"1.6" || return
    doctest(PublicAPI; manual = false)
end

end  # module


text:
abstract type PoissonBCs end

# PoissonBCs are named XYZ, where each of X, Y, and Z is either
# 'P' (for Periodic) or 'N' (for Neumann).
struct PPN <: PoissonBCs end
struct PNN <: PoissonBCs end

# Not yet supported:
#struct PPP <: PoissonBCs end
#struct NNN <: PoissonBCs end
#struct NPP <: PoissonBCs end
#struct PNP <: PoissonBCs end
#struct NPN <: PoissonBCs end
#struct NNP <: PoissonBCs end

poisson_bc_symbol(::BC) = :N
poisson_bc_symbol(::BC{<:Periodic}) = :P

"""
    PoissonBCs(bcs)

Returns the boundary conditions for the Poisson solver corresponding
to the model boundary conditions `bcs`.
"""
function PoissonBCs(bcs)
    # We assume that bounary conditions on all fields are
    # consistent with the boundary conditions on one side of ``u``.
    x = poisson_bc_symbol(bcs.u.x.left)
    y = poisson_bc_symbol(bcs.u.y.left)
    z = poisson_bc_symbol(bcs.u.z.left)

    return eval(Expr(:call, Symbol(x, y, z)))
end

PoissonBCs(model_bcs::ModelBoundaryConditions) = PoissonBCs(model_bcs.solution)

PoissonSolver(::CPU, pbcs::PoissonBCs, grid::AbstractGrid) = PoissonSolverCPU(pbcs, grid)
PoissonSolver(::GPU, pbcs::PoissonBCs, grid::AbstractGrid) = PoissonSolverGPU(pbcs, grid)

unpack_grid(grid::AbstractGrid) = grid.Nx, grid.Ny, grid.Nz, grid.Lx, grid.Ly, grid.Lz

"""
    ω(M, k)

Return the `M`th root of unity raised to the `k`th power.
"""
@inline ω(M, k) = exp(-2im*π*k/M)

"""
    λi(grid::AbstractGrid, ::PoissonBCs)

Return an Nx×1×1 array of eigenvalues satisfying the discrete form of Poisson's
equation with periodic boundary conditions in the x-dimension on `grid`.
"""
function λi(grid::AbstractGrid, ::PoissonBCs)
    Nx, Ny, Nz, Lx, Ly, Lz = unpack_grid(grid)
    is = reshape(1:Nx, Nx, 1, 1)
    @. (2sin((is-1)*π/Nx) / (Lx/Nx))^2
end

"""
    λj(grid::AbstractGrid, ::PPN)

Return an 1×Ny×1 array of eigenvalues satisfying the discrete form of Poisson's
equation with periodic boundary conditions in the y-dimension on `grid`.
"""
function λj(grid::AbstractGrid, ::PPN)
    Nx, Ny, Nz, Lx, Ly, Lz = unpack_grid(grid)
    js = reshape(1:Ny, 1, Ny, 1)
    @. (2sin((js-1)*π/Ny) / (Ly/Ny))^2
end

"""
    λj(grid::AbstractGrid, ::PNN)

Return an 1×Ny×1 array of eigenvalues satisfying the discrete form of Poisson's
equation with staggered Neumann boundary conditions in the y-dimension on `grid`.
"""
function λj(grid::AbstractGrid, ::PNN)
    Nx, Ny, Nz, Lx, Ly, Lz = unpack_grid(grid)
    js = reshape(1:Ny, 1, Ny, 1)
    @. (2sin((js-1)*π/(2Ny)) / (Ly/Ny))^2
end

"""
    λk(grid::AbstractGrid, ::PoissonBCs)

Return an 1×1×Nz array of eigenvalues satisfying the discrete form of Poisson's
equation with staggered Neumann boundary conditions in the y-dimension on `grid`.
"""
function λk(grid::AbstractGrid, ::PoissonBCs)
    Nx, Ny, Nz, Lx, Ly, Lz = unpack_grid(grid)
    ks = reshape(1:Nz, 1, 1, Nz)
    @. (2sin((ks-1)*π/(2Nz)) / (Lz/Nz))^2
end

"""
    plan_transforms(::PPN, A::Array; planner_flag=FFTW.PATIENT)

Return the transforms required to solve Poisson's equation with periodic
boundary conditions in the x,y dimensions and staggered Neumann boundary
conditions in the z-dimension.

Fast Fourier transforms (FFTs) are used in the x,y dimensions and real-to-real
discrete cosine transforms are used in the z dimension. Note that the DCT-II is
used for the DCT and the DCT-III for the IDCT which correspond to REDFT10 and
REDFT01 in FFTW.

They operatore on an array with the shape of `A`, which is needed to plan
efficient transforms. `A` will be mutated.

Note that the transforms returns operate on Arrays and so only work on CPUs.
"""
function plan_transforms(::PPN, A::Array, planner_flag=FFTW.PATIENT)
    FFT!  = FFTW.plan_fft!(A, [1, 2]; flags=planner_flag)
    IFFT! = FFTW.plan_ifft!(A, [1, 2]; flags=planner_flag)
    DCT!  = FFTW.plan_r2r!(A, FFTW.REDFT10, 3; flags=planner_flag)
    IDCT! = FFTW.plan_r2r!(A, FFTW.REDFT01, 3; flags=planner_flag)
    return FFT!, IFFT!, DCT!, IDCT!
end

"""
    plan_transforms(::PPN, A::Array; planer_flag=FFTW.PATIENT)

Similar to `plan_transforms(::PPN, ...)` but return the transforms required to
solve Poisson's equation with periodic boundary conditions in the x dimension
and staggered Neumann boundary conditions in the y,z dimensions.
"""
function plan_transforms(::PNN, A::Array, planner_flag=FFTW.PATIENT)
    FFT!  = FFTW.plan_fft!(A, 1; flags=planner_flag)
    IFFT! = FFTW.plan_ifft!(A, 1; flags=planner_flag)
    DCT!  = FFTW.plan_r2r!(A, FFTW.REDFT10, [2, 3]; flags=planner_flag)
    IDCT! = FFTW.plan_r2r!(A, FFTW.REDFT01, [2, 3]; flags=planner_flag)
    return FFT!, IFFT!, DCT!, IDCT!
end

struct PoissonSolverCPU{BC, AAR, AAC, FFTT, DCTT, IFFTT, IDCTT} <: AbstractPoissonSolver
        bcs :: BC
        kx² :: AAR
        ky² :: AAR
        kz² :: AAR
    storage :: AAC
       FFT! :: FFTT
       DCT! :: DCTT
      IFFT! :: IFFTT
      IDCT! :: IDCTT
end

function PoissonSolverCPU(pbcs::PoissonBCs, grid::AbstractGrid, planner_flag=FFTW.PATIENT)
    Nx, Ny, Nz, _ = unpack_grid(grid)

    # The eigenvalues of the discrete form of Poisson's equation correspond
    # to discrete wavenumbers so we call them k² to make the analogy with
    # solving Poisson's equation in spectral space.
    kx² = λi(grid, pbcs)
    ky² = λj(grid, pbcs)
    kz² = λk(grid, pbcs)

    # Storage for RHS and Fourier coefficients is hard-coded to be Float64
    # because of precision issues with Float32.
    # See https://github.com/climate-machine/Oceananigans.jl/issues/55
    storage = zeros(Complex{Float64}, grid.Nx, grid.Ny, grid.Nz)

    FFT!, IFFT!, DCT!, IDCT! = plan_transforms(pbcs, storage, planner_flag)

    PoissonSolverCPU(pbcs, kx², ky², kz², storage, FFT!, DCT!, IFFT!, IDCT!)
end

normalize_idct_output(::PPN, grid::AbstractGrid, ϕ::Array) = (@. ϕ = ϕ / (2grid.Nz))
normalize_idct_output(::PNN, grid::AbstractGrid, ϕ::Array) = (@. ϕ = ϕ / (4grid.Ny*grid.Nz))

"""
    solve_poisson_3d!(solver::PoissonSolverCPU, grid::RegularCartesianGrid)

Solve Poisson equation on a uniform staggered grid (Arakawa C-grid) with
appropriate boundary conditions as specified by `solver.bcs` using planned FFTs
and DCTs. The right-hand-side RHS is stored in solver.storage which the solver
mutates to produce the solution, so it will also be stored in solver.storage.

We should describe the algorithm in detail in the documentation.
"""
function solve_poisson_3d!(solver::PoissonSolverCPU, grid::RegularCartesianGrid)
    Nx, Ny, Nz, _ = unpack_grid(grid)
    kx², ky², kz² = solver.kx², solver.ky², solver.kz²

    # We can use the same storage for the RHS and the solution ϕ.
    RHS, ϕ = solver.storage, solver.storage

    solver.DCT! * RHS  # Calculate DCTᶻ(f) in place.
    solver.FFT! * RHS  # Calculate FFTˣʸ(f) in place.

    # Solve the discrete Poisson equation in spectral space. We are essentially
    # computing the Fourier coefficients of the solution from the Fourier
    # coefficients of the RHS.
    @. ϕ = -RHS / (kx² + ky² + kz²)

    # Setting DC component of the solution (the mean) to be zero. This is also
    # necessary because the source term to the Poisson equation has zero mean
    # and so the DC component comes out to be ∞.
    ϕ[1, 1, 1] = 0

    solver.IFFT! * ϕ  # Calculate IFFTˣʸ(ϕ) in place.
    solver.IDCT! * ϕ  # Calculate IDCTᶻ(ϕ) in place.

    # Must normalize by 2N for each dimension transformed via FFTW.REDFT.
    normalize_idct_output(solver.bcs, grid, ϕ)

    nothing
end

function plan_transforms(::PPN, A)
    FFT!      = plan_fft!(A, [1, 2])
    FFT_DCT!  = plan_fft!(A, 3)
    IFFT!     = plan_ifft!(A, [1, 2])
    IFFT_DCT! = plan_ifft!(A, 3)
    return FFT!, FFT_DCT!, IFFT!, IFFT_DCT!
end

function plan_transforms(::PNN, A)
    FFT!      = plan_fft!(A, 1)
    FFT_DCT!  = plan_fft!(A, [2, 3])
    IFFT!     = plan_ifft!(A, 1)
    IFFT_DCT! = plan_ifft!(A, [2, 3])
    return FFT!, FFT_DCT!, IFFT!, IFFT_DCT!
end

struct PoissonSolverGPU{BC, AAR, AAC, AAI, FFT, FFTD, IFFT, IFFTD} <: AbstractPoissonSolver
          bcs :: BC
          kx² :: AAR
          ky² :: AAR
          kz² :: AAR
       ω_4Ny⁺ :: AAC
       ω_4Ny⁻ :: AAC
       ω_4Nz⁺ :: AAC
       ω_4Nz⁻ :: AAC
     p_y_inds :: AAI
     p_z_inds :: AAI
     r_y_inds :: AAI
     r_z_inds :: AAI
         M_ky :: AAR
         M_kz :: AAR
      storage :: AAC
     storage2 :: AAC
         FFT! :: FFT
     FFT_DCT! :: FFTD
        IFFT! :: IFFT
    IFFT_DCT! :: IFFTD
end

function PoissonSolverGPU(pbcs::PoissonBCs, grid::AbstractGrid)
    Nx, Ny, Nz, _ = unpack_grid(grid)

    # The eigenvalues of the discrete form of Poisson's equation correspond
    # to discrete wavenumbers so we call them k² to make the analogy with
    # solving Poisson's equation in spectral space.
    kx² = λi(grid, pbcs) |> CuArray
    ky² = λj(grid, pbcs) |> CuArray
    kz² = λk(grid, pbcs) |> CuArray

    # Storage for RHS and Fourier coefficients is hard-coded to be Float64
    # because of precision issues with Float32.
    # See https://github.com/climate-machine/Oceananigans.jl/issues/55
    storage = zeros(Complex{Float64}, grid.Nx, grid.Ny, grid.Nz) |> CuArray

    # For solving the Poisson equation with PNN boundary conditions, we need a
    # second storage/buffer array to perform the 2D fast cosine transform. Maybe
    # we can get around this but for now we'll just use up a second array. It's
    # not needed for the PPN case so we just use the smallest 3D array we can
    # to keep the struct concretely typed.
    if pbcs == PPN()
        storage2 = zeros(Complex{Float64}, 1, 1, 1) |> CuArray
    elseif pbcs == PNN()
        storage2 = zeros(Complex{Float64}, grid.Nx, grid.Ny, grid.Nz) |> CuArray
    end

    ky⁺ = reshape(0:Ny-1,       1, Ny, 1)
    kz⁺ = reshape(0:Nz-1,       1, 1, Nz)
    ky⁻ = reshape(0:-1:-(Ny-1), 1, Ny, 1)
    kz⁻ = reshape(0:-1:-(Nz-1), 1, 1, Nz)

    ω_4Ny⁺ = ω.(4Ny, ky⁺) |> CuArray
    ω_4Nz⁺ = ω.(4Nz, kz⁺) |> CuArray
    ω_4Ny⁻ = ω.(4Ny, ky⁻) |> CuArray
    ω_4Nz⁻ = ω.(4Nz, kz⁻) |> CuArray

    # The zeroth coefficient of the IDCT (DCT-III or FFTW.REDFT01) is not
    # multiplied by 2. For some reason, we only need to account for this when
    # doing a 1D IDCT (for PPN boundary conditions) but not for a 2D IDCT. It's
    # possible that the masks are effectively doing this job.
    if pbcs == PPN()
        ω_4Nz⁻[1] *= 1/2
    end

    # Indices used when we need views to permuted arrays where the odd indices
    # are iterated over first followed by the even indices.
    p_y_inds = [1:2:Ny..., Ny:-2:2...] |> CuArray
    p_z_inds = [1:2:Nz..., Nz:-2:2...] |> CuArray

    # Indices used when we need views with reverse indexing but index N+1 should
    # return a 0. This can't be enforced using views so we just map N+1 to 1,
    # and use masks M_ky and M_kz to enforce that the value at N+1 is 0.
    r_y_inds = [1, collect(Ny:-1:2)...] |> CuArray
    r_z_inds = [1, collect(Nz:-1:2)...] |> CuArray

    M_ky = ones(1, Ny, 1) |> CuArray
    M_kz = ones(1, 1, Nz) |> CuArray

    M_ky[1] = 0
    M_kz[1] = 0

    FFT!, FFT_DCT!, IFFT!, IFFT_DCT! = plan_transforms(pbcs, storage)

    PoissonSolverGPU(pbcs, kx², ky², kz², ω_4Ny⁺, ω_4Ny⁻, ω_4Nz⁺, ω_4Nz⁻,
                     p_y_inds, p_z_inds, r_y_inds, r_z_inds, M_ky, M_kz,
                     storage, storage2, FFT!, FFT_DCT!, IFFT!, IFFT_DCT!)
end

"""
    solve_poisson_3d!(solver::PoissonSolverGPU, grid::RegularCartesianGrid)

Similar to solve_poisson_3d!(solver::PoissonSolverCPU, ...) except that since
the discrete cosine transform is not available through cuFFT, we perform our
own fast cosine transform (FCT) via an algorithm that utilizes the FFT.

Note that for the FCT algorithm to work, the input must have been permuted along
the dimension the FCT is to be calculated by ordering the odd elements first
followed by the even elements. For example,

    [a, b, c, d, e, f, g, h] -> [a, c, e, g, h, f, d, b]

The output will be permuted in this way and so the permutation must be undone.

We should describe the algorithm in detail in the documentation.
"""
function solve_poisson_3d!(solver::PoissonSolverGPU{<:PPN}, grid::RegularCartesianGrid)
    kx², ky², kz² = solver.kx², solver.ky², solver.kz²
    ω_4Ny⁺, ω_4Ny⁻, ω_4Nz⁺, ω_4Nz⁻ = solver.ω_4Ny⁺, solver.ω_4Ny⁻, solver.ω_4Nz⁺, solver.ω_4Nz⁻

    # We can use the same storage for the RHS and the solution ϕ.
    RHS, ϕ = solver.storage, solver.storage

    # Calculate DCTᶻ(f) in place using the FFT.
    solver.FFT_DCT! * RHS
    @. RHS = 2 * real(ω_4Nz⁺ * RHS)

    solver.FFT! * RHS  # Calculate FFTˣʸ(f) in place.

    # Solve the discrete Poisson equation in spectral space. We are essentially
    # computing the Fourier coefficients of the solution from the Fourier
    # coefficients of the RHS.
    @. ϕ = -RHS / (kx² + ky² + kz²)

    # Setting DC component of the solution (the mean) to be zero. This is also
    # necessary because the source term to the Poisson equation has zero mean
    # and so the DC component comes out to be ∞.
    ϕ[1, 1, 1] = 0

    solver.IFFT! * ϕ  # Calculate IFFTˣʸ(ϕ̂) in place.

    # Calculate IDCTᶻ(ϕ̂) in place using the FFT.
    ϕ .*= ω_4Nz⁻
    solver.IFFT_DCT! * ϕ

    nothing
end

function solve_poisson_3d!(solver::PoissonSolverGPU{<:PNN}, grid::RegularCartesianGrid)
    Nx, Ny, Nz, _ = unpack_grid(grid)
    kx², ky², kz² = solver.kx², solver.ky², solver.kz²
    ω_4Ny⁺, ω_4Ny⁻, ω_4Nz⁺, ω_4Nz⁻ = solver.ω_4Ny⁺, solver.ω_4Ny⁻, solver.ω_4Nz⁺, solver.ω_4Nz⁻
    r_y_inds, r_z_inds = solver.r_y_inds, solver.r_z_inds
    M_ky, M_kz = solver.M_ky, solver.M_kz

    # We can use the same storage for the RHS and the solution ϕ.
    RHS, ϕ = solver.storage, solver.storage

    B = solver.storage2  # Complex buffer storage.

    # Calculate DCTʸᶻ(f) in place using the FFT.
    solver.FFT_DCT! * RHS

    RHS⁻ = view(RHS, 1:Nx, r_y_inds, 1:Nz)
    @. B = 2 * real(ω_4Nz⁺ * (ω_4Ny⁺ * RHS + ω_4Ny⁻ * RHS⁻))

    solver.FFT! * B # Calculate FFTˣ(f) in place.

    @. B = -B / (kx² + ky² + kz²)

    B[1, 1, 1] = 0  # Setting DC component of the solution (the mean) to be zero.

    solver.IFFT! * B  # Calculate IFFTˣ(ϕ̂) in place.

    # Calculate IDCTʸᶻ(ϕ̂) in place using the FFT.
    B⁻⁺ = view(B, 1:Nx, r_y_inds, 1:Nz)
    B⁺⁻ = view(B, 1:Nx, 1:Ny, r_z_inds)
    B⁻⁻ = view(B, 1:Nx, r_y_inds, r_z_inds)

    @. ϕ = 1/4 *  ω_4Ny⁻ * ω_4Nz⁻ * ((B - M_ky * M_kz * B⁻⁻) - im*(M_kz * B⁺⁻ + M_ky * B⁻⁺))

    solver.IFFT_DCT! * ϕ

    nothing
end


text:
using Ruska
using Test

@testset "Ruska.jl" begin
    # Write your tests here.
end


text:
################# Parameters to interact with the web app #################

currencies_list = ["BTC", "LTC", "BCH", "ETH", "KNC", "LINK", "ETC", "BNB", "ADA", "XTZ",
              "EOS", "XRP", "XLM", "ZEC", "DASH", "XMR", "DOT", "UNI", "SOL", "MATIC",
              "THETA", "OMG", "ALGO", "GRT", "AAVE", "FIL", "BAT", "ZRX", "COMP"]

currencies = sort(currencies_list)
currencies_index = 1:length(currencies)

modes = ["Average price + Daily trade (AV)", "Candlestick + Volume (AV)",
         "Cumulative + Daily return (AV)", "Daily volatility (AV)", "MACD + Signal (AV)",
         "Linear regression channel (AV)", "FCAS data (AV)", "Developer + Community data (CG)",
          "Exchange volume data per currency (CG)", "Overall exchange volume data (CG)"]
modes_index = 1:length(modes)

durations = [7, 14, 30, 90, 180, 270, 365, 500, 750, 1000]
windows = [1, 5, 10, 30, 50, 75, 100]

################# CoinGecko URL #################

const URL = "https://api.coingecko.com/api/v3"

################# Cleanup function #################

function remove_old_files()
    # Cleanup data files from previous days
    try
        main_dir = pwd()
        cd("data")
        files = readdir()
        rx1 = "data"
        rx2 = "List"
        rx3 = ".csv"
        rx4 = ".txt"
        for file in files
            ts = Dates.unix2datetime(stat(file).mtime)
            file_date = Date(ts)
            if file_date != Dates.today() && (occursin(rx3, file) || occursin(rx4, file)) &&
                                             (occursin(rx1, file) || occursin(rx2, file))
                rm(file)
            end
        end
        cd(main_dir)
    catch
        @info "Unable to perform cleanup action"
    end
end

################# Run the app #################

"""
    run_app(port::Int64, key::String)

Start the app on the given port using the provided API key.

# Arguments
- `port::Int64` : Port number, for example 8056
- `key::String` : API key from Alpha Vantage

# Example
```julia-repl
julia> run_app(8056, "Your API key")
[ Info: data folder exists, cleanup action will be performed!
[ Info: Listening on: 0.0.0.0:8056
[ Info: Fetching BTC price/vol data from Alpha Vantage
```
"""
function run_app(port::Int64, key::String)

    # Check if "data" folder exists, if not, create a new one
    if isdir("data")
        @info "data folder exists, cleanup action will be performed!"
    else
        mkdir("data")
        @info "New data folder has been created"
    end

    # Perform cleanup
    remove_old_files()

    # Set API key
    AlphaVantage.global_key!(key)

    # UI part of the tool
    app = dash()

    app.layout = html_div() do

        # Title
        html_h1("Crypto Dashboard",
                style= (
                textAlign = "center",
                )
        ),

        # Subtitle
        html_div(style=(width="75%", margin="auto", textAlign="center"),

        ["Visualization of market data obtained from",
                html_a(" Alpha Vantage", href="https://www.alphavantage.co/documentation/"),

                html_a(" and"),

                html_a(" CoinGecko", href="https://www.coingecko.com/en/api/documentation"),

                html_a(" until $(Dates.today())"),

                html_p("                                                                 "),
                ]),

        # First row of drop down buttons
        html_div(className="row", [
                html_div(className="col-8",
                    html_table(style=(width="100%", textAlign="center"),
                        vcat(html_tr([html_th("Type of data", style=(width="50%", textAlign="center")),
                                      html_th("Cryptocurrency", style=(width="75%", textAlign="center"))]),
                             [html_tr([
                                     html_td(dcc_dropdown(
                                        id = "mode_ID",
                                        options = [
                                            (label = "$(modes[i])", value = i) for i in modes_index
                                        ],
                                        value = 1,
                                    )),

                                     html_td(dcc_dropdown(
                                        id = "pair_ID",
                                        options = [
                                            (label = "$(currencies[i])", value = i) for i in currencies_index
                                        ],
                                        value = 1,
                                    ))])]
                        )))]),

        # Add space between the two rows
        html_div(style=(width="75%", margin="auto", textAlign="center"),
        [html_p("                                                                 "), ]),

        # Second row of radio items
        html_div(className="row", [
            html_div(className="col-8",
                html_table(style=(width="100%", textAlign="center"),
                    vcat(html_tr([html_th("Moving average window", style=(width="50%", textAlign="center")),
                                    html_th("Historical data [days]", style=(width="75%", textAlign="center"))]),
                            [html_tr([
                                    html_td(dcc_radioitems(
                                        id = "window_ID",
                                        options = [(label = "$(i)-day", value = i) for i in windows],
                                        value = 1,
                                    )),

                                    html_td(dcc_radioitems(
                                        id = "duration_ID",
                                        options = [(label = "$(i)d", value = i) for i in durations],
                                        value = 7,
                                ))])]
                    )))]),

        # Plots
        html_div(style = (width="100%", display="inline-block", margin="0 5% 0 2%")) do
            dcc_graph(id = "graph")
        end
    end

    callback!(
        app,
        Output("graph", "figure"),
        Input("mode_ID", "value"),
        Input("pair_ID", "value"),
        Input("window_ID", "value"),
        Input("duration_ID", "value"),
    ) do mode_ID, pair_ID, window_ID, duration_ID

        if mode_ID == 1
            t1, t2, t3, t4, t5 = plot_price_ma_trade_data(pair_ID, duration_ID, window_ID)

            layout1 = Layout(;title="Daily average price data for $(currencies[pair_ID])",
                xaxis = attr(title="Time", showgrid=true, zeroline=true, linewidth=1.0),
                yaxis = attr(title="Price [euros]", showgrid=true, zeroline=true, linewidth=1.0),
                height = 500,
                width = 1000,
                paper_bgcolor="white"
            )
            layout2 = Layout(;title="Daily trade data (volume x price) for $(currencies[pair_ID])",
                xaxis = attr(title="Time", showgrid=true, zeroline=true, linewidth=1.0),
                yaxis = attr(title="Daily trade [euros]", showgrid=true, zeroline=true, linewidth=1.0),
                height = 500,
                width = 1000,
                paper_bgcolor="white"
            )
            P1 = Plot([t1, t2, t3, t4], layout1) # plots daily average price and three diferent moving averages
            P2 = Plot(t5, layout2)               # plots daily market cap
            return [P1 P2]

        elseif mode_ID == 2
            t1, t2 = plot_candle_vol_data(pair_ID, duration_ID)

            layout1 = Layout(;title="Candlestick data for $(currencies[pair_ID])",
                xaxis=attr(title="Time", showgrid=true, zeroline=true),
                yaxis=attr(title="Price [euros]", zeroline=true),
                height = 500,
                width = 1000,
            )
            layout2 = Layout(;title="Daily volume data for $(currencies[pair_ID])",
                xaxis=attr(title="Time", showgrid=true, zeroline=true),
                yaxis=attr(title="Volume [Number of coins]", zeroline=true),
                height = 100,
                width = 200
            )
            P1 = Plot(t1, layout1) # plots candlestick data
            P2 = Plot(t2, layout2) # plots daily volume
            return [P1 P2]

        elseif mode_ID == 3
            t1, t2_green, t2_red, _ , _ , _ = plot_cumul_daily_return_hist(pair_ID, duration_ID)

            layout1 = Layout(;title="Cumulative return for $(currencies[pair_ID])",
                xaxis=attr(title="Time", showgrid=true, zeroline=true),
                yaxis=attr(title="Return [in %]", zeroline=true),
                height = 500,
                width = 1000,
            )
            layout2 = Layout(;title="Daily return for $(currencies[pair_ID])",
                xaxis=attr(title="Time", showgrid=true, zeroline=true),
                yaxis=attr(title="Return [in %]", zeroline=true),
                height = 500,
                width = 1000,
                barmode = "group",
            )

            P1 = Plot(t1, layout1)                 # plots cumulative return %
            P2 = Plot([t2_green, t2_red], layout2) # plots daily change %
            return [P1 P2]

        elseif mode_ID == 4
            _ , _ , _ , t1 , σ, duration = plot_cumul_daily_return_hist(pair_ID, duration_ID)

            layout1 = Layout(;title="Distribution of daily price change for $(currencies[pair_ID]) over $(duration) days, 3σ = $(round(3*σ, digits = 2)) %",
                xaxis=attr(title="Change [in %]", showgrid=true, zeroline=true),
                yaxis=attr(title="Number of counts", zeroline=true),
                height = 500,
                width = 1000,
            )

            P1 = Plot(t1, layout1) # plots histogram of daily price change
            return [P1]

        elseif mode_ID == 5
            t1, t2, t3, t4, t5, t6_green, t6_red = plot_macd_signal(pair_ID, duration_ID)

            layout1 = Layout(;title="EMA and average price data for $(currencies[pair_ID])",
                xaxis = attr(title="Time", showgrid=true, zeroline=true, linewidth=1.0),
                yaxis = attr(title="Price [euros]", showgrid=true, zeroline=true, linewidth=1.0),
                height = 500,
                width = 1000,
            )
            layout2 = Layout(;title="MACD and signal for $(currencies[pair_ID])",
                xaxis = attr(title="Time", showgrid=true, zeroline=true, linewidth=1.0),
                yaxis = attr(title="Indicator [euros]", showgrid=true, zeroline=true, linewidth=1.0),
                height = 500,
                width = 1000,
            )

            P1 = Plot([t1, t2, t3], layout1)               # plots average price, EMA-12 and EMA-26
            P2 = Plot([t4, t5, t6_green, t6_red], layout2) # plots MACD, signal and their distance
            return [P1 P2]

        elseif mode_ID == 6
            t1, t2, t3, t4, R² = plot_linear_regression(pair_ID, duration_ID)

            layout1 = Layout(;title="Linear regression channel for $(currencies[pair_ID]) price over last $(duration_ID) days, R² = $(R²)",
                xaxis = attr(title="Time", showgrid=true, zeroline=true),
                yaxis = attr(title="Price [euros]", showgrid=true, zeroline=true),
                height = 500,
                width = 1000,
            )

            P1 = Plot([t1, t2, t3, t4], layout1) # plots linear regression channel
            return [P1]

        elseif mode_ID == 7
            t1, fr = plot_fcas_data(pair_ID)

            layout1 = Layout(;title="FCAS metrics data for $(currencies[pair_ID]), overall rating = $(fr)",
                xaxis = attr(title="Type of metric", showgrid=true, zeroline=true),
                yaxis = attr(title="Score", showgrid=true, zeroline=true),
                height = 500,
                width = 1000,
            )

            P1 = Plot(t1, layout1) # plots FCAS metrics
            return [P1]

        # From CoinGecko
        elseif mode_ID == 8
            t1, t2 = plot_dev_comm_data(pair_ID)

            layout1 = Layout(;title="Developer metrics for $(currencies[pair_ID])",
                xaxis = attr(title="", showgrid=true, zeroline=true, automargin=true),
                xaxis_tickangle = -22.5,
                yaxis = attr(title="Value", showgrid=true, zeroline=true),
                height = 500,
                width = 1000,
            )
            layout2 = Layout(;title="Community metrics for $(currencies[pair_ID])",
                xaxis = attr(title="", showgrid=true, zeroline=true, automargin=true),
                xaxis_tickangle = -22.5,
                yaxis = attr(title="Value", showgrid=true, zeroline=true),
                height = 500,
                width = 1000,
            )

            P1 = Plot(t1, layout1) # plots developer data
            P2 = Plot(t2, layout2) # plots community data
            return [P1 P2]

        elseif mode_ID == 9

            t1, t2 = plot_exchange_vol_data(pair_ID)

            layout1 = Layout(;title="Exchange volume data (24h) for $(currencies[pair_ID])",
                xaxis = attr(title="", showgrid=true, zeroline=true, automargin=true),
                xaxis_tickangle = -22.5,
                yaxis = attr(title="Number of coins", showgrid=true, zeroline=true),
                height = 500,
                width = 1000,
            )
            layout2 = Layout(;title="Exchange volume data (24h) for $(currencies[pair_ID])",
                xaxis = attr(title="", showgrid=true, zeroline=true, automargin=true),
                xaxis_tickangle = -22.5,
                yaxis = attr(title="Volume in USD", showgrid=true, zeroline=true),
                height = 500,
                width = 1000,
            )

            P1 = Plot(t1, layout1)      # plots coin volume data for a given currency
            P2 = Plot(t2, layout2)      # plots USD volume data for a given currency
            return [P1 P2]

        elseif mode_ID == 10

            t_all = plot_overall_vol_data(duration_ID)

            layout1 = Layout(;title="Overall historical volume data (24h) for top 10 exchanges",
                xaxis = attr(title="", showgrid=true, zeroline=true, automargin=true),
                xaxis_tickangle = 0.0,
                yaxis = attr(title="Volume in BTC", showgrid=true, zeroline=true),
                height = 500,
                width = 1000,
                barmode = "stack",
            )

            P1 = Plot(t_all, layout1)    # plots overall BTC volume data for all exchanges
            return P1

        end
    end

    # Allows access from a web browser, port can be changed to any valid entry
    run_server(app, "0.0.0.0", port)
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#     Functions to manage the I2C.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export init_i2c, close_i2c, i2c_get_funcs, i2c_slave, i2c_slave_force,
       i2c_smbus_read_byte, i2c_smbus_read_byte_data, i2c_smbus_read_word_data,
       i2c_smbus_read_i2c_block_data!, i2c_smbus_read_i2c_block_data,
       i2c_smbus_write_byte, i2c_smbus_write_byte_data,
       i2c_smbus_write_word_data, i2c_smbus_write_i2c_block_data

################################################################################
#                                Initialization
################################################################################

"""
    init_i2c(devices)

Initialize the I2C devices. `devices` can be a string with the path to `i2cdev`
or a vector of strings with a set of I2C devices that will be initialized.

"""
@inline init_i2c(device::String) = init_i2c([device])

function init_i2c(devices::AbstractVector{String})

    # Number of devices that the user wants to initialize.
    num_init_devices = length(devices)

    # If we have already initialized, then close all the devices first.
    objects.i2c_init && close_i2c()

    # Initialize each device.
    i2cdev = Vector{I2CDEV}(undef, 0)

    @inbounds for i = 1:num_init_devices
        try
            # Open the device
            # ==================================================================

            io = open(devices[i], "w+")

            # Get the available functionality
            # ==================================================================

            funcs = _i2c_get_funcs(io)

            # Pack values
            # ==================================================================
            i2cdev_i = I2CDEV(io, funcs)
            push!(i2cdev, i2cdev_i)
        catch
            error("Could not open `$(devices[i])`. Make sure you have the required permissions.")
        end
    end

    objects.i2cdev   = i2cdev
    objects.i2c_init = true

    return nothing
end

"""
    close_i2c()

Close all I2C connections.

"""
function close_i2c()
    if objects.i2c_init
        @inbounds for d in objects.i2cdev
            close(d.io)
        end
        objects.i2cdev = Vector{SPIDEV}(undef, 0)
        objects.i2c_init = false
    end

    return nothing
end

################################################################################
#                                   Commands
################################################################################

"""
    i2c_get_funcs(devid::Int)

Return a vector with the available I2C functions in the device `devid`.

`devid` is the ID of the I2C device considering the initialization order when
the function `init_i2c` was called.

"""
@inline function i2c_get_funcs(devid::Int)
    @assert objects.i2c_init "I2C not initialized. Run init_i2c()."
    @assert (0 < devid ≤ length(objects.i2cdev)) "I2C device ID is out of bounds."

    return objects.i2cdev[devid].funcs
end

################################################################################
#                                    SMBUS
################################################################################

"""
    i2c_slave(devid::Integer, address::Integer)

Select the slave device in address `address` using I2C device `devid`.

`devid` is the ID of the I2C device considering the initialization order when
the function `init_i2c` was called.

"""
function i2c_slave(devid::Integer, address::Integer)
    @assert objects.i2c_init "I2C not initialized. Run init_i2c()."
    @assert (0 < devid ≤ length(objects.i2cdev)) "I2C device ID is out of bounds."

    i2cdev = objects.i2cdev[devid]

    _ioctl(fd(i2cdev.io), I2C_SLAVE, __u8(address))

    return nothing
end

"""
    i2c_slave_force(devid::Integer, address::Integer)

Force selection of slave device in address `address` using I2C device `devid`.

`devid` is the ID of the I2C device considering the initialization order when
the function `init_i2c` was called.

"""
function i2c_slave_force(devid::Integer, address::Integer)
    @assert objects.i2c_init "I2C not initialized. Run init_i2c()."
    @assert (0 < devid ≤ length(objects.i2cdev)) "I2C device ID is out of bounds."

    i2cdev = objects.i2cdev[devid]

    _ioctl(fd(i2cdev.io), I2C_SLAVE_FORCE, __u8(address))

    return nothing
end

# Read
# ==============================================================================

"""
    i2c_smbus_read_byte(devid::Integer)

Perform a SMBUS read byte using the I2C device `devid`. This functions returns
the read byte.

`devid` is the ID of the I2C device considering the initialization order when
the function `init_i2c` was called.

"""
function i2c_smbus_read_byte(devid::Integer)
    @assert objects.i2c_init "I2C not initialized. Run init_i2c()."
    @assert (0 < devid ≤ length(objects.i2cdev)) "I2C device ID is out of bounds."

    i2cdev = objects.i2cdev[devid]

    # Create the SMBUS IOCTL structure to issue the ioctl command.
    smbus = struct_i2c_smbus_ioctl_data(I2C_SMBUS_READ,
                                        0,
                                        I2C_SMBUS_BYTE,
                                        pointer(i2cdev.bbyte))

    _ioctl(fd(i2cdev.io), I2C_SMBUS, Ref(smbus))

    return (@inbounds i2cdev.bbyte[1])
end

"""
    i2c_smbus_read_byte_data(devid::Integer, command::Integer)

Perform a SMBUS read byte with command `command` using the I2C device `devid`.
This function return the read byte.

`devid` is the ID of the I2C device considering the initialization order when
the function `init_i2c` was called.

"""
function i2c_smbus_read_byte_data(devid::Integer, command::Integer)
    @assert objects.i2c_init "I2C not initialized. Run init_i2c()."
    @assert (0 < devid ≤ length(objects.i2cdev)) "I2C device ID is out of bounds."

    i2cdev = objects.i2cdev[devid]

    # Create the SMBUS IOCTL structure to issue the ioctl command.
    smbus = struct_i2c_smbus_ioctl_data(I2C_SMBUS_READ,
                                        command,
                                        I2C_SMBUS_BYTE_DATA,
                                        pointer(i2cdev.bbyte))

    _ioctl(fd(i2cdev.io), I2C_SMBUS, Ref(smbus))

    return (@inbounds i2cdev.bbyte[1])
end

"""
    i2c_smbus_read_word_data(devid::Integer, command::Integer)

Perform a SMBUS read word with command `command` using the I2C device `devid`.
This function return the read word.

`devid` is the ID of the I2C device considering the initialization order when
the function `init_i2c` was called.

"""
function i2c_smbus_read_word_data(devid::Integer, command::Integer)
    @assert objects.i2c_init "I2C not initialized. Run init_i2c()."
    @assert (0 < devid ≤ length(objects.i2cdev)) "I2C device ID is out of bounds."

    i2cdev = objects.i2cdev[devid]

    # Create the SMBUS IOCTL structure to issue the ioctl command.
    smbus = struct_i2c_smbus_ioctl_data(I2C_SMBUS_READ,
                                        command,
                                        I2C_SMBUS_WORD_DATA,
                                        pointer(i2cdev.bword))

    _ioctl(fd(i2cdev.io), I2C_SMBUS, Ref(smbus))

    return (@inbounds i2cdev.bword[1])
end

"""
    i2c_smbus_read_i2c_block_data(devid::Integer, command::Integer, size::Integer)

Perform a SMBUS read block with command `command` and length `length` using the
I2C device `devid`. The read data will be returned in an array of `UInt8` with
size `size + 1`. This array is a view of the buffer that was allocated at the
initialization of the I2C device. Notice that the first element of the array is
always `size`.

`devid` is the ID of the I2C device considering the initialization order when
the function `init_i2c` was called.

!!! warning

    The returned array is a view of the buffer `bblock` created for each I2C
    device. This buffer is modified by other functions like
    `i2c_smbus_write_i2c_block_data`. Hence, before performing another I2C
    transfer that will modify this buffer, be sure to copy the returned value to
    another place.

!!! note "Allocations"

    Due to the view returned by this functions, it performs one allocation. If
    this is not wanted, then use the in-place version
    `i2c_smbus_read_i2c_block_data!`.

"""
function i2c_smbus_read_i2c_block_data(devid::Integer, command::Integer,
                                       size::Integer)

    @assert objects.i2c_init "I2C not initialized. Run init_i2c()."
    @assert (0 < devid ≤ length(objects.i2cdev)) "I2C device ID is out of bounds."
    @assert (size > 0) "The `size` must be larger than 0."
    @assert (size ≤ I2C_SMBUS_BLOCK_MAX) "I2C block data cannot read more than $(I2C_SMBUS_BLOCK_MAX) bytes."

    i2cdev = objects.i2cdev[devid]
    data = i2cdev.bblock

    # The first parameters is the length of the block we want to read.
    @inbounds data[1] = size

    # Create the SMBUS IOCTL structure to issue the ioctl command.
    smbus = struct_i2c_smbus_ioctl_data(I2C_SMBUS_READ,
                                        command,
                                        I2C_SMBUS_I2C_BLOCK_DATA,
                                        pointer(data))

    _ioctl(fd(i2cdev.io), I2C_SMBUS, Ref(smbus))

    return view(data, 1:size+1)
end

"""
    i2c_smbus_read_i2c_block_data!(devid::Integer, command::Integer, data::Vector{UInt8})

Perform a SMBUS read block with command `command` using the I2C device `devid`.
The read data will be written to `data`. Notice that the first element of the
array is always `size`. Hence, the number of read bytes will be equal to the
length of `data` minus 1.

`devid` is the ID of the I2C device considering the initialization order when
the function `init_i2c` was called.

"""
function i2c_smbus_read_i2c_block_data!(devid::Integer, command::Integer,
                                        data::Vector{UInt8})

    @assert objects.i2c_init "I2C not initialized. Run init_i2c()."
    @assert (0 < devid ≤ length(objects.i2cdev)) "I2C device ID is out of bounds."

    size = length(data)

    @assert (size ≤ I2C_SMBUS_BLOCK_MAX+1) "I2C block data cannot read more than $(I2C_SMBUS_BLOCK_MAX) bytes."
    @assert (size > 1) "The `data` array must have at least 2 elements."

    i2cdev = objects.i2cdev[devid]

    # The first parameters is the length of the block we want to read.
    @inbounds data[1] = size - 1

    # Create the SMBUS IOCTL structure to issue the ioctl command.
    smbus = struct_i2c_smbus_ioctl_data(I2C_SMBUS_READ,
                                        command,
                                        I2C_SMBUS_I2C_BLOCK_DATA,
                                        pointer(data))

    _ioctl(fd(i2cdev.io), I2C_SMBUS, Ref(smbus))

    return nothing
end

# Write
# ==============================================================================

"""
    i2c_smbus_write_byte(devid::Integer, value::Integer)

Perform a SMBUS write byte with value `value` using the I2C device `devid`.

`devid` is the ID of the I2C device considering the initialization order when
the function `init_i2c` was called.

"""
function i2c_smbus_write_byte(devid::Integer, value::Integer)
    @assert objects.i2c_init "I2C not initialized. Run init_i2c()."
    @assert (0 < devid ≤ length(objects.i2cdev)) "I2C device ID is out of bounds."

    i2cdev = objects.i2cdev[devid]

    # Create the SMBUS IOCTL structure to issue the ioctl command.
    smbus = struct_i2c_smbus_ioctl_data(I2C_SMBUS_WRITE,
                                        value,
                                        I2C_SMBUS_BYTE,
                                        C_NULL)

    _ioctl(fd(i2cdev.io), I2C_SMBUS, Ref(smbus))

    return nothing
end

"""
    i2c_smbus_write_byte_data(devid::Integer, command::Integer, value::Integer)

Perform a SMBUS write byte with command `command` and value `value` using the
I2C device `devid`.

`devid` is the ID of the I2C device considering the initialization order when
the function `init_i2c` was called.

"""
function i2c_smbus_write_byte_data(devid::Integer, command::Integer,
                                   value::Integer)
    @assert objects.i2c_init "I2C not initialized. Run init_i2c()."
    @assert (0 < devid ≤ length(objects.i2cdev)) "I2C device ID is out of bounds."

    i2cdev = objects.i2cdev[devid]

    # Create the SMBUS IOCTL structure to issue the ioctl command.
    data = i2cdev.bbyte
    @inbounds data[1] = value

    smbus = struct_i2c_smbus_ioctl_data(I2C_SMBUS_WRITE,
                                        command,
                                        I2C_SMBUS_BYTE_DATA,
                                        pointer(data))

    _ioctl(fd(i2cdev.io), I2C_SMBUS, Ref(smbus))

    return nothing
end

"""
    i2c_smbus_write_word_data(devid::Integer, command::Integer, value::Integer)

Perform a SMBUS write word with command `command` and value `value` using the
I2C device `devid`.

`devid` is the ID of the I2C device considering the initialization order when
the function `init_i2c` was called.

"""
function i2c_smbus_write_word_data(devid::Integer, command::Integer,
                                   value::Integer)
    @assert objects.i2c_init "I2C not initialized. Run init_i2c()."
    @assert (0 < devid ≤ length(objects.i2cdev)) "I2C device ID is out of bounds."

    i2cdev = objects.i2cdev[devid]

    # Create the SMBUS IOCTL structure to issue the ioctl command.
    data = i2cdev.bword
    @inbounds data[1] = value

    smbus = struct_i2c_smbus_ioctl_data(I2C_SMBUS_WRITE,
                                        command,
                                        I2C_SMBUS_WORD_DATA,
                                        pointer(data))

    _ioctl(fd(i2cdev.io), I2C_SMBUS, Ref(smbus))

    return nothing
end

"""
    i2c_smbus_write_i2c_block_data(devid::Integer, command::Integer, values::Vector{UInt8})

Perform a SMBUS write I2C block data with command `command` and data in `values`
using the I2C device `devid`.

`devid` is the ID of the I2C device considering the initialization order when
the function `init_i2c` was called.

"""
function i2c_smbus_write_i2c_block_data(devid::Integer, command::Integer,
                                        values::Vector{UInt8})

    @assert objects.i2c_init "I2C not initialized. Run init_i2c()."
    @assert (0 < devid ≤ length(objects.i2cdev)) "I2C device ID is out of bounds."

    size = length(values)

    @assert (size ≤ I2C_SMBUS_BLOCK_MAX) "I2C block data cannot write more than $(I2C_SMBUS_BLOCK_MAX) bytes."

    i2cdev = objects.i2cdev[devid]

    # Create copy the data to the buffer.
    data = i2cdev.bblock
    @inbounds data[1] = size

    # If we use a slice here, then one allocation will be performed because of
    # the required `view` that will be created.
    @inbounds for i = 1:size
        data[i+1] = values[i]
    end

    # Create the SMBUS IOCTL structure to issue the ioctl command.
    smbus = struct_i2c_smbus_ioctl_data(I2C_SMBUS_WRITE,
                                        command,
                                        I2C_SMBUS_I2C_BLOCK_DATA,
                                        pointer(data))

    _ioctl(fd(i2cdev.io), I2C_SMBUS, Ref(smbus))

    return nothing
end

################################################################################
#                                   Private
################################################################################

function _i2c_get_funcs(i2cdev_io)
    # Create the ioctl that will return the functionality mask.
    pfuncs = Ref(Culong(0))
    _ioctl(fd(i2cdev_io), I2C_FUNCS, pfuncs)
    funcs = pfuncs.x

    # Check the existing functionality.
    translated_funcs = Vector{Symbol}(undef, 0)
    (I2C_FUNC_I2C & funcs) > 0                    && push!(translated_funcs, :i2c)
    (I2C_FUNC_10BIT_ADDR & funcs) > 0             && push!(translated_funcs, :tenbit_addr)
    (I2C_FUNC_PROTOCOL_MANGLING & funcs) > 0      && push!(translated_funcs, :protocol_mangling)
    (I2C_FUNC_SMBUS_PEC & funcs) > 0              && push!(translated_funcs, :smbus_pec)
    (I2C_FUNC_NOSTART & funcs) > 0                && push!(translated_funcs, :nostart)
    (I2C_FUNC_SLAVE & funcs) > 0                  && push!(translated_funcs, :slave)
    (I2C_FUNC_SMBUS_BLOCK_PROC_CALL & funcs) > 0  && push!(translated_funcs, :smbus_block_proc_call)
    (I2C_FUNC_SMBUS_QUICK & funcs) > 0            && push!(translated_funcs, :smbus_quick)
    (I2C_FUNC_SMBUS_READ_BYTE & funcs) > 0        && push!(translated_funcs, :smbus_read_byte)
    (I2C_FUNC_SMBUS_WRITE_BYTE & funcs) > 0       && push!(translated_funcs, :smbus_write_byte)
    (I2C_FUNC_SMBUS_READ_BYTE_DATA & funcs) > 0   && push!(translated_funcs, :smbus_read_byte_data)
    (I2C_FUNC_SMBUS_WRITE_BYTE_DATA & funcs) > 0  && push!(translated_funcs, :smbus_write_byte_data)
    (I2C_FUNC_SMBUS_READ_WORD_DATA & funcs) > 0   && push!(translated_funcs, :smbus_read_word_data)
    (I2C_FUNC_SMBUS_WRITE_WORD_DATA & funcs) > 0  && push!(translated_funcs, :smbus_write_word_data)
    (I2C_FUNC_SMBUS_PROC_CALL & funcs) > 0        && push!(translated_funcs, :smbus_proc_call)
    (I2C_FUNC_SMBUS_READ_BLOCK_DATA & funcs) > 0  && push!(translated_funcs, :smbus_read_block_data)
    (I2C_FUNC_SMBUS_WRITE_BLOCK_DATA & funcs) > 0 && push!(translated_funcs, :smbus_write_block_data)
    (I2C_FUNC_SMBUS_READ_I2C_BLOCK & funcs) > 0   && push!(translated_funcs, :smbus_read_i2c_block)
    (I2C_FUNC_SMBUS_WRITE_I2C_BLOCK & funcs) > 0  && push!(translated_funcs, :smbus_write_i2c_block)
    (I2C_FUNC_SMBUS_HOST_NOTIFY & funcs) > 0      && push!(translated_funcs, :smbus_host_notify)
    (I2C_FUNC_SMBUS_BYTE & funcs) > 0             && push!(translated_funcs, :smbus_byte)
    (I2C_FUNC_SMBUS_BYTE_DATA & funcs) > 0        && push!(translated_funcs, :smbus_byte_data)
    (I2C_FUNC_SMBUS_WORD_DATA & funcs) > 0        && push!(translated_funcs, :smbus_word_data)
    (I2C_FUNC_SMBUS_BLOCK_DATA & funcs) > 0       && push!(translated_funcs, :smbus_block_data)
    (I2C_FUNC_SMBUS_I2C_BLOCK & funcs) > 0        && push!(translated_funcs, :smbus_i2c_block)
    (I2C_FUNC_SMBUS_EMUL & funcs) > 0             && push!(translated_funcs, :smbus_emul)

    return translated_funcs
end


text:
@test (MIDI.Note(60, 96, 0, 0) + 2).value == 62
@test (MIDI.Note(60, 96, 0, 0) - 2).value == 58
@test (MIDI.Note(60, 96, 0, 0) + 0).value == 60


text:
using HPL
using Test

@testset "HPL.jl" begin
    # Write your tests here.
end


text:
function test_shockcapturing_parallel()

  opts = Dict{String, Any}(
    "physics" => "Euler",
    "operator_type" => "SBPOmega",
    "dimensions" => 2,
    "run_type" => 5,
    "jac_method" => 2,
    "jac_type" => 3,
    "order" => 2,
    "IC_name" => "ICIsentropicVortex",
    "use_DG" => true,
    "volume_integral_type" => 2,
    "Volume_flux_name" => "IRFlux",
    "face_integral_type" => 2,
    "FaceElementIntegral_name" => "ESLFFaceIntegral",
    "Flux_name" => "IRFlux",
    "numBC" => 3,
    "BC1" => [0],
    "BC1_name" => "isentropicVortexBC",  # outlet
    "BC2" => [2],
    "BC2_name" => "isentropicVortexBC", # inlet
    "BC3" => [1, 3],
    "BC3_name" => "isentropicVortexBC",
    "aoa" => 0.0,
    "smb_name" => "SRCMESHES/vortex_3x3_p_.smb",
    "dmg_name" => ".null",
    "itermax" => 20,
    "res_abstol" => 1e-12,
    "res_reltol" => 1e-12,
    "do_postproc" => true,
    "exact_soln_func" => "ICIsentropicVortex",
    "solve" => false,
    )

  opts2 = copy(opts)
  opts2["operator_type"] = "SBPDiagonalE"
  opts2["face_integral_type"] = 1
  opts2["flux_name"] = "IRSLFFlux"
  opts2["use_lps"] = true

  mesh, sbp, eqn, opts = solvePDE(opts)
  mesh2, sbp2, eqn2, opts2 = solvePDE(opts2)

  @testset "Shock Capturing" begin
    test_shockmesh(mesh, sbp, eqn, opts)
    test_br2_parallelpart(mesh, sbp, eqn, opts)
    test_br2reduced_parallelpart(mesh2, sbp2, eqn2, opts2)
  end

end

function getEntireMesh(mesh, sbp, eqn::EulerData{Tsol, Tres}, opts) where {Tsol, Tres}

  degree = sbp.degree

  # construct the shock mesh with all elements in it that are fully interior
  shockmesh = EulerEquationMod.ShockedElements{Tres}(mesh)

  for i=1:mesh.numEl
    push!(shockmesh, i)
  end

  EulerEquationMod.completeShockElements(mesh, shockmesh)

  return shockmesh
end

function test_shockmesh(mesh, sbp, eqn, opts)

  shockmesh = getEntireMesh(mesh, sbp, eqn, opts)

  @test shockmesh.numShock == mesh.numEl
  @test shockmesh.numNeighbor == 0
  @test shockmesh.npeers == mesh.npeers

  for i=1:shockmesh.npeers
    @test shockmesh.peer_indices[i] == i
    @test shockmesh.numShared[i] == (mesh.shared_element_offsets[i+1] - mesh.shared_element_offsets[i])
    @test length(shockmesh.shared_els[i]) == shockmesh.numShared[i]
    @test shockmesh.numSharedInterfaces[i] == length(mesh.shared_interfaces[i])
  end

  return nothing
end

function test_br2_parallelpart(mesh, sbp, eqn::EulerData{Tsol, Tres},
                               _opts) where {Tsol, Tres}
# test consistency with the serial version

  opts = copy(_opts)
  sensor = EulerEquationMod.ShockSensorEverywhere{Tsol, Tres}(mesh, sbp, opts)
  capture = EulerEquationMod.SBPParabolicSC{Tsol, Tres}(mesh, sbp, eqn, opts, sensor)

  # solve the PDE to get a solution with non-zero jump between elements
  # that can be reproduced in parallel
  opts["solve"] = true
  solvePDE(mesh, sbp, eqn, opts)

  # the scheme is entropy stable when the dirichlet BC = 0
  for i=1:opts["numBC"]
    opts[string("BC", i, "_name")] = "zeroBC"
  end
  EulerEquationMod.getBCFunctors(mesh, sbp, eqn, opts)

  test_sc_parallelpart(mesh, sbp, eqn, opts, capture, "br2")

  return nothing
end

function test_br2reduced_parallelpart(mesh, sbp, eqn::EulerData{Tsol, Tres},
                                      _opts) where {Tsol, Tres}
# test consistency with the serial version

  opts = copy(_opts)
  sensor = EulerEquationMod.ShockSensorEverywhere{Tsol, Tres}(mesh, sbp, opts)
  capture = EulerEquationMod.SBPParabolicReducedSC{Tsol, Tres}(mesh, sbp, eqn, opts, sensor)

  # solve the PDE to get a solution with non-zero jump between elements
  # that can be reproduced in parallel
  opts["solve"] = true
  solvePDE(mesh, sbp, eqn, opts)

  test_sc_parallelpart(mesh, sbp, eqn, opts, capture, "br2reduced")

  return nothing
end

function test_sc_parallelpart(mesh, sbp, eqn::EulerData{Tsol, Tres},
                              opts, capture, prefix::String) where {Tsol, Tres}
  # need to communicate the final solution
  startSolutionExchange(mesh, sbp, eqn, opts, wait=true)

  w_vec = zeros(Tsol, mesh.numDof)
  copy!(w_vec, eqn.q_vec)
  EulerEquationMod.convertToIR(mesh, sbp, eqn, opts, w_vec)

  fill!(eqn.res, 0)
  EulerEquationMod.applyShockCapturing(mesh, sbp, eqn, opts, capture)

  val = dot(w_vec, eqn.res_vec)
  val = MPI.Allreduce(val, MPI.SUM, eqn.comm)

  println("val = ", val)
  @test val < 0

  # save this for parallel tests
  val2 = readdlm("$(prefix)_entropy_serial.dat")
  println("val2 = ", val2[1])
  println("diff = ", val - val2[1])

  @test abs(val - val2[1]) < 1e-8  # test the entropy dissipation is the same,
                                   # to the tolerance of the nonlinear solve

  saveSolutionToMesh(mesh, eqn.q_vec)
  writeVisFiles(mesh, "$(prefix)_parallel")

  return nothing
end

add_func1!(EulerTests, test_shockcapturing_parallel, [TAG_SHORTTEST, TAG_TMP])


text:
using Test, Random, AdvancedHMC
using Statistics: mean
include("common.jl")

@testset "`gdemo`" begin
    rng = MersenneTwister(0)

    n_samples = 5_000
    n_adapts = 1_000

    θ_init = randn(rng, 2)

    metric = DiagEuclideanMetric(2)
    h = Hamiltonian(metric, ℓπ_gdemo, ∂ℓπ∂θ_gdemo)
    init_eps = Leapfrog(0.1)
    prop = NUTS(init_eps)
    adaptor = StanHMCAdaptor(n_adapts, Preconditioner(metric), NesterovDualAveraging(0.8, prop.integrator))

    samples, _ = sample(rng, h, prop, θ_init, n_samples, adaptor, n_adapts)

    m_est = mean(map(_s -> [exp(_s[1]), _s[2]], samples[1000:end]))

    @test m_est ≈ [49 / 24, 7 / 6] atol=RNDATOL
end


text:
@doc """
weval(coef, wavename, res[, supp])

Evaluate the coefficients `coef` of the `wavename` scaling function
basis `supp` (depending on whether `coef` is a vector or a matrix).

The functions are evaluated at the dyadic rationals of resolution `res`.
"""->
function weval(coef, wavename::String, res::Integer, supp::Interval=Interval(-0.5, 0.5))

	supp_length = length(supp)
	isinteger(supp_length) || throw(DomainError())
	S = DaubSupport(0, supp_length)

	vm = van_moment(wavename)
	if vm == 1
		y = weval(coef, res, S)
	elseif vm >= 2
		y = weval(coef, vm, res, S)
	end

    if length(y) == 2
        return y[1] + left(supp), y[2]
    else
        return y
    end
end

@doc """
	weval(coef::Vector, res::Int, supp) -> x, y

Evaluate the coefficients `coef` of the Haar scaling function basis on the DaubSupport interval `supp`.
The functions are evaluated at the dyadic rationals of resolution `res`.
"""->
function weval(coef::AbstractVector, res::Integer, supp::DaubSupport)
	Ncoef = length(coef)
	lsupp = length(supp)
	Jpow = div(Ncoef, lsupp)
	Ncoef == lsupp * Jpow || throw(AssertionError())
	ispow2(Jpow) || throw(AssertionError())
	( J = Int(log2(Jpow)) ) <= res || throw(AssertionError())

	# Points in support
	x = dyadic_rationals(supp, res)

	dilation = 2.0^(-J)
	sqrt_dil = sqrt(dilation)
	y = zeros(x)

	translations = left(supp)*Jpow + (0:Ncoef-1)
	coef_idx = 0
	for k in translations
		# The indices of the x's in the support of phi_{J,k}
		start_idx = x2index( dilation*k, supp, res )
		end_idx = x2index( dilation*(k+1), supp, res )

		coef_idx+=1
		@inbounds for nx in start_idx:end_idx
			phi = HaarScaling(x[nx]/dilation - k) / sqrt_dil
			y[nx] += coef[coef_idx]*phi
		end
	end

	return x, y
end

@doc """
	weval(coef::Matrix, res::Int) -> Matrix

Evaluate the coefficients `coef` of the Haar scaling function basis on the DaubSupport interval `supp`.
The functions are evaluated at the dyadic rationals of resolution `res`.
"""->
function weval(coef::AbstractMatrix, res::Integer, supp::DaubSupport)
	(Ncoef = size(coef,1)) == size(coef,2) || throw(DimensionMismatch())
	lsupp = length(supp)
	Jpow = div(Ncoef, lsupp)
	Ncoef == lsupp * Jpow || throw(AssertionError())
	ispow2(Jpow) || throw(AssertionError())
	(J = Int(log2(Jpow)) ) <= res || throw(DomainError())

	x = dyadic_rationals(supp, res)

	Nx = length(x)
	y = zeros(Float64, Nx, Nx)
	dilation = 2.0^(-J)
	sqrt_dil = sqrt(dilation)

	translations = left(supp)*Jpow + (0:Ncoef-1)
	coefy_idx = 0
	for ky in translations
		starty_idx = x2index( dilation*ky, supp, res )
		endy_idx = x2index( dilation*(ky+1), supp, res )

		coefy_idx += 1
		coefx_idx = 0
		for kx in translations
			startx_idx = x2index( dilation*kx, supp, res )
			endx_idx = x2index( dilation*(kx+1), supp, res )

			coefx_idx += 1
			for nx in startx_idx:endx_idx
				@inbounds phix = HaarScaling(x[nx]/dilation - kx) / sqrt_dil

				@inbounds for ny in starty_idx:endy_idx
					phiy = HaarScaling(x[ny]/dilation - ky) / sqrt_dil
					y[ny,nx] += coef[coefy_idx,coefx_idx] * phix * phiy
				end
			end
		end
	end

	return y
end

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

function weval(coef::AbstractVector, p::Integer, R::Integer, supp::DaubSupport)
	Ncoef = length(coef)
	lsupp = length(supp)
	Jpow = div(Ncoef, lsupp)
	Ncoef == lsupp * Jpow || throw(AssertionError())
	ispow2(Jpow) || throw(AssertionError())
	2 <= p <= 8 || throw(AssertionError())
	log2( (2p-1)/lsupp ) <= ( J = Int(log2(Jpow)) ) <= R || throw(AssertionError())

	# Collect mother scaling functions
	PHI = allDaubScaling(p, R-J)
	scale!(PHI, sqrt(Jpow))

	# Allocate output and temporary array
	x = dyadic_rationals(supp, R)
	y = zeros(x)
	Nphi = size(PHI,1)
	phi = Array{Float64}(Nphi)

	# Translation with 1 is a fixed number of indices
	kinc = x2index( left(supp)+1/Jpow, supp, R ) - 1
	slice_idx = 1:Nphi

	sz = size(PHI)
	for k in 1:Ncoef
		unsafe_DaubScaling!(phi, k, PHI, Ncoef, sz, p)

		if p < k <= Ncoef - p + 1
			slice_idx += kinc
		end
		z = view(y, slice_idx)
		BLAS.axpy!( coef[k], phi, z )
	end

	return x, y
end

function weval(coef::AbstractMatrix, p::Integer, R::Integer, supp::DaubSupport)
	(Ncoef = size(coef,1)) == size(coef,2) || throw(DimensionMismatch())
	lsupp = length(supp)
	Jpow = div(Ncoef, lsupp)
	Ncoef == lsupp * Jpow || throw(AssertionError())
	ispow2(Jpow) || throw(AssertionError())
	2 <= p <= 8 || throw(AssertionError())
	log2( (2p-1)/lsupp ) <= ( J = Int(log2(Jpow)) ) <= R || throw(AssertionError())

	# Collect mother scaling functions
	PHI = allDaubScaling(p, R-J)
	scale!(PHI, sqrt(Jpow))

	# Allocate output and temporary arrays
	x = dyadic_rationals(supp, R)
	Nx = length(x)
	y = zeros(Float64, Nx, Nx)
	Nphi = size(PHI,1)
	phi = Array{Float64}(Nphi, Nphi)
	phix = Array{Float64}(Nphi)
	phiy = similar(phix)

	# Translation is a fixed number of indices
	kinc = x2index( left(supp)+1/Jpow, supp, R ) - 1

	slicex_idx = 1:Nphi
	sz = size(PHI)
	for kx in 1:Ncoef
		unsafe_DaubScaling!(phix, kx, PHI, Ncoef, sz, p)
		if p < kx <= Ncoef - p + 1
			slicex_idx += kinc
		end

		slicey_idx = 1:Nphi
		for ky in 1:Ncoef
			unsafe_DaubScaling!(phiy, ky, PHI, Ncoef, sz, p)
			A_mul_Bc!(phi, phiy, phix)

			if p < ky <= Ncoef - p + 1
				slicey_idx += kinc
			end

			z = view(y, slicey_idx, slicex_idx)
			# BLAS.axpy! is using *a lot* of memory
			for j in 1:Nphi, i in 1:Nphi
				@inbounds z[i,j] += coef[ky,kx] * phi[i,j]
			end
		end
	end

	return y
end

@doc """
Overwrite `phi` with the `m`'th scaling function from a matrix `Φ` where
the first `p` columns are the left scaling functions, the last `p` columns
are the right scaling functions and the middle column is the interior
scaling function.
"""->
function unsafe_DaubScaling!( phi::DenseVector, m::Integer, Φ::DenseMatrix,
	Ncoef::Integer, sz::Tuple=size(Φ), p::Integer=div(sz[2]-1,2) )

	if p < m <= Ncoef - p
		unsafe_copy!( phi, 1, Φ, sz[1]*p+1, sz[1] )
	elseif m <= p
		unsafe_copy!( phi, 1, Φ, sz[1]*(m-1)+1, sz[1] )
	else
		unsafe_copy!( phi, 1, Φ, sz[1]*(sz[2]-Ncoef+m-1)+1, sz[1] )
	end
end

@doc """
	allDaubScaling(p::Integer, R::Integer)

Collect the `p` left scaling functions, the interior scaling function
and the `p` right scaling functions as columns in a matrix (in that
order).

All functions are computed at resolution `R`.
"""->
function allDaubScaling(p::Integer, R::Integer)
	1 <= p <= 8 || throw(AssertionError())
	R >= 0 || throw(DomainError())

	lfilter = bfilter(p, 'L')
	int_filter = ifilter(p, true)
	lphi = DaubScaling(lfilter, int_filter, R)'

	iphi = DaubScaling(int_filter, R)

	rfilter = bfilter(p, 'R')
	rphi = DaubScaling(rfilter, int_filter, R)'

	hcat(lphi, iphi, flipdim(rphi,2))
end

@doc """
	IntervalScaling(p::Int, k::Int, J::Int, R::Int) -> x, y

On the interval [0,1] there are 2^J functions at scale `J` ordered such
that the first `p` are left boundary scaling functions, the last `p` are
right boundary scaling functions and the middle 2^J-2p are translated
interior scaling functions.

`IntervalScaling` returns the `k`'th function computed at the dyadic rationals of resolution `R`.

Note that the right boundary scaling functions are ordered differently
than the output of `DaubScaling`:
The output of

	IntervalScaling(p, 2^J, J, R)

is (a translated, dilated and scaled version of) the first column of `Y`
in

	x, Y = DaubScaling(p, 'R', R)
"""->
function IntervalScaling(p::Integer, k::Integer, J::Integer, R::Integer)
	2 <= p <= 8 || throw(AssertionError())
	log2(2*p-1) <= J <= R || throw(AssertionError("The scale is too small compared to p or too large compared to R"))
	0 <= k < (Ny = 2^J) || throw(AssertionError("This interval wavelet does not exist"))

	IF = ifilter(p)
	RmJ = R - J
	DS = DaubSupport(0,1)
	if 0 <= k < p
		BF = bfilter(p,'L')
		Y = DaubScaling(BF, IF, RmJ)
		y = vec( Y[k+1,:] )

		# The indices of the support
		start_idx = x2index(0, DS, R)
		end_idx = x2index( 2.0^-J*(2*p-1), DS, R)
	elseif p <= k < Ny-p
		y = DaubScaling(IF, RmJ)

		start_idx = x2index( 2.0^-J*(-p+k+1), DS, R )
		end_idx = x2index( 2.0^-J*(p+k), DS, R )
	else
		BF = bfilter(p,'R')
		Y = DaubScaling(BF, IF, RmJ)
		y = vec( Y[Ny-k,:] )

		start_idx = x2index( 1-2.0^-J*(2*p-1), DS, R )
		end_idx = x2index( 1, DS, R )
	end

	phi = zeros(Float64, 2^R+1)
	scale!(y, sqrt(Ny))
	phi[start_idx:end_idx] = y

	return phi
end

@doc """
	IntervalScaling(p::Int, J::Int, R::Int) -> matrix

Returns all 2^J scaling functions of scale `J` at resolution `R` on
[0,1] as a `2^R+1`-by-2^J`.

**Note:**
You should probably only use this function for small values of `J`.
"""->
function IntervalScaling(p::Integer, J::Integer, R::Integer)

	Y = Array{Float64}(2^R+1, 2^J)
	for k in 0:2^J-1
		Y[:,k+1] = IntervalScaling(p, k, J, R)
	end

	return Y
end