text:
mutable struct LossFunction{Δ, FT, ML, F, T, L, P}
          first_targets :: FT
  max_simulation_length :: ML
          field_weights :: F # scenario weights
            time_series :: T
                profile :: L
end

allsame(x) = all(y -> y ≈ first(x), x)
t_interval(data) = data.t[2:end] .- data.t[1:end-1]

function (loss::LossFunction)(simulation, observations, θ::Vector{<:FreeParameters})

    # iterate the model and record discrepancy summary in `time_series`
    evaluate!(loss, simulation, observations, θ)

    N_ens = ensemble_size(simulation.model)
    error = zeros((N_ens, 1))

    for ts in loss.time_series
        data_error = ts.analysis(ts.data) / N_ens
        error .+= data_error
    end

    return error
end

function LossFunction(simulation::Simulation{<:OneDimensionalEnsembleModel}, observations::OneDimensionalTimeSeriesBatch; data_weights=[1.0 for b in observations], relative_weights)

    @assert all([allsame(t_interval(data)) for data in observations]) "Simulation time steps are not uniformly spaced."
    @assert allsame([t_interval(data)[1] for data in observations]) "Time step differs between simulations."

    all_targets = getproperty.(observations, :targets)
    first_targets = getindex.(all_targets, 1)
    max_simulation_length = maximum(length.(all_targets))

    profile = ValueProfileAnalysis(simulation.model.grid, analysis = column_mean)

    field_weights = Dict(f => [] for f in [:u, :v, :b, :e])
    for (i, data) in enumerate(observations)
        data_fields = data.relevant_fields # e.g. (:b, :e)
        targets = all_targets[i]
        rw = [relative_weights[f] for f in data_fields]
        weights = estimate_weights(profile, data, rw) # e.g. (1.0, 0.5)

        for (j, field_name) in enumerate(data_fields)
            push!(field_weights[field_name], weights[j] * data_weights[i])
        end

        for field_name in keys(field_weights)
            field_name ∉ data_fields && push!(field_weights[field_name], 0)
        end
    end

    time_series = [EnsembleTimeSeriesAnalysis(observations[i].t[all_targets[i]], simulation.model.grid.Nx) for i in 1:length(observations)]

    return LossFunction(first_targets, max_simulation_length, field_weights, time_series, profile)
end

function calculate_value_discrepancy!(value, model_field, data_field)
    discrepancy = value.discrepancy

    interior(discrepancy) .= (interior(data_field) .- interior(model_field)) .^ 2

    return nothing
end

"""
    analyze_profile_discrepancy(value, model_field, data_field)

Calculates the discrepancy between model and data field values, and returns an
analysis of the discrepancy profile.
"""
function analyze_profile_discrepancy(value, model_field, data_field)
    calculate_value_discrepancy!(value, model_field, data_field) # MSE for each grid element
    return value.analysis(value.discrepancy) # e.g.column_mean of discrepancy field
end

function calculate_gradient_discrepancy!(prof, model_field, data_field)
    # Coarse grain the data
    ϵ = prof.ϵ
    set!(ϵ, data_field)

    # Calculate profients of both data and discrepancy
    ∇ϕ = prof.∇ϕ
    ∇ϵ = prof.∇ϵ
    ∂z!(∇ϵ, ϵ)
    ∂z!(∇ϕ, model_field)

    for i in eachindex(ϵ)
        @inbounds ϵ[i] = (ϵ[i] - model_field[i])^2
        @inbounds ∇ϵ[i] = (∇ϵ[i] - ∇ϕ[i])^2 # includes bottom boundary value, which will later be ignored.
    end

    # Top boundary contribution (ignored for now)
    #N = d.grid.N
    #@inbounds ∇d[N+1] = (∇d[N+1] - ∇ϕ[N+1])^2

    return nothing
end

"""
    analyze_profile_discrepancy(prof::GradientProfileAnalysis, model_field, data_field)

Calculates the discrepencies between both values and gradients of model and data fields,
and returns an analysis of the two discrepancy profiles.
"""
function analyze_profile_discrepancy(prof::GradientProfileAnalysis, model_field, data_field)
    calculate_gradient_discrepancy!(prof, model_field, data_field)

    # Calculate analysis on gradient, excluding boundary points.
    return prof.analysis(prof.ϵ) + prof.gradient_weight * prof.analysis(prof.∇ϵ.data[2:end-1])
end

#
# Loss function utils
#

@inline get_weight(::Nothing, field_index) = 1
@inline get_weight(weights, field_index) = @inbounds weights[field_index]

function new_field(field_name, field_data, grid)

    field_name == :u && return XFaceField(grid, field_data)
    field_name == :v && return YFaceField(grid, field_data)
    field_name == :b && return CenterField(grid, field_data)
    field_name == :e && return CenterField(grid, field_data)

end

function analyze_weighted_profile_discrepancy(loss::LossFunction, model::OneDimensionalEnsembleModel, observations::OneDimensionalTimeSeriesBatch, target)

    total_discrepancy = zeros(model.grid.Nx, model.grid.Ny, 1)

    for field_name in [:u, :v, :b, :e]

        model_field = get_model_field(model, field_name)

        # compensate for setting model time index 1 to to index `first_target` in data.
        data_indices = target .+ loss.first_targets .- 1

        field_data = column_ensemble_interior(observations, field_name, data_indices, model.grid.Nx)
        data_field = new_field(field_name, field_data, model_field.grid)

        # Calculate the per-field profile-based discrepancy
        field_discrepancy = analyze_profile_discrepancy(loss.profile, model_field, data_field)

        #=
        if any(isnan.(field_discrepancy))
            field_discrepency .= weight * remaining_time / stop_time
        end
        =#

        # Accumulate weighted profile-based discrepancies in the total discrepancyor
        total_discrepancy .+= loss.field_weights[field_name]' .* field_discrepancy # accumulate discrepancyor
    end

    return nan2inf.(total_discrepancy)
end

function evaluate!(loss::LossFunction, simulation::Simulation{<:OneDimensionalEnsembleModel}, observations::OneDimensionalTimeSeriesBatch, parameters)

    # Initialize
    initialize_forward_run!(simulation.model, observations, parameters, loss.first_targets)

    # this should be improved
    all_lengths = length.(getproperty.(observations, :t))
    longest_sim = observations[argmax(all_lengths)]

    # Calculate a loss function time-series
    for target in 1:loss.max_simulation_length

        simulation.stop_time = longest_sim.t[target]

        run!(simulation)

        discrepancy = analyze_weighted_profile_discrepancy(loss, simulation.model, observations, target)

        for (j, ts) in enumerate(loss.time_series)
            if target <= length(ts.time)
                # `ts.data` is N_ensemble x N_timesteps; `discrepancy` is N_ensemble x N_cases x 1
                ts.data[:, target] .= discrepancy[:, j, 1]
            end
        end

    end

    return nothing
end

#
# Miscellanea
#

function max_variance(data)
    fields = data.relevant_fields
    max_variances = zeros(length(fields))
    for (i, field) in enumerate(fields)
        max_variances[i] = get_weight(loss.weights, i) * max_variance(data, field)
    end
    return max_variances
end

function mean_variance(data)
    fields = data.relevant_fields
    mean_variance = zeros(length(fields))
    for (i, field) in enumerate(fields)
        mean_variance[i] = get_weight(loss.weights, i) * mean_variance(data, field)
    end
    return mean_variances
end


text:
# colors come from the ctioga2 project

goodColors = RGBA[RGBA(0.,0.,0.,1.0), # black
				RGBA(0.027,0.141,0.984,1.0), # blue
				RGBA(0.988,0.0352,0.106,1.0), # red
				RGBA(0.796,0.6118,0.988,1.0), # lilac
				RGBA(0.0745,0.4196,0.9843,1.0),  # bright blue
				RGBA(0.992,0.498,0.337,1.0), # coral
				RGBA(0.690,0.141,0.153,1.0)] # firebrick

ext: jl


text:
doc"""
`fzero(f,x1,x2,args...;tol=1e-12,maxsteps=100)`

finds the root $x_0$ of $f(x)$ in the interval $[x_1,x_2]$, such that
$f(x_0) = 0$ (up to specified tolerance `tol`).

`f` - a function f(x,args...) that returns Float64, where extra parameters
may be passed through to `f`.

throws an error if $f(x_1)$ and $f(x_2)$ have the same sign, or
if more than `maxsteps` iterations are taken.

uses the Dekker-Brent method (see `Numerical Recipes in C`, p. 361).
"""
@generated function fzero(f, x1::Float64, x2::Float64, args...; tol::Float64=1e-12, maxsteps::Int=100)
    quote
        local i,d,e,min1,min2,fc,p,q,r,s,tol1,xm
        a = x1
        b = x2
        c = x2
        fa = f(a,args...)::Float64
        fb = f(b,args...)::Float64
        δ = 3e-12
        ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) && error("f(a) and f(b) must have opposite signs")
        fc = fb
        for i in 1:maxsteps
            if (fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)
                c = a
                fc = fa
                e = d = b-a
            end
            if abs(fc) < abs(fb)
                a = b
                b = c
                c = a
                fa = fb
                fb = fc
                fc = fa
            end
            tol1 = 2*δ*abs(b) + 0.5*tol
            xm = 0.5*(c-b)
            if abs(xm) <= tol1 || fb == 0.0
                return b
            end
            if abs(e) >= tol1 && abs(fa) > abs(fb)
                s = fb/fa
                if a == c
                    p = 2.0*xm*s
                    q = 1.0-s
                else
                    q = fa/fc
                    r = fb/fc
                    p = s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0))
                    q = (q-1.0)*(r-1.0)*(s-1.0)
                end
                if p > 0.0
                    q = -q
                end
                p = abs(p)
                min1 = 3.0*xm*q - abs(tol1*q)
                min2 = abs(e*q)
                if 2.0*p < (min1 < min2 ? min1 : min2)
                    e = d
                    d = p/q
                else
                    d = xm
                    e = d
                end
            else
                d = xm
                e = d
            end
            a = b
            fa = fb
            if abs(d) > tol1
                b += d
            else
                b += abs(tol1)*sign(xm)
            end
            fb = f(b,args...)::Float64
        end
        error("maximum number of iterations exceeded")
    end
end

max_stars_repo_name: johncsnyder/Element


text:
import Base.convert

convert(::Type{Int}, s::String) = parse(Int, s)


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

using Test

using QuantumCircuits
using QuantumCircuits.QML
using QuantumCircuits.Execute

################################################################################
################################################################################
function test_findparam(n, expqc, params=nothing)
    expmat = tomatrix(expqc)

    qc = QCircuit(n)
    for i in 0:(n-1)
        qc.u3(i)
    end
    if !isnothing(params)
        setparameters!(qc, params)
    end
    params = getparameters(qc)

    @test expqc != qc
    @test expqc != simplify(qc)

    # fit the qc to expected
    findparam(expmat, qc)

    expqc != simplify(qc) && println(simplify(qc))
    expqc != simplify(qc) && println(params)
    @test expqc == simplify(qc)
end

################################################################################
n = 2
expqc = QCircuit(n)
expqc.x(0)

test_findparam(n, expqc)

# hard case

test_findparam(n, expqc, [3.266936951898874, 5.70062968255718, 4.413502891655148,
                    6.065806193650715, 4.586495424027681, 1.8946506980151154])

test_findparam(n, expqc, [3.146417882036646, 5.838659944632545, 0.7763087422764434,
                    3.890182287359087, 3.9144393001843962, 0.20934945580756525])

test_findparam(n, expqc, [6.027933728465509, 3.4939251107512956, 0.6586845450453815,
                    3.4399018540096193, 2.2051222189127144, 4.478774379878294])

################################################################################
n = 2
expqc = QCircuit(n)
expqc.x(0)
expqc.x(1)

test_findparam(n, expqc)

# hard case
test_findparam(n, expqc, [0.9904270500489234, 2.9438062464029695, 4.8575515480040545,
                    3.8392956982807447, 2.811189470141951, 5.58973178939339])

# hard case
test_findparam(n, expqc, [5.779528542008802, 4.8766335097778395, 3.6466973636363167,
                    2.9460221010727357, 2.924955576746307, 5.518914714482275])

################################################################################
n = 3
expqc = QCircuit(n)
expqc.z(0)
expqc.y(1)
expqc.x(2)

test_findparam(n, expqc)

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

n = 2
expqc = QCircuit(n)
expqc.h(0)
expqc.x(1)
expqc.cx(0, 1)
expqc.z(0)

expqc2 = QCircuit(n)
expqc2.h(0)
expqc2.cx(0, 1)
expqc2.z(0)
expqc2.x(1)

qc = QCircuit(n)
for i in 0:(n-1)
    qc.u3(i)
end
qc.cx(0, 1)
for i in 0:(n-1)
    qc.u3(i)
end

# test 1
expmat = tomatrix(expqc)
findparam(expmat, qc)
sqc = simplify(qc)
@test sqc == expqc || sqc == expqc2

# test 2
setparameters!(qc, getRandParameters(qc))
expmat = tomatrix(expqc2)
findparam(expmat, qc)
sqc = simplify(qc)
@test sqc == expqc || sqc == expqc2


text:
# class constructors
ZZ(::Type{T}) where T<:Integer = ZZ{T}

# construction
category_trait(::Type{<:ZZ}) = EuclidianDomainTrait
basetype(::Type{<:ZZ{T}}) where T = T
depth(::Type{<:ZZ}) = 1
_lcunit(a::ZZ) = sign(a.val)
issimpler(a::T, b::T) where T<:ZZ = abs(a.val) < abs(b.val)
iscoprime(a::T, b::T) where T<:ZZ = gcd(a.val, b.val) == 1
value(a::ZZ) = a.val

copy(a::ZZ) = typeof(a)(a.val)
ZZ(a::T) where T<:Integer = ZZ{T}(a)
ZZ{T}(a::ZZ) where T = ZZ{T}(T(a.val))

