text:
using Revise
using ADSeismic
using ADCME
using JLD2
using DelimitedFiles
ADCME.options.customop.verbose = false
mpi_init()
M = N = Int(sqrt(mpi_size()))
n = 100

param = MPIElasticPropagatorParams(NX= n * M, NY= n * N, n=n, NSTEP=2000, DELTAT=1e-4/2,
                        DELTAX=1.0, DELTAY=1.0)
compute_PML_Params!(param)
vp = 3300.
vs = 3300. / 1.732
rho = 2800.
λ, μ, ρ = compute_lame_parameters(param, vp, vs, rho)
source = Ricker(param, 50.0, 200.0)

srci = [div(param.NX,5)]
srcj = [div(param.NY,2)]
srctype = [0]
srcv = reshape(source, :, 1)
src = MPIElasticSource(param, srci, srcj, srctype, srcv)

rcvj = Array(20:n*param.M-19)
rcvi = ones(Int64, length(rcvj))*((n*param.N)÷5)
rcvtype = zeros(Int64, length(rcvj))
receiver =  MPIElasticReceiver(param, rcvi, rcvj, rcvtype)

xs = Float64[]
ys = Float64[]
for i = 1:n + 4
    for j = 1:n + 4
        x = (i - 2)*param.DELTAX + (param.II-1) * n * param.DELTAX
        y = (j - 2)*param.DELTAY + (param.JJ-1) * n * param.DELTAY
        push!(xs, x)
        push!(ys, y)
    end
end

init_guess = fc_init([2, 20, 20, 20, 1])
θ = placeholder(init_guess)
θ_ = mpi_bcast(θ)
c = fc([xs ys], [20, 20, 20, 1], θ_) * 1e10
λ = reshape(c, (n+4, n+4))

propagator = MPIElasticPropagatorSolver(param, src, ρ, λ, μ)
MPISimulatedObservation!(propagator, receiver)

@load "data/dat$(mpi_rank())-$(mpi_size()).jld2" rcvv

sess = Session(); init(sess)
@info "starting..."
rcvv, Vx = run(sess, [receiver.rcvv, propagator.vx])

local_loss = sum(propagator.vx[end,:]) * 1e-20
if !isnothing(rcvv)
    @info mpi_rank(), "value"
    global local_loss += sum((receiver.rcvv - rcvv )^2)
end

loss = mpi_sum(local_loss)
G = gradients(loss, θ)

sess = Session(); init(sess)

@info "running loss"
run(sess, loss)
@info "running gradients"
run(sess, G)

timer_fwd = []
timer_bwd = []
for i = 1:10
    @info "forward, $i"
    init(sess)
    d = @timed run(sess, loss)
    push!(timer_fwd, d[2])
end

for i = 1:10
    @info "backward, $i"
    init(sess)
    d = @timed run(sess, G)
    push!(timer_bwd, d[2])
end

timer_fwd = mean(timer_fwd)
timer_bwd = mean(timer_bwd)

if mpi_rank()==0
    open("timing.txt", "a") do io
        writedlm(io, [mpi_size() timer_fwd timer_bwd])
    end
end

mpi_finalize()


text:
# This file contains a type providing
# the main interface of the package

# Precompile the package
__precompile__()

"""
```text
A package for completing spectral correlation analysis of time series.

Important note: this package was developed primarily for
learning purposes and is not advised to be used in production.

Links:
- Repo: https://github.com/paveloom-j/Scats.jl
- Docs: https://paveloom-j.github.io/Scats.jl

Please, use the published docs instead of retrieving
information from these docstrings manually.
```
"""
baremodule Scats

"""
A module containing all inner parts of this package. Is not advisable to be used directly.

See documentation for more info about its contents.
"""
module Internal

# Include source code
include("Prec.jl")          # Precisions and formats of numbers (source code)
include("Extras/Extras.jl") # Extras (source code)
include("Input/Input.jl")   # Input data (source code)
include("Result/Result.jl") # Result data (source code)
include("Gen/Gen.jl")       # Generator (source code)

# Export contents of the modules into Internal
using .Prec   # Precisions and formats of numbers (module)
using .Extras # Extras (module)
using .Input  # Input data (module)
using .Result # Result data (module)
using .Gen    # Generator (module)

end

# Import custom structure types
using .Internal: InputStruct  # A structure type to contain input data
using .Internal: ResultStruct # A structure type to contain result data
using .Internal: GenStruct    # A structure type to contain generator parameters

"""
    API()

Instantiate an instance of Scats API to get access to the public interface.

# Types
- [`Input`](@ref Scats.Internal.Input.InputStruct): input data;
- [`Gen`](@ref Scats.Internal.Gen.GenStruct): generator;
- [`Result`](@ref Scats.Internal.Result.ResultStruct): result data.

# Methods
- for [`Input`](@ref Scats.Internal.Input):
    - [`read_input!`](@ref Scats.Internal.Input.read!)`(file::AbstractString)`:
        read input data from a file;
    - [`write_input`](@ref Scats.Internal.Input.write)`(file::AbstractString)`:
        write input data to a file;
    - [`input_example`](@ref Scats.Internal.Input.example)`(file::AbstractString)`:
        generate an example of the input/output file.

- for [`Gen`](@ref Scats.Internal.Gen):
    - [`read_gen!`](@ref Scats.Internal.Gen.read!)`(file::AbstractString)`:
        read generator parameters from a file;
    - [`gen_example`](@ref Scats.Internal.Gen.example)`(file::AbstractString)`:
        generate an example of a file containing the generator parameters;
    - [`gen!`](@ref Scats.Internal.Gen.gen!)`()`:
        generate time series.

# Usage

```julia
using Scats
s = Scats.API()
```
"""
mutable struct API

    # Input data
    Input::InputStruct      # A structure to contain input data
    read_input!::Function   # Read input data from a file
    write_input::Function   # Write input data to a file
    input_example::Function # Generate an example of the input/output file

    # Time series generator
    Gen::GenStruct        # A structure to contain generator parameters
    read_gen!::Function   # Read generator parameters from a file
    gen_example::Function # Generate an example of a file
                          # containing the generator parameters
    gen!::Function        # Generate time series

    # Result data
    Result::ResultStruct # A structure to contain result data

    # Reset all structures
    reset!::Function

    # Construct an object of this type
    function API()

        # Initialize new object
        this = new()

        # Initialize input data
        this.Input = InputStruct()               # A structure to contain input data
        this.read_input! = this.Input.read!      # Read input data from a file
        this.write_input = this.Input.write      # Write input data to a file
        this.input_example = this.Input.example  # Generate an example
                                                 # of the input/output file

        # Initialize result data
        this.Result = ResultStruct() # A structure to contain result data

        # Initialize time series generator
        this.Gen = GenStruct()              # A structure to contain generator parameters
        this.read_gen! = this.Gen.read!     # Read generator parameters from a file
        this.gen_example = this.Gen.example # Generate an example of a file
                                            # containing the generator parameters

        # Generate time series
        this.gen! = function ()
            this.Gen.gen!(this.Gen, this.Input)
        end

        # Reset all structures
        this.reset! = function ()
            this.Input.reset!()
            this.Result.reset!()
            this.Gen.reset!()
        end

        # Return constructed object
        return this

    end

end

end

ext: jl


text:
using Test

using MultilayerPerceptron

#===== Helpers =====#
function cross_entropy(ŷ::AbstractMatrix{Float64}, y::AbstractMatrix{Float64})::Float64
    n_batches = size(ŷ, 2)
    -sum(@. y * log(ŷ)) / n_batches
end

function binary_cross_entropy(ŷ::Float64, y::Bool)::Float64
    ϵ = eps(ŷ)
    -log(ϵ + (y ? ŷ : 1.0 - ŷ))
end

# https://math.stackexchange.com/a/2503773
function ∇binary_cross_entropy(ŷ::Float64, y::Bool)::Float64
    let y = Float64(y)
        (ŷ - y) / (ŷ - ŷ * ŷ)
    end
end

@testset "Activation" begin
    sz = (5, 3)
    y = rand(Float64, sz) * 10
    ŷ = rand(Float64, sz) * 10
    (e, _) = logit_cross_entropy(ŷ, y)
    @test e ≈ cross_entropy(softmax(ŷ), y)

    len = 5
    ŷ = rand(Float64, 1, len)
    y = rand(Bool, 1, len)
    (e, ∇) = logit_binary_cross_entropy(ŷ, y)
    logistic_ŷ = logistic.(ŷ)
    @test e ≈ sum(binary_cross_entropy.(logistic_ŷ, y)) / len
    @test ∇ ≈ @. ∇binary_cross_entropy(logistic_ŷ, y) * ∇logistic(ŷ, logistic_ŷ)
end

@testset "Xor" begin
    N = 1000

    inputs = unique(Iterators.product(0:1, 0:1))
    train_x = hcat(map(ab -> [Float64.(ab)...], inputs)...)
    train_y = reshape(map(((a, b),) -> Float64(a ⊻ b), inputs), 1, :)

    hidden_layer = Layer((2, 2); activation_fn=tanh, grad_fn=∇tanh)
    output_layer = Layer((1, 2))

    mlp = MLP([hidden_layer, output_layer])
    fit!(mlp, train_x, train_y; learning_rate=0.2, batch_size=2, n_epochs=N)

    test_x, test_y = train_x, train_y
    @test round.(Int, predict(mlp, Float64.(test_x))) == test_y
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
module CBinding
	using Libdl
	using Markdown
	using Scratch

	# TODO: something for Base?
	Base.iterate(::Type{T}, state::Type{<:Tuple} = T) where {T<:Tuple} = state === Tuple{} ? nothing : (Base.tuple_type_head(state), Base.tuple_type_tail(state))

	export @cc_cmd
	export @c_cmd, @c_str
	# export @cxx_cmd, @cxx_str
	export Cbool, Ccomplex, Cmem8, Cmem16, Cmem32, Cmem64, Cmem128, Clongdouble
	export Cptr, Cstruct, Cunion, Cenum, Carray
	export Cfunction, Cbinding, Cconst, Crestrict, Cvolatile
	export isqualifiedwith, unqualifiedtype, bitstype

	const Cbool    = Cuchar
	const Ccomplex = Complex
	primitive type Cmem8 <: Unsigned 8 end
	primitive type Cmem16 <: Unsigned 16 end
	primitive type Cmem32 <: Unsigned 32 end
	primitive type Cmem64 <: Unsigned 64 end
	primitive type Cmem128 <: Unsigned 128 end
	primitive type Clongdouble <: AbstractFloat 2*sizeof(Cdouble)*8 end

	primitive type Cptr{T} <: Ref{T} sizeof(Ptr)*8 end

	abstract type Caggregate end
	abstract type Cstruct <: Caggregate end
	abstract type Cunion <: Caggregate end

	abstract type Cenum end

	# TODO: whenever `incomplete type` becomes available, this can become Carray{T, N} with NTuple{N, T}
	struct Carray{T, N, S} <: AbstractVector{T}
		mem::NTuple{S, UInt8}

		Carray{T, N, S}(::UndefInitializer, mem::NTuple{S, UInt8}) where {T, N, S} = new{T, N, S}(mem)
	end

	struct Cfunction{ret, args<:Tuple, conv}
		let constructor = false end
	end

	abstract type Cbinding{T, lib, name} end
	Base.eltype(::Type{CB}) where {T, CB<:Cbinding{T}} = T

	struct Cvariadic
		let constructor = false end
	end

	abstract type Cqualifier{T} end
	Base.eltype(::Type{CQ}) where {T, CQ<:Cqualifier{T}} = T

	struct Cconst{T} <: Cqualifier{T}
		let constructor = false end
	end

	struct Crestrict{T} <: Cqualifier{T}
		let constructor = false end
	end

	struct Cvolatile{T} <: Cqualifier{T}
		let constructor = false end
	end

	const Cqualifiers{T} = Union{T, Cconst{T}, Crestrict{T}, Cvolatile{T}} where {T}
	const Cptrs{T} = Union{Cptr{T}, Cconst{Cptr{T}}} where {T}
	const Carrays{T, N} = Union{Carray{T, N}, Cconst{Carray{T, N}}} where {T, N}
	const Caggregates = Union{CA, Cconst{CA}} where {CA<:Caggregate}
	const Cenums = Union{CE, Cconst{CE}} where {CE<:Cenum}

	isqualifiedwith(::Type{T}, ::Type{CQ}) where {T, CQ<:Cqualifier} = T <: Cqualifier ? T <: CQ ? true : isqualifiedwith(eltype(T)) : false

	unqualifiedtype(::Type{T}) where {T} = T <: Cconst || T <: Cvolatile || T <: Crestrict ? unqualifiedtype(eltype(T)) : T

	function integraltype(::Type{T}) where {T}
		U = bitstype(T)
		isprimitivetype(U) || return Nothing
		U <: Cenum && return Base.Enums.basetype(U)
		I = sizeof(U) == 1 ? UInt8 : sizeof(U) == 2 ? UInt16 : sizeof(U) == 4 ? UInt32 : sizeof(U) == 8 ? UInt64 : UInt128
		return U <: Signed ? signed(I) : unsigned(I)
	end

	bitstype(::Type{T}) where {T} = isbitstype(T) ? T : error("Failure to obtain a bits-type")
	bitstype(::Type{CQ}) where {CQ<:Cqualifier} = bitstype(unqualifiedtype(CQ))
	bitstype(::Type{CA}) where {T, N, CA<:Carray{T, N}} = Carray{bitstype(T), N, sizeof(bitstype(T))*N}

	values(::Type{CE}) where {CE<:Cenum} = values(bitstype(CE))

	fields(::Type{CA}) where {CA<:Caggregate} = fields(bitstype(CA))

	include("longdouble.jl")
	include("arrays.jl")
	include("bindings.jl")
	include("pointers.jl")
	include("functions.jl")
	include("enums.jl")
	include("aggregates.jl")
	include("accessors.jl")

	include("libclang.jl")

	struct CodeLocation
		file::String
		line::Int
		column::Int
	end

	struct CodeBlock
		loc::LineNumberNode
		lines::UnitRange{Int64}
		flags::NamedTuple
		inlines::Vector{String}
	end

	mutable struct Context{lang}
		mod::Module
		args::Vector{String}
		ind::Union{CXIndex, Nothing}
		tu::Ref{CXTranslationUnit}
		line::Int
		libs::Vector{Pair}
		hdrs::Dict{String, String}
		macros::Dict{String, CXCursor}
		bindings::Vector{Symbol}
		blocks::Vector{CodeBlock}
		src::IOBuffer

		function Context{lang}(mod::Module, args...) where {lang}
			ctx = new{lang}(mod, map(String, collect(args)), nothing, Ref(CXTranslationUnit(C_NULL)), 0, Pair[], Dict{String, String}(), Dict{String, CXCursor}(), Symbol[], CodeBlock[], IOBuffer())
			finalizer(ctx) do x
				x.tu[] == C_NULL || clang_disposeTranslationUnit(x.tu[])
				isnothing(x.ind) || x.ind == C_NULL || clang_disposeIndex(x.ind)
				foreach(lib -> dlclose(lib.first), x.libs)
			end

			configure!(ctx)
			parse!(ctx)

			return ctx
		end
	end

	include("context.jl")
	include("context_c.jl")
	include("context_cxx.jl")
	include("comments.jl")