# promotion and conversion
_promote_rule(::Type{ZZ{T}}, ::Type{ZZ{S}}) where {S,T} = ZZ{promote_type(S,T)}
_promote_rule(::Type{ZZ{T}}, ::Type{QQ{S}}) where {S,T} = QQ{promote_type(S,T)}
promote_rule(::Type{ZZ{T}}, ::Type{S}) where {S<:Integer,T} = ZZ{promote_type(S,T)}
promote_rule(::Type{ZZ{T}}, ::Type{Rational{S}}) where {S<:Integer,T} = QQ{promote_type(S,T)}

# induced homomorphism
function (h::Hom{F,R,S})(p::Z) where {F,R,S,Z<:ZZ{<:R}}
    Z(h.f(p.val))
end

for op in (:+, :-, :/, :(==), :div, :rem, :divrem, :gcd, :gcdx, :pgcd, :pgcdx)
    @eval begin
        ($op)(a::ZZ{T}, b::Integer) where T = ($op)(promote(a, b)...)
        ($op)(a::Integer, b::ZZ{T}) where T = ($op)(promote(a, b)...)
    end
end

Base.isless(p::T, q::T) where T<:ZZ = isless(p.val, q.val)

# operations for ZZ
+(a::ZZ{T}, b::ZZ{T}) where T = ZZ(checked_add(a.val, b.val))
-(a::ZZ{T}, b::ZZ{T}) where T = ZZ(checked_sub(a.val, b.val))
-(a::ZZ{T}) where T = ZZ{T}(checked_neg(a.val))
*(a::ZZ{T}, b::ZZ{T}) where T = ZZ(checked_mul(a.val, b.val))
*(a::ZZ{T}, b::Integer) where T = ZZ(T(checked_mul(a.val, T(b))))
*(a::Integer, b::ZZ{T}) where T = ZZ(T(checked_mul(T(a), b.val)))
^(a::ZZ{BigInt}, b::Integer) = ZZ(a.val^b)
divrem(a::T, b::T) where T<:ZZ = T.(divrem(a.val, b.val))
div(a::T, b::T) where T<:ZZ = T(a.val ÷ b.val)
rem(a::T, b::T) where T<:ZZ = T(a.val % b.val)

isunit(a::ZZ) = abs(a.val) == 1
isone(a::ZZ) = isone(a.val)
iszero(a::ZZ) = iszero(a.val)
zero(::Type{ZZ{T}}) where T = ZZ(zero(T))
one(::Type{ZZ{T}}) where T = ZZ(one(T))
==(a::ZZ{T}, b::ZZ{T}) where T = a.val == b.val
hash(a::ZZ, h::UInt) = hash(a.val, h)

gcd(a::T, b::T) where T<: ZZ = T(gcd(a.val, b.val))
lcm(a::T, b::T) where T<:ZZ = T(lcm(a.val, b.val))

factor(a::Z) where Z <:ZZ = [Z(first(x)) => last(x) for x in Primes.factor(value(a))]
isirreducible(a::ZZ) = isirreducible(a.val)
isirreducible(a::Integer) = Primes.isprime(a)

Base.show(io::IO, z::ZZ) = print(io, z.val)


text:
Text provided under a Creative Commons Attribution license, CC-BY.  All code is made available under the FSF-approved BSD-3 license.  (c) Lorena A. Barba, Gilbert F. Forsyth 2017. Thanks to NSF for support via CAREER award #1149784.
[@LorenaABarba](https://twitter.com/LorenaABarba)

12 steps to Navier–Stokes
=====
***

Did you make it this far? This is the last step! How long did it take you to write your own Navier–Stokes solver in Python following this interactive module? Let us know!

Step 12: Channel Flow with Navier–Stokes
----
***

The only difference between this final step and Step 11 is that we are going to add a source term to the $u$-momentum equation, to mimic the effect of a pressure-driven channel flow. Here are our modified Navier–Stokes equations:

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=-\frac{1}{\rho}\frac{\partial p}{\partial x}+\nu\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right)+F$$

$$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}=-\frac{1}{\rho}\frac{\partial p}{\partial y}+\nu\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right)$$

$$\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2}=-\rho\left(\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial v}{\partial y}\right)
$$

### Discretized equations

With patience and care, we write the discretized form of the equations. It is highly recommended that you write these in your own hand, mentally following each term as you write it.

The $u$-momentum equation:

$$
\begin{split}
& \frac{u_{i,j}^{n+1}-u_{i,j}^{n}}{\Delta t}+u_{i,j}^{n}\frac{u_{i,j}^{n}-u_{i-1,j}^{n}}{\Delta x}+v_{i,j}^{n}\frac{u_{i,j}^{n}-u_{i,j-1}^{n}}{\Delta y} = \\
& \qquad -\frac{1}{\rho}\frac{p_{i+1,j}^{n}-p_{i-1,j}^{n}}{2\Delta x} \\
& \qquad +\nu\left(\frac{u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}}{\Delta x^2}+\frac{u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}}{\Delta y^2}\right)+F_{i,j}
\end{split}
$$

The $v$-momentum equation:

$$
\begin{split}
& \frac{v_{i,j}^{n+1}-v_{i,j}^{n}}{\Delta t}+u_{i,j}^{n}\frac{v_{i,j}^{n}-v_{i-1,j}^{n}}{\Delta x}+v_{i,j}^{n}\frac{v_{i,j}^{n}-v_{i,j-1}^{n}}{\Delta y} = \\
& \qquad -\frac{1}{\rho}\frac{p_{i,j+1}^{n}-p_{i,j-1}^{n}}{2\Delta y} \\
& \qquad +\nu\left(\frac{v_{i+1,j}^{n}-2v_{i,j}^{n}+v_{i-1,j}^{n}}{\Delta x^2}+\frac{v_{i,j+1}^{n}-2v_{i,j}^{n}+v_{i,j-1}^{n}}{\Delta y^2}\right)
\end{split}
$$

And the pressure equation:

$$
\begin{split}
& \frac{p_{i+1,j}^{n}-2p_{i,j}^{n}+p_{i-1,j}^{n}}{\Delta x^2} + \frac{p_{i,j+1}^{n}-2p_{i,j}^{n}+p_{i,j-1}^{n}}{\Delta y^2} = \\
& \qquad \rho\left[\frac{1}{\Delta t}\left(\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}+\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right) - \frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x} - 2\frac{u_{i,j+1}-u_{i,j-1}}{2\Delta y}\frac{v_{i+1,j}-v_{i-1,j}}{2\Delta x} - \frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right]
\end{split}
$$

As always, we need to re-arrange these equations to the form we need in the code to make the iterations proceed.

For the $u$- and $v$ momentum equations, we isolate the velocity at time step `n+1`:

$$
\begin{split}
u_{i,j}^{n+1} = u_{i,j}^{n} & - u_{i,j}^{n} \frac{\Delta t}{\Delta x} \left(u_{i,j}^{n}-u_{i-1,j}^{n}\right) - v_{i,j}^{n} \frac{\Delta t}{\Delta y} \left(u_{i,j}^{n}-u_{i,j-1}^{n}\right) \\
& - \frac{\Delta t}{\rho 2\Delta x} \left(p_{i+1,j}^{n}-p_{i-1,j}^{n}\right) \\
& + \nu\left[\frac{\Delta t}{\Delta x^2} \left(u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}\right) + \frac{\Delta t}{\Delta y^2} \left(u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}\right)\right] \\
& + \Delta t F
\end{split}
$$

$$
\begin{split}
v_{i,j}^{n+1} = v_{i,j}^{n} & - u_{i,j}^{n} \frac{\Delta t}{\Delta x} \left(v_{i,j}^{n}-v_{i-1,j}^{n}\right) - v_{i,j}^{n} \frac{\Delta t}{\Delta y} \left(v_{i,j}^{n}-v_{i,j-1}^{n}\right) \\
& - \frac{\Delta t}{\rho 2\Delta y} \left(p_{i,j+1}^{n}-p_{i,j-1}^{n}\right) \\
& + \nu\left[\frac{\Delta t}{\Delta x^2} \left(v_{i+1,j}^{n}-2v_{i,j}^{n}+v_{i-1,j}^{n}\right) + \frac{\Delta t}{\Delta y^2} \left(v_{i,j+1}^{n}-2v_{i,j}^{n}+v_{i,j-1}^{n}\right)\right]
\end{split}
$$

And for the pressure equation, we isolate the term $p_{i,j}^n$ to iterate in pseudo-time:

$$
\begin{split}
p_{i,j}^{n} = & \frac{\left(p_{i+1,j}^{n}+p_{i-1,j}^{n}\right) \Delta y^2 + \left(p_{i,j+1}^{n}+p_{i,j-1}^{n}\right) \Delta x^2}{2(\Delta x^2+\Delta y^2)} \\
& -\frac{\rho\Delta x^2\Delta y^2}{2\left(\Delta x^2+\Delta y^2\right)} \\
& \times \left[\frac{1}{\Delta t} \left(\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x} + \frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right) - \frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x} - 2\frac{u_{i,j+1}-u_{i,j-1}}{2\Delta y}\frac{v_{i+1,j}-v_{i-1,j}}{2\Delta x} - \frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right]
\end{split}
$$

The initial condition is $u, v, p=0$ everywhere, and at the boundary conditions are:

$u, v, p$ are periodic on $x=0,2$

$u, v =0$ at $y =0,2$

$\frac{\partial p}{\partial y}=0$ at $y =0,2$

$F=1$ everywhere.

Let's begin by importing our usual run of libraries:

```python
import numpy
from matplotlib import pyplot, cm
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
```

In step 11, we isolated a portion of our transposed equation to make it easier to parse and we're going to do the same thing here.  One thing to note is that we have periodic boundary conditions throughout this grid, so we need to explicitly calculate the values at the leading and trailing edge of our `u` vector.

```python
def build_up_b(rho, dt, dx, dy, u, v):
    b = numpy.zeros_like(u)
    b[1:-1, 1:-1] = (rho * (1 / dt * ((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx) +
                                      (v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy)) -
                            ((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx))**2 -
                            2 * ((u[2:, 1:-1] - u[0:-2, 1:-1]) / (2 * dy) *
                                 (v[1:-1, 2:] - v[1:-1, 0:-2]) / (2 * dx))-
                            ((v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy))**2))

    # Periodic BC Pressure @ x = 2
    b[1:-1, -1] = (rho * (1 / dt * ((u[1:-1, 0] - u[1:-1,-2]) / (2 * dx) +
                                    (v[2:, -1] - v[0:-2, -1]) / (2 * dy)) -
                          ((u[1:-1, 0] - u[1:-1, -2]) / (2 * dx))**2 -
                          2 * ((u[2:, -1] - u[0:-2, -1]) / (2 * dy) *
                               (v[1:-1, 0] - v[1:-1, -2]) / (2 * dx)) -
                          ((v[2:, -1] - v[0:-2, -1]) / (2 * dy))**2))

    # Periodic BC Pressure @ x = 0
    b[1:-1, 0] = (rho * (1 / dt * ((u[1:-1, 1] - u[1:-1, -1]) / (2 * dx) +
                                   (v[2:, 0] - v[0:-2, 0]) / (2 * dy)) -
                         ((u[1:-1, 1] - u[1:-1, -1]) / (2 * dx))**2 -
                         2 * ((u[2:, 0] - u[0:-2, 0]) / (2 * dy) *
                              (v[1:-1, 1] - v[1:-1, -1]) / (2 * dx))-
                         ((v[2:, 0] - v[0:-2, 0]) / (2 * dy))**2))

    return b
```

We'll also define a Pressure Poisson iterative function, again like we did in Step 11.  Once more, note that we have to include the periodic boundary conditions at the leading and trailing edge.  We also have to specify the boundary conditions at the top and bottom of our grid.

```python
def pressure_poisson_periodic(p, dx, dy):
    pn = numpy.empty_like(p)

    for q in range(nit):
        pn = p.copy()
        p[1:-1, 1:-1] = (((pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy**2 +
                          (pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx**2) /
                         (2 * (dx**2 + dy**2)) -
                         dx**2 * dy**2 / (2 * (dx**2 + dy**2)) * b[1:-1, 1:-1])

        # Periodic BC Pressure @ x = 2
        p[1:-1, -1] = (((pn[1:-1, 0] + pn[1:-1, -2])* dy**2 +
                        (pn[2:, -1] + pn[0:-2, -1]) * dx**2) /
                       (2 * (dx**2 + dy**2)) -
                       dx**2 * dy**2 / (2 * (dx**2 + dy**2)) * b[1:-1, -1])

        # Periodic BC Pressure @ x = 0
        p[1:-1, 0] = (((pn[1:-1, 1] + pn[1:-1, -1])* dy**2 +
                       (pn[2:, 0] + pn[0:-2, 0]) * dx**2) /
                      (2 * (dx**2 + dy**2)) -
                      dx**2 * dy**2 / (2 * (dx**2 + dy**2)) * b[1:-1, 0])

        # Wall boundary conditions, pressure
        p[-1, :] =p[-2, :]  # dp/dy = 0 at y = 2
        p[0, :] = p[1, :]  # dp/dy = 0 at y = 0

    return p
```

Now we have our familiar list of variables and initial conditions to declare before we start.

```python
##variable declarations
nx = 41
ny = 41
nt = 10
nit = 50
c = 1
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
x = numpy.linspace(0, 2, nx)
y = numpy.linspace(0, 2, ny)
X, Y = numpy.meshgrid(x, y)

##physical variables
rho = 1
nu = .1
F = 1
dt = .01

#initial conditions
u = numpy.zeros((ny, nx))
un = numpy.zeros((ny, nx))

v = numpy.zeros((ny, nx))
vn = numpy.zeros((ny, nx))

p = numpy.ones((ny, nx))
pn = numpy.ones((ny, nx))

b = numpy.zeros((ny, nx))
```

For the meat of our computation, we're going to reach back to a trick we used in Step 9 for Laplace's Equation.  We're interested in what our grid will look like once we've reached a near-steady state.  We can either specify a number of timesteps `nt` and increment it until we're satisfied with the results, or we can tell our code to run until the difference between two consecutive iterations is very small.