end


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using RDatasets, StatsPlots
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.23"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
iris = dataset("datasets", "iris");

# ╔═╡ 7c93ed28-250e-11eb-1983-7bf9991291ab
@df iris violin(string.(:Species), :SepalLength,
	fill=:blue, xlabel="Species", ylabel="Sepal Length", legend=false)

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.23"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═7c93ed28-250e-11eb-1983-7bf9991291ab
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
"""
    ConditionalCrayon

Sets conditional formatting for printing to the terminal. Each `ConditionalCrayon`
specifies a function `bad(x)` that returns `true` if the argument `x` is not acceptable,
printing with color `cbad`, and a function `good(x)` that returns `true` if `x` is
acceptable, printing with color `cgood`. If neither are true then `cdefault` is the
color. The `bad` function will be checked before `good`, so takes precedence if both
happen to return `true`. This is not checked.

All colors must be specified as a Crayon type from Crayons.jl

# Constructors
Usage will generally be through the following constructor:

    ConditionalCrayon(badfun::Function, goodfun::Function; [goodcolor, badcolor, defaultcolor])

The colors default to green for good, red for bad, and default color (`Crayon(reset=true)`)
for the default.

For convenience when working with ranges, the following constructor is also provided:

    ConditionalCrayon(lo, hi; [reverse, goodcolor, badcolor, defaultcolor])

Which will be good if the value is less than `lo` and bad if it's higher than `hi`.
This is reversed if `reverse=true`.

A default constructor

    ConditionalCrayon()

Is also provided, which always returns the default color.

# Usage
The `ConditionalCrayon` is a functor object that accepts a single argument of any type,
and returns the `Crayon` according to the output of the `good` and `bad` functions on that
input. It's the user's responsibility to make sure the input type is appropriate for the
functions (since these are all user-specified).
"""
struct ConditionalCrayon
    cgood::Crayon
    cdefault::Crayon
    cbad::Crayon
    bad::Function
    good::Function
end

ConditionalCrayon() = ConditionalCrayon(x->false, x->false)

function ConditionalCrayon(badfun::Function, goodfun::Function;
    goodcolor=crayon"green", badcolor=crayon"red", defaultcolor=Crayon(reset=true)
)
    ConditionalCrayon(goodcolor, defaultcolor, badcolor, badfun, goodfun)
end

function ConditionalCrayon(lo, hi; reverse::Bool=false, kwargs...)
    @assert lo <= hi
    if reverse
        badfun = x->x < lo
        goodfun = x->x >= hi
    else
        badfun = x->x > hi
        goodfun = x->x <= lo
    end
    ConditionalCrayon(badfun, goodfun; kwargs...)
end

function (c::ConditionalCrayon)(v)
    if c.bad(v)
        return c.cbad
    elseif c.good(v)
        return c.cgood
    else
        return c.cdefault
    end
end


text:
export filtersamplemap, filtersamplemap!

"""
    filtersamplemap(x; samples = nothing, experiments = nothing, colnames = nothing)

Filter the sample mapping `DataFrame` to the requested samples, experiments and column names.
`x` can either be a `MultiAssayExperiment` or its [`samplemap`](@ref).

If `samples` is `nothing`, it is not used for any filtering.
Otherwise, it may be a vector or set of strings specifying the samples to retain.
A single string may also be supplied.

If `experiments` is `nothing`, it is not used for any filtering.
Otherwise, it may be a vector or set of strings specifying the experiments to retain.
A single string may also be supplied.

If `colnames` is `nothing`, it is not used for any filtering.
Otherwise, it may be a vector or set of strings specifying the columns to retain.
A single string may also be supplied.

A row of the sample mapping is only retained if it passes all supplied filters.

# Examples
```jldoctest
julia> using MultiAssayExperiments

julia> x = exampleobject();

julia> filtersamplemap(samplemap(x); samples = ["Patient1", "Patient2"])
8×3 DataFrame
 Row │ sample    experiment  colname
     │ String    String      String
─────┼───────────────────────────────
   1 │ Patient1  foo         foo1
   2 │ Patient1  foo         foo2
   3 │ Patient1  foo         foo3
   4 │ Patient2  foo         foo4
   5 │ Patient2  foo         foo5
   6 │ Patient2  foo         foo6
   7 │ Patient2  bar         bar1
   8 │ Patient2  bar         bar2

julia> filtersamplemap(samplemap(x); experiments = "foo")
10×3 DataFrame
 Row │ sample    experiment  colname
     │ String    String      String
─────┼───────────────────────────────
   1 │ Patient1  foo         foo1
   2 │ Patient1  foo         foo2
   3 │ Patient1  foo         foo3
   4 │ Patient2  foo         foo4
   5 │ Patient2  foo         foo5
   6 │ Patient2  foo         foo6
   7 │ Patient3  foo         foo7
   8 │ Patient3  foo         foo8
   9 │ Patient3  foo         foo9
  10 │ Patient4  foo         foo10
```
"""
function filtersamplemap(x::DataFrames.DataFrame; samples = nothing, experiments = nothing, colnames = nothing)
    sampfun = create_filter_function(samples)
    expfun = create_filter_function(experiments)
    colfun = create_filter_function(colnames)
    return filter(row -> sampfun(row.sample) && expfun(row.experiment) && colfun(row.colname), x)
end

function create_filter_function(input)
    if input == nothing
        return x -> true
    end

    if isa(input, AbstractString)
        return x -> (x == input)
    end

    if isa(input, AbstractSet{<:AbstractString})
        return x -> (x in input)
    end

    inputset = Set(input)
    return x -> (x in inputset)
end

function filtersamplemap(x::MultiAssayExperiment; kwargs...)
    y = copy(x)
    filtersamplemap!(y; kwargs...)
    return y
end

"""
    filtersamplemap!(x; samples = nothing, experiments = nothing, colnames = nothing)

Modifies `samplemap(x)` in place by filtering based on [`filtersamplemap`](@ref).
A reference to the modified `x` is returned.

# Examples
```jldoctest
julia> using MultiAssayExperiments

julia> x = exampleobject();

julia> filtersamplemap!(x; samples = ["Patient1", "Patient2"]);

julia> samplemap(x)
8×3 DataFrame
 Row │ sample    experiment  colname
     │ String    String      String
─────┼───────────────────────────────
   1 │ Patient1  foo         foo1
   2 │ Patient1  foo         foo2
   3 │ Patient1  foo         foo3
   4 │ Patient2  foo         foo4
   5 │ Patient2  foo         foo5
   6 │ Patient2  foo         foo6
   7 │ Patient2  bar         bar1
   8 │ Patient2  bar         bar2
```
"""
function filtersamplemap!(x::MultiAssayExperiment; samples = nothing, experiments = nothing, colnames = nothing)
    x.samplemap = filtersamplemap(samplemap(x); samples = samples, experiments = experiments, colnames = colnames)
    return x
end


text:
function RouwenNormalize!(A)
    A[2:end-1, :] ./= 2
    nothing
end

function rouwenupdate!(Πold, Πnew, p, q)
    Nold = size(Πold, 1)
    Πnew[1:Nold, 1:Nold] = p*Πold
    Πnew[1:Nold, 2:(Nold+1)] .+= (1-p)*Πold
    Πnew[2:(Nold+1), 1:Nold] .+= (1-q)*Πold
    Πnew[2:(Nold+1), 2:(Nold+1)] .+= q*Πold
    nothing
end

function rouwenmat(N, ρ, p, q)
    initmat = [p 1-p;1-q q]
    if N==2
        ret = initmat
    else
        currmat = initmat
        for n = 3:N
            nextmat = fill(zero(ρ), (n, n))
            rouwenupdate!(currmat, nextmat, p, q)
            currmat = nextmat
            RouwenNormalize!(currmat)
        end
        ret = currmat
    end
    return(ret)
end

"""
    rouwenhorst(npts, ρ, σ)

Discretizes AR(1) process `y′ = ρ y +  σ ϵ`.

# Example

npts = 10
rho = 0.9
sig = 0.01
points, Π = rouwenhorst(npts, rho, sig)

"""
function rouwenhorst(npts, ρ, σ)
    ω = σ/sqrt(1-ρ^2) # long run std dev
    # println("long run variance: $ω")
    q = (1+ρ)/2

    points = range(-ω*sqrt(npts - 1),
                   ω*sqrt(npts-1),
                   length=npts)

    # points = GenGrid(uniform(npts, bds))
    Π = (npts > 1) ? rouwenmat(npts, ρ, q, q) : fill(one(ρ), (1,1))

    return points, Π
end

"Returns a MarkovChain with discretized AR(1) from applying `rouwenhorst`"
rouwenhorst_mc(npts, ρ, σ)=MarkovChain(rouwenhorst(npts, ρ, σ)...)


text:
# NOTE: This shouldn't require an additional implementation.
#       Just make sure that ϵᵢ== 0, and never update ϵ.


text:
"""
Perform one level of the maximal overlap discrete wavelet transform (MODWT) on
the the vector `v`, which is either the `j`th level scaling coefficients or,
if `j==1`, the raw signal.  The vectors `h` and `g` are the MODWT detail and
scaling filters.

Returns a tuple `(v, w)` of the scaling and detail coefficients at level `j+1`.
"""
function modwt_step(v::AbstractVector{T}, j::Integer, h::Array{S,1},
        g::Array{S,1}) where {T <: Number, S <: Number}
    N = length(v)
    L = length(h)
    v1 = zeros(T, N)
    w1 = zeros(T, N)
    for t in 1:N
        k = t
        w1[t] = h[1] * v[k]
        v1[t] = g[1] * v[k]
        for n in 2:L
            k -= 2^(j-1)
            if k <= 0
                k = mod1(k, N)
            end
            w1[t] += h[n] * v[k]
            v1[t] += g[n] * v[k]
        end
    end
    return v1, w1
end

"""
Perform a maximal overlap discrete wavelet transform (MODWT) of the array
`x`.  The wavelet type `wt` determines the transform type and the wavelet
class, see `wavelet`. (Note that the wavelet filter coefficients are scaled
by 1/√2 for the MODWT so that the transform maintains unit energy).

The number of transform levels `L` can be any number <= `maxmodwttransformlevels(x)`
(the default value).

Returns an `n × L+1` matrix (where `n` is the length of `x`) with the wavelet
coefficients for level j in column j.  The scaling coefficients are in the
last (L+1th) column.
"""
function modwt(x::AbstractVector{T}, wt::OrthoFilter,
        L::Integer=maxmodwttransformlevels(x)) where T <: Number
    L <= maxmodwttransformlevels(x) ||
        throw(ArgumentError("Too many transform levels (length(x) < 2^L)"))
    L >= 1 || throw(ArgumentError("L must be >= 1"))
    g, h = WT.makereverseqmfpair(wt)
    g /= sqrt(2)
    h /= sqrt(2)
    N = length(x)
    W = zeros(T, N, L)
    V = deepcopy(x)
    for j in 1:L
        V[:], W[:, j] = modwt_step(V, j, h, g)
    end
    return [W V]
end

"""
Perform one level of the inverse maximal overlap discrete wavelet transform
(MODWT). Takes the `j`th-level scaling and detail coefficients `v` and `w`
and returns a vector of the `j-1`th level scaling coefficients. The vectors
`h` and `g` are the MODWT detail and scaling filters.
"""
function imodwt_step(v::AbstractVector{T}, w::AbstractVector{T}, j::Integer,
        h::Array{S,1}, g::Array{S,1}) where {T<:Number, S<:Number}
    length(v) == length(w) ||
        throw(DimensionMismatch("Input array sizes must match"))
    length(h) == length(g) ||
        throw(DimensionMismatch("Filter sizes must match"))
    N = length(v)
    L = length(h)
    v0 = zeros(T, N)
    for t in 1:N
        k = t
        v0[t] = h[1] * w[k] + g[1] * v[k]
        for n in 2:L
            k += 2^(j-1)
            if k > N
                k = mod1(k, N)
            end
            v0[t] += h[n] * w[k] + g[n] * v[k]
        end
    end
    return v0
end

"""
Perform an inverse maximal overlap discrete wavelet transform (MODWT) of `xw`,
the inverse of `modwt(x, wt, L)`.
"""
function imodwt(xw::Array{T, 2}, wt::OrthoFilter) where T <: Number
    g, h = WT.makereverseqmfpair(wt)
    g /= sqrt(2)
    h /= sqrt(2)
    N, L = size(xw)
    x = deepcopy(xw[:, L])
    for j in L-1:-1:1
        x[:] = imodwt_step(x, xw[:, j], j, h, g)
    end
    return x
end


text:
abstract type AbstractLeakyDAFESN <: AbstractLeakyESN end

struct dafESN{T, S<:AbstractArray{T}, I, B, F1, F2, N} <: AbstractLeakyDAFESN
    res_size::I
    in_size::I
    out_size::I
    train_data::S
    alpha::T
    nla_type::N
    first_activation::F1
    second_activation::F2
    first_lambda::T
    second_lambda::T
    W::S
    W_in::S
    states::S
    extended_states::B
end

function dafESN(approx_res_size::Int,
        train_data::AbstractArray{T},
        degree::Int,
        radius::T,
        first_lambda::T,
        second_lambda::T;
        first_activation::Any = tanh,
        second_activation::Any = tanh,
        sigma::T = 0.1,
        alpha::T = 1.0,
        nla_type::NonLinearAlgorithm = NLADefault(),
        extended_states::Bool = false) where T<:AbstractFloat

    in_size = size(train_data, 1)
    out_size = size(train_data, 1)
    res_size = Int(floor(approx_res_size/in_size)*in_size)

    W = init_reservoir_givendeg(res_size, radius, degree)
    W_in = init_input_layer(res_size, in_size, sigma)
    states = daf_states_matrix(W, W_in, train_data, alpha,
    first_activation, second_activation, first_lambda, second_lambda, extended_states)

    return dafESN{T, typeof(train_data),
        typeof(res_size),
        typeof(extended_states),
        typeof(first_activation),
        typeof(second_activation),
        typeof(nla_type)}(res_size, in_size, out_size, train_data,
    alpha, nla_type, first_activation, second_activation, first_lambda,
    second_lambda, W, W_in, states, extended_states)
end

#reservoir matrix W given by the user
function dafESN(W::AbstractArray{T},
        train_data::AbstractArray{T},
        first_lambda::T,
        second_lambda::T;
        first_activation::Any = tanh,
        second_activation::Any = tanh,
        sigma::T = 0.1,
        alpha::T = 1.0,
        nla_type::NonLinearAlgorithm = NLADefault(),
        extended_states::Bool = false) where T<:AbstractFloat

    in_size = size(train_data, 1)
    out_size = size(train_data, 1)
    res_size = size(W, 1)

    W_in = init_input_layer(res_size, in_size, sigma)
    states = daf_states_matrix(W, W_in, train_data, alpha,
    first_activation, second_activation, first_lambda, second_lambda, extended_states)

    return dafESN{T, typeof(train_data),
        typeof(res_size),
        typeof(extended_states),
        typeof(first_activation),
        typeof(second_activation),
        typeof(nla_type)}(res_size, in_size, out_size, train_data,
    alpha, nla_type, first_activation, second_activation, first_lambda,
    second_lambda, W, W_in, states, extended_states)
end

#input layer W_in given by the user
function dafESN(approx_res_size::Int,
        train_data::AbstractArray{T},
        degree::Int,
        radius::T,
        first_lambda::T,
        second_lambda::T,
        W_in::AbstractArray{T};
        first_activation::Any = tanh,
        second_activation::Any = tanh,
        alpha::T = 1.0,
        nla_type::NonLinearAlgorithm = NLADefault(),
        extended_states::Bool = false) where T<:AbstractFloat

    in_size = size(train_data, 1)
    out_size = size(train_data, 1)
    res_size = Int(floor(approx_res_size/in_size)*in_size)
    W = init_reservoir_givendeg(res_size, radius, degree)

    if size(W_in, 1) != res_size
        throw(DimensionMismatch(W_in, "size(W_in, 1) must be equal to size(W, 1)"))
    elseif size(W_in, 2) != in_size
        throw(DimensionMismatch(W_in, "size(W_in, 2) must be equal to in_size"))
    end

    states = daf_states_matrix(W, W_in, train_data, alpha,
    first_activation, second_activation, first_lambda, second_lambda, extended_states)

    return dafESN{T, typeof(train_data),
        typeof(res_size),
        typeof(extended_states),
        typeof(first_activation),
        typeof(second_activation),
        typeof(nla_type)}(res_size, in_size, out_size, train_data,
    alpha, nla_type, first_activation, second_activation, first_lambda,
    second_lambda, W, W_in, states, extended_states)
end

#reservoir matrix W and input layer W_in given by the user
"""
    dafESN(W::AbstractArray{T}, train_data::AbstractArray{T}, first_lambda::T,
    second_lambda::T, W_in::AbstractArray{T} [, first_activation::Any,
    second_activation::Any, alpha::T, nla_type::NonLinearAlgorithm, extended_states::Bool])

Build a double activation function ESN struct given the input and reservoir matrices, as described in [1].

[1]: Lun, Shu-Xian, et al. "A novel model of leaky integrator echo state network for time-series prediction." Neurocomputing 159 (2015): 58-66.
"""
function dafESN(W::AbstractArray{T},
        train_data::AbstractArray{T},
        first_lambda::T,
        second_lambda::T,
        W_in::AbstractArray{T};
        first_activation::Any = tanh,
        second_activation::Any = tanh,
        alpha::T = 1.0,
        nla_type::NonLinearAlgorithm = NLADefault(),
        extended_states::Bool = false) where T<:AbstractFloat

    in_size = size(train_data, 1)
    out_size = size(train_data, 1)
    res_size = size(W, 1)

    if size(W_in, 1) != res_size
        throw(DimensionMismatch(W_in, "size(W_in, 1) must be equal to size(W, 1)"))
    elseif size(W_in, 2) != in_size
        throw(DimensionMismatch(W_in, "size(W_in, 2) must be equal to in_size"))
    end

    states = daf_states_matrix(W, W_in, train_data, alpha,
    first_activation, second_activation, first_lambda, second_lambda, extended_states)

    return dafESN{T, typeof(train_data),
        typeof(res_size),
        typeof(extended_states),
        typeof(first_activation),
        typeof(second_activation),
        typeof(nla_type)}(res_size, in_size, out_size, train_data,
    alpha, nla_type, first_activation, second_activation, first_lambda,
    second_lambda, W, W_in, states, extended_states)
end

function daf_states_matrix(W::AbstractArray{Float64},
        W_in::AbstractArray{Float64},
        train_data::AbstractArray{Float64},
        alpha::Float64,
        first_activation::Function,
        second_activation::Function,
        first_lambda::Float64,
        second_lambda::Float64,
        extended_states::Bool)

    train_len = size(train_data, 2)
    res_size = size(W, 1)
    in_size = size(train_data, 1)
    states = zeros(Float64, res_size, train_len)
    for i=1:train_len-1
        states[:, i+1] = double_leaky_fixed_rnn(alpha, first_activation, second_activation, first_lambda, second_lambda, W, W_in, states[:, i], train_data[:, i])
    end

    if extended_states == true
        ext_states = vcat(states, hcat(zeros(Float64, in_size), train_data[:,1:end-1]))
        return ext_states
    else
        return states
    end

    return states
end

"""
    dafESNpredict(esn::AbstractLeakyESN, predict_len::Int, W_out::AbstractArray{Float64})

Return the prediction for a given length of the constructed dafESN struct.
"""
function dafESNpredict(esn::AbstractLeakyDAFESN,
    predict_len::Int,
    W_out::AbstractArray{Float64})

    output = zeros(Float64, esn.in_size, predict_len)
    x = esn.states[:, end]

    if esn.extended_states == false
        for i=1:predict_len
            x_new = nla(esn.nla_type, x)
            out = (W_out*x_new)
            output[:, i] = out
            x = double_leaky_fixed_rnn(esn.alpha, esn.first_activation, esn.second_activation, esn.first_lambda, esn.second_lambda, esn.W, esn.W_in, x, out)
        end
    elseif esn.extended_states == true
        for i=1:predict_len
            x_new = nla(esn.nla_type, x)
            out = (W_out*x_new)
            output[:, i] = out
            x = vcat(double_leaky_fixed_rnn(esn.alpha, esn.first_activation, esn.second_activation, esn.first_lambda, esn.second_lambda, esn.W, esn.W_in, x[1:esn.res_size], out), out)
        end
    end

    return output
end

"""
    dafESNpredict_h_steps(esn::AbstractLeakyESN, predict_len::Int, h_steps::Int,
    test_data::AbstractArray{Float64}, W_out::AbstractArray{Float64})

Return the prediction for h steps ahead of the constructed ESN struct.
"""
function dafESNpredict_h_steps(esn::AbstractLeakyESN,
    predict_len::Int,
    h_steps::Int,
    test_data::AbstractArray{Float64},
    W_out::AbstractArray{Float64})

    output = zeros(Float64, esn.in_size, predict_len)
    x = esn.states[:, end]

    if esn.extended_states == false
        for i=1:predict_len
            x_new = nla(esn.nla_type, x)
            out = (W_out*x_new)
            output[:, i] = out
            if mod(i, h_steps) == 0
                x = double_leaky_fixed_rnn(esn.alpha, esn.first_activation, esn.second_activation, esn.first_lambda, esn.second_lambda, esn.W, esn.W_in, x, test_data[:,i])
            else
                x = double_leaky_fixed_rnn(esn.alpha, esn.first_activation, esn.second_activation, esn.first_lambda, esn.second_lambda, esn.W, esn.W_in, x, out)
            end
        end
    elseif esn.extended_states == true
        for i=1:predict_len
            x_new = nla(esn.nla_type, x)
            out = (W_out*x_new)
            output[:, i] = out
            if mod(i, h_steps) == 0
                x = vcat(double_leaky_fixed_rnn(esn.alpha, esn.first_activation, esn.second_activation, esn.first_lambda, esn.second_lambda, esn.W, esn.W_in, x[1:esn.res_size], test_data[:,i]), test_data[:,i])
            else
                x = vcat(double_leaky_fixed_rnn(esn.alpha, esn.first_activation, esn.second_activation, esn.first_lambda, esn.second_lambda, esn.W, esn.W_in, x[1:esn.res_size], out), out)
            end
        end
    end
    return output
end


text:
# rule_generation.jl

struct Rule
    p::Array{Int32,1}
    q::Int32
    supp::Float64
    conf::Float64
    lift::Float64
end

function update_support_cnt!(supp_dict::Dict, nd::Node)
    supp_dict[nd.item_ids] = nd.supp
end

function grow_support_dict!(supp_cnt::Dict{Array{Int32,1}, Int}, node::Node)
    if has_children(node)
        for nd in node.children
            update_support_cnt!(supp_cnt, nd)
            grow_support_dict!(supp_cnt, nd)
        end
    end
end

# This function generates a dictionary whose keys are the frequent
# itemsets (their integer represenations, actually), and whose values
# are the support count for the given itemset. This function is used
# for computing support, confidence, and lift of association rules.
function gen_support_dict(root::Node, num_transacts)
    supp_cnt = Dict{Array{Int32, 1}, Int}()
    supp_cnt[Int32[]] = num_transacts
    grow_support_dict!(supp_cnt, root)
    return supp_cnt
end

# Given a single node in a frequent item tree, this function generates all the
# rules for that node. This does not include rules for the node's children.
function gen_node_rules(node::Node, supp_dict::Dict{Array{Int32,1}, Int}, k, num_transacts, minconf)
    lhs_keep = trues(k)
    rules = Array{Rule, 1}(undef, 0)
    for i = 1:k
        lhs_keep[i] = false
        if i > 1
            lhs_keep[i-1] = true
        end
        p = node.item_ids[lhs_keep]
        supp = node.supp/num_transacts
        # for debugging
        if supp == 0.0
            warn("zero supp for: ", node.item_ids)
            return rules            # NOTE: returning early if zero support
        else
            conf = supp/((supp_dict[node.item_ids[lhs_keep]])/num_transacts)
            if conf ≥ minconf
                rhs_keep = .!lhs_keep
                q = first(node.item_ids[rhs_keep])
                lift = conf/((supp_dict[node.item_ids[rhs_keep]])/num_transacts)

                rule = Rule(p, q, supp, conf, lift)
                push!(rules, rule)
            end
        end
    end
    rules
end

function gen_rules!(rules::Array{Rule, 1}, node::Node, supp_dict::Dict{Array{Int32, 1}, Int}, k, num_transacts, minconf)
    m = length(node.children)

    for child in node.children
        rules_tmp = gen_node_rules(child, supp_dict, k, num_transacts, minconf)
        append!(rules, rules_tmp)
        if !isempty(child.children)
            gen_rules!(rules, child, supp_dict, k+1, num_transacts, minconf)
        end
    end