We also have to manage **8** separate boundary conditions for each iteration.  The code below writes each of them out explicitly.  If you're interested in a challenge, you can try to write a function which can handle some or all of these boundary conditions.  If you're interested in tackling that, you should probably read up on Python [dictionaries](http://docs.python.org/2/tutorial/datastructures.html#dictionaries).

```python
udiff = 1
stepcount = 0

while udiff > .001:
    un = u.copy()
    vn = v.copy()

    b = build_up_b(rho, dt, dx, dy, u, v)
    p = pressure_poisson_periodic(p, dx, dy)

    u[1:-1, 1:-1] = (un[1:-1, 1:-1] -
                     un[1:-1, 1:-1] * dt / dx *
                    (un[1:-1, 1:-1] - un[1:-1, 0:-2]) -
                     vn[1:-1, 1:-1] * dt / dy *
                    (un[1:-1, 1:-1] - un[0:-2, 1:-1]) -
                     dt / (2 * rho * dx) *
                    (p[1:-1, 2:] - p[1:-1, 0:-2]) +
                     nu * (dt / dx**2 *
                    (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +
                     dt / dy**2 *
                    (un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1])) +
                     F * dt)

    v[1:-1, 1:-1] = (vn[1:-1, 1:-1] -
                     un[1:-1, 1:-1] * dt / dx *
                    (vn[1:-1, 1:-1] - vn[1:-1, 0:-2]) -
                     vn[1:-1, 1:-1] * dt / dy *
                    (vn[1:-1, 1:-1] - vn[0:-2, 1:-1]) -
                     dt / (2 * rho * dy) *
                    (p[2:, 1:-1] - p[0:-2, 1:-1]) +
                     nu * (dt / dx**2 *
                    (vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, 0:-2]) +
                     dt / dy**2 *
                    (vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[0:-2, 1:-1])))

    # Periodic BC u @ x = 2
    u[1:-1, -1] = (un[1:-1, -1] - un[1:-1, -1] * dt / dx *
                  (un[1:-1, -1] - un[1:-1, -2]) -
                   vn[1:-1, -1] * dt / dy *
                  (un[1:-1, -1] - un[0:-2, -1]) -
                   dt / (2 * rho * dx) *
                  (p[1:-1, 0] - p[1:-1, -2]) +
                   nu * (dt / dx**2 *
                  (un[1:-1, 0] - 2 * un[1:-1,-1] + un[1:-1, -2]) +
                   dt / dy**2 *
                  (un[2:, -1] - 2 * un[1:-1, -1] + un[0:-2, -1])) + F * dt)

    # Periodic BC u @ x = 0
    u[1:-1, 0] = (un[1:-1, 0] - un[1:-1, 0] * dt / dx *
                 (un[1:-1, 0] - un[1:-1, -1]) -
                  vn[1:-1, 0] * dt / dy *
                 (un[1:-1, 0] - un[0:-2, 0]) -
                  dt / (2 * rho * dx) *
                 (p[1:-1, 1] - p[1:-1, -1]) +
                  nu * (dt / dx**2 *
                 (un[1:-1, 1] - 2 * un[1:-1, 0] + un[1:-1, -1]) +
                  dt / dy**2 *
                 (un[2:, 0] - 2 * un[1:-1, 0] + un[0:-2, 0])) + F * dt)

    # Periodic BC v @ x = 2
    v[1:-1, -1] = (vn[1:-1, -1] - un[1:-1, -1] * dt / dx *
                  (vn[1:-1, -1] - vn[1:-1, -2]) -
                   vn[1:-1, -1] * dt / dy *
                  (vn[1:-1, -1] - vn[0:-2, -1]) -
                   dt / (2 * rho * dy) *
                  (p[2:, -1] - p[0:-2, -1]) +
                   nu * (dt / dx**2 *
                  (vn[1:-1, 0] - 2 * vn[1:-1, -1] + vn[1:-1, -2]) +
                   dt / dy**2 *
                  (vn[2:, -1] - 2 * vn[1:-1, -1] + vn[0:-2, -1])))

    # Periodic BC v @ x = 0
    v[1:-1, 0] = (vn[1:-1, 0] - un[1:-1, 0] * dt / dx *
                 (vn[1:-1, 0] - vn[1:-1, -1]) -
                  vn[1:-1, 0] * dt / dy *
                 (vn[1:-1, 0] - vn[0:-2, 0]) -
                  dt / (2 * rho * dy) *
                 (p[2:, 0] - p[0:-2, 0]) +
                  nu * (dt / dx**2 *
                 (vn[1:-1, 1] - 2 * vn[1:-1, 0] + vn[1:-1, -1]) +
                  dt / dy**2 *
                 (vn[2:, 0] - 2 * vn[1:-1, 0] + vn[0:-2, 0])))

    # Wall BC: u,v = 0 @ y = 0,2
    u[0, :] = 0
    u[-1, :] = 0
    v[0, :] = 0
    v[-1, :]=0

    udiff = (numpy.sum(u) - numpy.sum(un)) / numpy.sum(u)
    stepcount += 1
```

You can see that we've also included a variable `stepcount` to see how many iterations our loop went through before our stop condition was met.

```python
print(stepcount)
```

    499

If you want to see how the number of iterations increases as our `udiff` condition gets smaller and smaller, try defining a function to perform the `while` loop written above that takes an input `udiff` and outputs the number of iterations that the function runs.

For now, let's look at our results.  We've used the quiver function to look at the cavity flow results and it works well for channel flow, too.

```python
fig = pyplot.figure(figsize = (11,7), dpi=100)
pyplot.quiver(X[::3, ::3], Y[::3, ::3], u[::3, ::3], v[::3, ::3]);
```

![png](output_23_0.png)

The structures in the `quiver` command that look like `[::3, ::3]` are useful when dealing with large amounts of data that you want to visualize.  The one used above tells `matplotlib` to only plot every 3rd data point.  If we leave it out, you can see that the results can appear a little crowded.

```python
fig = pyplot.figure(figsize = (11,7), dpi=100)
pyplot.quiver(X, Y, u, v);
```

![png](output_25_0.png)

## Learn more
***

##### What is the meaning of the $F$ term?

Step 12 is an exercise demonstrating the problem of flow in a channel or pipe. If you recall from your fluid mechanics class, a specified pressure gradient is what drives Poisseulle flow.

Recall the $x$-momentum equation:

$$\frac{\partial u}{\partial t}+u \cdot \nabla u = -\frac{\partial p}{\partial x}+\nu \nabla^2 u$$

What we actually do in Step 12 is split the pressure into steady and unsteady components $p=P+p'$. The applied steady pressure gradient is the constant $-\frac{\partial P}{\partial x}=F$ (interpreted as a source term), and the unsteady component is $\frac{\partial p'}{\partial x}$. So the pressure that we solve for in Step 12 is actually $p'$, which for a steady flow is in fact equal to zero everywhere.

Why did we do this?

Note that we use periodic boundary conditions for this flow. For a flow with a constant pressure gradient, the value of pressure on the left edge of the domain must be different from the pressure at the right edge. So we cannot apply periodic boundary conditions on the pressure directly. It is easier to fix the gradient and then solve for the perturbations in pressure.

Shouldn't we always expect a uniform/constant $p'$ then?

That's true only in the case of steady laminar flows. At high Reynolds numbers, flows in channels can become turbulent, and we will see unsteady fluctuations in the pressure, which will result in non-zero values for $p'$.

In step 12, note that the pressure field itself is not constant, but it's the pressure perturbation field that is. The pressure field varies linearly along the channel with slope equal to the pressure gradient. Also, for incompressible flows, the absolute value of the pressure is inconsequential.

##### And explore more CFD materials online

The interactive module **12 steps to Navier–Stokes** is one of several components of the Computational Fluid Dynamics class taught by Prof. Lorena A. Barba in Boston University between 2009 and 2013.