end

function gen_rules(root::Node, supp_dict::Dict{Array{Int32, 1}, Int}, num_transacts, minconf)
    rules = Array{Rule, 1}(undef, 0)
    n_kids = length(root.children)
    if n_kids > 0
        for i = 1:n_kids
            gen_rules!(rules, root.children[i], supp_dict, 2, num_transacts, minconf)
        end
    end
    rules
end

function rules_to_dataframe(rules::Array{Rule, 1}, item_lkup::Dict{T, S}; join_str = " | ") where {T <: Integer, S}
    n_rules = length(rules)
    dt = DataFrame(lhs = fill("", n_rules),
                   rhs = fill("", n_rules),
                   supp = zeros(n_rules),
                   conf = zeros(n_rules),
                   lift = zeros(n_rules))
    for i = 1:n_rules
        lhs_items = map(x -> string.(item_lkup[x]), rules[i].p)

        lhs_string = "{" * join(lhs_items, join_str) * "}"
        dt[i, :lhs] = lhs_string
        dt[i, :rhs] = string.(item_lkup[rules[i].q])
        dt[i, :supp] = rules[i].supp
        dt[i, :conf] = rules[i].conf
        dt[i, :lift] = rules[i].lift
    end
    dt
end

"""
    apriori(transactions; supp, conf, maxlen)
Given an array of transactions (a vector of string vectors), this function runs the a-priori
algorithm for generating frequent item sets. These frequent items are then used to generate
association rules. The `supp` argument allows us to stipulate the minimum support
required for an itemset to be considered frequent. The `conf` argument allows us to exclude
association rules without at least `conf` level of confidence. The `maxlen` argument stipulates
the maximum length of an association rule (i.e., total items on left- and right-hand sides)
"""
function apriori(transactions::Array{Array{S, 1}, 1}; supp::Float64 = 0.01, conf = 0.8, maxlen::Int = 5) where S
    n = length(transactions)
    uniq_items = unique_items(transactions)
    item_lkup = Dict{Int32, S}()
    for (i, itm) in enumerate(uniq_items)
        item_lkup[i] = itm
    end
    minsupp = floor(Int, supp * n)
    if minsupp == 0
        minsupp += 1
    end
    freq_tree = frequent_item_tree(transactions, uniq_items, minsupp, maxlen)
    supp_lkup = gen_support_dict(freq_tree, n)
    rules = gen_rules(freq_tree, supp_lkup, n, conf)
    rules_dt = rules_to_dataframe(rules, item_lkup)
    return rules_dt
end

"""
apriori(occurrences, item_lkup; supp, conf, maxlen)

Given an boolean occurrence matrix of transactions (rows are transactions, columns are items) and
a lookup dictionary of column-index to items-string, this function runs the a-priori
algorithm for generating frequent item sets. These frequent items are then used to generate
association rules. The `supp` argument allows us to stipulate the minimum support
required for an itemset to be considered frequent. The `conf` argument allows us to exclude
association rules without at least `conf` level of confidence. The `maxlen` argument stipulates
the maximum length of an association rule (i.e., total items on left- and right-hand sides)
"""
function apriori(occurrences::BitArray{2}; supp::Float64 = 0.01, conf = 0.8, maxlen::Int = 5)
    n = size(occurrences, 1)
    minsupp = floor(Int, supp * n)
    if minsupp == 0
        minsupp = 1
    end
    freq_tree = frequent_item_tree(occurrences, minsupp, maxlen)
    supp_lkup = gen_support_dict(freq_tree, n)
    rules = gen_rules(freq_tree, supp_lkup, n, conf)

    return rules
end


text:
# Routines for observations of light sources in particular images,
# rather than for sources in the abstract, in physical units,
# and rather than for images alone (that's image_model.jl).

"""
Attributes of the patch of sky surrounding a single
celestial object in a single image.

Attributes:
  - center: The approximate source location in world coordinates
  - radius_pix: The width of the influence of the object in pixel coordinates
  - psf: The point spread function in this region of the sky
  - wcs_jacobian: The jacobian of the WCS transform in this region of the
                  sky for each band
  - pixel_center: The pixel location of center in each band.
"""
immutable SkyPatch
    center::Vector{Float64}
    radius_pix::Float64

    psf::Vector{PsfComponent}
    wcs_jacobian::Matrix{Float64}
    pixel_center::Vector{Float64}
end

"""
Initialize a SkyPatch from an existing SkyPatch and a new PSF.
"""
SkyPatch(patch::SkyPatch, psf::Vector{PsfComponent}) =
    SkyPatch(patch.center, patch.radius_pix, psf, patch.wcs_jacobian,
             patch.pixel_center)

"""Centers of patches in pixel coordinates"""
patch_ctrs_pix(patches::Vector{SkyPatch}) = [p.pixel_center for p in patches]

"""Radii of patches in pixel coordinates"""
patch_radii_pix(patches::Vector{SkyPatch}) = [p.radius_pix for p in patches]

"""
A fast function to determine which sources might belong to which tiles.

Args:
  - tiles: A TiledImage
  - patch_ctrs: Vector of length-2 vectors, giving pixel center of each patch.
  - patch_radii_px: Radius of each patch.

Returns:
  - An array (over tiles) of a vector of candidate
    source patches.  If a patch is a candidate, it may be within the patch
    radius of a point in the tile, though it might not.
"""
function local_source_candidates(tile::ImageTile,
                                 patch_ctrs::Vector{Vector{Float64}},
                                 patch_radii_px::Vector{Float64})
    @assert length(patch_ctrs) == length(patch_radii_px)

    ret = Int[]

    # Find the patches that are less than the radius plus diagonal from the
    # center of the tile.  These are candidates for having some
    # overlap with the tile.
    tile_center = (mean(tile.h_range), mean(tile.w_range))
    h_width, w_width = size(tile.pixels)
    tile_diag = (0.5 ^ 2) * (h_width ^ 2 + w_width ^ 2)

    for s in 1:length(patch_ctrs)
        patch_dist = (tile_center[1] - patch_ctrs[s][1])^2
                    + (tile_center[2] - patch_ctrs[s][2])^2
        if patch_dist <= (tile_diag + patch_radii_px[s])^2
            push!(ret, s)
        end
    end

    return ret
end

"""
Args:
  - tile: An ImageTile (containing tile coordinates)
  - patch_ctrs: Vector of length-2 vectors, giving pixel center of each patch.
  - patch_radii_px: Radius of each patch.

Returns:
  - A vector of source ids (from 1 to length(patches)) that influence
    pixels in the tile.  A patch influences a tile if
    there is any overlap in their squares of influence.
"""
function get_local_sources(tile::ImageTile,
                           patch_ctrs::Vector{Vector{Float64}},
                           patch_radii_px::Vector{Float64})
    @assert length(patch_ctrs) == length(patch_radii_px)
    tile_source_map = Int[]
    tile_ctr = (mean(tile.h_range), mean(tile.w_range))
    h_width, w_width = size(tile.pixels)

    for i in eachindex(patch_ctrs)
        patch_ctr = patch_ctrs[i]
        patch_r = patch_radii_px[i]

        # This is a "ball" in the infinity norm.
        if ((abs(tile_ctr[1] - patch_ctr[1]) < patch_r + 0.5 * h_width) &&
            (abs(tile_ctr[2] - patch_ctr[2]) < patch_r + 0.5 * w_width))
            push!(tile_source_map, i)
        end
    end

    tile_source_map
end

"""
Get the sources associated with each tile in a TiledImage.

Args:
  - tiled_image: A TiledImage
  - patch_ctrs: Vector of length-2 vectors, giving pixel center of each patch.
  - patch_radii_px: Radius of each patch.
Returns:
  - An array (same dimensions as the tiles) of vectors of indices
    into patches indicating which patches are affected by any pixels
    in the tiles.
"""
function get_sources_per_tile(image_tiles::Matrix{ImageTile},
                              patch_ctrs::Vector{Vector{Float64}},
                              patch_radii_px::Vector{Float64})
    out = similar(image_tiles, Vector{Int})

    HH, WW = size(image_tiles)
    for ww in 1:WW, hh in 1:HH
        cands = local_source_candidates(image_tiles[hh, ww],
                                        patch_ctrs,
                                        patch_radii_px)
        # get indicies in cands that truly overlap the tile.
        idx = get_local_sources(image_tiles[hh, ww],
                                patch_ctrs[cands],
                                patch_radii_px[cands])
        out[hh, ww] = cands[idx]
    end

    return out
end

function choose_patch_radius(ce::CatalogEntry,
                             b::Int64,
                             psf_width::AbstractFloat,
                             epsilon::AbstractFloat;
                             width_scale=1.0)
    # The galaxy scale is the point with half the light -- if the light
    # were entirely in a univariate normal, this would be at 0.67 standard
    # deviations.  We are being a bit conservative here.
    obj_width =
      ce.is_star ? psf_width: width_scale * ce.gal_scale / 0.67 + psf_width

    flux = ce.is_star ? ce.star_fluxes[b] : ce.gal_fluxes[b]
    @assert flux > 0.

    # Choose enough pixels that the light is either 90% of the light
    # would be captured from a 1d gaussian or 5% of the sky noise,
    # whichever is a larger radius.
    pdf_90 = exp(-0.5 * (1.64)^2) / (sqrt(2pi) * obj_width)
    pdf_target = min(pdf_90, epsilon / (20 * flux))
    rhs = log(pdf_target) + 0.5 * log(2pi) + log(obj_width)
    sqrt(-2 * (obj_width ^ 2) * rhs)
end

"""
Choose a reasonable patch radius based on the catalog.
TODO: Select this by rendering the object and solving an optimization
problem.

Args:
  - pixel_center: The pixel location of the object
  - ce: The catalog entry for the object
  - psf: The psf at this location
  - img: The image
  - width_scale: A multiple of standard deviations to use
  - max_radius: The maximum radius in pixels.

Returns:
  - A radius in pixels chosen from the catalog entry.
"""
function choose_patch_radius(
            pixel_center::Vector{Float64},
            ce::CatalogEntry,
            psf::Array{PsfComponent},
            img::TiledImage;
            width_scale=1.0,
            max_radius=100)
    psf_width = Model.get_psf_width(psf, width_scale=width_scale)

    # Get the average sky noise around a source
    close_tile = get_containing_tile(pixel_center, img)
    epsilon = mean(close_tile.epsilon_mat)

    radius_req = choose_patch_radius(ce, img.b, psf_width, epsilon;
                                        width_scale=width_scale)
    min(radius_req, max_radius)
end


text:
using GPAR: is_closed_downwards

@testset "util" begin
    @testset "is_closed_downwards" begin
        @test is_closed_downwards(randn(10, 2))

        X = fill!(Matrix{Union{Float64, Missing}}(undef, 10, 3), 5.0)
        @test is_closed_downwards(X)

        X[1, end] = missing
        @test is_closed_downwards(X)

        X[1, end-1] = missing
        @test is_closed_downwards(X)

        X[1, end-2] = missing
        @test is_closed_downwards(X)

        X[2, 2] = missing
        @test !is_closed_downwards(X)
    end
end


text:
################################################################################
function init(::Type{Custom}, config::ConfigCustom, objective::Objective, difftune)
    return Custom(ModelWrappers.ℓDensityResult(objective), config.proposal)
end

function infer(
    _rng::Random.AbstractRNG,
    diagnostics::Type{AbstractDiagnostics},
    kernel::Custom,
    model::ModelWrapper,
    data::D,
) where {D}
    return DiagnosticsCustom
end

################################################################################
# Export
export infer, init


text:
struct Weight <: EdgeModifier
    w
end

struct Start <: EdgeModifier
    s
end

struct Observed <: NodeModifier end
struct Ordinal <: NodeModifier end

@testset "Multiplication of Edges" begin
    @test Edge(:a, :b) * Weight(1) == ModifiedEdge(Edge(:a, :b), Weight(1))
    @test ModifiedEdge(Edge(:a, :b), Weight(1)) == Weight(1) * Edge(:a, :b)
    @test Edge(:a, :b) * Weight(1) == Weight(1) * Edge(:a, :b)
end

@testset "Multiplication of Nodes" begin
    @test ModifyingNode(Node(:a), Weight(2)) == Weight(2) * Node(:a)
    @test Node(:a) * Weight(3) == Weight(3) * Node(:a)
    @test ModifiedEdge(Edge(:a, :b), Weight(1)) == Edge(:a * Weight(1), :b)
    @test Edge(:a * Weight(1), :b) == Edge(Weight(1) * :a, :b)
    @test ModifiedEdge(Edge(:a, :b), [Weight(1), Start(1)]) == Edge(Weight(1) * :a, :b * Start(1))
end

@testset "Multiplication of Modifiers" begin
    @test (ModifiedEdge(Edge(:a, :b), [Weight(1) Start(1)])) ==
        (Edge(:a, Weight(1) * Start(1) * :b)) ==
        (Edge(:a, Weight(1) * :b * Start(1))) ==
        (Edge(:a, :b * Weight(1) * Start(1))) ==
        (Edge(:a, :b) * Weight(1) * Start(1)) ==
        (Weight(1) * Edge(:a, :b) * Start(1)) ==
        (Edge(:a, :b) * (Weight(1) * Start(1)))

    @test (ModifyingNode(:b, [Weight(1) Start(1)])) ==
        (Weight(1) * Start(1) * :b) ==
        (Weight(1) * :b * Start(1)) ==
        (:b * Weight(1) * Start(1)) ==
        (Node(:b) * Weight(1) * Start(1)) ==
        (Weight(1) * Node(:b) * Start(1)) ==
        (Node(:b) * (Weight(1) * Start(1)))

    @test Start(3)*Start(4)*Weight(5) ==
        [Start(3), Start(4), Weight(5)] ==
        (Start(3)*Start(4))*Weight(5) ==
        Start(3)*(Start(4)*Weight(5))
end

@testset "ModifiedNode" begin
    @test Node(:a)^Observed() == Observed()^Node(:a) == ModifiedNode(Node(:a), Observed())
    @test Node(:a)^Observed()^Ordinal() == Node(:a)^[Observed(), Ordinal()] == ModifiedNode(Node(:a), [Observed(), Ordinal()])
    @test Node(:a)^Observed()^Ordinal()^Ordinal() == ModifiedNode(Node(:a), [Observed(), Ordinal()])
    @test Observed()^Observed()^Ordinal() ==
        [Observed(), Observed(), Ordinal()] ==
        (Observed()^Observed())^Ordinal() ==
        Observed()^(Observed()^Ordinal())
end

@testset "Modification of Arrow" begin
    @test unarrow((:a → :b) * Start(1)) == unarrow(:a → Start(1) * :b)
    @test unarrow(([:a :b] → [:c :d]) * Start(1)) == unarrow([:a :b] → [:c :d] .* Ref(Start(1)))
end


text:
"""
    IntervalMesh <: AbstractMesh

A 1D mesh on an `IntervalDomain`.

# Constuctors

    IntervalMesh(domain::IntervalDomain, faces::AbstractVector)

Construct a 1D mesh with face locations at `faces`.

    IntervalMesh(domain::IntervalDomain[, stetching=Uniform()]; nelems=)

Constuct a 1D mesh on `domain` with `nelems` elements, using `stretching`. Possible values of `stretching` are:

- [`Uniform()`](@ref)
- [`ExponentialStretching(H)`](@ref)
"""
struct IntervalMesh{I <: IntervalDomain, V <: AbstractVector} <: AbstractMesh1D
    domain::I
    faces::V
end

domain(mesh::IntervalMesh) = mesh.domain
nelements(mesh::IntervalMesh) = length(mesh.faces) - 1
elements(mesh::IntervalMesh) = Base.OneTo(nelements(mesh))

function Base.show(io::IO, mesh::IntervalMesh)
    print(io, nelements(mesh), "-element IntervalMesh of ")
    print(io, mesh.domain)
end

coordinates(mesh::IntervalMesh, elem::Integer, vert::Integer) =
    mesh.faces[elem + vert - 1]

function coordinates(
    mesh::IntervalMesh,
    elem::Integer,
    (ξ,)::Tuple{T},
) where {T}
    ca = mesh.faces[elem]
    cb = mesh.faces[elem + 1]
    Geometry.linear_interpolate((ca, cb), ξ)
end

function containing_element(mesh::IntervalMesh, coord)
    i = min(searchsortedlast(mesh.faces, coord), nelements(mesh))
    lo = Geometry.component(mesh.faces[i], 1)
    hi = Geometry.component(mesh.faces[i + 1], 1)
    val = Geometry.component(coord, 1)
    ξ = ((val - lo) + (val - hi)) / (hi - lo)
    return i, (ξ,)
end
function is_boundary_face(mesh::IntervalMesh, elem::Integer, face)
    !Domains.isperiodic(mesh.domain) &&
        ((elem == 1 && face == 1) || (elem == nelements(mesh) && face == 2))
end

function boundary_face_name(mesh::IntervalMesh, elem::Integer, face)
    if !Domains.isperiodic(mesh.domain)
        if elem == 1 && face == 1
            return mesh.domain.boundary_names[1]
        elseif elem == nelements(mesh) && face == 2
            return mesh.domain.boundary_names[2]
        end
    end
    return nothing
end

abstract type StretchingRule end

"""
    Uniform()

Use uniformly-sized elements.
"""
struct Uniform <: StretchingRule end

function IntervalMesh(
    domain::IntervalDomain{CT},
    ::Uniform = Uniform();
    nelems,
) where {CT <: Geometry.Abstract1DPoint{FT}} where {FT}
    faces = range(domain.coord_min, domain.coord_max; length = nelems + 1)
    IntervalMesh(domain, faces)
end

"""
    ExponentialStretching(H)

Apply exponential stretching to the domain when constructing elements. `H` is
the scale height (a typical atmospheric scale height `H ≈ 7.5`km).

For an interval ``[z_0,z_1]``, this makes the elements uniformally spaced in
``\\zeta``, where
```math
\\zeta = \\frac{1 - e^{-\\eta/h}}{1-e^{-1/h}},
```
where ``\\eta = \\frac{z - z_0}{z_1-z_0}``, and ``h = \\frac{H}{z_1-z_0}`` is
the non-dimensional scale height.
"""
struct ExponentialStretching{FT} <: StretchingRule
    H::FT
end

function IntervalMesh(
    domain::IntervalDomain{CT},
    stretch::ExponentialStretching;
    nelems,
) where {CT <: Geometry.Abstract1DPoint{FT}} where {FT}
    cmin = Geometry.component(domain.coord_min, 1)
    cmax = Geometry.component(domain.coord_max, 1)
    R = cmax - cmin
    h = stretch.H / R
    η(ζ) = -h * log1p((expm1(-1 / h)) * ζ)
    faces =
        [CT(cmin + R * η(ζ)) for ζ in range(FT(0), FT(1); length = nelems + 1)]
    IntervalMesh(domain, faces)
end


text:
# ################################################################
# Quadratic solver abstract configuration
# ################################################################
#
abstract type Abstract_BC_QuadSolver_Conf end


text:
function fetchsolution(app,master,dmd)

nt = length(app.dt);
if nt==1
    Uout = getsolution("dataout/out",dmd,master.npe);
else
    if app.soltime != []
        app.soltime = 1:nt;
    end
    fn = "dataout/out_t" * string(app.soltime[1]);
    tmp = getsolution(fn,dmd,master.npe);
    if app.wave==1
        w = getsolution("dataout/out_wdg_t" * string(app.soltime[1]),dmd,master.npe);
        tmp = hcat(tmp,w);
    end
    Uout = zeros((size(tmp,1), size(tmp,2), size(tmp,3), length(app.soltime)));
    Uout[:,:,:,1] = tmp;
    for i = 2:length(app.soltime)
        fn = "dataout/out_t" * string(app.soltime[i]);
        tmp = getsolution(fn,dmd,master.npe);
        if app.wave==1
            w = getsolution("dataout/out_wdg_t" * string(app.soltime[i]),dmd,master.npe);
            tmp = hcat(tmp,w);
        end
        Uout[:,:,:,i] = tmp;
    end
end