For a sample of what the othe components of this class are, you can explore the **Resources** section of the Spring 2013 version of [the course's Piazza site](https://piazza.com/bu/spring2013/me702/resources).

***

```python
from IPython.core.display import HTML
def css_styling():
    styles = open("../styles/custom.css", "r").read()
    return HTML(styles)
css_styling()
```

<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>

(The cell above executes the style for this notebook.)

max_stars_repo_name: mkitti/CFDJulia

max_issues_repo_name: mkitti/CFDJulia

max_forks_repo_name: mkitti/CFDJulia


text:
############################################################################################
#Import
include("diagnostics.jl")
include("tune.jl")

############################################################################################
#export


text:
# Stata Files

## ----------  File reading

function stata_dir(ctry :: Country)
    return joinpath(baseDir, "MAFE_data/$(ctry.name)/Stata");
end

function stata_file_name(ctry:: Country, area :: String)
    return "$(area).dta"
end

function stata_file_path(ctry :: Country, area :: String)
    return joinpath(stata_dir(ctry), stata_file_name(ctry, area))
end

function read_stata_file(ctry :: Country, area :: String)
    return DataFrame(load(stata_file_path(ctry, area)))
end

"""
    $(SIGNATURES)

Read person (general) file and rename variables.
"""
function read_person_file(ctry :: Country)
    # This is one row per person
    df = read_stata_file(ctry, "general");
    rename!(df, :q1 => :sex, :q1a => :birthYear);
    println("No of persons: $(length(df.ident))");

    # Filter
    validV = (df.birthYear .>= birthYearMin) .& (df.birthYear .<= birthYearMax);
    println("N after age filter: $(sum(validV .== true))")
    vIdxV = findall(validV);
    println("Final N: $(length(vIdxV))");

    return df[vIdxV, [:ident, :sex, :birthYear]]
end

"""
	$(SIGNATURES)

Given a file with spells, clean up and sort.
"""
function clean_spells(df :: DataFrame)
    # Ongoing spells. Use future year as end year
    df.endYear = coalesce.(df.endYear,  2020);

    dropmissing!(df, [:startYear, :endYear],  disallowmissing = true);

    # Types of years are sometimes strange in the Stata files
    df.startYear = convert(Vector{Int}, df.startYear);
    df.endYear = convert(Vector{Int}, df.endYear);

    # Fix a typo for Congo
    df.endYear[df.endYear .== 1695] .= 1965;

    validV = (df.startYear .> 1900)  .&  (df.endYear .< 2050);
    dfOut = df[findall(validV .== true),  :];

    # Sort by start year of spell (for each ident)
    sort!(dfOut, [:ident, :startYear, :endYear]);
    return dfOut
end

"""
	$(SIGNATURES)

Validate start and end years in a DataFrame
"""
function validate_start_end_years(df :: DataFrame)
    isValid = true;
    if !all(df.endYear .>= df.startYear)
        @warn "End years before start years"
        isValid = false;
    end

    nRows = size(df, 1);
    for i1 = 2 : nRows
        currId = df.ident[i1];
        if currId == df.ident[i1-1]
            if df.startYear[i1] < df.endYear[i1-1]
                @warn "$currId: Start year before previous end year"
                isValid = false;
            end
            if df.startYear[i1] < df.startYear[i1-1]
                @warn "$currId: Start years not increasing"
                isValid = false;
            end
            if df.endYear[i1] < df.endYear[i1-1]
                @warn "$currId: End years not increasing"
                isValid = false;
            end
        end
    end

    return isValid
end

"""
	$(SIGNATURES)

Read activity file. Rename variables.
Replace missing end years with interview year.
Missing startYear should indicate a person without any spells.
Activity spells are not consecutive.
"""
function read_activity_file(ctry; skipValidation :: Bool = false)
    df = read_stata_file(ctry, "activity");
    rename!(df, :q401d => :startYear,
        :q401f => :endYear,
        :q1a => :birthYear,
        :q408 => :income,
        :q408_dev => :currency);

    dfOut = clean_spells(df);

    dropmissing!(dfOut, [:income, :currency],  disallowmissing = true);

    validV = (dfOut.income .> 0.0);
    dfOut = dfOut[findall(validV .== true),  :];
    nRows = size(dfOut, 1);
    println("No of valid rows: $nRows")

    if skipValidation
        @assert validate_start_end_years(dfOut)  "Invalid activity file"
    end

    return dfOut
end

"""
	$(SIGNATURES)

Read migration file. Rename variables.
Replace missing end years with interview year.
"""
function read_migration_file(ctry; skipValidation :: Bool = false)
    df = read_stata_file(ctry, "migration");
    rename!(df, :q601d => :startYear,
        :q601f => :endYear,
        :q1a => :birthYear,
        :q602 => :countryCode,
        :q606 => :kindOfStay);

    dropmissing!(df, [:kindOfStay, :countryCode],  disallowmissing = true);

    dfOut = clean_spells(df);

    # Delete <= 1 year spells where kindOfStay != 3 (vactations, transit)
    # This is a little risky, but needed to avoid stays within longer stays
    dfOut.duration = dfOut.endYear .- dfOut.startYear;
    # validV = findall((dfOut.duration .> 1) .| (dfOut.kindOfStay .= 3));
    dropV = (dfOut.duration .< 2) .& (dfOut.kindOfStay .!= 3);

    # This person claims to stay in two places at the same time (invalid history)
    dropV[dfOut.ident .== "N013000"] .= true;

    deleterows!(dfOut, findall(dropV .== true));

    if skipValidation
        @assert validate_start_end_years(dfOut)  "Invalid migration file"
    end

    return dfOut
end

## -----------  Make job histories

"""
    $(SIGNATURES)

Make job histories for all persons. Also collects fixed info, such as schooling.

test this +++++
"""
function make_job_histories(ctry :: Country)
    println("Making job histories for $(ctry.name)");

    dfInd = read_person_file(ctry);
    dfAct = read_activity_file(ctry);
    dfMig = read_migration_file(ctry);

    nInd = length(dfInd.ident);
    println("Number of persons: $nInd")
    phV = Vector{PersonHistory}(undef, nInd);
    for i1 = 1 : nInd
        p = make_person_entry(ctry, i1, dfInd);
        # Make a country history for this person
        mh = make_migration_history(dfInd.ident[i1], dfMig);
        jh = make_job_history(dfInd.ident[i1], dfAct);
        # Merge migration info into job history
        merge_countries_to_jobs!(jh, mh);
        phV[i1] = PersonHistory(p, jh.spellV);
    end

    return phV
end

"""
	$(SIGNATURES)

Make entry for one person from individual (general) file
"""
function make_person_entry(ctry :: Country,
    i1 :: Int, dfInd :: DataFrame)

    return Person(dfInd.ident[i1],
        make_gender(dfInd.sex[i1]),
        dfInd.birthYear[i1],
        ctry)
end

"""
	$(SIGNATURES)

Make migration history for one person. Person may not be in the migration file.
Then return empty migration history.
"""
function make_migration_history(ident :: String, dfMig :: DataFrame)
    iMigV = findall(dfMig.ident .== ident);
    nMig = length(iMigV);
    migV = Vector{Migration}(undef, nMig);
    if !isempty(iMigV)
        for i1 = 1 : nMig
            ir = iMigV[i1];
            migV[i1] = Migration(dfMig.startYear[ir], dfMig.endYear[ir], Int(dfMig.countryCode[ir]));
        end
    end
    mh = MigrationHistory(migV);
    return mh
end

"""
	$(SIGNATURES)

Make one person's job history.
Contains each job with its country. Country is still unknown.
All persons have at least one spell in the original files, but not in the filtered ones returned by `read_activity_file`.
"""
function make_job_history(ident :: String, dfAct :: DataFrame)
    # Rows in activity file
    iActV = findall(dfAct.ident .== ident);

    nSpells = length(iActV);
    jobV = Vector{Job}(undef, nSpells);
    if nSpells >= 1
        for j = 1 : nSpells
            iRow = iActV[j];
            jobV[j] = Job(dfAct.startYear[iRow], dfAct.endYear[iRow],
                -1, dfAct.income[iRow]);

            @assert jobV[j].endYear >= jobV[j].startYear
        end
    end
    return JobHistory(jobV)
end

"""
	$(SIGNATURES)

Find country where a person was in a given year.
Returns MAFE country code for foreign country or 0 for home country.
Return negative for years that span spells.
"""
function find_country(mh :: MigrationHistory, startYear :: Int, endYear :: Int)
    nMig = n_spells(mh);
    if nMig == 0
        return homeCountry
    elseif endYear <= start_year(mh, 1)
        return homeCountry
    elseif startYear >= end_year(mh, nMig)
        return homeCountry
    else
        idx = spell_from_year_range(mh, startYear, endYear);
        if idx < 0
            return unknownCountry;
        else
            return mh.spellV[idx].ctry
        end
    end
end

"""
	$(SIGNATURES)

Merge countries from migration histories into job histories

test this +++++
"""
function merge_countries_to_jobs!(jh :: JobHistory, mh :: MigrationHistory)
    nJobs = n_spells(jh);
    if nJobs > 0
        for iJob = 1 : nJobs
            job = get_spell(jh, iJob);
            job.ctry = find_country(mh, job.startYear, job.endYear);
        end
    end
    return nothing
end

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


text:
# Quelques fonctions ou macros utiles en mode interactif
# (dont des alias courts pour les types ou les méthodes)

# Raccourci pour la gestion des options (Args)
as(sym::Symbol, val) = Args.set(sym, val)
ag(sym::Symbol, val) = Args.set(sym, val)
ag(sym::Symbol)      = Args.get(sym)
ag()                 = Args.show_args()

# Raccourci pour la gestion du level (pour les messages d'erreur)
lv(L) = Args.set(:level, L)
lv() = Args.get(:level)

p = println
pn = print
# inc=include

data = dirname(dirname(@__FILE__)) * "/data"

# Les chemins d'instance (String)
p01="$data/alp_01_p10.alp"
p02="$data/alp_02_p15.alp"
p03="$data/alp_03_p20.alp"
p04="$data/alp_04_p20.alp"
p05="$data/alp_05_p20.alp"
p06="$data/alp_06_p30.alp"
p07="$data/alp_07_p44.alp"
p08="$data/alp_08_p50.alp"
p09="$data/alp_09_p100.alp"
p10="$data/alp_10_p150.alp"
p11="$data/alp_11_p200.alp"
p12="$data/alp_12_p250.alp"
p13="$data/alp_13_p500.alp"

i01=Instance(p01)
i02=Instance(p02)
i03=Instance(p03)
i04=Instance(p04)
i05=Instance(p05)
i06=Instance(p06)
i07=Instance(p07)
i08=Instance(p08)
i09=Instance(p09)
i10=Instance(p10)
i11=Instance(p11)
i12=Instance(p12)
i13=Instance(p13)

"""
    macro @i(pats...)

Pour simplifier les "include" en interactif

Exemple d'itilisation :
   @i "04" "05"
Recherche tous les fichiers correspondant à la gpatternes *04* et *05*
puis les recharges par include (dans l'ordre)
Recherche d'abord dans le répertoire racine, puis dans src/ puis dans test/

ASTUCE :
Si les arguments sont purement alphanumériques, alors on peut éviter des
guillemets
  @i inst       ok
  => charge tous les fichiers contenant "inst" dans src/ ou dans test/
  @i  src/inst   KO
  @i "src/inst"  ok
  @i  04         KO (car est équivalent à @i "4" qui couvre plus large)
  @i "04"        ok
"""
macro i(pats...)
    files = get_files_from_pats(pats...)

    if length(files) == 0
        println("Aucun fichier ne correspond.")
    else
        for absfile in files
            relfiles = replace(absfile, "$APPDIR/" => "")
            println("### include file: \"", relfiles, "\" ...")
            include(absfile)
            println("### include file: \"", relfiles, "\" FAIT")
            # println("include done.")
        end
    end
end

"""
    macro f (find)

Affiche les fichiers correspondants aux patternes passées (avec date de modif)
"""
macro f(pats...)
    files = get_files_from_pats(pats...)
    wd = pwd(); cd(APPDIR)
    # On rend le chemin des fichiers relatifs par rapport au projet pour
    # alléger le listing (quelque soit le répertoire courant)
    files = replace.(files, "$APPDIR/" => "")
    for file in files
        run(`ls -l $file`)
    end
    cd(wd)
end

function get_files_from_pats(pats...)
    # @show typeof(pats)
    # @show pats
    files = []
    for pat in pats
        append!(files, glob("$pat", "$APPDIR") )
        append!(files, glob("*$pat*", "$APPDIR") )
        append!(files, glob("*$pat*", "$APPDIR/src") )
        append!(files, glob("*$pat*", "$APPDIR/test") )
    end
    return files
end
#./


text:
# Since we'll be working with economic varible which are all positive
# We'll modify this algorithm to limit search for root only in the postive quadrant.
using Calculus

function NewtonRoot(g, init_guess)

    xn1 = init_guess;
    tol = 10^(-8);
    distance = 1;
    iter = 1;

    while distance >= tol
        xn = xn1;
        xn1 = xn - g(xn) / derivative(g, xn);
        xn1 = max(xn1, 0.0001);

        distance = abs(xn1 - xn);
        iter = iter + 1;
    end

    return xn1
end

## Checks
# f1(x) = log(x)
# NewtonRoot(f1, 10)


text:
function try_sat_pure(model,V0,f!,T,result,error_val,method = LineSearch(Newton()))
    if !isfinite(V0[1]) | !isfinite(V0[2])
        return false
    end
    try
        res = sat_pure(model,V0,f!,T,method)
        result[] = res
    catch e #normally, failures occur near the critical point
        error_val[] = e
        return false
    end

    (P_sat,V_l,V_v) = result[]
    _,dpdvl = p∂p∂V(model,V_l,T,SA[1.0])
    _,dpdvv = p∂p∂V(model,V_v,T,SA[1.0])
    (dpdvl > 0) | (dpdvv > 0) && return false
    ε = abs(V_l-V_v)/(eps(typeof(V_l-V_v)))
    ups = (eps(typeof(V_l-V_v)))
    #@show abs(V_l-V_v),ups
    if (ε > 5e9)
        #if ΔV > ε then Vl and Vv are different values
        return true
    end

    return false
end

"""
    saturation_pressure(model::EoSModel, T, V0 = x0_sat_pure(model,T))

Performs a single component saturation equilibrium calculation, at the specified temperature `T`, of one mol of pure sustance specified by `model`

Returns `(p₀, Vₗ, Vᵥ)` where `p₀` is the saturation pressure (in Pa), `Vₗ` is the liquid saturation volume (in m³) and `Vᵥ` is the vapour saturation volume (in m³).

If the calculation fails, returns  `(NaN, NaN, NaN)`

`V0` is `[log10(Vₗ₀),log10(Vᵥ₀)]` , where `Vₗ₀`  and `Vᵥ₀` are initial guesses for the liquid and vapour volumes.
"""
function saturation_pressure(model::EoSModel, T, V0 = x0_sat_pure(model,T))
    T = T*T/T
    V_lb = lb_volume(model,SA[1.0])
    TYPE = promote_type(typeof(T),typeof(V_lb))
    nan = zero(TYPE)/zero(TYPE)
    scales = scale_sat_pure(model)
    f! = (F,x) -> Obj_Sat(model, F, T, exp10(x[1]), exp10(x[2]),scales)
    res0 = (nan,nan,nan)
    result = Ref(res0)
    error_val = Ref{Any}(nothing)
    converged = try_sat_pure(model,V0,f!,T,result,error_val)
    #did not converge, but didnt error.
    if converged
        return result[]
    end
    (T_c, p_c, V_c) = crit_pure(model)
    if abs(T_c-T) < eps(typeof(T))
        return (p_c,V_c,V_c)
    end
    if T_c < T
        @error "initial temperature $T greater than critical temperature $T_c. returning NaN"
    else
        V0 = x0_sat_pure_crit(model,T,T_c,p_c,V_c)
        converged = try_sat_pure(model,V0,f!,T,result,error_val)
        if converged
            return result[]
        end
    end
    #not converged, even trying with better critical aprox.
    return res0
end

#remove when Clapeyron is at 0.3.0
sat_pure(model::EoSModel, T, V0 = x0_sat_pure(model,T)) = saturation_pressure(model,T,V0)
Base.@deprecate sat_pure(model,T) saturation_pressure(model,T)

#with the critical point, we can perform a
#corresponding states approximation with the
#propane reference equation of state
function x0_sat_pure_crit(model,T,T_c,P_c,V_c)
    h = V_c*5000
    T0 = 369.89*T/T_c
    Vl0 = (1.0/_propaneref_rholsat(T0))*h
    Vv0 = (1.0/_propaneref_rhovsat(T0))*h
    _1 = SA[1.0]
    #μ_l = only(VT_chemical_potential(model,Vl0,T,_1))
    #μ_v = only(VT_chemical_potential(model,Vv0,T,_1))
    #@show (μ_l < μ_v,T/T_c)
    #if μ_l < μ_v
      #@show μ_l,μ_v
    #end
    # _,dpdvv = p∂p∂V(model,Vv0,T,SA[1.0])
    # @show dpdvv*Vv0
    # _,dpdvv = p∂p∂V(model,2*Vv0,T,SA[1.0])
    # @show dpdvv*Vv0
    return [log10(Vl0),log10(Vv0)]
end
function sat_pure(model::EoSModel,V0,f!,T,method =LineSearch(Newton()) )
    r = Solvers.nlsolve(f!, V0 ,method )
    Vsol = Solvers.x_sol(r)
    V_l = exp10(Vsol[1])
    V_v = exp10(Vsol[2])
    P_sat = pressure(model,V_v,T)
    return (P_sat,V_l,V_v)
end

function Obj_Sat(model::EoSModel, F, T, V_l, V_v,scales)
    fun(x) = eos(model, x[2], T,SA[x[1]])
    df(x)  = ForwardDiff.gradient(fun,x)
    df_l = df(SA[one(V_l*T),V_l*one(T)])
    df_v = df(SA[one(V_v),V_v*one(T)])
    (p_scale,μ_scale) = scales
    #T̄ = T/T_scale(model)
    F[1] = (df_l[2]-df_v[2])*p_scale
    F[2] = (df_l[1]-df_v[1])*μ_scale
    return F

    #=
    This is the original algorithm, but it seems
    that evaluating both pressure and chemical potential
    at the same time is faster.

    #GERG2008(["water"]) - ForwardDiff one pass - 74.100 μs
    #GERG2008(["water"]) - original - 108.000 μs
    _1 = SA[1.0]
    μ_l = only(VT_chemical_potential(model,V_l,T,_1))
    μ_v = only(VT_chemical_potential(model,V_v,T,_1))
    p_l = -pressure(model,V_l,T,_1)
    p_v = -pressure(model,V_v,T,_1)
    F[1] = (μ_l-μ_v)*p_scale
    F[2] = (p_l-p_v)*μ_scale
    return F
    =#
end

function saturation_temperature(model::EoSModel,p)
    f(z) = Obj_sat_pure_T(model,z,p)
    Tc,pc,vc = crit_pure(model)
    if p>pc
        return (NaN,NaN,NaN)
    else
        T = Roots.find_zero(f,(0.3*Tc,0.99*Tc))
        p,v_l,v_v = saturation_pressure(model,T)
        return T,v_l,v_v
    end
end

function Obj_sat_pure_T(model,T,p)
    p̃,v_l,v_v = saturation_pressure(model,T)
    return p̃-p
end

#=
function Obj_Sat(model::ABCubicModel, F, T, V_l, V_v,V_lb)
    #components = model.components
    ar(V) = a_res(model,V,T)
    function dar(V)
        a,da = Solvers.f∂f(ar,V)
        @show Vval = ForwardDiff.value(V)
        @show aval = ForwardDiff.value(da)
        @show rp = pressure(model,Vval,T)

        p = rp
        @show pval = ForwardDiff.value(p)

         Z = p*V/(R̄*T)
        ϕ = a + (Z-1) - log(Z)
        @show phival = ForwardDiff.value(ϕ)

        return p,exp(ϕ)
    end
    pl,ϕl = dar(V_l)
    pv,ϕv = dar(V_v)
    (p_scale,μ_scale) = scale_sat_pure(model)
    #T̄ = T/T_scale(model)
    F[1] = (pl-pv)*p_scale#*exp(5e-10*(V_l-V_v)^-2)*exp(1e-7*(V_l-V_lb)^-2)
    F[2] = (ϕl-ϕv)*μ_scale#*exp(5e-10*(V_l-V_v)^-2)*exp(1e-7*(V_l-V_lb)^-2)
    return F
end
=#

"""
    crit_pure(model::EoSModel,x0=nothing)

Calculates the critical point of a single component modelled by `model`.

Returns `(Tc, pc, Vc)` where `Tc` is the critical temperature (in K), `pc` is the critical pressure (in Pa) and `Vc` is the critical volume (in  m³)
"""
function crit_pure(model::EoSModel,x0=nothing)
    T̄  = T_scale(model)
    f! = (F,x) -> obj_crit(model, F, x[1]*T̄, exp10(x[2]))
    if x0 === nothing
        x0 = x0_crit_pure(model)
    end

    solver_res = Solvers.nlsolve(f!, x0)
    #print(solver_res)
    r  = Solvers.x_sol(solver_res)
    T_c = r[1]*T̄
    V_c = exp10(r[2])
    p_c = pressure(model, V_c, T_c)
    return (T_c, p_c, V_c)
end

function obj_crit(model::EoSModel, F, T_c, V_c)
    _1 = one(T_c+V_c)
    ∂²A∂V², ∂³A∂V³ = ∂²³f(model, V_c, T_c, SA[_1])
    F[1] = -∂²A∂V²
    F[2] = -∂³A∂V³
    return F
end
"""
    enthalpy_vap(model::EoSModel, T)

Calculates `ΔH`, the difference between saturated vapour and liquid enthalpies at temperature `T`, in J
"""
function enthalpy_vap(model::EoSModel, T)
    (P_sat,V_l,V_v) = saturation_pressure(model,T)
    H_v = VT_enthalpy(model,V_v,T)
    H_l = VT_enthalpy(model,V_l,T)
    H_vap=H_v -H_l
    return H_vap
end

"""
    acentric_factor(model::EoSModel)

calculates the acentric factor using its definition:

    ω = -log10(psatᵣ) -1, at Tᵣ = 0.7
To do so, it calculates the critical temperature (using `crit_pure`) and performs a saturation calculation (with `sat_pure`)

"""
function acentric_factor(model)
    T_c,p_c,_ = crit_pure(model)
    T = 0.7*T_c
    p = first(saturation_pressure(model,T))
    p_r = p/p_c
    return -log10(p_r) - 1.0
end

# function p_sat_t0(model,p)
#     #ln(pr) = h*(1-1/tr)
#     #1/(1-ln(pr)/h)*tc = t
#     tc = T_scale(model)
#     pc = p_scale(model)
#     scale = 0.7
#     t7 = scale*tc
#     p7 = first(sat_pure(model,t7))
#     h = 2.3333333333333335*log(pc/p7)
#     #T/Tc = 1(1-log(p/pc)/h)
#     return 1/(1-log(p/pc)/h)*tc
# end

# function sat_pure_p(model,p)
#     #=
#     model_pred = single_sat_aprox(model)
#     t_pred = p_sat_t0(model,p) #prediction temperature
#     Vv = V_zeros(_pt,model,p,t_pred)

#     V1 = first(Vv)
#     #(V1 == V2) && return V1
#     V2 = last(Vv)

#     =#
#     T0 = p_sat_t0(model,p)
#     px = p

#     _A(V, t) = eos(model, V, t,SA[1.0])
#     V1old = zero(p)
#     V2old = zero(p)
#     Told = zero(p)
#     V1,V2 = x0_sat_pure(model,T0)
#     Tx = T0
#     for i = 1:20
#         if i > 1
#             V1old = V1
#             V2old = V2
#         end
#         A = z -> _A(z, Tx)
#         _V1 = Threads.@spawn volume(model,$p,$Tx;phase=:liquid)
#         _V2 = Threads.@spawn volume(model,$p,$Tx;phase=:gas)
#         V1 = fetch(_V1)
#         V2 = fetch(_V2)

#         if abs(V1 - V1old) / V1 < 1e-15 && i > 1
#             #println("V1 condition")
#             break
#         elseif abs(V2 - V2old) / V2 < 1e-15 && i > 1
#             #println("V2 condition")
#             break
#         end

#         _px(T) = (_A(V1, T) - _A(V2, T)) / (V2 - V1) - p
#         #_px(T) = exp(_A(V1,T)-_A(V2,T)) - one(Tx)
#         Told = Tx
#         Tx = Solvers.ad_newton(_px, Tx)
#         if abs(Told - Tx) < 1e-10 * Tx
#             #println("P condition")
#             break
#         end
#     end

#     return (Tx,V1,V2)
# end

# function naive_sat_pure_p(model,p)
#     T0 = p_sat_t0(model,p)
#     ft0(t) = log(first(sat_pure(model,t))/p)
#     tt0 =T0
#     res = Solvers.ad_newton(ft0,tt0)
#     T = res
#     Vl = volume(model,p,T,phase=:l)
#     Vv = volume(model,p,T,phase=:v)
#     return (T,Vl,Vv)
# end

export saturation_pressure, saturation_temperature, crit_pure, enthalpy_vap

# function spinodals(model,T,Vx = nothing)
#     T7 = (0.7)*one(T)*T_scale(model)
#     #@show T/T_scale(model)
#     #f0(_V) = 1/last(p∂p∂V(model,_V,T))

#     lb_v = lb_volume(model)
#     f0(k) = last(p∂p∂V(model,lb_v + exp(k),T))

#     if Vx === nothing
#     if T > T7
#         sp_l7 = volume_compress(model,zero(T),T7)
#         Tc,pc,Vc = crit_pure(model)
#         V00 = log(sp_l7 - lb_v)
#         sp_l = Roots.find_zero(f0,V00)

#         sp_l = lb_v + exp(sp_l)
#     else
#         #zero pressure approach, the volume_compress function tends to stop at a spinodal at the liquid side
#         sp_l = volume_compress(model,zero(T),T)
#     end
#     else
#         V00 = log(Vx - lb_v)
#         sp_l = Roots.find_zero(f0,V00)
#         sp_l = lb_v + exp(sp_l)

#     end
#     #now calculate gas spinodal
#     #@show Pmin = pressure(model,sp_l,T)
#     #p could be negative, and a volume search fails here.
#     #@show volume_virial(model,Pmin,T)
#     return sp_l#sp_v)
# end


text:
# Provides a structure for turning a mdp into a pomdp where observations are the states of the original mdp

struct FullyObservablePOMDP{S, A} <: POMDP{S,A,S}
    mdp::MDP{S, A}
end

# observations are the state of the MDP
# The observation distribution is modeled by a SparseCat distribution with only one element
POMDPs.observations(pomdp::FullyObservablePOMDP) = states(pomdp.mdp)
POMDPs.n_observations(pomdp::FullyObservablePOMDP) = n_states(pomdp.mdp)
POMDPs.obs_index{S, A}(pomdp::FullyObservablePOMDP{S, A}, o::S) = state_index(pomdp.mdp, o)

POMDPs.convert_o(T::Type{V}, o, pomdp::FullyObservablePOMDP) where {V<:AbstractArray} = convert_s(T, s, pomdp.mdp)
POMDPs.convert_o(T::Type{S}, vec::V, pomdp::FullyObservablePOMDP) where {S,V<:AbstractArray} = convert_s(T, vec, pomdp.mdp)

function POMDPs.generate_o(pomdp::FullyObservablePOMDP, s, a, rng::AbstractRNG)
    return s
end

function POMDPs.observation(pomdp::FullyObservablePOMDP, s, a)
    return SparseCat((s,), (1.,))
end

function POMDPs.observation(pomdp::FullyObservablePOMDP, s, a, sp)
    return SparseCat((sp,), (1.,))
end

POMDPs.isterminal_obs{S,A}(problem::FullyObservablePOMDP{S,A}, o::S) = isterminal(pomdp.mdp, o)

# inherit other function from the MDP type

POMDPs.states(pomdp::FullyObservablePOMDP) = states(pomdp.mdp)
POMDPs.actions(pomdp::FullyObservablePOMDP) = actions(pomdp.mdp)
POMDPs.transition{S,A}(pomdp::FullyObservablePOMDP{S,A}, s::S, a::A) = transition(pomdp.mdp, s, a)
POMDPs.initial_state_distribution(pomdp::FullyObservablePOMDP) = initial_state_distribution(pomdp.mdp)
POMDPs.initial_state(pomdp::FullyObservablePOMDP, rng::AbstractRNG) = initial_state(pomdp.mdp, rng)
POMDPs.generate_s(pomdp::FullyObservablePOMDP, s, a, rng::AbstractRNG) = generate_s(pomdp.mdp, s, a, rng)
POMDPs.generate_sr(pomdp::FullyObservablePOMDP, s, a, rng::AbstractRNG) = generate_sr(pomdp.mdp, s, a, rng)
POMDPs.reward{S,A}(pomdp::FullyObservablePOMDP{S, A}, s::S, a::A) = reward(pomdp.mdp, s, a)
POMDPs.isterminal(pomdp::FullyObservablePOMDP, s) = isterminal(pomdp.mdp, s)
POMDPs.discount(pomdp::FullyObservablePOMDP) = discount(pomdp.mdp)
POMDPs.n_states(pomdp::FullyObservablePOMDP) = n_states(pomdp.mdp)
POMDPs.n_actions(pomdp::FullyObservablePOMDP) = n_actions(pomdp.mdp)
POMDPs.state_index{S,A}(pomdp::FullyObservablePOMDP{S,A}, s::S) = state_index(pomdp.mdp, s)
POMDPs.action_index{S, A}(pomdp::FullyObservablePOMDP{S, A}, a::A) = action_index(pomdp.mdp, a)
POMDPs.convert_s(T::Type{V}, s, pomdp::FullyObservablePOMDP) where V<:AbstractArray = convert_s(T, s, pomdp.mdp)
POMDPs.convert_s(T::Type{S}, vec::V, pomdp::FullyObservablePOMDP) where {S,V<:AbstractArray} = convert_s(T, vec, pomdp.mdp)
POMDPs.convert_a(T::Type{V}, a, pomdp::FullyObservablePOMDP) where V<:AbstractArray = convert_a(T, a, pomdp.mdp)
POMDPs.convert_a(T::Type{A}, vec::V, pomdp::FullyObservablePOMDP) where {A,V<:AbstractArray} = convert_a(T, vec, pomdp.mdp)


text:
export SchedulerHaltedException
export sch_handle, halt!, exec!, get_dag_ids, add_thunk!

"Identifies a thunk by its ID, and preserves the thunk in the scheduler."
struct ThunkID
    id::Int
    ref::Union{DRef,Nothing}
end
ThunkID(id::Int) = ThunkID(id, nothing)

"A handle to the scheduler, used by dynamic thunks."
struct SchedulerHandle
    thunk_id::ThunkID
    out_chan::RemoteChannel
    inp_chan::RemoteChannel
end

"Gets the scheduler handle for the currently-executing thunk."
sch_handle() = task_local_storage(:_dagger_sch_handle)::SchedulerHandle

"Thrown when the scheduler halts before finishing processing the DAG."
struct SchedulerHaltedException <: Exception end

"Thrown when a dynamic thunk encounters an exception in Dagger's utilities."
struct DynamicThunkException <: Exception
    reason::String
end

struct RescheduleSignal end

function safepoint(state)
    if state.halt.set
        # Force dynamic thunks and listeners to terminate
        for (inp_chan,out_chan) in values(state.worker_chans)
            close(inp_chan)
            close(out_chan)
        end
        # Throw out of scheduler
        throw(SchedulerHaltedException())
    end
end

"Processes dynamic messages from worker-executing thunks."
function dynamic_listener!(ctx, state, wid)
    task = current_task() # The scheduler's main task
    inp_chan, out_chan = state.worker_chans[wid]
    listener_task = @async begin
        while isopen(inp_chan) && !state.halt.set
            tid, f, data = try
                take!(inp_chan)
            catch err
                if !(unwrap_nested_exception(err) isa Union{SchedulerHaltedException,
                                                            ProcessExitedException,
                                                            InvalidStateException})
                    iob = IOContext(IOBuffer(), :color=>true)
                    println(iob, "Error in sending dynamic request:")
                    Base.showerror(iob, err)
                    Base.show_backtrace(iob, catch_backtrace())
                    println(iob)
                    seek(iob.io, 0)
                    write(stderr, iob)
                end
                break
            end
            res = try
                (false, lock(state.lock) do
                    Base.invokelatest(f, ctx, state, task, tid, data)
                end)
            catch err
                (true, RemoteException(CapturedException(err,catch_backtrace())))
            end
            try
                put!(out_chan, res)
            catch err
                if !(unwrap_nested_exception(err) isa Union{SchedulerHaltedException,
                                                            ProcessExitedException,
                                                            InvalidStateException})
                    iob = IOContext(IOBuffer(), :color=>true)
                    println(iob, "Error in sending dynamic result from $f:")
                    Base.showerror(iob, err)
                    Base.show_backtrace(iob, catch_backtrace())
                    println(iob)
                    seek(iob.io, 0)
                    write(stderr, iob)
                    println(stderr, res)
                end
            end
        end
    end
    @async begin
        wait(state.halt)
        # TODO: Not sure why we need the @async here, but otherwise we
        # don't stop all the listener tasks
        @async Base.throwto(listener_task, SchedulerHaltedException())
    end
end

## Worker-side methods for dynamic communication

const DYNAMIC_EXEC_LOCK = Threads.ReentrantLock()

"Executes an arbitrary function within the scheduler, returning the result."
function exec!(f, h::SchedulerHandle, args...)
    failed, res = lock(DYNAMIC_EXEC_LOCK) do
        put!(h.out_chan, (h.thunk_id.id, f, args))
        take!(h.inp_chan)
    end
    failed && throw(res)
    res
end

"Commands the scheduler to halt execution immediately."
halt!(h::SchedulerHandle) = exec!(_halt, h, nothing)
function _halt(ctx, state, task, tid, _)
    notify(state.halt)
    put!(state.chan, (1, nothing, SchedulerHaltedException(), nothing))
    Base.throwto(task, SchedulerHaltedException())
end

"Waits on a thunk to complete, and fetches its result."
function Base.fetch(h::SchedulerHandle, id::ThunkID)
    future = ThunkFuture(Future(1))
    exec!(_register_future!, h, future, id, true)
    fetch(future; proc=thunk_processor())
end
"""
Waits on a thunk to complete, and fetches its result. If `check` is set to
`true` (the default), then a domination check will occur to ensure that the
future isn't being registered on a thunk dominated by the calling thunk.
"""
register_future!(h::SchedulerHandle, id::ThunkID, future::ThunkFuture, check::Bool=true) =
        exec!(_register_future!, h, future, id, check)
function _register_future!(ctx, state, task, tid, (future, id, check)::Tuple{ThunkFuture,ThunkID,Bool})
    tid != id.id || throw(DynamicThunkException("Cannot fetch own result"))
    GC.@preserve id begin
        thunk = unwrap_weak_checked(state.thunk_dict[id.id])
        ownthunk = unwrap_weak_checked(state.thunk_dict[tid])
        function dominates(target, t)
            t == target && return true
            seen = Set{Thunk}()
            to_visit = Thunk[t]
            while !isempty(to_visit)
                t = pop!(to_visit)
                if t == target
                    return true
                end
                for input in t.inputs
                    # N.B. Skips expired tasks
                    input = Dagger.unwrap_weak(input)
                    istask(input) || continue
                    input in seen && continue
                    push!(seen, input)
                    push!(to_visit, input)
                end
            end
            return false
        end
        if check && dominates(ownthunk, thunk)
            throw(DynamicThunkException("Cannot fetch result of dominated thunk"))
        end
        # TODO: Assert that future will be fulfilled
        if haskey(state.cache, thunk)
            put!(future, state.cache[thunk]; error=state.errored[thunk])
        else
            futures = get!(()->ThunkFuture[], state.futures, thunk)
            push!(futures, future)
        end
    end
    nothing
end

# TODO: Optimize wait() to not serialize a Chunk
"Waits on a thunk to complete."
function Base.wait(h::SchedulerHandle, id::ThunkID; future=ThunkFuture(1))
    register_future!(h, id, future)
    wait(future)
end

"Returns all Thunks IDs as a Dict, mapping a Thunk to its downstream dependents."
get_dag_ids(h::SchedulerHandle) =
    exec!(_get_dag_ids, h, nothing)::Dict{ThunkID,Set{ThunkID}}
function _get_dag_ids(ctx, state, task, tid, _)
    deps = Dict{ThunkID,Set{ThunkID}}()
    for (id,thunk) in state.thunk_dict
        thunk = unwrap_weak_checked(thunk)
        # TODO: Get at `thunk_ref` for `thunk_id.ref`
        thunk_id = ThunkID(id, nothing)
        if haskey(state.waiting_data, thunk)
            deps[thunk_id] = Set(map(t->ThunkID(t.id, nothing), collect(state.waiting_data[thunk])))
        else
            deps[thunk_id] = Set{ThunkID}()
        end
    end
    deps
end

"Adds a new Thunk to the DAG."
add_thunk!(f, h::SchedulerHandle, args...; future=nothing, ref=nothing, kwargs...) =
    exec!(_add_thunk!, h, f, args, kwargs, future, ref)
function _add_thunk!(ctx, state, task, tid, (f, args, kwargs, future, ref))
    timespan_start(ctx, :add_thunk, tid, 0)
    _args = map(arg->arg isa ThunkID ? state.thunk_dict[arg.id] : arg, args)
    GC.@preserve _args begin
        thunk = Thunk(f, _args...; kwargs...)
        # Create a `DRef` to `thunk` so that the caller can preserve it
        thunk_ref = poolset(thunk; size=64)
        thunk_id = ThunkID(thunk.id, thunk_ref)
        state.thunk_dict[thunk.id] = WeakThunk(thunk)
        reschedule_inputs!(state, thunk)
        if future !== nothing
            # Ensure we attach a future before the thunk is scheduled
            _register_future!(ctx, state, task, tid, (future, thunk_id, false))
        end
        if ref !== nothing
            # Preserve the `EagerThunkFinalizer` through `thunk`
            thunk.eager_ref = ref
        end
        state.valid[thunk] = nothing
        put!(state.chan, RescheduleSignal())
        timespan_finish(ctx, :add_thunk, tid, 0)
        return thunk_id
    end
end


text:
function cube_morphing(D; vis=Visualizer(), fps=30, rot=0.00, vis_truth::Bool=true,
		vis_learned::Bool=true, translate::Bool=true, cam_pos=[0,-10,4.], zoom=1.8,
		color_truth=RGBA(1.0, 153.0 / 255.0, 51.0 / 255.0, 1.0),
		color_learn=RGBA(51.0 / 255.0, 1.0, 1.0, 1.0),
		background=true, alt=0.0, b0=0.2, b1=1.5)
	setvisible!(vis["/Background"], background)
	setprop!(vis["/Background"], "top_color", RGBA(1.0, 1.0, 1.0, 1.0))
	setprop!(vis["/Background"], "bottom_color", RGBA(1.0, 1.0, 1.0, 1.0))
	setvisible!(vis["/Axes"], false)
	setvisible!(vis["/Grid"], false)
	set_camera!(vis, zoom=zoom, cam_pos=cam_pos)

	anim = MeshCat.Animation(fps)
	b0 = Int(floor(b0*fps))
	b1 = Int(floor(b1*fps))
	D = [fill(D[1], b0); D; fill(D[end] .+ 0.005*rand(length(D[end])), 2b1)]
	setobject!(vis[:cube],
		Rect(Vec(-1,-1,-1),Vec(2,2,2.0)),
		MeshPhongMaterial(color=color_truth))

	for (i,d) in enumerate(D)
		v = vrep([d[1 + (i-1)*3 .+ (1:3)] for i = 1:8])
		p = polyhedron(v)
		Polyhedra.Mesh(p)
		setobject!(vis[Symbol("polyhedra$i")],
			Polyhedra.Mesh(p),
			MeshPhongMaterial(color=color_learn))
	end

	for (i,d) in enumerate(D)
		atframe(anim, i) do
			MeshCat.setvisible!(vis[:cube], vis_truth)
			for (j,d) in enumerate(D)
				MeshCat.setvisible!(vis[Symbol("polyhedra$j")], i==j && vis_learned)
			end
			settransform!(vis[:cube], MeshCat.compose(
				MeshCat.Translation(-2.0*translate,0,1+alt),
				MeshCat.LinearMap(RotZ(rot*min(i,length(D)-b1))),
				))

			settransform!(vis[Symbol("polyhedra$i")], MeshCat.compose(
				MeshCat.Translation(2.0*translate,0,1+alt),
				MeshCat.LinearMap(RotZ(rot*min(i,length(D)-b1))),
				))

		end
	end
	setanimation!(vis, anim)
	return vis, anim
end

function cone_morphing(D; vis=Visualizer(), fps=30, rot=0.00, vis_truth::Bool=true,
		vis_learned::Bool=true, translate::Bool=true, cam_pos=[0,-10,4.], zoom=1.8,
		color_truth=RGBA(1.0, 153.0 / 255.0, 51.0 / 255.0, 1.0),
		color_learn=RGBA(51.0 / 255.0, 1.0, 1.0, 1.0),
		background=true, alt=0.0, b0=0.2, b1=1.5)
	setvisible!(vis["/Background"], background)
	setprop!(vis["/Background"], "top_color", RGBA(1.0, 1.0, 1.0, 1.0))
	setprop!(vis["/Background"], "bottom_color", RGBA(1.0, 1.0, 1.0, 1.0))
	setvisible!(vis["/Axes"], false)
	setvisible!(vis["/Grid"], false)
	set_camera!(vis, cam_pos=cam_pos, zoom=zoom)

	anim = MeshCat.Animation(fps)
	b0 = Int(floor(b0*fps))
	b1 = Int(floor(b1*fps))
	D = [fill(D[1], b0); D; fill(D[end] .+ 0.005*rand(length(D[end])), 2b1)]
	friction_coefficient_truth = 0.21
	setobject!(vis[:cone],
		MeshCat.Cone(Point(0,0,1.0),Point(0,0,0.0),friction_coefficient_truth),
		MeshPhongMaterial(color=color_truth))

	for (i,d) in enumerate(D)
		setobject!(vis[Symbol("cone$i")],
			MeshCat.Cone(Point(0,0,1.0),Point(0,0,0.0),d[1]),
			MeshPhongMaterial(color=color_learn))
	end

	for (i,d) in enumerate(D)
		atframe(anim, i) do
			MeshCat.setvisible!(vis[:cone], vis_truth)
			for (j,d) in enumerate(D)
				MeshCat.setvisible!(vis[Symbol("cone$j")], i==j && vis_learned)
			end
			settransform!(vis[:cone], MeshCat.compose(
				MeshCat.Translation(-0.5*translate,0,1+alt),
				MeshCat.LinearMap(RotZ(rot*min(i,length(D)-b1))),
				))

			settransform!(vis[Symbol("cone$i")], MeshCat.compose(
				MeshCat.Translation(0.5*translate,0,1+alt),
				MeshCat.LinearMap(RotZ(rot*min(i,length(D)-b1))),
				))

		end
	end
	setanimation!(vis, anim)
	return vis, anim
end

function cube_sim_v_truth(d, traj_truth::Storage{T,HT}, traj_sim::Storage{T,HS};
		vis=Visualizer(), fps=30, zoom=1.8,
		transparency_truth=1.0,
		transparency_learn=1.0,
		color_truth=RGBA(1.0, 153.0 / 255.0, 51.0 / 255.0, transparency_truth),
		color_learn=RGBA(51.0 / 255.0, 1.0, 1.0, transparency_learn),
		background=true, b0=0.2, b1=1.5) where {T,HT,HS}
	setvisible!(vis["/Background"], background)
	setprop!(vis["/Background"], "top_color", RGBA(1.0, 1.0, 1.0, 1.0))
	setprop!(vis["/Background"], "bottom_color", RGBA(1.0, 1.0, 1.0, 1.0))
	setvisible!(vis["/Axes"], false)
	# setvisible!(vis["/Grid"], false)
	set_camera!(vis, zoom=zoom)
	setprop!(vis["/Lights/AmbientLight/<object>"], "intensity", 0.80)

	anim = MeshCat.Animation(fps)
	b0 = Int(floor(b0*fps))
	b1 = Int(floor(b1*fps))
	H = max(HT,HS)

	# Create Objects
	setobject!(vis[:truth],
		Rect(Vec(-1,-1,-1),Vec(2,2,2.0)),
		MeshPhongMaterial(color=color_truth))

	v = vrep([d[1 + (i-1)*3 .+ (1:3)] for i = 1:8])
	p = polyhedron(v)
	Polyhedra.Mesh(p)
	setobject!(vis[:sim],
		Polyhedra.Mesh(p),
		MeshPhongMaterial(color=color_learn))

	# Animate
	for i = 1:b0
		atframe(anim, i) do
			settransform!(vis[:truth], MeshCat.compose(
				MeshCat.Translation(traj_truth.x[1][1]),
				MeshCat.LinearMap(traj_truth.q[1][1]),
				))
			settransform!(vis[:sim], MeshCat.compose(
				MeshCat.Translation(traj_sim.x[1][1]),
				MeshCat.LinearMap(traj_sim.q[1][1]),
				))
		end
	end
	for i = 1:H
		atframe(anim, b0+i) do
			settransform!(vis[:truth], MeshCat.compose(
				MeshCat.Translation(traj_truth.x[1][min(i,HT)]),
				MeshCat.LinearMap(traj_truth.q[1][min(i,HT)]),
				))
			settransform!(vis[:sim], MeshCat.compose(
				MeshCat.Translation(traj_sim.x[1][min(i,HS)]),
				MeshCat.LinearMap(traj_sim.q[1][min(i,HS)]),
				))
		end
	end
	for i = 1:b1
		atframe(anim, b0+H+i) do
			settransform!(vis[:truth], MeshCat.compose(
				MeshCat.Translation(traj_truth.x[1][end]),
				MeshCat.LinearMap(traj_truth.q[1][end]),
				))
			settransform!(vis[:sim], MeshCat.compose(
				MeshCat.Translation(traj_sim.x[1][end]),
				MeshCat.LinearMap(traj_sim.q[1][end]),
				))
		end
	end
	setanimation!(vis, anim)
	return vis, anim
end

function cube_ghost_sim_v_truth(d, traj_truth::Storage{T,HT}, traj_sim::Storage{T,HS};
		vis=Visualizer(), zoom=1.8,
		transparency_truth=1.0,
		transparency_learn=1.0,
		color_truth=RGBA(1.0, 153.0 / 255.0, 51.0 / 255.0, transparency_truth),
		color_learn=RGBA(51.0 / 255.0, 1.0, 1.0, transparency_learn),
		background=true) where {T,HT,HS}
	setvisible!(vis["/Background"], background)
	setprop!(vis["/Background"], "top_color", RGBA(1.0, 1.0, 1.0, 1.0))
	setprop!(vis["/Background"], "bottom_color", RGBA(1.0, 1.0, 1.0, 1.0))
	setvisible!(vis["/Axes"], false)
	# setvisible!(vis["/Grid"], false)
	set_camera!(vis, zoom=zoom)
	setprop!(vis["/Lights/AmbientLight/<object>"], "intensity", 0.80)

	H = max(HT,HS)
	frames = [5, 30, H]
	α = [0.3, 0.5,1.0]
	# Create Objects
	for (j,i) in enumerate(frames)
		setobject!(vis["truth$(i)"],
			Rect(Vec(-1,-1,-1),Vec(2,2,2.0)),
			MeshPhongMaterial(color=color_truth))

		v = vrep([d[1 + (i-1)*3 .+ (1:3)] for i = 1:8])
		p = polyhedron(v)
		Polyhedra.Mesh(p)
		setobject!(vis["sim$(i)"],
			Polyhedra.Mesh(p),
			MeshPhongMaterial(color=color_learn))
	end

	# Animate
	for i in frames
		settransform!(vis["truth$(i)"], MeshCat.compose(
			MeshCat.Translation(traj_truth.x[1][min(i,HT)]),
			MeshCat.LinearMap(traj_truth.q[1][min(i,HT)]),
			))
		settransform!(vis["sim$(i)"], MeshCat.compose(
			MeshCat.Translation(traj_sim.x[1][min(i,HS)]),
			MeshCat.LinearMap(traj_sim.q[1][min(i,HS)]),
			))
	end
	return vis
end


text:
export setsrc!

function  setsrc! end

function get_src_comp(sim::Simulation{Dim},  dir::Direction{D}) where {D,Dim}
	if D == 1 return sim.J_x end
	if D == 2 return sim.J_y end
	if D == 3 return sim.J_z end
end

function (sim::Simulation)(::typeof(setsrc!),dir::Direction,src::AbstractArray)
  get_src_comp(sim,  dir) = src
end

###################
function (sim::Simulation)(::typeof(setsrc!),dir::Direction,val::Function, reg::Function = (x -> true) , gridtype::GridType= p̂)
  sim(set!, get_src_comp(sim,  dir) , val,reg,gridtype)
end

include("PlaneWave.jl")


text:
module SpikingRNNs
using LinearAlgebra,Statistics,StatsBase,Random,Distributions
using SparseArrays
using Colors # to save spike rasters as png
using DSP # this is to have convolutions to analyze spike trains
using ExtendableSparse # to optimize adding elements in sparse matrices

# all type declarations should go here :-/

abstract type NeuronType end
abstract type SynapticKernel end
abstract type NTConductance <:NeuronType end
abstract type SpikeGenFunction end

# Base elements

# neuron parameterns and behavior
Base.broadcastable(t::NeuronType)=Ref(t)

# Dynamical state and memory allocation for a particular group of neurons
# n is always the number of neurons, and neurontype is the neuron type
abstract type PopulationState{NT<:NeuronType} end
nneurons(ps::PopulationState) = ps.n

# Everything related to connection (including plasticity, etc)
abstract type Connection{N} end
abstract type PlasticityRule end

struct NoPlasticity <: PlasticityRule end
reset!(::NoPlasticity) = nothing

struct FakeConnection{N,PL<:NTuple{N,PlasticityRule}} <: Connection{N}
  weights::SparseMatrixCSC{Float64,Int64}
  plasticities::PL
  function FakeConnection()
    w = sparse(fill(NaN,1,1))
    return new{0,NTuple{0,NoPlasticity}}(w,())
  end
end

struct BaseConnection{N,PL<:NTuple{N,PlasticityRule}} <: Connection{N}
  weights::SparseMatrixCSC{Float64,Int64}
  plasticities::PL
end
function BaseConnection(w::SparseMatrixCSC)
  return BaseConnection{0,NTuple{0,NoPlasticity}}(w,())
end

# weights evolve according to plasticity rules and pre-post activity
# BUT neurons do not exchange signals of any kind
# this works well in tandem with input neurons with specific spike times
# to test plasticity rules
struct ConnectionPlasticityTest{N,PL<:NTuple{N,PlasticityRule}} <: Connection{N}
  weights::SparseMatrixCSC{Float64,Int64}
  plasticities::PL
end
function ConnectionPlasticityTest(weights::SparseMatrixCSC,
    (plasticities::PlasticityRule)...)
  return ConnectionPlasticityTest(weights,plasticities)
end
function reset!(conn::ConnectionPlasticityTest)
  reset!.(conn.plasticities)
  return nothing
end

@inline function n_plasticity_rules(::Connection{N}) where N
  return N
end

struct PSSimpleInput{In} <: PopulationState{In}
  neurontype::In # not really a neuron, but I keep the name for consistency
  n::Int64
  function PSSimpleInput(in::N) where N<:NeuronType
    return new{N}(in,1)
  end
end

function rand_pop_label()
  return Symbol(randstring(3))
end

abstract type AbstractPopulation{N,PS} end

# this is for populations without incoming synapses (so no presynaptic
# population and no connections). Useful for input units.
struct UnconnectedPopulation{N,PS} <: AbstractPopulation{N,PS}
  label::Symbol
  state::PS
end
function UnconnectedPopulation(ps::PS;
    label::Union{Nothing,String}=nothing) where PS<:PopulationState
  label = isnothing(label) ? rand_pop_label() : Symbol(label)
  return UnconnectedPopulation{0,PS}(label,ps)
end

struct Population{N,PS<:PopulationState,
    TC<:NTuple{N,Connection},
    TP<:NTuple{N,PopulationState}} <:AbstractPopulation{N,PS}
  label::Symbol # I use symbols because it might be a DataFrame column name
  state::PS
  connections::TC
  pre_states::TP
end
nneurons(p::Population) = nneurons(p.state)

function Population(state::PopulationState,
    (conn_pre::Tuple{C,PS} where {C<:Connection,PS<:PopulationState})... ;
    label::Union{Nothing,String}=nothing)
  connections = Tuple(getindex.(conn_pre,1))
  pre_states = Tuple(getindex.(conn_pre,2))
  label = isnothing(label) ? rand_pop_label() : Symbol(label)
  return Population(label,state,connections,pre_states)
end

@inline function n_prepops(p::Population{N,PS,TC,TP}) where {N,PS,TC,TP}
  return N
end
@inline function n_prepops(p::AbstractPopulation{N,PS}) where {N,PS}
  return N
end

abstract type AbstractNetwork end
struct RecurrentNetwork{N,TP<:NTuple{N,AbstractPopulation}} <: AbstractNetwork
  dt::Float64
  populations::TP
end

function RecurrentNetwork(dt,(pops::P where P<:AbstractPopulation)...)
  RecurrentNetwork(dt,pops)
end

# fallback functions

# updates the state of the population locally, depending on
# some input gathered in previous steps
# usually it's simply the vector p.input
function local_update!(tnow::Real,dt::Real,p::PopulationState)
  return nothing
end
function local_update!(tnow::Real,dt::Real,p::AbstractPopulation)
  return local_update!(tnow,dt,p.state)
end

# reset the input gathered in the previous steps
# usually it's simply the vector p.input
@inline function reset_input!(ps::PopulationState)
  fill!(ps.input,0.0)
  return nothing
end
function reset_input!(p::AbstractPopulation)
  return reset_input!(p.state)
end
# or sometimes there are spikes
# stored in the BitVector called `isfiring`
@inline function reset_spikes!(spk::BitArray{1})
  fill!(spk,false)
	return nothing
end
@inline function reset_spikes!(popst::PopulationState)
	return reset_spikes!(popst.isfiring)
end

# Computes the input, considering each connection and each presynaptic population
function forward_signal!(tnow::Real,dt::Real,p_post::PopulationState,c::Connection,p_pre::PopulationState)
  return nothing
end
function forward_signal!(tnow::Real,dt::Real,p::AbstractPopulation)
  for i in 1:n_prepops(p)
    forward_signal!(tnow,dt,p.state,p.connections[i],p.pre_states[i])
  end
  return nothing
end
function plasticity_update!(tnow::Real,dt::Real,p::AbstractPopulation)
  for i in 1:n_prepops(p)
    for r in 1:n_plasticity_rules(p.connections[i])
      plasticity_update!(tnow,dt,p.state,p.connections[i],
        p.pre_states[i],p.connections[i].plasticities[r])
    end
  end
  return nothing
end

# this is the network iteration
function dynamics_step!(t_now::Float64,ntw::RecurrentNetwork)
  reset_input!.(ntw.populations)
  forward_signal!.(t_now,ntw.dt,ntw.populations)
  local_update!.(t_now,ntw.dt,ntw.populations)
  plasticity_update!.(t_now,ntw.dt,ntw.populations)
  return nothing
end

dynamics_step!(ntw::RecurrentNetwork) = dynamics_step!(NaN,ntw)

##

@inline hardbounds(x::R,low::R,high::R) where R = min(high,max(x,low))

##

include("firingneurons_shared.jl")
include("firingneurons_plasticity_shared.jl")
include("firingneurons_heterosynaptic_plasticity.jl")
include("firingneurons_structural_plasticity.jl")
include("connectivity_shared.jl")
include("rate_models.jl")
include("lif_current.jl")
include("lif_exponential.jl")
include("lif_conductance.jl")
include("inputs_shared.jl")
include("recorders_shared.jl")

#=
include("./hawkes.jl")
=#

end # of SpikingRNNs module


text:
## ExoplanetsSysSim/src/koi_table.jl
## (c) 2015 Eric B. Ford

# Note this file is currently not used by SysSim.
# This functionality is now in kepler_catalog.jl

module KoiTable
using ExoplanetsSysSim
#using DataArrays
using DataFrames
using CSV

export setup_koi_table, koi_table, num_koi_for_kepid

df = DataFrame()
usable = Array{Int64}(0)

default_koi_symbols_to_keep = [ :kepoi_name, :koi_vet_stat, :koi_pdisposition, :koi_period, :koi_time0bk, :koi_duration, :koi_ingress, :koi_depth, :koi_ror, :koi_prad, :koi_srad, :koi_smass, :koi_steff, :koi_slogg, :koi_smet ]

function setup(sim_param::SimParam; force_reread::Bool = false, symbols_to_keep::Vector{Symbol} = default_koi_symbols_to_keep )
  global df, usable
  if haskey(sim_param,"read_koi_catalog") && !force_reread
     return df
  end
  planetary_catalog = joinpath(dirname(pathof(ExoplanetsSysSim)),"..", "data", get(sim_param,"planetary_catalog","q1_q17_dr25_koi.csv") )
  add_param_fixed(sim_param,"read_koi_catalog",true)
  try
    #df = readtable(planetary_catalog)
    #df = CSV.read(planetary_catalog,nullable=true)
    df = CSV.read(planetary_catalog,allowmissing=:all)
  catch
    error(string("# Failed to read koi catalog >",planetary_catalog,"<."))
  end

  has_planet = ! (isna(df[:koi_period]) | isna(df[:koi_time0bk]) | isna(df[:koi_duration]) | isna(:koi_depth) )
  has_star = ! ( isna(:koi_srad) )
  is_usable = has_planet & has_star

  delete!(df, [~(x in symbols_to_keep) for x in names(df)])    # delete columns that we won't be using anyway
  usable = find(is_usable)
  df = df[usable, symbols_to_keep]
end

setup_koi_table(sim_param::SimParam) = setup(sim_param::SimParam)

function kepids_w_kois()
  unique(df[:,:kepid])
end

function df_for_kepid(kepid::Integer)
  df[df[:kepid].==kepid,:]
end

function num_koi(kepid::Integer)
  sum(df[:kepid].==kepid)
end

function koi_by_kepid(kepid::Integer, plid::Integer, sym::Symbol)
  kepid_idx = df[:kepid].==kepid
  per_perm = sortperm(df[kepid_idx,:koi_period])
  @assert( 1<= plid <= length(per_perm) )
  df[kepid_idx,sym][per_perm[plid]]
end

function num_usable()
  global usable
  length(usable)
end

num_usable_in_koi_table() = num_usable()

function idx(i::Integer)
  global usable
  @assert( 1<=i<=length(usable) )
  usable[i]
end

function koi_table(i::Integer, sym::Symbol)
  global df, usable
  @assert( 1<=i<=length(usable) )
  return df[i,sym]
  #return df[usable[i],sym]
end

function koi_table(i::Integer)
  global data
  return df[i,:]
  #return df[usable[i],:]
end

function koi_table(i::Integer, sym::Vector{Symbol})
  global df, usable
  @assert( 1<=i<=length(usable) )
  return df[i,sym]
  #return df[usable[i],sym]
end

function koi_table(i::Vector{Integer}, sym::Symbol)
  global df, usable
  return df[i,sym]
  #return df[usable[i],sym]
end

function koi_table(i::Vector{Integer}, sym::Vector{Symbol})
  global df, usable
  return df[i,sym]
  #return df[usable[i],sym]
end

end # module KoiTable

using ExoplanetsSysSim.KoiTable


text:
"""
    mpiexec(fn)

A wrapper function for the MPI launcher executable. Calls `fn(cmd)`, where `cmd` is a `Cmd` object of the MPI launcher.

# Usage

```jldoctest
julia> mpiexec(cmd -> run(`\$cmd -n 3 echo hello world`));
hello world
hello world
hello world
```
"""
mpiexec

# Administrative functions
function _warn_if_wrong_mpi()
    # warn if we have only one process but environment variables
    # suggest we should have more
    if Comm_rank(MPI.COMM_WORLD) == 0 && Comm_size(MPI.COMM_WORLD) == 1
        known_envvars = ("MPI_LOCALNRANKS", #MPICH
                         "OMPI_COMM_WORLD_SIZE" # OpenMPI
                         )
        if any(v -> haskey(ENV, v) && parse(Int, ENV[v]) > 1, known_envvars)
            @warn """
                You appear to have run julia under a different `mpiexec` than the one used by MPI.jl.
                See the documentation for details.
            """
        end
    end
end

"""
    MPI.free(obj)

Free the MPI object handle `obj`. This is typically used as the finalizer, and so need not be called directly unless otherwise noted.
"""
function free
end

const mpi_init_hooks = Any[]
add_init_hook!(f) = push!(mpi_init_hooks, f)
function run_init_hooks()
    for f in mpi_init_hooks
        f()
    end
    empty!(mpi_init_hooks)
    return nothing
end

"""
    Init(;threadlevel=:serialized, finalize_atexit=true, errors_return=true)

Initialize MPI in the current process. The keyword options:

- `threadlevel`: either `:single`, `:funneled`, `:serialized` (default),
  `:multiple`, or an instance of [`ThreadLevel`](@ref).
- `finalize_atexit`: if `true` (default), adds an `atexit` hook to call
  [`MPI.Finalize`](@ref) if it hasn't already been called.
- `errors_return`: if `true` (default), will set the default error handlers for
  [`MPI.COMM_SELF`](@ref) and [`MPI.COMM_WORLD`](@ref) to be
  `MPI.ERRORS_RETURN`. MPI errors will then appear as Julia exceptions.

It will return the [`ThreadLevel`](@ref) value which MPI is initialized at.

All MPI programs must call this function at least once before calling any other
MPI operations: the only MPI functions that may be called before `MPI.Init` are
[`MPI.Initialized`](@ref) and [`MPI.Finalized`](@ref).

It is safe to call `MPI.Init` multiple times, however it is not valid to call
it after calling [`MPI.Finalize`](@ref).

# External links
$(_doc_external("MPI_Init"))
$(_doc_external("MPI_Init_thread"))
"""
function Init(;threadlevel=:serialized, finalize_atexit=true, errors_return=true)
    if threadlevel isa Symbol
        threadlevel = ThreadLevel(threadlevel)
    end
    if MPI.Finalized()
        error("MPI cannot be initialized after MPI.Finalize has been called.")
    end
    if MPI.Initialized()
        provided = Query_thread()
        if provided < threadlevel
            @warn "MPI already initialized with thread level $provided, requested = $threadlevel"
        end
    else
        provided = _init_thread(threadlevel)
        if provided < threadlevel
            @warn "MPI thread level requested = $threadlevel, provided = $provided"
        end

        if finalize_atexit
            atexit() do
                # MPI_Finalize is a collective and can act like a barrier (this may be implementation
                # specific). If we are terminating due to a Julia exception, we shouldn't call
                # MPI_Finalize. We thus peek at the current exception, and only if that field is
                # nothing do we terminate.
                if !Finalized() && ccall(:jl_current_exception, Any, ()) === nothing
                    Finalize()
                end
            end
        end

        run_init_hooks()
        if errors_return
            set_default_error_handler_return()
        end
        _warn_if_wrong_mpi()
    end
    return provided
end

"""
    ThreadLevel

An Enum denoting the level of threading support in the current process:

 - `MPI.THREAD_SINGLE`: Only one thread will execute.

 - `MPI.THREAD_FUNNELED`: The process may be multi-threaded, but the application must
   ensure that only the main thread makes MPI calls. See [`Is_thread_main`](@ref).

 - `MPI.THREAD_SERIALIZED`: The process may be multi-threaded, and multiple threads may
   make MPI calls, but only one at a time (i.e. all MPI calls are serialized).

 - `MPI.THREAD_MULTIPLE`: Multiple threads may call MPI, with no restrictions.

# See also

- [`Init`](@ref)
- [`Query_thread`](@ref)
"""
mutable struct ThreadLevel
    val::Cint
end
const THREAD_SINGLE     = ThreadLevel(Consts.MPI_THREAD_SINGLE[])
const THREAD_FUNNELED   = ThreadLevel(Consts.MPI_THREAD_FUNNELED[])
const THREAD_SERIALIZED = ThreadLevel(Consts.MPI_THREAD_SERIALIZED[])
const THREAD_MULTIPLE   = ThreadLevel(Consts.MPI_THREAD_MULTIPLE[])
add_load_time_hook!(() -> THREAD_SINGLE.val     = Consts.MPI_THREAD_SINGLE[]    )
add_load_time_hook!(() -> THREAD_FUNNELED.val   = Consts.MPI_THREAD_FUNNELED[]  )
add_load_time_hook!(() -> THREAD_SERIALIZED.val = Consts.MPI_THREAD_SERIALIZED[])
add_load_time_hook!(() -> THREAD_MULTIPLE.val   = Consts.MPI_THREAD_MULTIPLE[]  )
ThreadLevel(threadlevel::Symbol) =
    threadlevel == :single ? THREAD_SINGLE :
    threadlevel == :funneled ? THREAD_FUNNELED :
    threadlevel == :serialized ? THREAD_SERIALIZED :
    threadlevel == :multiple ? THREAD_MULTIPLE :
    error("Invalid threadlevel: must be one of :single, :funneled, :serialized, or :multiple")
Base.:(==)(tl1::ThreadLevel, tl2::ThreadLevel) = tl1.val == tl2.val
Base.isless(tl1::ThreadLevel, tl2::ThreadLevel) = tl1.val < tl2.val

function _init_thread(required::ThreadLevel)
    r_provided = Ref{Cint}()
    # int MPI_Init_thread(int *argc, char ***argv, int required, int *provided)
    @mpichk ccall((:MPI_Init_thread, libmpi), Cint,
                    (Ptr{Cint},Ptr{Cvoid}, Cint, Ref{Cint}),
                    C_NULL, C_NULL, required.val, r_provided)
    return ThreadLevel(r_provided[])
end

"""
    Query_thread()

Query the level of threading support in the current process.
Returns a [`ThreadLevel`](@ref) value denoting

# External links
$(_doc_external("MPI_Query_thread"))
"""
function Query_thread()
    r_provided = Ref{Cint}()

    # int MPI_Query_thread(int *provided)
    @mpichk ccall((:MPI_Query_thread, libmpi), Cint,
                  (Ref{Cint},), r_provided)
    return ThreadLevel(r_provided[])
end

"""
    Is_thread_main()

Queries whether the current thread is the main thread according to MPI. This can be called
by any thread, and is useful for the  `THREAD_FUNNELED` [`ThreadLevel`](@ref).

# External links
$(_doc_external("MPI_Is_thread_main"))
"""
function Is_thread_main()
    r_flag = Ref{Cint}()
    # int MPI_Is_thread_main(int *flag)
    @mpichk ccall((:MPI_Is_thread_main, libmpi), Cint,
                  (Ref{Cint},), r_flag)
    return r_flag[] != 0
end

"""
    Finalize()

Marks MPI state for cleanup. This should be called after [`MPI.Init`](@ref), and
can be called at most once. No further MPI calls (other than
[`Initialized`](@ref) or [`Finalized`](@ref)) should be made after it is called.

[`MPI.Init`](@ref) will automatically insert a hook to call this function when
Julia exits, if it hasn't already been called.

# External links
$(_doc_external("MPI_Finalize"))
"""
function Finalize()
    if !MPI.Finalized()
        @mpichk ccall((:MPI_Finalize, libmpi), Cint, ())
    end
    return nothing
end

"""
    Abort(comm::Comm, errcode::Integer)

Make a “best attempt” to abort all tasks in the group of `comm`. This function does not
require that the invoking environment take any action with the error code. However, a Unix
or POSIX environment should handle this as a return errorcode from the main program.

# External links
$(_doc_external("MPI_Abort"))
"""
function Abort(comm::Comm, errcode::Integer)
    @mpichk ccall((:MPI_Abort, libmpi), Cint, (MPI_Comm, Cint), comm, errcode)
end

"""
    Initialized()

Returns `true` if [`MPI.Init`](@ref) has been called, `false` otherwise.

It is unaffected by [`MPI.Finalize`](@ref), and is one of the few functions that may be
called before [`MPI.Init`](@ref).

# External links
$(_doc_external("MPI_Intialized"))
"""
function Initialized()
    flag = Ref{Cint}()
    @mpichk ccall((:MPI_Initialized, libmpi), Cint, (Ptr{Cint},), flag)
    flag[] != 0
end

"""
    Finalized()

Returns `true` if [`MPI.Finalize`](@ref) has completed, `false` otherwise.

It is safe to call before [`MPI.Init`](@ref) and after [`MPI.Finalize`](@ref).

# External links
$(_doc_external("MPI_Finalized"))
"""
function Finalized()
    flag = Ref{Cint}()
    @mpichk ccall((:MPI_Finalized, libmpi), Cint, (Ptr{Cint},), flag)
    flag[] != 0
end

function Wtick()
    @mpicall ccall((:MPI_Wtick, libmpi), Cdouble, ())
end

function Wtime()
    @mpicall ccall((:MPI_Wtime, libmpi), Cdouble, ())
end

"""
    MPI.has_cuda()

Check if the MPI implementation is known to have CUDA support. Currently only Open MPI
provides a mechanism to check, so it will return `false` with other implementations
(unless overriden).

This can be overriden by setting the `JULIA_MPI_HAS_CUDA` environment variable to `true`
or `false`.
"""
function has_cuda()
    flag = get(ENV, "JULIA_MPI_HAS_CUDA", nothing)
    if flag === nothing
        # Only Open MPI provides a function to check CUDA support
        @static if MPI_LIBRARY == OpenMPI
            # int MPIX_Query_cuda_support(void)
            return 0 != ccall((:MPIX_Query_cuda_support, libmpi), Cint, ())
        elseif MPI_LIBRARY == IBMSpectrumMPI
            return true
        else
            return false
        end
    else
        return parse(Bool, flag)
    end
end


text:
using KitBase, FluxReconstruction, OrdinaryDiffEq, LinearAlgebra, Plots
using ProgressMeter: @showprogress

begin
    x0 = 0
    x1 = 1
    ncell = 100
    dx = (x1 - x0) / ncell
    deg = 2 # polynomial degree
    nsp = deg + 1
    γ = 5 / 3
    cfl = 0.05
    dt = cfl * dx
    t = 0.0
end
ps = FRPSpace1D(x0, x1, ncell, deg)

u = zeros(ncell, nsp, 3)
for i = 1:ncell, ppp1 = 1:nsp
    if ps.x[i] <= 0.5
        prim = [1.0, 0.0, 0.5]
    else
        prim = [0.3, 0.0, 0.625]
    end
    #prim = [1 + 0.1*sin(2π * ps.xpg[i, ppp1]), 1.0, 1.0]

    u[i, ppp1, :] .= prim_conserve(prim, γ)
end

function dudt!(du, u, p, t)
    du .= 0.0
    nx, nsp, J, ll, lr, lpdm, dgl, dgr, γ = p

    ncell = size(u, 1)
    nsp = size(u, 2)

    f = zeros(ncell, nsp, 3)
    for i = 1:ncell, j = 1:nsp
        f[i, j, :] .= euler_flux(u[i, j, :], γ)[1] ./ J[i]
    end

    u_face = zeros(ncell, 3, 2)
    f_face = zeros(ncell, 3, 2)
    for i = 1:ncell, j = 1:3
        # right face of element i
        u_face[i, j, 1] = dot(u[i, :, j], lr)
        f_face[i, j, 1] = dot(f[i, :, j], lr)

        # left face of element i
        u_face[i, j, 2] = dot(u[i, :, j], ll)
        f_face[i, j, 2] = dot(f[i, :, j], ll)
    end

    f_interaction = zeros(nx + 1, 3)
    for i = 2:nx
        fw = @view f_interaction[i, :]
        flux_hll!(fw, u_face[i-1, :, 1], u_face[i, :, 2], γ, 1.0)
    end
    fw = @view f_interaction[1, :]
    flux_hll!(fw, u_face[nx, :, 1], u_face[1, :, 2], γ, 1.0)
    fw = @view f_interaction[nx+1, :]
    flux_hll!(fw, u_face[nx, :, 1], u_face[1, :, 2], γ, 1.0)

    rhs1 = zeros(ncell, nsp, 3)
    for i = 1:ncell, ppp1 = 1:nsp, k = 1:3
        rhs1[i, ppp1, k] = dot(f[i, :, k], lpdm[ppp1, :])
    end

    idx = 2:ncell-1 # ending points are Dirichlet
    for i in idx, ppp1 = 1:nsp, k = 1:3
        du[i, ppp1, k] =
            -(
                rhs1[i, ppp1, k] +
                (f_interaction[i, k] / J[i] - f_face[i, k, 2]) * dgl[ppp1] +
                (f_interaction[i+1, k] / J[i] - f_face[i, k, 1]) * dgr[ppp1]
            )
    end
end

tspan = (0.0, 0.15)
p = (ps.nx, ps.deg + 1, ps.J, ps.ll, ps.lr, ps.dl, ps.dhl, ps.dhr, γ)
prob = ODEProblem(dudt!, u, tspan, p)
nt = tspan[2] ÷ dt |> Int
itg = init(prob, Midpoint(), saveat = tspan[2], adaptive = false, dt = dt)

@showprogress for iter = 1:nt
    step!(itg)
end

sol = zero(itg.u)
for i in axes(sol, 1), j in axes(sol, 2)
    sol[i, j, :] .= conserve_prim(itg.u[i, j, :], γ)
    sol[i, j, end] = 1 / sol[i, j, end]
end
plot(ps.x, sol[:, 2, :])


text:
module FESpacesTestsAll

using Test

@testset "FESpaces" begin include("FESpacesTests.jl") end

@testset "FESpaceConstructors" begin include("FESpaceConstructorsTests.jl") end

@testset "ConformingFESpaces" begin include("ConformingFESpacesTests.jl") end

@testset "CLagrangianFESpaces" begin include("CLagrangianFESpacesTests.jl") end

@testset "DLagrangianFESpaces" begin include("DLagrangianFESpacesTests.jl") end

@testset "DiscFESpaces" begin include("DiscFESpacesTests.jl") end

@testset "RTFESpacesTests" begin include("RTFESpacesTests.jl") end

@testset "NedelecFESpacesTests" begin include("NedelecFESpacesTests.jl") end

@testset "ConstrainedFESpaces" begin include("ConstrainedFESpacesTests.jl") end

@testset "ZeroMeanFESpaces" begin include("ZeroMeanFESpacesTests.jl") end

@testset "FEFunctions" begin include("FEFunctionsTests.jl") end

@testset "FEBases" begin include("FEBasesTests.jl") end

@testset "Assemblers" begin include("AssemblersTests.jl") end

@testset "FETerms" begin include("FETermsTests.jl") end

@testset "FEOperators" begin include("FEOperatorsTests.jl") end

@testset "DGFEOperators" begin include("DGFEOperatorsTests.jl") end

@testset "VectorValuedFEOperators" begin include("VectorValuedFEOperatorsTests.jl") end

@testset "NonLinearFEOperators" begin include("NonLinearFEOperatorsTests.jl") end

@testset "VectorValuedNonLinearFEOperators" begin include("VectorValuedNonLinearFEOperatorsTests.jl") end

end # module