return Uout

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
function Base.map!(
    op::UnaryUnion, C::GBArray, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    mask, accum = _handlenothings(mask, accum)
    desc === nothing && (desc = DEFAULTDESC)
    op = getoperator(op, eltype(A))
    accum = getaccum(accum, eltype(C))
    A, desc = _handletranspose(A, desc)
    if C isa GBVector && A isa GBVector
        libgb.GrB_Vector_apply(C, mask, accum, op, A, desc)
    elseif C isa GBMatrix && A isa GBMatrix
        libgb.GrB_Matrix_apply(C, mask, accum, op, A, desc)
    end
    return C
end

Base.map!(op::Function, C::GBArray, A::GBArray; mask = nothing, accum = nothing, desc = nothing) =
    map!(UnaryOp(op), C, A; mask, accum, desc)

function Base.map!(
    op::UnaryUnion, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    return map!(op, A, A; mask, accum, desc)
end

Base.map!(op::Function, A::GBArray; mask = nothing, accum = nothing, desc = nothing) =
    map!(UnaryOp(op), A; mask, accum, desc)

function Base.map(
    op::UnaryUnion, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    t = inferoutputtype(A, op)
    return map!(op, similar(A, t), A; mask, accum, desc)
end
Base.map(op::Function, A::GBArray; mask = nothing, accum = nothing, desc = nothing) =
    map(UnaryOp(op), A; mask, accum, desc)

function Base.map!(
    op::BinaryUnion, C::GBArray, x, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    mask, accum = _handlenothings(mask, accum)
    desc === nothing && (desc = DEFAULTDESC)
    op = getoperator(op, optype(eltype(A), typeof(x)))
    accum = getaccum(accum, eltype(C))
    _, desc, A = _handletranspose(nothing, desc, A)
    if C isa GBVector && A isa GBVector
        libgb.scalarvecapply1st[optype(typeof(x), eltype(A))](C, mask, accum, op, x, A, desc)
    elseif C isa GBMatrix && A isa GBMatrix
        libgb.scalarmatapply1st[optype(typeof(x), eltype(A))](C, mask, accum, op, x, A, desc)
    end
    return C
end

function Base.map!(
    op::Function, C::GBArray, x, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    map!(BinaryOps.BinaryOp(op), C, x, A; mask, accum, desc)
end

function Base.map!(
    op::BinaryUnion, x, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    return map!(op, A, x, A; mask, accum, desc)
end

function Base.map!(
    op::Function, x, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    map!(BinaryOps.BinaryOp(op), x, A; mask, accum, desc)
end

function Base.map(
    op::BinaryUnion, x, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    t = inferoutputtype(A, op)
    return map!(op, similar(A, t), x, A; mask, accum, desc)
end

function Base.map(
    op::Function, x, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    map(BinaryOps.BinaryOp(op), x, A; mask, accum, desc)
end

function Base.map!(
    op::BinaryUnion, C::GBArray, A::GBArray, x;
    mask = nothing, accum = nothing, desc = nothing
)
    mask, accum = _handlenothings(mask, accum)
    desc === nothing && (desc = DEFAULTDESC)
    op = getoperator(op, optype(eltype(A), typeof(x)))
    accum = getaccum(accum, eltype(C))
    A, desc, _ = _handletranspose(A, desc)
    if C isa GBVector && A isa GBVector
        libgb.scalarvecapply2nd[optype(typeof(x), eltype(A))](C, mask, accum, op, A, x, desc)
    elseif C isa GBMatrix && A isa GBMatrix
        libgb.scalarmatapply2nd[optype(typeof(x), eltype(A))](C, mask, accum, op, A, x, desc)
    end
    return C
end

function Base.map!(
    op::Function, C::GBArray, A::GBArray, x;
    mask = nothing, accum = nothing, desc = nothing
)
    map!(BinaryOps.BinaryOp(op), C, A, x; mask, accum, desc)
end

function Base.map!(
    op::BinaryUnion, A::GBArray, x;
    mask = nothing, accum = nothing, desc = nothing
)
    return map!(op, A, A, x; mask, accum, desc)
end

function Base.map!(
    op::Function, A::GBArray, x;
    mask = nothing, accum = nothing, desc = nothing
)
    map!(BinaryOps.BinaryOp(op), A, A, x; mask, accum, desc)
end

function Base.map(
    op::BinaryUnion, A::GBArray, x;
    mask = nothing, accum = nothing, desc = nothing
)
    t = inferoutputtype(A, op)
    return map!(op, similar(A, t), A, x; mask, accum, desc)
end

function Base.map(
    op::Function, A::GBArray, x;
    mask = nothing, accum = nothing, desc = nothing
)
    map(BinaryOps.BinaryOp(op), A, x; mask, accum, desc)
end

Base.:+(x::valid_union, u::GBArray; mask = nothing, accum = nothing, desc = nothing) =
    map(BinaryOps.PLUS, x, u; mask, accum, desc)
Base.:+(u::GBArray, x::valid_union; mask = nothing, accum = nothing, desc = nothing) =
    map(BinaryOps.PLUS, u, x; mask, accum, desc)

Base.:*(x::valid_union, u::GBArray; mask = nothing, accum = nothing, desc = nothing) =
    map(BinaryOps.TIMES, x, u; mask, accum, desc)
Base.:*(u::GBArray, x::valid_union; mask = nothing, accum = nothing, desc = nothing) =
    map(BinaryOps.TIMES, u, x; mask, accum, desc)

Base.:-(x::valid_union, u::GBArray; mask = nothing, accum = nothing, desc = nothing) =
    map(BinaryOps.MINUS, x, u; mask, accum, desc)
Base.:-(u::GBArray, x::valid_union; mask = nothing, accum = nothing, desc = nothing) =
    map(BinaryOps.MINUS, u, x; mask, accum, desc)

"""
    map(op::UnaryOp, A::GBArray; kwargs...)::GBArray
    map(op::BinaryOp, A::GBArray, x; kwargs...)::GBArray
    map(op::BinaryOp, x, A::GBArray, kwargs...)::GBArray

Transform a GBArray by applying `op` to each element.

UnaryOps apply elementwise in the usual fashion.
BinaryOps require the additional argument `x` which is substituted as the first or second
argument of `op` depending on its position.

# Arguments
- `op::MonoidBinaryOrRig = BinaryOps.PLUS`: the binary operation which is applied such that
    `C[i,j] = op(A[i,j], B[i,j])` for all `i,j` present in either `A` and/or `B`.
- `A::GBArray`: `GBVector` or optionally transposed `GBMatrix`.
- `x`: Position dependent argument to binary operators.

# Keywords
- `mask::Union{Nothing, GBMatrix} = nothing`: optional mask.
- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation
    where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied.
- `desc = nothing`
"""
function Base.map(
    op::AbstractOp, A::GBArray, x = nothing;
    mask = nothing, accum = nothing, desc = nothing
)
end

"""
    mask!(C::GBArray, A::GBArray, mask::GBArray)

Apply a mask to matrix `A`, storing the results in C.

"""
function mask!(C::GBArray, A::GBArray, mask::GBArray; structural = false, complement = false)
    desc = DEFAULTDESC
    structural && (desc = desc + S)
    complement && (desc = desc + C)
    map!(identity, C, A; mask, desc)
    return C
end

"""
    mask(A::GBArray, mask::GBArray)

Apply a mask to matrix `A`.
"""
function mask(A::GBArray, mask::GBArray; structural = false, complement = false)
    return mask!(similar(A), A, mask; structural, complement)
end


text:
########### Utilitaires ###########

function distance(p1, p2)
    # Retourne la distance (euclidienne) entre les deux points
    dx = abs(p1.x - p2.x)
    dy = abs(p1.y - p2.y)
    return sqrt(dx^2 + dy^2)
end

function creationDistancier(nbPoints, listePoints)
    D = Array{Float64,2}(undef, nbPoints, nbPoints)
    for i in 1:nbPoints
        for j in 1:nbPoints
            D[i,j] = distance(listePoints[i], listePoints[j])
        end
    end
    return D
end

function coutTour(ordreVisite, instance)
    # Retourne le coût du tour passé en paramètre avec le distancier
    distance = 0
    for i in 1:size(ordreVisite,1)-1
        distance += instance.D[ordreVisite[i], ordreVisite[i+1]]
    end
    return (distance / instance.vitesseCamion)
end

########### Contrôles de saisie utilisateur ###########

function choixListe(borneSup, message=" --> Votre choix : ")
    # Permet de demander à l'utilisateur de faire un choix entre plusieurs items [[1, borneSup]]
    liste = []
    for i in 1:borneSup
        push!(liste, string(i))
    end

    choix = "-1"
    tentative = 1
    while !(choix in liste)
        # Tant que l'utilisateur ne saisit pas un numéro correct, on redemande
        if tentative == 1 # Première tentative
            print("\n")
            print(message)
        else
            print(message)
        end
        choix = readline()
        tentative += 1
    end
    return choix
end

function choixBinaire(message)
    # Permet de demander à l'utilisateur de faire un choix binaire (oui ou non)
    choix = "-1"
    tentative = 1
    while choix != "o" && choix != "n"
        # Tant que l'utilisateur ne saisit pas une lettre correcte, on redemande
        if tentative == 1 # Première tentative
            print(message)
        else
            print(message)
        end
        choix = readline()
        tentative += 1
    end
    return choix
end

########### Affichages ###########

function caracteristiquesInstance(instance)
    # Permet d'afficher les différents paramètres d'une instance
    println("\n========== Caractéristiques de l'instance ==========")
    println("\n Nombre de points (dépôt inclus) : " * string(instance.nbPoints))
    println(" Vitesse du camion : " * string(instance.vitesseCamion))
    println(" Vitesse du drone : " * string(instance.vitesseDrone) * "\n")
    choix = choixBinaire(" --> Souhaitez-vous afficher la liste des points (o/n) ? ")
    if choix == "o"
        println("\n  Liste des points :")
        println("  ------------------")
        for i in 1:instance.nbPoints
            println("    ", i, ") x = ", instance.listePoints[i].x, " | y = ", instance.listePoints[i].y)
        end
        println("\n")
    end

    choix = choixBinaire(" --> Souhaitez-vous afficher le distancier (o/n) ? ")
    if choix == "o"
        println("\n  Distancier :")
        print("  ------------")
        for i in 1:instance.nbPoints
            print("\n    ")
            for j in 1:instance.nbPoints
                print(round(instance.D[i,j], digits=2), "  ")
            end
        end
        println("")
    end
    println("\n========== Caractéristiques de l'instance ==========")
end

function affichageSolution(instance)
    println("\n========== Détail de la solution ==========\n")
    println("  Temps de parcours total : " * string(instance.solution.tempsParcours))
    println("  Temps d'attente total : " * string(instance.solution.tempsAttente))
    println("\n  Ordre de visite du camion : ")
    println("  --------------------------- \n")
    if length(instance.solution.ordreVisiteCamion) <= 10
        print("  ")
        for i in 1:length(instance.solution.ordreVisiteCamion)-1
            print(string(instance.solution.ordreVisiteCamion[i]) * " --> ")
        end
        println(string(instance.solution.ordreVisiteCamion[end]))
    else
        print("  ")
        for i in 1:9
            print(string(instance.solution.ordreVisiteCamion[i]) * " --> ")
        end
        print("... --> ")
        println(string(instance.solution.ordreVisiteCamion[end-1]) * " --> " * string(instance.solution.ordreVisiteCamion[end]))
        println("\n /!\\ Affichage partiel (trop de clients) /!\\")
    end

    if length(instance.solution.operationsDrone) != 0
        println("\n  Liste des opérations du drone : ")
        println("  ------------------------------- \n")
        op = instance.solution.operationsDrone
        if length(op) <= 10
            for i in 1:length(op)
                #println(" $i) Décollage : $(op[i][1]); Livraison : $(op[i][2]); Atterrissage : $(op[i][3])")
                println(" $i) $(op[i][1]) --drone--> $(op[i][2]) --drone--> $(op[i][3])")
            end
        else
            for i in 1:9
                #println(" $i) Décollage : $(op[i][1]); Livraison : $(op[i][2]); Atterrissage : $(op[i][3])")
                println(" $i) $(op[i][1]) --drone--> $(op[i][2]) --drone--> $(op[i][3])")
            end
            println(" ...")
            println(" ...")
            println(" ...")
            println(" $(length(op))) $(op[(length(op))][1]) --drone--> $(op[(length(op))][2]) --drone--> $(op[(length(op))][3])")
            println("\n /!\\ Affichage partiel (trop d'opérations) /!\\")
        end
    end
    println("\n========== Détail de la solution ==========")
end

function initGraphique(instance)
    graphique = plt.figure() # Sortie graphique unique
    affichageGraphique(instance, graphique)
    return graphique
end

function majGraphique(instance, graphique)
    clf() # On supprime l'ancienne figure pour afficher la nouvelle
    affichageGraphique(instance, graphique)
end

function affichageGraphique(instance, fig=plt.figure(), animation=false, col="blue")
    # Affiche une solution de TSP, éventuellement juste les points si ordreVisite vide
    ax = fig.gca()
    if instance.codeResolution == 0
        #title("Localisation du dépôt et des clients")
        title(L"Customers and depot locations ($speed_{t} = $"*string(instance.vitesseCamion)*L", $speed_{d} = $"*string(instance.vitesseDrone)*")")
    elseif instance.codeResolution == 1
        #title("Solution optimale TSP - Concorde")
        title(L"TSP optimal solution - Concorde ($speed_{t} = $"*string(instance.vitesseCamion)*L", $speed_{d} = $"*string(instance.vitesseDrone)*")")
    elseif instance.codeResolution == 2
        #title("Solution admissible TSP - Lin-Kernighan")
        title(L"TSP solution - Lin-Kernighan ($speed_{t} = $"*string(instance.vitesseCamion)*L", $speed_{d} = $"*string(instance.vitesseDrone)*")")
    elseif instance.codeResolution == 3
        #title(L"Solution admissible TSP-D - AEP ($vitesse_{t} = $"*string(instance.vitesseCamion)*L", $vitesse_{d} = $"*string(instance.vitesseDrone)*")")
        title(L"TSP-D solution - A1 ($speed_{t} = $"*string(instance.vitesseCamion)*L", $speed_{d} = $"*string(instance.vitesseDrone)*")")
    elseif instance.codeResolution == 4
        #title(L"Solution admissible TSP-D - EP ($vitesse{t} = $"*string(instance.vitesseCamion)*L", $vitesse{d} = $"*string(instance.vitesseDrone)*")")
        title(L"TSP-D solution - A2 ($speed_{t} = $"*string(instance.vitesseCamion)*L", $speed_{d} = $"*string(instance.vitesseDrone)*")")
    elseif instance.codeResolution == 5
        #title(L"Solution admissible TSP-D - EP ($vitesse{t} = $"*string(instance.vitesseCamion)*L", $vitesse{d} = $"*string(instance.vitesseDrone)*")")
        title(L"TSP-D solution - 2A1 ($speed_{t} = $"*string(instance.vitesseCamion)*L", $speed_{d} = $"*string(instance.vitesseDrone)*")")
    end
    xlabel(L"x ($z_1 = $"*string(round(instance.solution.tempsParcours, digits=2))*L", $z_2 = $"*string(round(instance.solution.tempsAttente, digits=2))*")")
    ylabel("y")
    grid(true)
    xCoordClients = []
    yCoordClients = []
    #ax.plot(instance.listePoints[1].x, instance.listePoints[1].y, ls="", marker="x", color="green", label="Dépôt")
    ax.plot(instance.listePoints[1].x, instance.listePoints[1].y, ls="", marker="x", color="green", label="Depot", picker=true)
    for i in 2:instance.nbPoints
        push!(xCoordClients, instance.listePoints[i].x)
        push!(yCoordClients, instance.listePoints[i].y)
    end
    #ax.plot(xCoordClients, yCoordClients, ls="", marker="x", color="red", label="Clients")
    ax.plot(xCoordClients, yCoordClients, ls="", marker="x", color="red", label="Customers", picker=true)

    for i in 1:instance.nbPoints
        ax.text(instance.listePoints[i].x+1.5, instance.listePoints[i].y+1.5, instance.listePoints[i].id)
    end

    if instance.solution.ordreVisiteCamion != []
        plot(0, 0, linewidth=0.7, color="blue", label="Truck") # Pour forcer l'affichage de la légende
        for i in 2:size(instance.solution.ordreVisiteCamion,1)
            p1 = instance.listePoints[instance.solution.ordreVisiteCamion[i-1]]
            p2 = instance.listePoints[instance.solution.ordreVisiteCamion[i]]
            connecterPointsCamion(p1, p2, col)
            if animation == true
                pause(0.5)
            end
        end
    end

    if instance.solution.operationsDrone != []
        plot(0, 0, "--", linewidth=0.7, color="orange", label="Drone") # Pour forcer l'affichage de la légende
        for i in 1:size(instance.solution.operationsDrone,1)
            p1 = instance.listePoints[instance.solution.operationsDrone[i][1]]
            p2 = instance.listePoints[instance.solution.operationsDrone[i][2]]
            p3 = instance.listePoints[instance.solution.operationsDrone[i][3]]
            connecterPointsDrone(p1, p2)
            if p3 != p1 # Si p3 == p1 alors il s'agit d'un aller-retour, on ne trace qu'un seul trait en pointillé
                connecterPointsDrone(p2, p3)
            end
        end
    end
    legend(loc="best")
end

function connecterPointsCamion(p1, p2, col="blue")
    # Trace une ligne en trait plein entre deux points (clients) pour le camion
    plot([p1.x, p2.x], [p1.y, p2.y], linewidth=0.7, color=col)
end

function connecterPointsDrone(p1, p2)
    # Trace une ligne en pointillé entre deux points pour le drone
    plot([p1.x, p2.x], [p1.y, p2.y], "--", linewidth=0.7, color="orange")
end

########### Reset ###########

function reset()
    close() # On ferme la figure courante
    run(`clear`) # On efface le contenu de la console
end

########### Etude vitesse ###########

max_stars_repo_name: xgandibleux/tspd-solver

max_issues_repo_name: xgandibleux/tspd-solver

max_forks_repo_name: xgandibleux/tspd-solver


text:
# big integers and floats: pass as strings
function mlput(ml::MLink,n::BigInt)
    mlputnext(ml,MLTKINT)
    mlput(ml,string(n))
end
mlget(ml::MLink,::Type{BigInt}) = BigInt(mlget(ml,String))

function mlput(ml::MLink,x::BigFloat)
    mlputnext(ml,MLTKREAL)
    mlput(ml,string(x))
end

# a wrapper around a string for storing Mathematica arbitrary precision floats
type MLFloatStr
    str::ASCIIString
end

function mlput(ml::MLink,x::MLFloatStr)
    mlputnext(ml,MLTKREAL)
    mlput(ml,x.str)
end

function convert(::Type{BigFloat}, x::MLFloatStr)
    exp_ind = search(x.str,'*')
    prc_ind = search(x.str,'`')
    if exp_ind == 0
        bstr = prc_ind == 0 ? x.str : x.str[1:prc_ind-1]
    else
        sstr = prc_ind == 0 ? x.str[1:exp_ind-1] : x.str[1:prc_ind-1]
        estr = x.str[exp_ind+2:end]
        bstr = sstr*"e"*estr
    end
    BigFloat(bstr)
end

mlget(ml::MLink,::Type{MLFloatStr}) = MLFloatStr(mlget(ml,ASCIIString))
mlget(ml::MLink,::Type{BigFloat}) = convert(BigFloat,mlget(ml,MLFloatStr))

# default representation
mlget(ml::MLink,::Type{String}) = mlget(ml,UTF8String)
mlget(ml::MLink,::Type{Integer}) = mlget(ml,BigInt)
mlget(ml::MLink,::Type{FloatingPoint}) = mlget(ml,BigFloat)

# automatically handle types
mlget(ml::MLink) = mlget(ml, token_type[mlgetrawtype(ml)])

const symbolsub = {
             :+ => :Plus,
             :- => :Subtratc,
             :* => :Times,
             :/ => :Divide,
             :^ => :Power,
             :> => :Greater,
             :>= => :GreaterEqual,
             :< => :Less,
             :<= => :LessEqual,
             }

# Expressions
function mlput(ml, ex::Expr)
    if ex.head == :call
        fn = ex.args[1]
        mlput(ml, MLFunction(get(symbolsub, fn, fn), length(ex.args)-1))
        for arg = ex.args[2:end]
            mlput(ml, arg)
        end
    elseif ex.head == :comparison
        fn = ex.args[2]
        mlput(ml, MLFunction(get(symbolsub, fn, fn), (length(ex.args) + 1) >> 1))
        mlput(ml, ex.args[1])
        for i = 2:2:length(ex.args)
            if ex.args[i] == fn
                mlput(ml, ex.args[i+1])
            else
                error("comparison operators should all be the same")
            end
        end
    elseif ex.head == :vcat || ex.head == :cell1d
        mlput(ml, MLFunction(:List, length(ex.args)))
        for arg = ex.args
            mlput(ml, arg)
        end
    elseif ex.head == :(=)
        mlput(ml, MLFunction(:Set, length(ex.args)))
        for arg = ex.args
            mlput(ml, arg)
        end
    elseif ex.head == :(=>)
        mlput(ml, MLFunction(:Rule, length(ex.args)))
        for arg = ex.args
            mlput(ml, arg)
        end
    end
end

# automatically construct expression
function mlget(ml, ::Type{Expr})
    y = mlget(ml)
    if isa(y, MLFunction)
        args = [mlget(ml, Expr) for i = 1:y.nargs]
        return Expr(:call,y.name,args...)
    else
        return y
    end
end


text:
using Mamba

## Data
leuk = Dict{Symbol, Any}(
  :t_obs =>
    [1, 1, 2, 2, 3, 4, 4, 5, 5, 8, 8, 8, 8, 11, 11, 12, 12, 15, 17, 22, 23, 6,
     6, 6, 6, 7, 9, 10, 10, 11, 13, 16, 17, 19, 20, 22, 23, 25, 32, 32, 34, 35],
  :fail =>
    [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
     1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0],
  :Z =>
    [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
     0.5, 0.5, 0.5, 0.5, 0.5, 0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
     -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
     -0.5, -0.5],
  :t => [1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 13, 15, 16, 17, 22, 23, 35]
)
leuk[:N] = N = length(leuk[:t_obs])
leuk[:T] = T = length(leuk[:t]) - 1

leuk[:Y] = Array{Int}(N, T)
leuk[:dN] = Array{Int}(N, T)
for i in 1:N
  for j in 1:T
    leuk[:dN][i, j] = leuk[:fail][i] * (leuk[:t_obs][i] == leuk[:t][j])
    leuk[:Y][i, j] = Int(leuk[:t_obs][i] >= leuk[:t][j])
  end
end

leuk[:c] = 0.001
leuk[:r] = 0.1

## Model Specification
model = Model(

  dN = Stochastic(2,
    (Y, beta, Z, dL0, N, T) ->
      UnivariateDistribution[
        Y[i, j] > 0 ? Poisson(exp(beta * Z[i]) * dL0[j]) : Flat()
        for i in 1:N, j in 1:T
      ],
    false
  ),

  mu = Logical(1,
    (c, r, t) -> c * r * (t[2:end] - t[1:(end - 1)]),
    false
  ),

  dL0 = Stochastic(1,
    (mu, c, T) -> UnivariateDistribution[Gamma(mu[j], 1 / c) for j in 1:T],
    false
  ),

  beta = Stochastic(
    () -> Normal(0, 1000)
  ),

  S0 = Logical(1,
    dL0 -> exp(-cumsum(dL0)),
    false
  ),

  S_treat = Logical(1,
    (S0, beta) -> S0.^exp(-0.5 * beta)
  ),

  S_placebo = Logical(1,
    (S0, beta) -> S0.^exp(0.5 * beta)
  )

)

## Initial Values
inits = [
  Dict(:dN => leuk[:dN], :beta => 0, :dL0 => fill(1, leuk[:T])),
  Dict(:dN => leuk[:dN], :beta => 1, :dL0 => fill(2, leuk[:T]))
]

## Sampling Scheme
scheme = [AMWG(:dL0, 0.1),
          Slice(:beta, 3.0)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, leuk, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule BLST_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("BLST")
JLLWrappers.@generate_main_file("BLST", UUID("07413bf9-d2c6-5469-9e83-945f30a2c5ad"))
end  # module BLST_jll


text:
using Ipopt
using PowerModelsPrivacyPreserving
using PowerModels
using Distributions
using Random

const PMs = PowerModels
const PMPP = PowerModelsPrivacyPreserving

ipopt = Ipopt.Optimizer
file =  "test/data/matpower/case5.m"
data_unpert = parse_file(file)
data_min_loss = deepcopy(data_unpert)
data_min_cost = deepcopy(data_unpert)

"this is the canonical OPF problem with generation cost minimization"
result_unpert_cost = PMPP.run_ac_opf_cost(data_unpert, ipopt)
PMs.print_summary(result_unpert_cost["solution"])
PMPP.calculate_losses!(result_unpert_cost, data_unpert)
"store faithfulness info"
data_min_cost["cost"] = Dict()
data_min_cost["cost"]["value"] = result_unpert_cost["objective"]
data_min_cost["cost"]["beta"] = 1

"this variant of the OPF problem minimizes grid losses instead of generation cost"
result_unpert_loss = PMPP.run_ac_opf_loss(data_unpert, ipopt)
PMs.print_summary(result_unpert_loss["solution"])
PMPP.calculate_losses!(result_unpert_loss, data_unpert)
"store faithfulness info"
data_min_loss["loss"] = Dict()
data_min_loss["loss"]["value"] = result_unpert_loss["totalloss"]
data_min_loss["loss"]["beta"] = 1

# Set constant parameters
α = 0.01
ϵ = 1
λ = 50

# Add impedance perturbation to both data dictionaries.
data_pert_min_loss = PMPP.create_impedance_perturbation(data_min_loss, α, ϵ, λ)
data_pert_min_cost = PMPP.create_impedance_perturbation(data_min_cost, α, ϵ, λ)

# Calculate the result
result_pert_loss = PMPP.run_opf_variable_impedance_loss(data_pert_min_loss, ipopt)
PMs.print_summary(result_pert_loss["solution"])
PMPP.calculate_losses!(result_pert_loss, data_pert_min_loss)

result_pert_cost = PMPP.run_opf_variable_impedance_cost(data_pert_min_cost, ipopt)
PMs.print_summary(result_pert_cost["solution"])
PMPP.calculate_losses!(result_pert_cost, data_pert_min_cost)


text:
import PyPlot
plt = PyPlot
import InsarTimeseries
include("./testl1.jl")
include("./point_analysis.jl")

rowcol2 = [245, 189]  # small/looked uplift
# rowcol2 = [1224, 944]  # full uplift
# rowcol3 = [368, 131]  # Small/looked subs
# rowcol3 = [2229, 2236]  # TESTING for diffs
rowcol3 = [312, 135]  # Subs line that disappears with long baseline
# rowcol3 = [920, 329]  # Subs line that disappears with long baseline (big map)
rowcol4 = [32, 352]  # Subs around the decorrelated top

# rowcol = [1632, 687]  # Test near pecos, disappears with too much filtering
# rowcol2 = [1478, 743]  # bowl of strong subsidence that appears in 400-day soln, not in 600+
# rowcol3 = [1649, 699]  # patch of subsidence that appears in 400-day soln, not in 600+

# geolistig, intlist500ig, valid_igram_indices500ig = load_geolist_intlist("unw_stack.h5", "geolist_ignore.txt", 500)
geolist, intlist400, valid_igram_indices400 =
    load_geolist_intlist("unw_stack.h5", "geolist_ignore.txt", 400)
_, intlist500, valid_igram_indices500 =
    load_geolist_intlist("unw_stack.h5", "geolist_ignore.txt", 500)
_, intlist600, valid_igram_indices600 =
    load_geolist_intlist("unw_stack.h5", "geolist_ignore.txt", 600)
_, intlist700, valid_igram_indices700 =
    load_geolist_intlist("unw_stack.h5", "geolist_ignore.txt", 700)
_, intlist800, valid_igram_indices800 =
    load_geolist_intlist("unw_stack.h5", "geolist_ignore.txt", 800)
_, intlistall, valid_igram_indicesall =
    load_geolist_intlist("unw_stack.h5", "geolist_ignore.txt", 1500)
# geolist, intlist400, valid_igram_indices400 = load_geolist_intlist("unw_stack.h5", nothing, 400)

B400 = InsarTimeseries.build_B_matrix(geolist, intlist400);
Blin400 = sum(B400, dims = 2);
B500 = InsarTimeseries.build_B_matrix(geolist, intlist500);
Blin500 = sum(B500, dims = 2);
Blin600 = sum(InsarTimeseries.build_B_matrix(geolist, intlist600), dims = 2);
Blin700 = sum(InsarTimeseries.build_B_matrix(geolist, intlist700), dims = 2);
Blin800 = sum(InsarTimeseries.build_B_matrix(geolist, intlist800), dims = 2);
Blinall = sum(InsarTimeseries.build_B_matrix(geolist, intlistall), dims = 2);

println("point 0")

unw_vals500_uplift = get_stack_vals(
    "unw_stack.h5",
    rowcol2...,
    1,
    "stack_flat_shifted",
    valid_igram_indices500,
    reference_station = nothing,
);
unw_valsall_uplift = get_stack_vals(
    "unw_stack.h5",
    rowcol2...,
    1,
    "stack_flat_shifted",
    valid_igram_indicesall,
);
cc500_uplift = get_stack_vals("cc_stack.h5", rowcol2..., 1, "stack", valid_igram_indices500);

println("point 1")

unw_vals500_subs = get_stack_vals(
    "unw_stack.h5",
    rowcol3...,
    1,
    "stack_flat_shifted",
    valid_igram_indices500,
)
cc500_subs = get_stack_vals("cc_stack.h5", rowcol3..., 1, "stack", valid_igram_indices500);
unw_valsall_subs = get_stack_vals(
    "unw_stack.h5",
    rowcol3...,
    1,
    "stack_flat_shifted",
    valid_igram_indicesall,
);
ccall_subs = get_stack_vals("cc_stack.h5", rowcol3..., 1, "stack", valid_igram_indicesall);

println("point 2")
unw_valsall_txoz =
    get_stack_vals("unw_stack.h5", "TXOZ", 1, "stack_flat_shifted", valid_igram_indicesall);
unw_vals500_txoz = get_stack_vals(
    "unw_stack.h5",
    "TXOZ",
    1,
    "stack_flat_shifted",
    valid_igram_indices500,
    reference_station = nothing,
);
cc500_txoz = get_stack_vals("cc_stack.h5", "TXOZ", 1, "stack", valid_igram_indices500);
unw_vals500_nmhb = get_stack_vals(
    "unw_stack.h5",
    "NMHB",
    1,
    "stack_flat_shifted",
    valid_igram_indices500,
    reference_station = nothing,
);
cc500_nmhb = get_stack_vals("cc_stack.h5", "NMHB", 1, "stack", valid_igram_indices500);

println("point 3")
unw_vals500_txmh = get_stack_vals(
    "unw_stack.h5",
    "TXMH",
    1,
    "stack_flat_shifted",
    valid_igram_indices500,
    reference_station = nothing,
);
cc500_txmh = get_stack_vals("cc_stack.h5", "TXMH", 1, "stack", valid_igram_indices500);
cc500_txkm = get_stack_vals("cc_stack.h5", "TXKM", 1, "stack", valid_igram_indices500);
#
println("done ")

_, intlist600, valid_igram_indices600 =
    load_geolist_intlist("unw_stack.h5", "geolist_ignore.txt", 600)
B600 = InsarTimeseries.build_B_matrix(geolist, intlist600);
Blin600 = sum(B600, dims = 2);
unw_vals600_uplift = get_stack_vals(
    "unw_stack.h5",
    rowcol2...,
    1,
    "stack_flat_shifted",
    valid_igram_indices600,
    reference_station = nothing,
);
cc600_uplift = get_stack_vals(
    "cc_stack.h5",
    rowcol2...,
    1,
    "stack",
    valid_igram_indices600,
    reference_station = nothing,
);
unw_vals600_subs = get_stack_vals(
    "unw_stack.h5",
    rowcol3...,
    1,
    "stack_flat_shifted",
    valid_igram_indices600,
    reference_station = nothing,
);
cc600_subs = get_stack_vals(
    "cc_stack.h5",
    rowcol3...,
    1,
    "stack",
    valid_igram_indices600,
    reference_station = nothing,
);
unw_vals600_txoz = get_stack_vals(
    "unw_stack.h5",
    "TXOZ",
    1,
    "stack_flat_shifted",
    valid_igram_indices600,
    reference_station = nothing,
);
cc600_txoz = get_stack_vals(
    "cc_stack.h5",
    "TXOZ",
    1,
    "stack",
    valid_igram_indices600,
    reference_station = nothing,
);
unw_vals600_nmhb = get_stack_vals(
    "unw_stack.h5",
    "NMHB",
    1,
    "stack_flat_shifted",
    valid_igram_indices600,
    reference_station = nothing,
);
cc600_nmhb = get_stack_vals(
    "cc_stack.h5",
    "NMHB",
    1,
    "stack",
    valid_igram_indices600,
    reference_station = nothing,
);

function get_vals(days; station = nothing, rowcol = nothing, ignore = true)
    gi = ignore ? "geolist_ignore.txt" : nothing

    gl, il, vv = load_geolist_intlist("unw_stack.h5", gi, days)
    Bl = sum(InsarTimeseries.build_B_matrix(gl, il), dims = 2)
    p = isnothing(station) ? rowcol : [station]
    println(p)

    u = get_stack_vals("unw_stack.h5", p..., 1, "stack_flat_shifted", vv)
    c = get_stack_vals("cc_stack.h5", p..., 1, "stack", vv)
    return gl, il, vv, Bl, u, c
end

#### Points for looked at the diffs:
# Looked, path78 asc
# rowcol_asc = [226, 40]
rowcol_asc = [259, 82]  # Uplift looks stronger in stack than not
# # Looked, path85 desc
# rowcol_desc = [106, 400]
#
#
# rowcol_top = [25, 47]  # Uplift bowl in top left
rowcol_top = [14, 66]  # Uplift Edge of top left (random point)

unw_vals500_asc = get_stack_vals(
    "unw_stack.h5",
    rowcol_asc...,
    1,
    "stack_flat_shifted",
    valid_igram_indices500,
    reference_station = nothing,
);
cc500_asc = get_stack_vals(
    "cc_stack.h5",
    rowcol_asc...,
    1,
    "stack",
    valid_igram_indices500,
    reference_station = nothing,
);
unw_vals600_asc = get_stack_vals(
    "unw_stack.h5",
    rowcol_asc...,
    1,
    "stack_flat_shifted",
    valid_igram_indices600,
    reference_station = nothing,
);
# unw_vals500_desc = get_stack_vals("unw_stack.h5", rowcol_desc..., 1, "stack_flat_shifted", valid_igram_indices500, reference_station=nothing);

unw_vals500_top = get_stack_vals(
    "unw_stack.h5",
    rowcol_top...,
    1,
    "stack_flat_shifted",
    valid_igram_indices500,
);
cc500_top = get_stack_vals("cc_stack.h5", rowcol_top..., 1, "stack", valid_igram_indices500);

# g35, i35, v35 = load_geolist_intlist("unw_stack.h5", "geolist_ignore_35pct.txt", 500)
# Blin35 = sum(InsarTimeseries.build_B_matrix(g35, i35), dims=2)
# g50, i50, v50 = load_geolist_intlist("unw_stack.h5", "geolist_ignore_50pct.txt", 500)
# Blin50 = sum(InsarTimeseries.build_B_matrix(g50, i50), dims=2)
# unw_vals_up1 = get_stack_vals("unw_stack.h5", rowcol..., 1, "stack_flat_shifted", valid_igram_indices500);
# unw_vals_up35 = get_stack_vals("unw_stack.h5", rowcol..., 1, "stack_flat_shifted", v35);
# unw_vals_up50 = get_stack_vals("unw_stack.h5", rowcol..., 1, "stack_flat_shifted", v50);

####
# # Splitting by 2017 stuff
# geolist1, intlist1, valid_igram_indices1 = load_geolist_intlist("unw_stack.h5", "geolist_ignore.txt", 500, max_date=Date(2017,1,1))
# geolist2, intlist2, valid_igram_indices2 = load_geolist_intlist("unw_stack.h5", "geolist_ignore.txt", 500, min_date=Date(2017,1,1))
#
# B1 = InsarTimeseries.build_B_matrix(geolist1, intlist1)
# Blin1 = sum(B1, dims=2);
# B2 = InsarTimeseries.build_B_matrix(geolist2, intlist2)
# Blin2 = sum(B2, dims=2);
#
# unw_vals1_txoz = get_stack_vals("unw_stack.h5", "TXOZ", 1, "stack_flat_shifted", valid_igram_indices1, reference_station=nothing);
# cc1_txoz = get_stack_vals("cc_stack.h5", "TXOZ", 1, "stack", valid_igram_indices1, reference_station=nothing);
#
# unw_vals2_txoz = get_stack_vals("unw_stack.h5", "TXOZ", 1, "stack_flat_shifted", valid_igram_indices2, reference_station=nothing);
# cc2_txoz = get_stack_vals("cc_stack.h5", "TXOZ", 1, "stack", valid_igram_indices2, reference_station=nothing);
#
# # Good area: txmh
# unw_vals1_txmh = get_stack_vals("unw_stack.h5", "TXMH", 1, "stack_flat_shifted", valid_igram_indices1, reference_station=nothing);
# cc1_txmh = get_stack_vals("cc_stack.h5", "TXMH", 1, "stack", valid_igram_indices1, reference_station=nothing);
#
# unw_vals2_txmh = get_stack_vals("unw_stack.h5", "TXMH", 1, "stack_flat_shifted", valid_igram_indices2, reference_station=nothing);
# cc2_txmh = get_stack_vals("cc_stack.h5", "TXMH", 1, "stack", valid_igram_indices2, reference_station=nothing);
#
# # TXS3: worst for 2017 only
# unw_vals1_txs3 = get_stack_vals("unw_stack.h5", "TXS3", 1, "stack_flat_shifted", valid_igram_indices1, reference_station=nothing);
# cc1_txs3 = get_stack_vals("cc_stack.h5", "TXS3", 1, "stack", valid_igram_indices1, reference_station=nothing);
#
# unw_vals2_txs3 = get_stack_vals("unw_stack.h5", "TXS3", 1, "stack_flat_shifted", valid_igram_indices2, reference_station=nothing);
# cc2_txs3 = get_stack_vals("cc_stack.h5", "TXS3", 1, "stack", valid_igram_indices2, reference_station=nothing);

function plot_unregs(; sigma = 3, max_temp = 700, shrink = false, max_date = nothing)
    constant = false
    g, i, igram_idxs = load_geolist_intlist(
        "unw_stack.h5",
        "geolist_ignore.txt",
        max_temp,
        max_date = max_date,
    )
    rowcols = [
        [(414, 678), "noise? -8/yr for sigma=4+, 0 on sigma 3"],
        [(410, 667), "No signal, lots of stuff removed"],
        [(363, 594), "Bot right (real) bowl"],
        [(301, 301), "middle of map, vague red"],
        [(257, 81), "Bowl with eastward near pecos"],
        [(365, 128), "Deepest part of pecos wing"],
    ]
    plt.figure()
    for (rowcol, label) in rowcols
        row, col = rowcol
        unw_vals =
            get_stack_vals("unw_stack.h5", row, col, 1, "stack_flat_shifted", igram_idxs)

        geo, ints, unws = prune_igrams(g, i, unw_vals, sigma, nothing, nothing, false)
        ints, unws = shrink ? shrink_baseline(geo, ints, unws) : (ints, unws)
        B = InsarTimeseries.prepB(geo, ints, constant)
        soln =
            p2c .*
            InsarTimeseries.integrate_velocities(B \ unws, InsarTimeseries.day_diffs(geo))

        Blin = sum(B, dims = 2)
        soln_lin = (Blin\unws)[1] * p2mm
        println("$label: $soln_lin cm/yr")

        plt.plot(
            geo,
            soln,
            ".-",
            label = "($row, $col): $label, Linear: $(round(soln_lin)) mm/yr",
        )
    end
    plt.legend()
    plt.title("Sigma = $sigma, max baseline = $max_temp, shrink = $shrink")
end


text:
# Test that Richard's equation agrees with solution from Bonan's book,
# simulation 8.2
using MPI
using OrderedCollections
using StaticArrays
using Statistics
using Dierckx
using Test
using Pkg.Artifacts
using DelimitedFiles

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using ClimateMachine
using ClimateMachine.Land
using ClimateMachine.Land.SoilWaterParameterizations
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.DGMethods: BalanceLaw, LocalGeometry
using ClimateMachine.MPIStateArrays
using ClimateMachine.GenericCallbacks
using ClimateMachine.SystemSolvers
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates
using ClimateMachine.SingleStackUtils
using ClimateMachine.BalanceLaws:
    BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux, vars_state
using ClimateMachine.ArtifactWrappers

haverkamp_dataset = ArtifactWrapper(
    joinpath("test", "Land", "Model", "Artifacts.toml"),
    "richards",
    ArtifactFile[ArtifactFile(
        url = "https://caltech.box.com/shared/static/dfijf07io7h5dk1k87saaewgsg9apq8d.csv",
        filename = "bonan_haverkamp_data.csv",
    ),],
)
haverkamp_dataset_path = get_data_folder(haverkamp_dataset)

@testset "Richard's equation - Haverkamp test" begin
    ClimateMachine.init()
    FT = Float64

    function init_soil_water!(land, state, aux, coordinates, time)
        myfloat = eltype(aux)
        state.soil.water.ϑ_l = myfloat(land.soil.water.initialϑ_l(aux))
        state.soil.water.θ_i = myfloat(land.soil.water.initialθ_i(aux))
    end

    soil_heat_model = PrescribedTemperatureModel()

    soil_param_functions = SoilParamFunctions{FT}(
        porosity = 0.495,
        Ksat = 0.0443 / (3600 * 100),
        S_s = 1e-3,
    )
    # Keep in mind that what is passed is aux⁻
    # Fluxes are multiplied by ẑ (normal to the surface, -normal to the bottom,
    # where normal point outs of the domain.)
    surface_value = FT(0.494)
    bottom_flux_multiplier = FT(1.0)
    initial_moisture = FT(0.24)

    surface_state = (aux, t) -> surface_value
    bottom_flux = (aux, t) -> aux.soil.water.K * bottom_flux_multiplier
    ϑ_l0 = (aux) -> initial_moisture

    soil_water_model = SoilWaterModel(
        FT;
        moisture_factor = MoistureDependent{FT}(),
        hydraulics = Haverkamp{FT}(),
        initialϑ_l = ϑ_l0,
        dirichlet_bc = Dirichlet(
            surface_state = surface_state,
            bottom_state = nothing,
        ),
        neumann_bc = Neumann(surface_flux = nothing, bottom_flux = bottom_flux),
    )

    m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model)
    sources = ()
    m = LandModel(
        param_set,
        m_soil;
        source = sources,
        init_state_prognostic = init_soil_water!,
    )

    N_poly = 5
    nelem_vert = 10

    # Specify the domain boundaries
    zmax = FT(0)
    zmin = FT(-1)

    driver_config = ClimateMachine.SingleStackConfiguration(
        "LandModel",
        N_poly,
        nelem_vert,
        zmax,
        param_set,
        m;
        zmin = zmin,
        numerical_flux_first_order = CentralNumericalFluxFirstOrder(),
    )

    t0 = FT(0)
    timeend = FT(60 * 60 * 24)

    dt = FT(200)

    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_dt = dt,
    )

    #################### Change the ode_solver

    dg = solver_config.dg
    Q = solver_config.Q

    vdg = DGModel(
        driver_config.bl,
        driver_config.grid,
        driver_config.numerical_flux_first_order,
        driver_config.numerical_flux_second_order,
        driver_config.numerical_flux_gradient,
        state_auxiliary = dg.state_auxiliary,
        direction = VerticalDirection(),
    )

    linearsolver = BatchedGeneralizedMinimalResidual(
        dg,
        Q;
        max_subspace_size = 30,
        atol = -1.0,
        rtol = 1e-9,
    )

    """
    N(q)(Q) = Qhat  => F(Q) = N(q)(Q) - Qhat

    F(Q) == 0
    ||F(Q^i) || / ||F(Q^0) || < tol

    """
    nonlinearsolver =
        JacobianFreeNewtonKrylovSolver(Q, linearsolver; tol = 1e-9)

    ode_solver = ARK548L2SA2KennedyCarpenter(
        dg,
        vdg,
        NonLinearBackwardEulerSolver(
            nonlinearsolver;
            isadjustable = true,
            preconditioner_update_freq = 100,
        ),
        Q;
        dt = dt,
        t0 = 0,
        split_explicit_implicit = false,
        variant = NaiveVariant(),
    )

    solver_config.solver = ode_solver

    #######################################

    mygrid = solver_config.dg.grid
    aux = solver_config.dg.state_auxiliary
    ClimateMachine.invoke!(solver_config)
    ϑ_l_ind = varsindex(vars_state(m, Prognostic(), FT), :soil, :water, :ϑ_l)
    ϑ_l = Array(Q[:, ϑ_l_ind, :][:])
    z_ind = varsindex(vars_state(m, Auxiliary(), FT), :z)
    z = Array(aux[:, z_ind, :][:])

    # Compare with Bonan simulation data at 1 day.
    data = joinpath(haverkamp_dataset_path, "bonan_haverkamp_data.csv")
    ds_bonan = readdlm(data, ',')
    bonan_moisture = reverse(ds_bonan[:, 1])
    bonan_z = reverse(ds_bonan[:, 2]) ./ 100.0

    # Create an interpolation from the Bonan data
    bonan_moisture_continuous = Spline1D(bonan_z, bonan_moisture)
    bonan_at_clima_z = bonan_moisture_continuous.(z)
    MSE = mean((bonan_at_clima_z .- ϑ_l) .^ 2.0)
    @test MSE < 1e-5
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA