text:
P_lwr_reflecting = zeros(size(am.P_lwr))
P_lwr_reflecting[:,rates(am).>0.0] .= 1.0
P_lwr_reflecting = P_lwr_reflecting./sum(P_lwr_reflecting,dims=2)
P_upr_reflecting = zeros(size(am.P_upr))
P_upr_reflecting[:,rates(am).<0.0] .= 1.0
P_upr_reflecting = P_upr_reflecting./sum(P_upr_reflecting,dims=2)
refl_model = BoundedFluidQueue(am.T,am.S,P_lwr_reflecting,P_upr_reflecting,nodes[end])

rng = StableRNGs.StableRNG(94)
P_lwr_rand = rand(rng,size(am.P_lwr)...)
P_lwr_rand[:,(rates(am).==0.0).&DiscretisedFluidQueues._has_right_boundary(am.S)] .= 0.0
P_lwr_rand = P_lwr_rand./sum(P_lwr_rand,dims=2)
P_upr_rand = rand(rng,size(am.P_upr)...)
P_upr_rand[:,(rates(am).==0.0).&DiscretisedFluidQueues._has_left_boundary(am.S)] .= 0.0
P_upr_rand = P_upr_rand./sum(P_upr_rand,dims=2)
rand_model = BoundedFluidQueue(am.T,am.S,P_lwr_rand,P_upr_rand,nodes[end])

@testset "Generators" begin
    @testset "augment_model" begin
        @testset "Lazy" begin
            for i in (:dgmesh,:frapmesh), tmp_model in (:am,:refl_model,:rand_model)
                msh = @eval $i
                (i==:dgmesh) && include("test_data/"*string(typeof(msh.nodes))[1:4]*"/"*string(tmp_model)*"/test_DG_B_data.jl")
                (i==:frapmesh) && include("test_data/"*string(typeof(msh.nodes))[1:4]*"/"*string(tmp_model)*"/test_FRAP_B_data.jl")
                println((tmp_model,i))
                dq = @eval DiscretisedFluidQueue($tmp_model,$i)
                B = build_lazy_generator(dq)
                Matrix(B)
                # types
                @test typeof(B)<:LazyGenerator
                idx = DiscretisedFluidQueues.qbd_idx(dq)
                @test all(isapprox.(B,B_data[idx,idx],atol=1e-4))
                # multiplcation (values)
                @test all(DiscretisedFluidQueues.fast_mul(B,Matrix{Float64}(I(size(B,1)))) .≈ B)
                @test all(DiscretisedFluidQueues.fast_mul(Matrix{Float64}(I(size(B,1))),B) .≈ B)
                # row sums
                @test sum(abs.(sum(B,dims=2)))≈0.0 atol=√eps()
                @test all(isapprox.(B*B,B_data[idx,idx]*B_data[idx,idx],atol=1e-3))
                # multiplication (types)
                @test typeof(DiscretisedFluidQueues.fast_mul(B,Matrix{Float64}(I(size(B,1)))))==Array{Float64,2}
                @test typeof(DiscretisedFluidQueues.fast_mul(Matrix{Float64}(I(size(B,1))),B))==Array{Float64,2}
                @test typeof(DiscretisedFluidQueues.fast_mul(B,SparseMatrixCSC{Float64,Int}(I(size(B,1)))))==SparseMatrixCSC{Float64,Int}
                @test typeof(DiscretisedFluidQueues.fast_mul(SparseMatrixCSC{Float64,Int}(I(size(B,1))),B))==SparseMatrixCSC{Float64,Int}
                # size
                @test size(B) == (40,40)
                @test size(B,1) == 40
                @test size(B,2) == 40
                # getindex
                @testset "getindex" begin
                    sz = size(B,1)
                    ind = true
                    for i in 1:sz, j in 1:sz
                        ei = zeros(1,sz)
                        ei[i] = 1
                        ej = zeros(sz)[:,:]
                        ej[j] = 1
                        !(B[i,j] ≈ (ei*B*ej)[1]) && (ind = false)
                    end
                    @test ind
                end
            end
        end

        @testset "Full" begin
            for i in (:dgmesh,:frapmesh,:fvmesh), tmp_model in (:am,:refl_model,:rand_model)
                println((tmp_model,i))
                msh = @eval $i
                (i==:dgmesh) && include("test_data/"*string(typeof(msh.nodes))[1:4]*"/"*string(tmp_model)*"/test_DG_B_data.jl")
                (i==:frapmesh) && include("test_data/"*string(typeof(msh.nodes))[1:4]*"/"*string(tmp_model)*"/test_FRAP_B_data.jl")
                ((i==:fvmesh)&&(tmp_model==:am)) && include("test_data/"*string(typeof(msh.nodes))[1:4]*"/"*string(tmp_model)*"/test_FV_B_data.jl")
                if (!(i==:fvmesh)||(tmp_model==:am))
                    @eval begin
                        dq = DiscretisedFluidQueue($tmp_model,$i)
                        B_Full = build_full_generator(dq)
                        if !(typeof($i)<:FVMesh)
                            B = build_lazy_generator(dq)
                            @test build_full_generator(B)==B_Full
                            @test all(isapprox.(B_Full*B_Full,B*B,atol=√eps()))
                            # size
                            @test size(B_Full) == (40,40)
                            @test size(B_Full,1) == 40
                            @test size(B_Full,2) == 40
                            #row sums
                            @test all(isapprox.(sum(B_Full,dims=2),0,atol=√eps()))
                        else
                            @test size(B_Full) == (16,16)
                            @test size(B_Full,1) == 16
                            @test size(B_Full,2) == 16
                        end
                        @test typeof(B_Full.B)==SparseMatrixCSC{Float64,Int}
                        # types
                        B_data = B_data[DiscretisedFluidQueues.qbd_idx(dq),DiscretisedFluidQueues.qbd_idx(dq)]
                        @test all(isapprox.(B_Full,B_data,atol=1e-4))
                        # multiplcation (values)
                        @test B_Full*SparseMatrixCSC{Float64,Int}(I(size(B_Full,1)))==B_Full
                        @test B_Full==SparseMatrixCSC{Float64,Int}(I(size(B_Full,1)))*B_Full
                        # row sums
                        @test all(isapprox.(B_Full*B_Full,B_data*B_data,atol=1e-3))
                    end
                end
            end
        end
    end
    @testset "normal model" begin
        @testset "Lazy" begin
            for i in (:dgmesh,:frapmesh)
                println(i)
                msh = @eval $i
                dq = @eval DiscretisedFluidQueue(model,$i)
                B = build_lazy_generator(dq)
                # types
                @test typeof(B)<:LazyGenerator
                # multiplcation (values)
                @test all(DiscretisedFluidQueues.fast_mul(B,Matrix{Float64}(I(size(B,1)))) .≈ B)
                @test all(DiscretisedFluidQueues.fast_mul(Matrix{Float64}(I(size(B,1))),B) .≈ B)
                # row sums
                @test all(isapprox.(sum(B,dims=2),0,atol=√eps()))
                # multiplication (types)
                @test typeof(DiscretisedFluidQueues.fast_mul(B,Matrix{Float64}(I(size(B,1)))))==Array{Float64,2}
                @test typeof(DiscretisedFluidQueues.fast_mul(Matrix{Float64}(I(size(B,1))),B))==Array{Float64,2}
                @test typeof(DiscretisedFluidQueues.fast_mul(B,SparseMatrixCSC{Float64,Int}(I(size(B,1)))))==SparseMatrixCSC{Float64,Int}
                @test typeof(DiscretisedFluidQueues.fast_mul(SparseMatrixCSC{Float64,Int}(I(size(B,1))),B))==SparseMatrixCSC{Float64,Int}
                # size
                @test size(B) == (31,31)
                @test size(B,1) == 31
                @test size(B,2) == 31
                # getindex
                @testset "getindex" begin
                    sz = size(B,1)
                    getindex_does_not_match_mul = true
                    for i in 1:sz, j in 1:sz
                        ei = zeros(1,sz)
                        ei[i] = 1
                        ej = zeros(sz)[:,:]
                        ej[j] = 1
                        !(B[i,j] .≈ (ei*B*ej)[1]) && (getindex_does_not_match_mul = false)
                    end
                    @test getindex_does_not_match_mul
                end
            end
        end

        @testset "Full" begin
            for i in (:dgmesh,:frapmesh,:fvmesh)
                println(i)
                msh = @eval $i
                @eval begin
                    dq = DiscretisedFluidQueue(model,$i)
                    B_Full = build_full_generator(dq)
                    if !(typeof($i)<:FVMesh)
                        B = build_lazy_generator(dq)
                        @test build_full_generator(B)==B_Full
                        @test all(isapprox.(B_Full*B_Full,B*B,atol=√eps()))
                        # size
                        @test size(B_Full) == (31,31)
                        @test size(B_Full,1) == 31
                        @test size(B_Full,2) == 31
                        #row sums
                        @test all(isapprox.(sum(B_Full,dims=2),0,atol=√eps()))
                    else
                        @test size(B_Full) == (13,13)
                        @test size(B_Full,1) == 13
                        @test size(B_Full,2) == 13
                    end
                    @test typeof(B_Full.B)==SparseMatrixCSC{Float64,Int}
                    # types
                    # multiplcation (values)
                    @test B_Full*SparseMatrixCSC{Float64,Int}(I(size(B_Full,1)))==B_Full
                    @test B_Full==SparseMatrixCSC{Float64,Int}(I(size(B_Full,1)))*B_Full
                end
            end
        end
    end
end
@testset "bnded qs default" begin
    P_lwr = zeros(sum(rates(model.S).<0.0),n_phases(model))
    P_lwr[end] = 1.0
    P_upr = zeros(sum(rates(model.S).>0.0),n_phases(model))
    P_upr[2] = 1.0
    model_bnd = BoundedFluidQueue(T,S,P_lwr,P_upr,nodes[end])

    dq = DiscretisedFluidQueue(model,dgmesh)
    dq_bnd = DiscretisedFluidQueue(model_bnd,dgmesh)
    B = build_lazy_generator(dq)
    B_bnd = build_lazy_generator(dq_bnd)
    @test B==B_bnd
    @test B_bnd*B_bnd≈B*B
    @test Matrix{Float64}(I(size(B_bnd,1)))*B_bnd≈Matrix{Float64}(I(size(B_bnd,1)))*B
end

max_stars_repo_name: angus-lewis/StochasticFluidQueues

max_issues_repo_name: angus-lewis/StochasticFluidQueues

max_forks_repo_name: angus-lewis/StochasticFluidQueues


text:
"""
    PlotStochasticProcesses

This module implements methods for plotting realizations of various stochastic processes.
"""
module PlotStochasticProcesses
using StochasticProcesses
using Plots
using GraphRecipes
using GraphPlot
using LightGraphs: AbstractGraph, nv, vertices
using Distributions: Sampleable
using StochasticProcesses.RandomWalks: RandomWalk
using StochasticProcesses.PoissonProcesses

export

# methods
plot

# include scripts
include("plotrandomwalks.jl")
include("plotpoissonprocesses.jl")
end

ext: jl

max_stars_repo_name: LoganDGraham/StochasticProcesses

max_issues_repo_name: LoganDGraham/StochasticProcesses

max_forks_repo_name: LoganDGraham/StochasticProcesses


text:
@defcomp TotalForcing begin
    f_CO2forcing = Parameter(index=[time], unit="W/m2")
    f_CH4forcing = Parameter(index=[time], unit="W/m2")
    f_N2Oforcing = Parameter(index=[time], unit="W/m2")
    f_lineargasforcing = Parameter(index=[time], unit="W/m2")
    exf_excessforcing = Parameter(index=[time], unit="W/m2")
    fs_sulfateforcing = Parameter(index=[time, region], unit="W/m2")

    area = Parameter(index=[region], unit="km2")
    area_e_eartharea = Parameter(unit="km2", default=5.1e8)

    # Total anthropogenic forcing
    ft_totalforcing = Variable(index=[time], unit="W/m2")
    fant_anthroforcing = Variable(index=[time], unit="W/m2")

    function run_timestep(p, v, d, tt)
        # From equation 16 of Hope (2006)
        v.ft_totalforcing[tt] = p.f_CO2forcing[tt] + p.f_CH4forcing[tt] + p.f_N2Oforcing[tt] + p.f_lineargasforcing[tt] + p.exf_excessforcing[tt]

        # Grand total forcing
        v.fant_anthroforcing[tt] = v.ft_totalforcing[tt] + sum(p.area .* p.fs_sulfateforcing[tt, :]) / p.area_e_eartharea
    end
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
"""
`IndividualsSelector` that implements a "trivial geography" similar to Spector and Kline (2006)
by first sampling an individual randomly and then selecting additional
individuals for the same tournament within a certain deme of limited size (`radius`)
for the sub-sequent individuals in the population.

The version we implement here is from:
    I. Harvey, "The Microbial Genetic Algorithm", in Advances in Artificial Life
    Darwin Meets von Neumann, Springer, 2011.
The original paper is:
    Spector, L., and J. Klein. 2005. Trivial Geography in Genetic Programming.
    In Genetic Programming Theory and Practice III, edited by T. Yu, R.L. Riolo,
    and B. Worzel, pp. 109-124. Boston, MA: Kluwer Academic Publishers.
    http://faculty.hampshire.edu/lspector/pubs/trivial-geography-toappear.pdf
"""
mutable struct RadiusLimitedSelector <: IndividualsSelector
    radius::Int
end

function select(sel::RadiusLimitedSelector, population, n::Integer)
    # The radius must be at least as big as the number of samples + 2 so that
    # there is something to sample from.
    radius = max(sel.radius, n+2)
    psize = popsize(population)
    deme_start = rand(1:psize)
    ixs = sample(deme_start:(deme_start+radius-1), n, replace=false, ordered=false)
    # Ensure they are not out of bounds by wrapping over at the end.
    ixs .= mod1.(ixs, psize)
end


text:
struct AWAP <: RasterDataSource end

layers(::Type{AWAP}) = (:solar, :rainfall, :vprpress09, :vprpress15, :tmin, :tmax)

@doc """
AWAP <: RasterDataSource

Daily weather data from the Australian Water Availability Project, developed by CSIRO.

See: [www.csiro.au/awap](http://www.csiro.au/awap/)

The available layers are: `$(layers(AWAP))`.
""" AWAP

const AWAP_PATHSEGMENTS = (
    solar = ("solar", "solarave", "daily"),
    rainfall = ("rainfall", "totals", "daily"),
    vprpress09 = ("vprp", "vprph09", "daily"),
    vprpress15 = ("vprp", "vprph15", "daily"),
    tmin = ("temperature", "minave", "daily"),
    tmax = ("temperature", "maxave", "daily"),
)
# Add ndvi monthly?  ndvi, ndviave, month

"""
    getraster(T::Type{AWAP}, layer::Symbol; date) => String
    getraster(T::Type{AWAP}, layer::Symbol, date) => String

Download data from the AWAP weather dataset, for `layer` in `$(layers(AWAP))`,
and `date` as a `DateTime` or iterable of `DateTime`.

AWAP is available on a daily timestep. If no `layer` is specified,
all layers will be downloaded, and a `Tuple` of `Vector{String}` will be returned.

## Example

Rainfall for the first month of 2001:

```julia
getraster(AWAP, :rainfall; date=Date(2001, 1, 1):Day(1):Date(2001, 1, 31))
```
"""
getraster(T::Type{AWAP}, layer::Symbol; date) = getraster(T, layer, date)
function getraster(T::Type{AWAP}, layer::Symbol, dates::Tuple)
    getraster(T, layer, _date_sequence(dates, Day(1)))
end
function getraster(T::Type{AWAP}, layer::Symbol, date::Dates.TimeType)
    _check_layer(T, layer)
    mkpath(_rasterpath(T, layer))
    raster_path = rasterpath(T, layer; date=date)
    if !isfile(raster_path)
        zip_path = zippath(T, layer; date=date)
        _maybe_download(zipurl(T, layer; date=date), zip_path)
        run(`uncompress $zip_path -f`)
    end
    return raster_path
end
function getraster(T::Type{AWAP}, layer::Symbol, dates::AbstractArray)
    getraster.(T, layer, dates)
end

rasterpath(T::Type{AWAP}) = joinpath(rasterpath(), "AWAP")
rasterpath(T::Type{AWAP}, layer; date::Dates.AbstractTime) =
    joinpath(_rasterpath(T, layer), rastername(T, layer; date))
_rasterpath(T::Type{AWAP}, layer) = joinpath(rasterpath(T), AWAP_PATHSEGMENTS[layer][1:2]...)
rastername(T::Type{AWAP}, layer; date::Dates.AbstractTime) =
    joinpath(_date2string(T, date) * ".grid")

function zipurl(T::Type{AWAP}, layer; date)
    s = AWAP_PATHSEGMENTS[layer]
    d = _date2string(T, date)
    # The actual zip name has the date twice, which is weird.
    # So we getraster in to a different name as there no output
    # name flages for `uncompress`. It's ancient.
    uri = URI(scheme="http", host="www.bom.gov.au", path="/web03/ncc/www/awap")
    joinpath(uri, s..., "grid/0.05/history/nat/$d$d.grid.Z")
end
zipname(T::Type{AWAP}, layer; date) = _date2string(T, date) * ".grid.Z"
zippath(T::Type{AWAP}, layer; date) =
    joinpath(_rasterpath(T, layer), zipname(T, layer; date))

_dateformat(::Type{AWAP}) = DateFormat("yyyymmdd")


text:
# This file holds the definition of the Retirement type. This type is used to
#   describe the particulars of retirement due to age/end of career.

# The Retirement type requires no additional types.
requiredTypes = [ ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

export Retirement
"""
This type defines a retirement scheme for mandatory retirement of personnel
members based on age or tenure.

The type contains the following fields:
* `maxCareerLength::Float64`: the maximum length of a personnel member's career.
  Set this to 0 to ignore this criterion.
* `retireAge::Float64`: the mandatory retirement age. Set this to 0 to ignore
  this criterion.
* `retireFreq::Float64`: the length of the time interval between two retirement
  checks.
* `retireOffset::Float64`: the offset of the retirement schedule with respect to
  the start of the simulation.
* 'isEither::Bool': a flag indicating whether either of the criteria (age or
  tenure) must be satisfied for retirement, or both.
"""
type Retirement

    maxCareerLength::Float64
    retireAge::Float64
    retireFreq::Float64
    retireOffset::Float64
    isEither::Bool

    function Retirement( ; freq::T1 = 1.0, offset::T2 = 0.0,
        maxCareer::T3 = 0.0, retireAge::T4 = 0.0, isEither::Bool = true ) where T1 <: Real where T2 <: Real where T3 <: Real where T4 <: Real

        if freq <= 0.0
            error( "Retirement cycle length must be > 0.0." )
        end  # if freq < 0.0

        if maxCareer < 0.0
            error( "Maximal career length must be ⩾ 0.0." )
        end  # if maxCareer < 0.0

        if retireAge < 0.0
            error( "Mandatory retirement age must be ⩾ 0.0." )
        end  # if maxCareer < 0.0

        newRet = new()
        newRet.maxCareerLength = maxCareer
        newRet.retireAge = retireAge
        newRet.retireFreq = freq
        newRet.retireOffset = freq > 0.0 ?
            ( offset % freq + ( offset < 0.0 ? freq : 0.0 ) ) : 0.0
        newRet.isEither = isEither
        return newRet

    end  # Retirement( freq, offset, maxCareer, retireAge )

end  # type Retirement


text:
module Ripser
export ripser

using SparseArrays: AbstractSparseMatrix, findnz

const depsfile = joinpath(dirname(@__DIR__), "deps", "deps.jl")
if isfile(depsfile)
    include("../deps/deps.jl")
else
    error("Ripser in not properly installed. Please run Pkg.build(\"Ripser\").")
end
__init__() = check_deps()

# The value_t type from ripser source code.
const Cvalue_t = Cfloat

# RawResult contains pointers to all output arrays required by ripser.
struct RawResult{T}
    dim_max        ::Int
    n_intervals    ::Ref{Ptr{Cint}}
    births_deaths  ::Ref{Ptr{Cvalue_t}}
    cocycle_length ::Ref{Ptr{Cint}}
    cocycles       ::Ref{Ptr{Cint}}
end

RawResult{T}(dim_max) where T =
    RawResult{T}(dim_max,
                 Ref{Ptr{Cint}}(),
                 Ref{Ptr{Cvalue_t}}(),
                 Ref{Ptr{Cint}}(),
                 Ref{Ptr{Cint}}())

# Converts a Matrix to Vector{Cvalue_t} since ripser expects a flat array as input.
function flatten_distmat(dists)
    dists_flat = Cvalue_t[]
    for i in 1:size(dists, 1)-1, j in i+1:size(dists, 1)
        push!(dists_flat, Cvalue_t(dists[j, i]))
    end
    dists_flat
end

function isprime(n)
    if iseven(n) || n < 2
        n == 2
    else
        p = 3
        q = n / p
        while p ≤ q
            iszero(n % p) && return false
            p += 2
            q = n / p
        end
        true
    end
end

function check_args(dists, modulus, dim_max, threshold)
    size(dists, 1) == size(dists, 2) || throw(ArgumentError("distance matrix must be square"))
    isprime(modulus) || throw(ArgumentError("modulus must be a prime number"))
    dim_max ≥ 0      || throw(ArgumentError("dim_max must be non-negative"))
    threshold > 0    || throw(ArgumentError("threshold must be positive"))
end

function split_cocycle(cocycle, dim)
    c = Tuple{Vector{Int}, Int}[]
    for i in 1:dim+2:length(cocycle)-dim-1
        push!(c, (cocycle[i:i+dim] .+ 1, cocycle[i+dim+1]))
    end
    c
end

# Unpack RawResult{T} to barcode and cocycles (if return_cocycles is true)
function unpack_results(raw::RawResult{T}, return_cocycles) where T
    dim_max = raw.dim_max

    n_intervals = unsafe_wrap(Vector{Cint}, raw.n_intervals[], raw.dim_max + 1,
                              own = true)
    intervals = unsafe_wrap(Matrix{Cvalue_t}, raw.births_deaths[], (2, sum(n_intervals)),
                            own = true)
    cocycle_length = unsafe_wrap(Vector{Cint}, raw.cocycle_length[], sum(n_intervals),
                                 own = true)
    if sum(cocycle_length) > 0
        cocycles_flat = unsafe_wrap(Vector{Cint}, raw.cocycles[], sum(cocycle_length),
                                    own = true)
    else
        cocycles_flat = Cint[]
    end

    if !return_cocycles
        barcodes = Matrix{T}[]
        start = 0
        for int in n_intervals
            push!(barcodes, T.(intervals[:, start+1:start+int]))
            start += int
        end

        map(barcodes) do bc
            collect(vec(reinterpret(Tuple{T, T}, bc)))
        end
    else
        barcodes = Matrix{T}[]
        cocycles = Vector{Vector{Tuple{Vector{Int}, Int}}}[]

        start_bc = 0
        start_cc = 0
        for int in n_intervals
            push!(barcodes, T.(intervals[:, start_bc+1:start_bc+int]))
            push!(cocycles, Int[])
            for i in 1:int
                len = cocycle_length[start_bc + i]
                push!(cocycles[end],
                      split_cocycle(cocycles_flat[start_cc+1:start_cc+len],
                                    length(cocycles)-1))
                start_cc += len
            end
            start_bc += int
        end

        map(barcodes) do bc
            collect(vec(reinterpret(Tuple{T, T}, bc)))
        end, cocycles
    end
end

"""
    ripser(dists; modulus = 2, dim_max = 1, threshold = Inf, cocycles = false)
"""
function ripser(dists     ::AbstractMatrix{T};
                modulus   ::Integer = 2,
                dim_max   ::Integer = 1,
                threshold ::Real = Inf,
                cocycles  ::Bool = false) where T<:AbstractFloat
    check_args(dists, modulus, dim_max, threshold)

    res = RawResult{T}(dim_max)
    dists_flat = flatten_distmat(dists)

    ripser_fptr = Libdl.dlsym(Libdl.dlopen(libripser), :c_rips_dm)
    n_edges = ccall(ripser_fptr,
                    Cint,
                    (Ptr{Ptr{Cint}}, Ptr{Ptr{Cvalue_t}}, Ptr{Ptr{Cint}}, Ptr{Ptr{Cint}},
                     Ptr{Cvalue_t}, Cint,
                     Cint, Cint, Cvalue_t, Cint),
                    res.n_intervals, res.births_deaths, res.cocycle_length, res.cocycles,
                    dists_flat, length(dists_flat),
                    modulus, dim_max, threshold, cocycles)

    unpack_results(res, cocycles)
end

function ripser(dists     ::AbstractSparseMatrix{T};
                modulus   ::Integer = 2,
                dim_max   ::Integer = 1,
                threshold ::Real = Inf,
                cocycles  ::Bool = false) where T<:AbstractFloat
    check_args(dists, modulus, dim_max, threshold)

    J, I, V = findnz(dists)
    I .-= 1
    J .-= 1
    res = RawResult{T}(dim_max)

    ripser_fptr = Libdl.dlsym(Libdl.dlopen(libripser), :c_rips_dm_sparse)
    n_edges = ccall(ripser_fptr,
                    Cint,
                    (Ptr{Ptr{Cint}}, Ptr{Ptr{Cvalue_t}}, Ptr{Ptr{Cint}}, Ptr{Ptr{Cint}},
                     Ptr{Cint}, Ptr{Cint}, Ptr{Cvalue_t}, Cint, Cint,
                     Cint, Cint, Cvalue_t, Cint),
                    res.n_intervals, res.births_deaths, res.cocycle_length, res.cocycles,
                    Cint.(I), Cint.(J), Cvalue_t.(V), length(I), size(dists, 1),
                    modulus, dim_max, threshold, cocycles)

    unpack_results(res, cocycles)
end
end


text:
@testset "Relativistic SGNHT" begin

    #----------------------
    # Gaussian sanity check.
    #----------------------
    dm = GaussianMixtureModel([0.],[1.],[1.])

    gradient = DataModel.getgrad(dm)
    s = RelSGNHTState(zeros(1),0.1)

    x_avg = [0.]
    x2_avg = [0.]
    for i in 1:50000
        sample!(s,gradient)
        x_avg =(i-1)/i*x_avg + s.x/i
        x2_avg = (i-1)/i*x2_avg + s.x.^2/i
    end

    @test abs(x_avg-0)[1] < 0.1
    @test abs(x2_avg-1)[1] < 0.2
end


text:
export PoissonTTSolver,solve!
using LinearAlgebra

struct PoissonTTSolver
    Mx::Array{Float64,2}
    My::Array{Float64,2}
    Mxm1::Array{Float64,2}
    Mym1::Array{Float64,2}
    Dx::Array{Float64,1}
    Dy::Array{Float64,1}
    t1::Array{Float64,2}
    t2::Array{Float64,2}

    function PoissonTTSolver(lp2D)
        # Lx=l1dpαI_dirichlet(nx,dx,0.0)
        Lx=SymTridiagonal(lp2D.lpx)
        nx=size(Lx,1)
        Ly=SymTridiagonal(lp2D.lpy)
        ny=size(Ly,1)

        # Ly=l1dpαI_dirichlet(ny,dy,0.0)

        Ex=eigen(Lx)
        Mx=Ex.vectors
        Dx=Ex.values
        Mxm1=inv(Mx)

        Ey=eigen(Ly)
        My=Ey.vectors
        Dy=Ey.values
        Mym1=inv(My)
        t1=zeros(nx,ny)
        t2=zeros(nx,ny)
        new(Mx,My,Mxm1,Mym1,Dx,Dy,t1,t2)
    end
end

function solve!(Xxy,Bxy,ps::PoissonTTSolver)
    (Mx,My,Mxm1,Mym1,Dx,Dy,t1,t2)=(ps.Mx,ps.My,ps.Mxm1,ps.Mym1,ps.Dx,ps.Dy,ps.t1,ps.t2)
    (nx,ny)=size(Xxy)
    @assert size(Xxy)==size(Bxy)==(size(Mx,1),size(My,2))

    LinearAlgebra.mul!(t1,Mym1,Bxy)
    LinearAlgebra.mul!(t2,Mxm1,transpose(t1))

    for j=1:ny
        @simd for i=1:nx
            @inbounds t2[i,j]/=(Dx[i]+Dy[j])
        end
    end

    LinearAlgebra.mul!(t1,My,t2)
    LinearAlgebra.mul!(Xxy,Mx,transpose(t1))
    nothing
    # Xxy
end


text:
module Boltzmann

export BernoulliRBM,
       GRBM,
       fit,
       transform,
       generate,
       components

include("rbm.jl")

end


text:
# Autogenerated wrapper script for vtkfig_jll for x86_64-linux-gnu-cxx11
export libvtkfig

using VTKMinimalQT_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `libvtkfig`
const libvtkfig_splitpath = ["lib", "libvtkfig.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libvtkfig_path = ""

# libvtkfig-specific global declaration
# This will be filled out by __init__()
libvtkfig_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libvtkfig = "libvtkfig.so"

# Inform that the wrapper is available for this platform
wrapper_available = true

"""
Open all libraries
"""
function __init__()
    # This either calls `@artifact_str()`, or returns a constant string if we're overridden.
    global artifact_dir = find_artifact_dir()

    global PATH_list, LIBPATH_list
    # Initialize PATH and LIBPATH environment variable listings
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (VTKMinimalQT_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (VTKMinimalQT_jll.LIBPATH_list,))

    global libvtkfig_path = normpath(joinpath(artifact_dir, libvtkfig_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libvtkfig_handle = dlopen(libvtkfig_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libvtkfig_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')

end  # __init__()


text:
# computes overapproximation of Φ * set(F[k-1]) with a zonotope
# this operations adds n generators, hence we use an order reduction
# function
function reach_homog_ASB07!(F::Vector{ReachSet{N, Zonotope{N, VN, MN}}},
                            Ω0::Zonotope{N, VN, MN},
                            Φ::AbstractMatrix,
                            NSTEPS::Integer,
                            δ::N,
                            max_order::Integer,
                            X::Universe,
                            recursive::Val{true},
                            reduction_method::AbstractReductionMethod,
                            Δt0::TN) where {N, TN, VN, MN}
    # initial reach set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = ReachSet(Ω0, Δt)

    # split the interval matrix into center and radius
    Φc, Φs = _split(Φ)

    k = 1
    @inbounds while k <= NSTEPS - 1
        Zk = set(F[k])
        ck = Zk.center
        Gk = Zk.generators

        Zₖ₊₁ = _overapproximate_interval_linear_map(Φc, Φs, ck, Gk)
        Zₖ₊₁ʳ = _reduce_order(Zₖ₊₁, max_order, reduction_method)

        k += 1
        Δt += δ
        F[k] = ReachSet(Zₖ₊₁ʳ, Δt)
    end
    return F
end

# non-recursive implementation; to get more accurate interval matrix powers Φ^k
# we use the IntervalMatrices.IntervalMatrixPower interface
function reach_homog_ASB07!(F::Vector{ReachSet{N, Zonotope{N, VN, MN}}},
                            Ω0::Zonotope{N, VN, MN},
                            Φ::AbstractMatrix,
                            NSTEPS::Integer,
                            δ::N,
                            max_order::Integer,
                            X::Universe,
                            recursive::Val{false},
                            reduction_method::AbstractReductionMethod,
                            Δt0::TN) where {N, TN, VN, MN}
    # initial reach set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = ReachSet(Ω0, Δt)
    Z0 = Ω0
    c0 = Z0.center
    G0 = Z0.generators

    Φpow = IntervalMatrixPower(Φ) # lazy interval matrix power

    k = 2
    @inbounds while k <= NSTEPS
        Φ_power_k = get(Φpow)
        Φc, Φs = _split(Φ_power_k)

        Zₖ = _overapproximate_interval_linear_map(Φc, Φs, c0, G0)
        Zₖʳ = _reduce_order(Zₖ, max_order, reduction_method)

        Δt += δ
        F[k] = ReachSet(Zₖʳ, Δt)
        increment!(Φpow)
        k += 1
    end
    return F
end

# case with an invariant
function reach_homog_ASB07!(F::Vector{ReachSet{N, Zonotope{N, VN, MN}}},
                            Ω0::Zonotope{N, VN, MN},
                            Φ::AbstractMatrix,
                            NSTEPS::Integer,
                            δ::N,
                            max_order::Integer,
                            X::LazySet,
                            recursive::Val{true},
                            reduction_method::AbstractReductionMethod,
                            Δt0::TN) where {N, TN, VN, MN}
    # initial reach set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = ReachSet(Ω0, Δt)

    # split the interval matrix into center and radius
    Φc, Φs = _split(Φ)

    k = 2
    @inbounds while k <= NSTEPS
        Zₖ₋₁ = set(F[k-1])
        cₖ₋₁ = Zₖ₋₁.center
        Gₖ₋₁ = Zₖ₋₁.generators

        Zₖ = _overapproximate_interval_linear_map(Φc, Φs, cₖ₋₁, Gₖ₋₁)
        Zₖʳ = _reduce_order(Zₖ, max_order, reduction_method)
        _is_intersection_empty(X, Zₖʳ) && break

        k += 1
        Δt += δ
        F[k] = ReachSet(Zₖʳ, Δt)
    end
    if k < NSTEPS + 1
        resize!(F, k-1)
    end
    return F
end


text:
struct IntervalLabels <: AbstractData end

Base.@kwdef mutable struct LabelledInterval <: AbstractModel
    label::String = ""
    nested_interval::NestedInterval
end

Base.@kwdef mutable struct LabelledIntervals{T₁ <: Dict{String, LabelledInterval}}  <: AbstractModel
    description::String = ""
    labelled_intervals::T₁
end

function get_description(lis::LabelledIntervals)
    lis.description
end

function set_description!(lis::LabelledIntervals, description::String)
    lis.description = description
end

function get_labelled_intervals(lis::LabelledIntervals)
    lis.labelled_intervals
end

function set_labelled_intervals!(lis::LabelledIntervals, labelled_intervals::Dict{String, LabelledInterval})
    lis.labelled_intervals = labelled_intervals
end

function get_label(li::LabelledInterval)
    li.label
end

function set_label!(li::LabelledInterval, label)
    li.label = label
end

function get_nested_interval(li::LabelledInterval)
    li.nested_interval
end

function set_nested_interval!(li::LabelledInterval, nested_interval)
    li.nested_interval = nested_interval
end


text:
using DiffEqBayesStan
using OrdinaryDiffEq
using ParameterizedFunctions
using RecursiveArrayTools
using StatsPlots, Statistics
using StanSample

ProjDir = @__DIR__

println("Four parameter case (using StanLanguage)")

f1 = @ode_def begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d
u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob1 = ODEProblem(f1,u0,tspan,p)
sol = solve(prob1,Tsit5())
t = collect(range(1,stop=10,length=50))
randomized = VectorOfArray([(sol(t[i]) + .3randn(2)) for i in 1:length(t)])
data_array = convert(Array, randomized)'

stan_pes = "
functions {
  vector sho(real t,vector internal_var___u,
    real internal_var___p_1, real internal_var___p_2,
    real internal_var___p_3, real internal_var___p_4) {

    vector[2] internal_var___du;
    internal_var___du[1] = internal_var___p_1 * internal_var___u[1] -
      internal_var___p_2 * internal_var___u[1] * internal_var___u[2];
    internal_var___du[2] = internal_var___u[2] * -internal_var___p_3 +
      internal_var___p_4 * internal_var___u[1] * internal_var___u[2];

    return internal_var___du;
  }
}
data {
  vector[2] u0;
  int<lower=1> T;
  real internal_var___u[T,2];
  real t0;
  real ts[T];
}
parameters {
  row_vector<lower=0>[2] sigma1;
  real<lower=0.5,upper=2.0> theta_1;real<lower=0.5,upper=1.5> theta_2;real<lower=0.5,upper=4.0> theta_3;real<lower=0.5,upper=2.0> theta_4;
}
model{
  vector[2] u_hat[T];
  sigma1 ~ inv_gamma(4.0, 1.0);
  theta_1 ~ normal(1.1, 0.5) T[0.5,2.0];
  theta_2 ~ normal(1.5, 0.5) T[0.5,1.5];
  theta_3 ~ normal(2.0, 0.5) T[0.5,4.0];
  theta_4 ~ normal(1.5, 0.5) T[0.5,2.0];
  u_hat = ode_rk45_tol(sho, u0, t0, ts, 0.001, 1.0e-6, 100000,
    theta_1,theta_2,theta_3,theta_4);
  for (t in 1:T){
    internal_var___u[t,:] ~ normal(u_hat[t,1:2],sigma1);
    }
}";

pes4 = SampleModel("pes4", stan_pes, joinpath(ProjDir, "tmp"))
data = (u0=u0, T=size(t, 1), internal_var___u=data_array, t0=0.0, ts=t)
rc = stan_sample(pes4; data)

if success(rc)
  df = read_samples(pes4, :dataframe)

  println("\ntheta means:\n",
    [mean(df.theta_1), mean(df.theta_2), mean(df.theta_3), mean(df.theta_4)])
  println("\ntheta vars:\n",
    [var(df.theta_1), var(df.theta_2), var(df.theta_3), var(df.theta_4)])
  println()

  dfs = read_samples(pes4, :dataframes)
  for j in 1:4
    for i in 1:4
      if i == 1
        plot(dfs[i][:, "theta_$j"], lab="chain[$j]")
      else
        plot!(dfs[i][:, "theta_$j"], lab="chain[$j]")
      end
    end
    savefig(joinpath(ProjDir, "tmp", "theta_$j.png"))
  end

end


text:
using Revise
using AutomotiveDrivingModels
using Parameters
using AutoViz

includet("environment.jl")

env = MergingEnvironment(main_lane_angle = 0.0, merge_lane_angle = pi/12)

AutoViz.render(Scene(),env.roadway, [LaneOverlay(main_lane(env), colorant"blue")])

struct LaneOverlay <: SceneOverlay
    lane::Lane
    color::Colorant
end
function AutoViz.render!(rendermodel::RenderModel, overlay::LaneOverlay, scene::Scene, roadway::Roadway)
    render!(rendermodel, overlay.lane, roadway, color_asphalt=overlay.color)
    return rendermodel
end

# test modulo roadway

state1 = VehicleState(VecSE2(45.0, 0.0, 0.0), env.roadway, 10.0)
veh1 = Vehicle(state1, VehicleDef(), 1)

s = Scene()
push!(s, veh1)

AutoViz.render(s, env.roadway, cam=FitToContentCamera(0.0))

veh1 = Vehicle(propagate(veh1, LaneFollowingAccel(0.0), env.roadway, 1.0), VehicleDef(), 1)

s = Scene()
push!(s, veh1)
AutoViz.render(s, env.roadway, cam=FitToContentCamera(0.0))


text:
import Adapt

using Oceananigans: short_show
using Oceananigans.Utils: user_function_arguments
using Oceananigans.Operators: assumed_field_location, index_and_interp_dependencies
using Oceananigans.Fields: show_location
using Oceananigans.Utils: tupleit

"""
    ContinuousForcing{X, Y, Z, P, F, D, I}

A callable object that implements a "continuous form" forcing function
on a field at the location `X, Y, Z` with optional parameters.
"""
struct ContinuousForcing{X, Y, Z, P, F, D, I, ℑ}
                          func :: F
                    parameters :: P
            field_dependencies :: D
    field_dependencies_indices :: I
     field_dependencies_interp :: ℑ

    # Non-public "temporary" constructor that stores func, parameters, and field_dependencies
    # for later regularization
    function ContinuousForcing(func, parameters, field_dependencies)
        field_dependencies = tupleit(field_dependencies)

        return new{Nothing, Nothing, Nothing,
                   typeof(parameters),
                   typeof(func),
                   typeof(field_dependencies),
                   Nothing,
                   Nothing}(func, parameters, field_dependencies, nothing, nothing)
    end

    # Non-public "final" constructor.
    function ContinuousForcing{X, Y, Z}(func, parameters=nothing, field_dependencies=(),
                                        field_dependencies_indices=(), field_dependencies_interp=()) where {X, Y, Z}
        return new{X, Y, Z,
                   typeof(parameters),
                   typeof(func),
                   typeof(field_dependencies),
                   typeof(field_dependencies_indices),
                   typeof(field_dependencies_interp)}(func, parameters, field_dependencies,
                                                      field_dependencies_indices, field_dependencies_interp)
    end
end

"""
    ContinuousForcing(func; parameters=nothing, field_dependencies=())

Construct a "continuous form" forcing with optional `parameters` and optional
`field_dependencies` on other fields in a model.

If neither `parameters` nor `field_dependencies` are provided, then `func` must be
callable with the signature

    `func(x, y, z, t)`

where `x, y, z` are the east-west, north-south, and vertical spatial coordinates, and `t` is time.

If `field_dependencies` are provided, the signature of `func` must include them.
For example, if `field_dependencies=(:u, :S)` (and `parameters` are _not_ provided), then
`func` must be callable with the signature

    `func(x, y, z, t, u, S)`

where `u` is assumed to be the `u`-velocity component, and `S` is a tracer. Note that any field
which does not have the name `u`, `v`, or `w` is assumed to be a tracer and must be present
in `model.tracers`.

If `parameters` are provided, then the _last_ argument to `func` must be `parameters`.
For example, if `func` has no `field_dependencies` but does depend on `parameters`, then
it must be callable with the signature

    `func(x, y, z, t, parameters)`

With `field_dependencies=(:u, :v, :w, :c)` and `parameters`, then `func` must be
callable with the signature

    `func(x, y, z, t, u, v, w, c, parameters)`

"""
ContinuousForcing(func; parameters=nothing, field_dependencies=()) =
    ContinuousForcing(func, parameters, field_dependencies)

"""
    regularize_forcing(forcing::ContinuousForcing, field, field_name, model_field_names)

Regularize `forcing::ContinuousForcing` by determining the indices of `forcing.field_dependencies`
in `model_field_names`, and associated interpolation functions so `forcing` can be used during
time-stepping `IncompressibleModel`.
"""
function regularize_forcing(forcing::ContinuousForcing, field, field_name, model_field_names)

    X, Y, Z = location(field)

    indices, interps = index_and_interp_dependencies(X, Y, Z,
                                                     forcing.field_dependencies,
                                                     model_field_names)

    return ContinuousForcing{X, Y, Z}(forcing.func, forcing.parameters, forcing.field_dependencies,
                                      indices, interps)
end

#####
##### Functions for calling ContinuousForcing in a time-stepping kernel
#####

@inline function (forcing::ContinuousForcing{X, Y, Z, F})(i, j, k, grid, clock, model_fields) where {X, Y, Z, F}

    args = user_function_arguments(i, j, k, grid, model_fields, forcing.parameters, forcing)

    return @inbounds forcing.func(xnode(X, i, grid), ynode(Y, j, grid), znode(Z, k, grid), clock.time, args...)
end

"""Show the innards of a `ContinuousForcing` in the REPL."""
Base.show(io::IO, forcing::ContinuousForcing{X, Y, Z, P}) where {X, Y, Z, P} =
    print(io, "ContinuousForcing{$P} at ", show_location(X, Y, Z), '\n',
        "├── func: $(short_show(forcing.func))", '\n',
        "├── parameters: $(forcing.parameters)", '\n',
        "└── field dependencies: $(forcing.field_dependencies)")

"""Show the innards of an "non-regularized" `ContinuousForcing` in the REPL."""
Base.show(io::IO, forcing::ContinuousForcing{Nothing, Nothing, Nothing, P}) where P =
    print(io, "ContinuousForcing{$P}", '\n',
        "├── func: $(short_show(forcing.func))", '\n',
        "├── parameters: $(forcing.parameters)", '\n',
        "└── field dependencies: $(forcing.field_dependencies)")

Adapt.adapt_structure(to, forcing::ContinuousForcing{X, Y, Z}) where {X, Y, Z} =
    ContinuousForcing{X, Y, Z}(Adapt.adapt(to, forcing.func),
                               Adapt.adapt(to, forcing.parameters),
                               nothing,
                               Adapt.adapt(to, forcing.field_dependencies_indices),
                               Adapt.adapt(to, forcing.field_dependencies_interp))


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""
    AcceleratedParametricSurface{T,N,S} <: ParametricSurface{T,N}

Wrapper class for [`ParametricSurface`](@ref)s where analytical intersection isn't feasible (e.g. [`ZernikeSurface`](@ref), [`ChebyshevSurface`](@ref)).
The surface is instead triangulated and an iterative (newton raphson) process carried out to determine precise ray intersection points.
`S` is the type of the ParametricSurface being wrapped.

```julia
AcceleratedParametricSurface(surf::ParametricSurface{T,N}, numsamples::Int = 17; interface::NullOrFresnel{T} = nullinterface(T))
```
"""
struct AcceleratedParametricSurface{T,N,S<:ParametricSurface{T,N}} <: ParametricSurface{T,N}
    surface::S
    triangles::Vector{Triangle{T}}
    sidelengths::Vector{T}
    triangles_bbox::BoundingBox{T}
    interface::NullOrFresnel{T}

    function AcceleratedParametricSurface(surf::S, numsamples::Int = 17; interface::NullOrFresnel{T} = NullInterface(T)) where {T<:Real,N,S<:ParametricSurface{T,N}}
        a = AcceleratedParametricSurface(surf, triangulate(surf, numsamples, true, true), interface = interface)
        emptytrianglepool!(T)
        return a
    end

    function AcceleratedParametricSurface(surf::S, triangles::Vector{Triangle{T}}; interface::NullOrFresnel{T} = NullInterface(T)) where {T<:Real,N,S<:ParametricSurface{T,N}}
        sidelengths = Vector{T}(undef, length(triangles))
        @inbounds @simd for i in 1:length(triangles)
            t = triangles[i]
            sidelengths[i] = max(norm(t.BA), norm(t.CA))
        end
        return new{T,N,S}(surf, copy(triangles), sidelengths, BoundingBox(triangles), interface)
    end
end
export AcceleratedParametricSurface

uvrange(::Type{AcceleratedParametricSurface{T,N,S}}) where {T<:Real,N,S<:ParametricSurface{T,N}} = uvrange(S)
uvrange(surf::AcceleratedParametricSurface{T,N,S}) where {T<:Real,N,S<:ParametricSurface{T,N}} = uvrange(surf.surface)
point(surf::AcceleratedParametricSurface{T,N,S}, u::T, v::T) where {T<:Real,N,S<:ParametricSurface{T,N}} = point(surf.surface, u, v)
partials(surf::AcceleratedParametricSurface{T,N,S}, u::T, v::T) where {T<:Real,N,S<:ParametricSurface{T,N}} = partials(surf.surface, u, v)
normal(surf::AcceleratedParametricSurface{T,N,S}, u::T, v::T) where {T<:Real,N,S<:ParametricSurface{T,N}} = normal(surf.surface, u, v)
inside(surf::AcceleratedParametricSurface{T,N,S}, x::T, y::T, z::T) where {T<:Real,N,S<:ParametricSurface{T,N}} = inside(surf.surface, x, y, z)
onsurface(surf::AcceleratedParametricSurface{T,N,S}, x::T, y::T, z::T) where {T<:Real,N,S<:ParametricSurface{T,N}} = onsurface(surf.surface, x, y, z)
"""
    interface(surf::Surface{T}) -> OpticalInterface{T}

Return the [`OpticalInterface`](@ref) associated with `surf`.
"""
interface(a::AcceleratedParametricSurface{T}) where {T<:Real} = a.interface
uv(surf::AcceleratedParametricSurface{T,3}, x::T, y::T, z::T) where {T<:Real} = uv(surf.surface, x, y, z)

function BoundingBox(surf::AcceleratedParametricSurface{T,3}) where {T<:Real}
    # the stored bbox is a tight one around the triangulated surface only, in reality the half-space created is an infinite prism capped by the triangulated surface
    # this bounding box is not accurate to the actual surface, but it is accurate to the triangulated surface
    # if a ray hits the accurate bounding box but not the triangulated bounding box then it would miss the triangulated surface
    # anyway (erroneously) and not return an intersection - so we don't lose anything by making this approximation here
    tbox = surf.triangles_bbox
    return BoundingBox(tbox.xmin, tbox.xmax, tbox.ymin, tbox.ymax, typemin(T), tbox.zmax)
end

# function makemesh(surface::AcceleratedParametricSurface{S,N,T}, ::Int) where {S<:Real,N,T<:ParametricSurface{S,N}}
#     # If we've already triangulated the surface then use this directly as it will be the most
#     # accurate representation of the surface being intersected with
#     return TriangleMesh(surface.triangles)
# end

# "computes all the intersections of the ray with the surface. Each intersection is represented as an Interval that encloses the interior of the object. If the ray is entering the surface, dot(ray,normal) < 0, then the interval is (intersection,infinity) where α is the parametric value of the line surface intersection. If the ray is leaving the surface, dot(ray,normal) > 0, then the interval is (rayorigin, α). If there are two intersections with the surface then returns an Interval(intersection,intersection). If there are multiple segments enclosing the object then will return a list that will contain some combination of the above types of Interval."
# fallback method for any accelerated surface which doesn't have an overridden method

"""
    surfaceintersection(surf::Surface{T}, r::AbstractRay{T}) where {T}

Calculates the intersection of `r` with a surface of any type, `surf`.
Note that some surfaces cannot be intersected analytically so must be wrapped in an [`AcceleratedParametricSurface`](@ref) in order to be intersected.

Returns an [`EmptyInterval`](@ref) if there is no [`Intersection`](@ref), an [`Interval`](@ref) if there is one or two intersections and a [`DisjointUnion`](@ref) if there are more than two intersections.
"""
function surfaceintersection(surf::AcceleratedParametricSurface{T,N}, r::AbstractRay{T,N}) where {T,N}
    if doesintersect(surf.triangles_bbox, r)
        i = triangulatedintersection(surf, r)
        if i isa EmptyInterval{T}
            if inside(surf, origin(r))
                return rayorigininterval(Infinity(T))
            else
                return EmptyInterval(T)
            end
        else
            return i
        end
    else
        if inside(surf, origin(r))
            return rayorigininterval(Infinity(T))
        else
            return EmptyInterval(T)
        end
    end
end

"""
    triangulatedintersection(surf::AcceleratedParametricSurface{T,N,S}, r::AbstractRay{T,N})

Intersection of a ray, `r`, with a triangulated surface, `surf`, no concept of inside so never returns a [`RayOrigin`](@ref) [`Interval`](@ref).
"""
function triangulatedintersection(surf::AcceleratedParametricSurface{T,N,S}, r::AbstractRay{T,N}) where {T<:Real,N,S}
    result = newinintervalpool!(T)
    offset_ray = Ray(origin(r) - ACCEL_SURF_RAY_OFFSET * direction(r), direction(r))
    for (i, tri) in enumerate(surf.triangles)
        # test if the ray is within the maximum side length from one of the points, if not then skip it
        w = vertex(tri, 1) - origin(r)
        dist = norm(w - dot(w, direction(r)) * direction(r))
        if dist > surf.sidelengths[i]
            continue
        end
        # find the triangle intersection - it's quite possible that for rays starting near the surface we will we miss
        # an intersection because of imprecise triangulation, so introduce an offset
        intv = surfaceintersection(tri, offset_ray)
        if !(intv isa EmptyInterval{T})
            ints = halfspaceintersection(intv)
            if α(ints) <= ACCEL_SURF_RAY_OFFSET
                # skip any intersection behind the true start point (with some wiggle room)
                continue
            end
            let dup = false
                # check if a point with the same alpha had already been included, if so ignore this one
                # doing this check now avoids executing newton many times unnecessarily
                for eintv in result
                    eints = halfspaceintersection(eintv)
                    if samepoint(α(eints), α(ints) + ACCEL_SURF_RAY_OFFSET)
                        dup = true
                        break
                    end
                end
                if !dup
                    # evaluate the precise point of the intersection with the surface
                    intuv = newton(surf.surface, r, uv(ints))
                    # check that the evaluated uv point lies on the (correct side of the) ray
                    if intuv !== nothing
                        u, v = intuv
                        surfpt = point(surf, u, v)
                        t = α(r, surfpt) # this is the 'true' α again (i.e. not offset)
                        if t > zero(T)
                            # check again that this point hasn't already been included (α might be slightly different, i.e. more accurate, now)
                            for eintv in result
                                eints = halfspaceintersection(eintv)
                                if samepoint(α(eints), t)
                                    dup = true
                                    break
                                end
                            end
                            if !dup && samepoint(point(r, t), surfpt)
                                # intersection is good so calc the normal and add to list
                                surfnormal = normal(surf, u, v)
                                intsct = Intersection(t, surfpt, surfnormal, u, v, interface(surf))
                                if dot(direction(r), surfnormal) < zero(T)
                                    intvl = positivehalfspace(intsct)
                                else
                                    intvl = rayorigininterval(intsct)
                                end
                                push!(result, intvl)
                            end
                        end
                    end
                end
            end
        end
    end

    nr = length(result)
    if nr === 0
        return EmptyInterval(T)
    elseif nr === 1
        return result[1]
    end

    # sort!(result, lt = (a, b) -> α(halfspaceintersection(a)) < α(halfspaceintersection(b))) #this will sort the intervals by α
    # alloc free sort of result
    @inbounds for i in 2:nr
        value = result[i]
        j = i - 1
        while j > 0 && α(halfspaceintersection(result[j])) > α(halfspaceintersection(value))
            result[j + 1] = result[j]
            j = j - 1
        end
        result[j + 1] = value
    end

    let start = 1
        temp = newinintervalpool!(T)

        if lower(result[1]) isa RayOrigin{T}
            # Starting "inside" the surface, i.e., n̂⋅d̂, > 0.
            # Can't truly define inside for Bezier patches because they don't define a half-space.
            start = 2
            push!(temp, result[1])
        end

        @inbounds for i in start:2:(nr - 1)
            intvlintsct = intervalintersection(result[i], result[i + 1])
            if !(intvlintsct isa EmptyInterval{T})
                push!(temp, intvlintsct)
            else
                # this can happen in some edge cases, seems like the best thing to do is juse ignore it...
                @warn "Funny behavior in triangulated intersection" maxlog = 1
                # throw(ErrorException("Triangulated intersection error - this should never happen."))
            end
        end

        if mod(nr - (start - 1), 2) === 1 #one unpaired interval
            push!(temp, last(result))
        end

        return DisjointUnion(temp)
    end
end

"""
    jacobian(surf::ParametricSurface{T,N}, u::T, v::T, P1::SVector{M,T}, P2::SVector{M,T})

Computes Jacobian of `f(t,u,v) = ( dot(P1,[surf(u,v),1],P2,[surf(u,v),1]) )`.
`P1`, `P2` are orthogonal planes that pass through the ray.
`J = [ ∂f1/∂u ∂f1/∂v ; ∂f2/∂u ∂f2/∂v]`
"""
function jacobian(surf::ParametricSurface{T,N}, u::T, v::T, P1::SVector{M,T}, P2::SVector{M,T}) where {T<:Real,N,M}
    (du, dv) = partials(surf, u, v)
    return SMatrix{2,2,T}(dot(view(P1, 1:(M - 1)), du), dot(view(P2, 1:(M - 1)), du), dot(view(P1, 1:(M - 1)), dv), dot(view(P2, 1:(M - 1)), dv))
end

"""
    newton(surf::ParametricSurface{T,N}, r::AbstractRay{T,N}, startingpoint::SVector{2,T})

Newton iteration to find the precise intersection of a parametric surface with a ray given a starting point (in uv space) on the surface.
"""
function newton(surf::ParametricSurface{T,N}, r::AbstractRay{T,N}, startingpoint::SVector{2,T}) where {T,N}
    tolerance = 1e-12
    maxiterations = 6

    urange, vrange = uvrange(surf)

    dx, dy, dz = direction(r)
    if abs(dx) > abs(dy) && abs(dx) > abs(dz)
        N1 = SVector{3,T}(dy, -dx, zero(T))
    else
        N1 = SVector{3,T}(zero(T), dz, -dy)
    end

    N2 = cross(N1, direction(r))
    P1 = SVector{4,T}(N1[1], N1[2], N1[3], -dot(N1, origin(r)))
    P2 = SVector{4,T}(N2[1], N2[2], N2[3], -dot(N2, origin(r)))

    @inline function f(uv::SVector{2,T})::SVector{2,T}
        psurf = point(surf, uv[1], uv[2])
        res1 = zero(T)
        res2 = zero(T)
        @inbounds for i in 1:3
            res1 += P1[i] * psurf[i]
            res2 += P2[i] * psurf[i]
        end
        res1 += P1[4]
        res2 += P2[4]
        # return [dot(P1, vcat(psurf, 1)), dot(P2, vcat(psurf, 1))]
        return SVector{2,T}(res1, res2)
    end

    xilast = startingpoint
    xi = zeros(SVector{2,T})
    error = typemax(T)
    i = 0
    finalpass = false

    while i < maxiterations
        if any(isnan.(xilast))
            # if something is NaN then just return the starting point and hope that the surface point matches
            xilast = startingpoint
            error = zero(T)
            break
        end

        j = jacobian(surf, xilast[1], xilast[2], P1, P2) # xi has values u,v

        if det(j) == zero(T)
            # if the jacobian is singular then break at the current iteration and return whatever we have
            error = zero(T)
            break
        end

        lstf = f(xilast)
        xi = xilast - (j \ lstf)
        error = sum(abs.(xi - xilast))

        if !finalpass && error < tolerance
            i = maxiterations - 3 # iterate a few more times after the xi and xi-1 are the same value, at least to 64 bit precision. Seems to give more accurate solutions. Perhaps because of higher precision for intermediate results?
            finalpass = true
        end

        i += 1

        xilast = xi
    end

    if error > tolerance
        return nothing
    else
        if xi[1] < urange[1] || xi[1] > urange[2] || xi[2] < vrange[1] || xi[2] > vrange[2]
            return nothing
        else
            return xilast
        end
    end
end

# Only show the underlying Bezier surface control points, not the approximating triangle mesh which is generally huge
Base.show(io::IO, a::AcceleratedParametricSurface{S,N,T}) where {S<:Real,N,T<:ParametricSurface{S,N}} = println(io, "Accelerated$(string(a.surface))")


text:
export Vandermonde

immutable Vandermonde{T} <: AbstractArray{T, 2}
	c :: Vector{T}
end

#XXX Inefficient but works
getindex(V::Vandermonde, i, j) = getindex(full(V), i, j)
isassigned(V::Vandermonde, i, j) = isassigned(full(V), i, j)

size(V::Vandermonde, r::Int) = (r==1 || r==2) ? length(V.c) :
    throw(ArgumentError("Invalid dimension $r"))
size(V::Vandermonde) = length(V.c), length(V.c)

function full{T}(V::Vandermonde{T})
	n=size(V, 1)
	M=Array(T, n, n)
	for i=1:n
		M[:,i] = V.c.^(i-1)
	end
	M
end


text:
# We have different sources for problem functions:
#  S1 = CEC 2013 competition on large-scale optimization
#  S2 = JADE paper http://150.214.190.154/EAMHCO/pdf/JADE.pdf
#  S3 = "Test Suite for the Special Issue of Soft Computing on Scalability of
#        Evolutionary Algorithms and other Metaheuristics for Large Scale
#        Continuous Optimization Problems", http://sci2s.ugr.es/eamhco/functions1-19.pdf
# Our primary focus is to implement all the problems from S1 since our
# focus is on large-scale optimization but these problems also can be used
# in lower dimensions.
include(joinpath(dirname(@__FILE__()), "single_objective_base_functions.jl"))
Shekel10 = fixeddim_problem(shekel10; range = (0.0, 10.0), dims = 4, name = "Shekel10", fmins = [-10.5364])
Shekel7 = fixeddim_problem(shekel7; range = (0.0, 10.0), dims = 4, name = "Shekel7", fmins = [-10.4029])
Shekel5 = fixeddim_problem(shekel5; range = (0.0, 10.0), dims = 4, name = "Shekel5", fmins = [-10.1532])
Hartman6 = fixeddim_problem(hartman6; range = (0.0, 1.0), dims = 6, name = "Hartman6", fmins = [-3.32237])
Hartman3 = fixeddim_problem(hartman3; range = (0.0, 1.0), dims = 3, name = "Hartman3", fmins = [-3.860038442])

# We skip (for now) f12 and f13 in the JADE paper since they are penalized
# functions which are quite nonstandard. We also skip f8 since we are unsure
# about its proper implementation.
JadeFunctionSet = {
  1   => anydim_problem("Sphere",        sphere,        (-100.0, 100.0), 0.0),
  2   => anydim_problem("Schwefel2.22",  schwefel2_22,  ( -10.0,  10.0), 0.0),
  3   => anydim_problem("Schwefel1.2",   schwefel1_2,   (-100.0, 100.0), 0.0),
  4   => anydim_problem("Schwefel2.21",  schwefel2_21,  (-100.0, 100.0), 0.0),
  5   => anydim_problem("Rosenbrock",    rosenbrock,    ( -30.0,  30.0), 0.0),
  6   => anydim_problem("Step",          s2_step,       (-100.0, 100.0), 0.0),
  7   => anydim_problem("Noisy quartic", noisy_quartic, ( -30.0,  30.0)),
#  8   => anydim_problem("Schwefel2.26",  schwefel2_26,  (-500.0, 500.0)),
  9   => anydim_problem("Rastrigin",     rastrigin,     ( -5.12,  5.12), 0.0),
  10  => anydim_problem("Ackley",        ackley,        ( -32.0,  32.0), 0.0),
  11  => anydim_problem("Griewank",      griewank,      (-600.0, 600.0), 0.0)
}

#####################################################################
# S3 Base functions.
#####################################################################

#####################################################################
# S3 Transformations
#####################################################################

# A TransformedProblem just makes a few changes in a sub-problem but refers
# most func calls to it. Concrete types must implement a sub_problem func.
abstract TransformedProblem <: OptimizationProblem
search_space(tp::TransformedProblem) = search_space(sub_problem(tp))
is_fixed_dimensional(tp::TransformedProblem) = is_fixed_dimensional(sub_problem(tp))
numfuncs(tp::TransformedProblem) = numfuncs(sub_problem(tp))
numdims(tp::TransformedProblem) = numdims(sub_problem(tp))
fmins(tp::TransformedProblem) = fmins(sub_problem(tp))
name(tp::TransformedProblem) = name(sub_problem(tp))

# A ShiftedAndBiasedProblem shifts the minimum value and biases the returned
# function values.
type ShiftedAndBiasedProblem <: TransformedProblem
  xshift::Array{Float64, 1}
  funcshift::Float64
  subp::OptimizationProblem

  ShiftedAndBiasedProblem(sub_problem::OptimizationProblem;
    xshift = false, funcshift = 0.0) = begin
    xshift = (xshift != false) ? xshift : rand_individual(search_space(sub_problem))
    new(xshift[:], funcshift, sub_problem)
  end
end

sub_problem(sp::ShiftedAndBiasedProblem) = sp.subp

is_fixed_dimensional(p::ShiftedAndBiasedProblem) = is_fixed_dimensional(sub_problem(p))

# Evaluate by first shifting x and then biasing the returned function value.
evalfunc(x, i::Int64, sp::ShiftedAndBiasedProblem) = begin
  ofunc(sub_problem(sp), i)(x - sp.xshift) + sp.funcshift
end

shifted(p::OptimizationProblem; funcshift = 0.0) = ShiftedAndBiasedProblem(p;
  funcshift = funcshift)

#####################################################################
# S1 Base functions. Typically slightly transformed to break symmetry
#   and introduce irregularities.
#####################################################################
s1_sphere = sphere

function s1_elliptic(x)
  xt = t_irreg(x)
  elliptic(xt)
end

function s1_rastrigin(x)
  xt = t_diag(t_asy(t_irreg(x), 0.2), 10)
  rastrigin(xt)
end

function s1_ackley(x)
  xt = t_diag(t_asy(t_irreg(x), 0.2), 10)
  ackley(xt)
end

function s1_schwefel(x)
  xt = t_asy(t_irreg(x), 0.2)
  schwefel1_2(xt)
end

s1_rosenbrock = rosenbrock

#####################################################################
# S1 Transformations
#####################################################################

# This transformation function is used to break the symmetry of symmetric
# functions.
function t_asy(f, beta)
  D = length(f)
  g = copy(f)
  temp = beta * linspace(0, 1, D)
  ind = collect(1:D)[f .> 0]
  t = f[ind] .^ (1 + temp[ind] .* sqrt(f[ind]))
  setindex!(g, t, ind)
  g
end

# This transformation is used to create the ill-conditioning effect.
function t_diag(f, alpha)
  D = length(f)
  scales = sqrt(alpha) .^ linspace(0, 1, D)
  scales .* f
end

# This transformation is used to create smooth local irregularities.
function t_irreg(f)
   a = 0.1
   g = copy(f)
   indices = collect(1:length(f))

   idxp = indices[f .> 0]
   t = log(f[idxp])/a
   r = exp(t + 0.49*(sin(t) + sin(0.79*t))).^a
   setindex!(g, r, idxp)

   idxn = indices[f .< 0]
   t = log(-f[idxn])/a
   r = -exp(t + 0.49*(sin(0.55*t) + sin(0.31*t))).^a
   setindex!(g, r, idxn)

   g
end

function xshifted(n, f)
  move = 10.0 * randn(n, 1)
  transformed_f(x) = f(x .- move)
end

function xrotatedandshifted(n, f, shiftAmplitude = 1.0, rotateAmplitude = 1.0)
  shift = shiftAmplitude * randn(n, 1)
  rotmatrix = rotateAmplitude * rand(n, n)
  transformed_f(x) = f(rotmatrix * (x .- shift))
end

example_problems = {
  "Sphere" => JadeFunctionSet[1],
  "Rosenbrock" => JadeFunctionSet[5],
  "Schwefel2.22" => JadeFunctionSet[2],
  "Schwefel1.2" => JadeFunctionSet[3],
  "Schwefel2.21" => JadeFunctionSet[4],
  "Step" => JadeFunctionSet[6],
  "Rastrigin" => JadeFunctionSet[9],
  "Ackley" => JadeFunctionSet[10],
  "Griewank" => JadeFunctionSet[11],
  "Ellipsoid" => anydim_problem("Ellipsoid", ellipsoid, (-65.536, 65.536), 0.0),
  "Cigar" => anydim_problem("Cigar", cigar, (-100.0, 100.0), 0.0),
  "DeceptiveCuccu2011_15_2" => anydim_problem("DeceptiveCuccu2011_15_2", deceptive_cuccu2011_15_2, (-100.0, 100.0), 0.0),
  "Shekel10" => Shekel10,
  "Shekel7" => Shekel7,
  "Shekel5" => Shekel5,
  "Hartman6" => Hartman6,
  "Hartman3" => Hartman3
}


text:
#
# using Mux

include( "server_libs.jl")

const DEFAULT_PORT=8054
const DEFAULT_SERVER="http://localhost:$DEFAULT_PORT/"

port = DEFAULT_PORT
if length(ARGS) > 0
   port = parse(Int, ARGS[1])
end

# the wait() bit here should allow 'clean breaks' may not be needed
# in a live server.
# see note in: https://github.com/JuliaWeb/Mux.jl/blob/master/src/server.jl
# this doesn't clear the port on shutdown
# use: fuser 8054/tcp to kill everything
# see: https://stackoverflow.com/questions/750604/freeing-up-a-tcp-ip-port
# wait(serve( stbdemo, port))
host = "0.0.0.0"
HTTP.serve(Mux.http_handler(stbdemo), host, port )

#=
while true # FIXME better way?
   @debug "new STB Server; main loop; server running on port $port"
   sleep( 60 )
end
=#


text:
module _XK2En
export XK2En
using Statistics: std, mean
    """
        XK2, Ci = XK2En(Sig)

    Returns the cross-Kolmogorov entropy estimates (`XK2`) and the correlation
    integrals (`Ci`) for m = [1,2] estimated between the data sequences
    contained in `Sig` using the default parameters:
    embedding dimension = 2, time delay = 1, distance threshold (r) = 0.2*SD(Sig),
    logarithm = natural

        XK2, Ci = XK2En(Sig::AbstractArray{T,2} where T<:Real; m::Int=2, tau::Int=1, r::Real=0.2*std(Sig,corrected=false), Logx::Real=exp(1))

    Returns the cross-Kolmogorov entropy estimates (`XK2`) estimated between
    the data sequences contained in `Sig` using the specified 'keyword' arguments:

    # Arguments:
    `m`     - Embedding Dimension, a positive integer [default: 2]  \n
    `tau`   - Time Delay, a positive integer         [default: 1]  \n
    `r`     - Radius Distance Threshold, a positive scalar  [default: 0.2*SD(`Sig`)]  \n
    `Logx`  - Logarithm base, a positive scalar      [default: natural]  \n

    # See also `XSampEn`, `XFuzzEn`, `XApEn`, `K2En`, `XMSEn`, `XDistEn`

    # References:
        [1]  Matthew W. Flood,
             "XK2En - EntropyHub Project"
             (2021) https://github.com/MattWillFlood/EntropyHub

    """
    function XK2En(Sig::AbstractArray{T,2} where T<:Real; m::Int=2,
        tau::Int=1, r::Real=0.2*std(Sig,corrected=false), Logx::Real=exp(1))

    size(Sig,2) > size(Sig,1) ? Sig = transpose(Sig) : nothing

    N = size(Sig,1)
    (N>10 && size(Sig,2)==2) ? nothing :  error("Sig:   must be a 2-columns matrix")
    (m > 0) ? nothing : error("m:     must be an integer > 0")
    (tau>0) ? nothing : error("tau:   must be an integer > 0")
    (Logx>0) ? nothing : error("Logx: must be a positive number > 0")
    (r>0) ? nothing : error("r:  must be 2 element tuple of positive values")

    m   += 1
    Zm1 = zeros(N,m)
    Zm2 = zeros(N,m)
    Ci = zeros(m)
    for n = 1:m
        N2 = N-(n-1)*tau
        Zm1[1:N2,n] = Sig[(n-1)*tau + 1:N,1]
        Zm2[1:N2,n] = Sig[(n-1)*tau + 1:N,2]
        Norm = zeros(N2,N2)
        for k = 1:N2
            Temp = repeat(transpose(Zm1[k,1:n]),outer=N2) .- Zm2[1:N2,1:n]
            Norm[k,:] = sqrt.(sum(Temp.*Temp,dims=2))
        end
        Ci[n] = mean(Norm .< r)
    end

    XK2 = log.(Logx, Ci[1:m-1]./Ci[2:m])/tau
    XK2[isinf.(XK2)] .= NaN

    return XK2, Ci

    end
end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
@testset "Simple tests" begin

    model(x, θ) = θ[1] * exp(x[1] * θ[2])

    gmodel!(g, x, θ) = begin

        g[1] = exp(x[1] * θ[2])
        g[2] = x[1] * θ[1] * exp(x[1] * θ[2])

    end

    @testset "Basic usage" begin

        data = [-1.0  3.2974425414002564;
                -0.75 2.9099828292364025;
                -0.5   2.568050833375483;
                -0.25 2.2662969061336526;
                0.0                  2.0;
                0.25   1.764993805169191;
                0.5   1.5576015661428098;
                0.75  1.5745785575819442; #noise
                1.0   1.2130613194252668;
                1.25  1.0705228570379806;
                1.5   0.9447331054820294;
                1.75  0.8337240393570168;
                2.0   0.7357588823428847;
                2.25  0.6493049347166995;
                2.5   0.5730095937203802;
                2.75  0.5056791916094929;
                3.0  0.44626032029685964;
                3.25  0.5938233504083881; #noise
                3.5   0.3475478869008902;
                3.75 0.30670993368985694;
                4.0   0.5706705664732254; #noise
                ]

        answer = [2.0, -0.5]

        θ = [0.0, 0.0]

        rout = lmlovo(model, θ, data, 2, 18)

        @test rout.status == 1
        @test rout.solution ≈ answer atol=1.0e-5
        @test rout.p == 18
        @test rout.nf >= rout.iter
        @test rout.nj >= rout.iter

        θ = [0.0, 0.0]

        # Test with small p
        rout = lmlovo(model, θ, data, 2, 3)

        @test rout.status == 1
        @test rout.p == 3

        θ = [0.0, 0.0]

        rout = raff(model, data, 2)

        @test rout.f ≈ 0.0 atol=1.0e-5
        @test rout.solution ≈ answer atol=1.0e-5
        @test rout.p == 18
        @test rout.iter >= size(data)[1]
        @test rout.nf >= 1
        @test rout.nj >= 1

        @test_throws AssertionError lmlovo(model, θ, data, 0, 1)
        @test_throws AssertionError lmlovo(model, θ, data, 2, -1)

        rout = lmlovo(model, θ, data, 2, 0)

        @test rout.status == 1
        @test rout.iter == 0
        @test rout.f == 0
        @test rout.outliers == [1:size(data)[1];]
        @test rout.solution == θ
        @test rout.nf == 0
        @test rout.nj == 0

        # lmlovo with function and gradient

        θ = [0.0, 0.0]

        rout = lmlovo(model, gmodel!, θ, data, 2, 18)

        @test rout.status == 1
        @test rout.solution ≈ answer atol=1.0e-5
        @test rout.p == 18

        θ = [0.0, 0.0]

        rout = raff(model, gmodel!, data, 2)

        @test rout.f ≈ 0.0 atol=1.0e-5
        @test rout.solution ≈ answer atol=1.0e-5
        @test rout.p == 18
        @test rout.iter >= size(data)[1]
        @test rout.nf >= 1
        @test rout.nj >= 1

        @test_throws AssertionError lmlovo(model, gmodel!, θ, data, 0, 1)
        @test_throws AssertionError lmlovo(model, gmodel!, θ, data, 2, -1)

        rout = lmlovo(model, gmodel!, θ, data, 2, 0)

        @test rout.status == 1
        @test rout.iter == 0
        @test rout.f == 0
        @test rout.outliers == [1:size(data)[1];]
        @test rout.solution == θ
        @test rout.nf == 0
        @test rout.nj == 0

    end

    # Test to check Issue #1

    @testset "Error in printing" begin

        m(x, θ) = θ[1] * x[1]^2 + θ[2]

        A = [ -2.0  5.00;
              -1.5  3.25;
              -1.0  2.00;
              -0.5  1.25;
              0.0  1.00;
              0.5  1.25;
              1.0  2.00;
              1.5  3.25;
              2.0  5.00 ]

        θ = [0.0, 0.0]

        # Changes log just for this test
        rout = with_logger(Logging.NullLogger()) do

            lmlovo(m, θ, A, 2, 4)

        end

        @test rout.status == 1
        @test rout.p == 4

    end

    @testset "Test parameters" begin

        data = [-1.0   3.2974425414002564;
                -0.75  2.9099828292364025;
                -0.5    2.568050833375483;
                -0.25  2.2662969061336526;
                 0.0                  2.0;
                 0.25   1.764993805169191;
                 0.5   1.5576015661428098;
                 0.75  1.5745785575819442; #noise
                 1.0   1.2130613194252668;
                 1.25  1.0705228570379806;
                 1.5   0.9447331054820294;
                 1.75  0.8337240393570168;
                 2.0   0.7357588823428847;
                 2.25  0.6493049347166995;
                 2.5   0.5730095937203802;
                 2.75  0.5056791916094929;
                 3.0  0.44626032029685964;
                 3.25  0.5938233504083881; #noise
                 3.5   0.3475478869008902;
                 3.75 0.30670993368985694;
                 4.0   0.5706705664732254; #noise
                ]

        answer = [2.0, -0.5]

        rout = raff(model, gmodel!, data, 2; noutliers=0)

        @test rout.p == 21

        rout = raff(model, data, 2; noutliers=5)

        @test rout.f ≈ 0.0 atol=1.0e-5
        @test rout.solution ≈ answer atol=1.0e-5
        @test rout.p == 18

        rout = raff(model, data, 2; ftrusted=(21 - 5)/21)

        @test rout.f ≈ 0.0 atol=1.0e-5
        @test rout.solution ≈ answer atol=1.0e-5
        @test rout.p == 18

        rout = raff(model, data, 2; ftrusted=(18/21, 18/21))

        @test rout.f ≈ 0.0 atol=1.0e-5
        @test rout.solution ≈ answer atol=1.0e-5
        @test rout.p == 18

        @test raff(model, data, 2; ftrusted=(0.5, 1.1)) == RAFFOutput()
        @test raff(model, data, 2; ftrusted=-0.1) == RAFFOutput()

    end

    # Tests for RAFFOutput

    @testset "RAFFOutput tests" begin

        nullOut = RAFFOutput(0, [], -1, 0, Inf, -1, -1, [])

        @test RAFFOutput() == nullOut

        @test RAFFOutput(0) == nullOut

        # Check if deprecated version is creating `nf` and `nj`
        @test RAFFOutput(0, [], -1, 0, Inf, []) == nullOut

        nullPOut = RAFFOutput(0, [], -1, 10, Inf, -1, -1, [])

        @test nullPOut == RAFFOutput(10)

        # Test output

        raff_output = RAFFOutput(1, ones(5), 2, 6, - 1.0, 10, 20, ones(Int, 6))

        io = IOBuffer()

        print(io, raff_output)

        s = String(take!(io))

        rx = Regex("\\(\\.status\\) = " * string(raff_output.status))

        @test match(rx, s) !== nothing

        svec = replace(string(raff_output.solution), r"([\[\]])"=>s"\\\1")

        rx = Regex("\\(\\.solution\\) = " * svec)

        @test match(rx, s) !== nothing

        rx = Regex("\\(\\.iter\\) = " * string(raff_output.iter))

        @test match(rx, s) !== nothing

        rx = Regex("\\(\\.p\\) = " * string(raff_output.p))

        @test match(rx, s) !== nothing

        rx = Regex("\\(\\.f\\) = " * string(raff_output.f))

        @test match(rx, s) !== nothing

        rx = Regex("\\(\\.nf\\) = " * string(raff_output.nf))

        @test match(rx, s) !== nothing

        rx = Regex("\\(\\.nj\\) = " * string(raff_output.nj))

        @test match(rx, s) !== nothing

        svec = replace(string(raff_output.outliers), r"([\[\]])"=>s"\\\1")

        rx = Regex("\\(\\.outliers\\) = " * svec)

        @test match(rx, s) !== nothing

    end

end


text:
#
# Graph type as wrapper around the c++ mglGraph class
# There are other types of graphs (qt, wx, glut, ...) which are not
# supported right now
#

type Graph
    # Pointer to the c++ mglGraph class
    ptr::Ptr{Void}

    function Graph(width::Integer=600, height::Integer=400)
        graph = new(cmgl.create_graph(width, height))
        finalizer(graph, x -> cmgl.delete_graph(x.ptr))
        return graph
    end
end

#
# Data type as wrapper around the standard mglData class for real floating
# point numbers
#

type Data
    # Pointer to the c++ mglData class
    ptr::Ptr{Void}

    function Data()
        data = new(cmgl.create_data())
        finalizer(data, x -> cmgl.delete_data(x.ptr))
        return data
    end
end

# Create Mgl data class with link to a julia array
function Data(values::Array{cmgl.Float})
    @assert ndims(values) <= 3 "Dimension of the input array is $(ndims(values))." +
                               "The MathGl data class only supports up to 3 dimensions."
    data = Data()
    cmgl.data_link(data.ptr, values, size(values, 1), size(values, 2), size(values, 3))
    return data
end

# Link or Relink Mgl data class to a julia array
function Data(data::Data, values::Array{cmgl.Float})
    @assert ndims(values) <= 3 "Dimension of the input array is $(ndims(values))." +
                               "The MathGl data class only supports up to 3 dimensions."
    cmgl.data_link(data.ptr, values, size(values, 1), size(values, 2), size(values, 3))
    return data
end

#
# DataC type as wrapper around the standard mglData class for complex floating
# point numbers
#

type DataC
    # Pointer to the c++ mglDataC class
    ptr::Ptr{Void}

    function DataC()
        data = new(cmgl.create_datac())
        finalizer(data, x -> cmgl.delete_datac(x.ptr))
        return data
    end
end

# Create Mgl dataC class with link to a julia array
function DataC(values::Array{cmgl.Dual})
    @assert ndims(values) <= 3 "Dimension of the input array is $(ndims(values))." +
                               "The MathGl data class only supports up to 3 dimensions."
    data = DataC()
    cmgl.datac_link(data.ptr, values, size(values, 1), size(values, 2), size(values, 3))
    return data
end

# Link or relink Mgl dataC class to a julia array
function DataC(data::DataC, values::Array{cmgl.Dual})
    @assert ndims(values) <= 3 "Dimension of the input array is $(ndims(values))." +
                               "The MathGl data class only supports up to 3 dimensions."
    cmgl.datac_link(data.ptr, values, size(values, 1), size(values, 2), size(values, 3))
    return data
end


text:
# Quelques méthodes de manipulation de tableau/vecteur

"""
Déplacement d'un élément dans un tableau

Args:
    v: vecteur à modifier
    idx1: indice de l'élément à décaller
    idx2: indice de l'émément après déplacement

Le tableau passé en paramètre est modifié.

Returns:
    Le tableau complet

Example:
    ```
    v = collect(1:7)

    shift!(v, 7, 1)'

    => 1×7 LinearAlgebra.Adjoint{Int64,Array{Int64,1}}:

    7  1  2  3  4  5  6
    ```

See Also:

    - permutate()
    - circshift()

"""
function shift!(v::Vector{T}, idx1::Int, idx2::Int) where {T}
    if idx1<idx2
        # si idx1=2 et idx2=6
        # alors v[2,3,4,5 , 6] devient v[3,4,5,6 , 2]
        v[[idx1:idx2-1;idx2]] = v[[idx1+1:idx2;idx1]]
    else
        v[[idx2+1:idx1;idx2]] = v[[idx2:idx1-1;idx1]]
    end
    return v
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VpnClientConnectionHealthDetail <: SwaggerModel
    vpnConnectionId::Any # spec type: Union{ Nothing, String } # spec name: vpnConnectionId
    vpnConnectionDuration::Any # spec type: Union{ Nothing, Int64 } # spec name: vpnConnectionDuration
    vpnConnectionTime::Any # spec type: Union{ Nothing, String } # spec name: vpnConnectionTime
    publicIpAddress::Any # spec type: Union{ Nothing, String } # spec name: publicIpAddress
    privateIpAddress::Any # spec type: Union{ Nothing, String } # spec name: privateIpAddress
    vpnUserName::Any # spec type: Union{ Nothing, String } # spec name: vpnUserName
    maxBandwidth::Any # spec type: Union{ Nothing, Int64 } # spec name: maxBandwidth
    egressPacketsTransferred::Any # spec type: Union{ Nothing, Int64 } # spec name: egressPacketsTransferred
    egressBytesTransferred::Any # spec type: Union{ Nothing, Int64 } # spec name: egressBytesTransferred
    ingressPacketsTransferred::Any # spec type: Union{ Nothing, Int64 } # spec name: ingressPacketsTransferred
    ingressBytesTransferred::Any # spec type: Union{ Nothing, Int64 } # spec name: ingressBytesTransferred
    maxPacketsPerSecond::Any # spec type: Union{ Nothing, Int64 } # spec name: maxPacketsPerSecond

    function VpnClientConnectionHealthDetail(;vpnConnectionId=nothing, vpnConnectionDuration=nothing, vpnConnectionTime=nothing, publicIpAddress=nothing, privateIpAddress=nothing, vpnUserName=nothing, maxBandwidth=nothing, egressPacketsTransferred=nothing, egressBytesTransferred=nothing, ingressPacketsTransferred=nothing, ingressBytesTransferred=nothing, maxPacketsPerSecond=nothing)
        o = new()
        validate_property(VpnClientConnectionHealthDetail, Symbol("vpnConnectionId"), vpnConnectionId)
        setfield!(o, Symbol("vpnConnectionId"), vpnConnectionId)
        validate_property(VpnClientConnectionHealthDetail, Symbol("vpnConnectionDuration"), vpnConnectionDuration)
        setfield!(o, Symbol("vpnConnectionDuration"), vpnConnectionDuration)
        validate_property(VpnClientConnectionHealthDetail, Symbol("vpnConnectionTime"), vpnConnectionTime)
        setfield!(o, Symbol("vpnConnectionTime"), vpnConnectionTime)
        validate_property(VpnClientConnectionHealthDetail, Symbol("publicIpAddress"), publicIpAddress)
        setfield!(o, Symbol("publicIpAddress"), publicIpAddress)
        validate_property(VpnClientConnectionHealthDetail, Symbol("privateIpAddress"), privateIpAddress)
        setfield!(o, Symbol("privateIpAddress"), privateIpAddress)
        validate_property(VpnClientConnectionHealthDetail, Symbol("vpnUserName"), vpnUserName)
        setfield!(o, Symbol("vpnUserName"), vpnUserName)
        validate_property(VpnClientConnectionHealthDetail, Symbol("maxBandwidth"), maxBandwidth)
        setfield!(o, Symbol("maxBandwidth"), maxBandwidth)
        validate_property(VpnClientConnectionHealthDetail, Symbol("egressPacketsTransferred"), egressPacketsTransferred)
        setfield!(o, Symbol("egressPacketsTransferred"), egressPacketsTransferred)
        validate_property(VpnClientConnectionHealthDetail, Symbol("egressBytesTransferred"), egressBytesTransferred)
        setfield!(o, Symbol("egressBytesTransferred"), egressBytesTransferred)
        validate_property(VpnClientConnectionHealthDetail, Symbol("ingressPacketsTransferred"), ingressPacketsTransferred)
        setfield!(o, Symbol("ingressPacketsTransferred"), ingressPacketsTransferred)
        validate_property(VpnClientConnectionHealthDetail, Symbol("ingressBytesTransferred"), ingressBytesTransferred)
        setfield!(o, Symbol("ingressBytesTransferred"), ingressBytesTransferred)
        validate_property(VpnClientConnectionHealthDetail, Symbol("maxPacketsPerSecond"), maxPacketsPerSecond)
        setfield!(o, Symbol("maxPacketsPerSecond"), maxPacketsPerSecond)
        o
    end
end # type VpnClientConnectionHealthDetail

const _property_map_VpnClientConnectionHealthDetail = Dict{Symbol,Symbol}(Symbol("vpnConnectionId")=>Symbol("vpnConnectionId"), Symbol("vpnConnectionDuration")=>Symbol("vpnConnectionDuration"), Symbol("vpnConnectionTime")=>Symbol("vpnConnectionTime"), Symbol("publicIpAddress")=>Symbol("publicIpAddress"), Symbol("privateIpAddress")=>Symbol("privateIpAddress"), Symbol("vpnUserName")=>Symbol("vpnUserName"), Symbol("maxBandwidth")=>Symbol("maxBandwidth"), Symbol("egressPacketsTransferred")=>Symbol("egressPacketsTransferred"), Symbol("egressBytesTransferred")=>Symbol("egressBytesTransferred"), Symbol("ingressPacketsTransferred")=>Symbol("ingressPacketsTransferred"), Symbol("ingressBytesTransferred")=>Symbol("ingressBytesTransferred"), Symbol("maxPacketsPerSecond")=>Symbol("maxPacketsPerSecond"))
const _property_types_VpnClientConnectionHealthDetail = Dict{Symbol,String}(Symbol("vpnConnectionId")=>"String", Symbol("vpnConnectionDuration")=>"Int64", Symbol("vpnConnectionTime")=>"String", Symbol("publicIpAddress")=>"String", Symbol("privateIpAddress")=>"String", Symbol("vpnUserName")=>"String", Symbol("maxBandwidth")=>"Int64", Symbol("egressPacketsTransferred")=>"Int64", Symbol("egressBytesTransferred")=>"Int64", Symbol("ingressPacketsTransferred")=>"Int64", Symbol("ingressBytesTransferred")=>"Int64", Symbol("maxPacketsPerSecond")=>"Int64")
Base.propertynames(::Type{ VpnClientConnectionHealthDetail }) = collect(keys(_property_map_VpnClientConnectionHealthDetail))
Swagger.property_type(::Type{ VpnClientConnectionHealthDetail }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VpnClientConnectionHealthDetail[name]))}
Swagger.field_name(::Type{ VpnClientConnectionHealthDetail }, property_name::Symbol) =  _property_map_VpnClientConnectionHealthDetail[property_name]

function check_required(o::VpnClientConnectionHealthDetail)
    true
end

function validate_property(::Type{ VpnClientConnectionHealthDetail }, name::Symbol, val)
end


text:
using DataDeps

register(DataDep(
  "cnae-9",
  "http://archive.ics.uci.edu/ml/datasets/CNAE-9",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/00233/CNAE-9.data",
  "fec557bade36db8d35b29596d73eab7fd88c26f75f0ba07ee69636614c5d6055",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=1,
        feature_indices=2:857,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
# # An equation of forth order in one dimmension
#
# Consider the following problem. Given \\(f \in C([0, 1])\\), find a function \\(u\\) satisfying
# u'''' = f in (0, 1)
# u(0) = u(1) = 0
# u'(0) = pi; u'(1) = pi

using SymFEL
using SymPy
using LinearAlgebra
using SparseArrays
using PyPlot
close("all")

# discretization parameters
N = 101; # number of nodes
dx = 1 / (N - 1); # discretization step
nodes = range(0, stop=1, length=N);
bound_nodes = [1, N];

# exact solution
u_exact = sin.(pi*nodes);
ud_exact = pi * cos.(pi*nodes)
# right hand
f = pi^4 * sin.(pi * nodes);
fd = pi^5 * cos.(pi * nodes);
ff = zeros(2*N);
ff[1:2:2*N-1] = f;
ff[2:2:2*N] = fd;

# elementary matrices
elem_K = SymFEL.get_hermite_em(3, 2, 2);
elem_M = SymFEL.get_em(3, 3, 0, 0; fe1="Hermite", fe2="Hermite")
elem_K_dx = convert(Matrix{Float64}, elem_K.subs(h, dx));
elem_M_dx = convert(Matrix{Float64}, elem_M.subs(h, dx));

K = SymFEL.assemble_1d_FE_matrix(elem_K_dx, N,
                                 intNodes1=0, intNodes2=0,
                                 dof1=2, dof2=2);
M = SymFEL.assemble_1d_FE_matrix(elem_M_dx, N,
                                 intNodes1=0, intNodes2=0,
                                 dof1=2, dof2=2);

F = M * ff;

# boundary conditions
tgv = 1e100
KB = copy(K)
KB[2*bound_nodes .- 1, 2*bound_nodes .- 1] += tgv*sparse(Matrix{Float64}(I, 2, 2));
KB[2*bound_nodes, 2*bound_nodes] += tgv*sparse(Matrix{Float64}(I, 2, 2));

F[2*bound_nodes .- 1] = zeros(2);
F[2] = pi*tgv
F[2*N] = -pi*tgv

u = KB \ F;

plot(nodes, u_exact);
plot(nodes, u[1:2:2*N-1]);
legend(["Exact solution", "Approximate solution"]);

U_exact = zeros(2*N)
U_exact[1:2:2*N-1] = u_exact
U_exact[2:2:2*N] = ud_exact
err = u - U_exact;
println("H1 error = ", sqrt(err' * K * err))


text:
###############################################################################
#
# Whole plant hydraulics --- Element and whole-tree conductance ratio
#
###############################################################################
"""
    plant_conductances!(tree::TreeSimple{FT}) where {FT<:AbstractFloat}

Update plant total hydraulic information in each element and whole plant, given
- `tree` [`TreeSimple`](@ref) type struct
"""
function plant_conductances!(tree::TreeSimple{FT}) where {FT<:AbstractFloat}
    # Calculate the resistances
    r_m_root = 1 / (tree.root).k_max;
    r_m_stem = 1 / (tree.stem).k_max;
    r_m_leaf = 1 / ((tree.leaf).k_sla * (tree.leaf).area);
    r_m_tree = r_m_root + r_m_stem + r_m_leaf

    # Calculate the relative resistances
    r_root = FT(0);
    r_stem = FT(0);
    r_leaf = FT(0);
    for i in 1:tree.root.N
        r_root += (tree.root).k_history[i];
    end
    for i in 1:tree.stem.N
        r_stem += (tree.stem).k_history[i];
    end
    for i in 1:tree.leaf.N
        r_leaf += (tree.leaf).k_history[i];
    end
    r_root /= tree.root.N;
    r_stem /= tree.stem.N;
    r_leaf /= tree.leaf.N;
    r_tree  = r_root * r_m_root + r_stem * r_m_stem + r_leaf * r_m_leaf;

    # Update the kr info
    tree.krs[1] = 1 / r_root;
    tree.krs[2] = 1 / r_stem;
    tree.krs[3] = 1 / r_leaf;
    tree.krs[4] = r_m_tree / r_tree;

    return nothing
end

###############################################################################
#
# Root pressure and conductance
#
###############################################################################
"""
    root_pk(root::RootHydraulics{FT},
            flow::FT
    ) where {FT<:AbstractFloat}
    root_pk(root::RootHydraulics{FT},
            q_in::FT,
            flow::Array{FT,1}
    ) where {FT<:AbstractFloat}

Return root xylem end pressure and root hydraulic conductance (reverse of
    summed resistance), given
- `root` [`RootHydraulics`](@ref) struct
- `flow` Given flow rate(s) in the root layer, array for non-steady state with
    capacitance enabled
- `q_in` Flow rate into the root
"""
function root_pk(
            root::RootHydraulics{FT},
            flow::FT
) where {FT<:AbstractFloat}
    @unpack f_st, f_vis, k_element, k_history, k_rhiz, p_gravity, p_history,
            p_osm, p_ups, sh, T_sap, vc = root;

    # make sure that p_ups is not p_25 and then convert
    p_end::FT = p_ups;
    p_25 ::FT = p_end / f_st;
    r_all::FT = FT(0);

    # compute pressure drop along rhizosphere, using p_25 for Θ
    _dr = 1 / k_rhiz * f_vis / 10;
    _dp = flow * _dr;
    for i in 1:10
        # No idea why soil_k_ratio_p25 results in unnecessary allocations
        # _f   = soil_k_ratio_p25(sh, p_25);
        _rwc   = soil_rwc(sh, p_25);
        _f     = soil_k_ratio_rwc(sh, _rwc);
        p_25  -= _dp / _f;
        r_all += _dr / _f;
    end
    p_end = p_25 * f_st + p_osm * T_sap / T_25(FT);

    # compute k from temperature and history, then update pressure
    for (_k, _kh, _pg, _ph) in zip(k_element, k_history, p_gravity, p_history)
        p_25 = p_end / f_st;
        if p_25 < _ph
            k = xylem_k_ratio(vc, p_25, f_vis) * _k;
        else
            k = _kh * _k;
        end
        p_end -= flow / k + _pg;
        r_all += 1 / k;
    end

    k_all = 1 / r_all;

    return p_end, k_all
end

function root_pk(
            root::RootHydraulics{FT},
            q_in::FT,
            flow::Array{FT,1}
) where {FT<:AbstractFloat}
    @unpack f_st, f_vis, k_element, k_history, k_rhiz, p_gravity, p_history,
            p_osm, p_ups, sh, T_sap, vc = root;

    # make sure that p_ups is not p_25 and then convert
    p_end::FT = p_ups;
    p_25 ::FT = p_end / f_st;
    r_all::FT = FT(0);

    # compute pressure drop along rhizosphere, using p_25 for Θ
    _dr = 1 / k_rhiz * f_vis / 10;
    _dp = q_in * _dr;
    for i in 1:10
        # No idea why soil_k_ratio_p25 results in unnecessary allocations
        # _f   = soil_k_ratio_p25(sh, p_25);
        _rwc   = soil_rwc(sh, p_25);
        _f     = soil_k_ratio_rwc(sh, _rwc);
        p_25  -= _dp / _f;
        r_all += _dr / _f;
    end
    p_end = p_25 * f_st + p_osm * T_sap / T_25(FT);

    # compute k from temperature and history, then update pressure
    for (_k, _kh, _pg, _ph, _fl) in zip(k_element, k_history, p_gravity,
                                        p_history, flow)
        p_25 = p_end / f_st;
        if p_25 < _ph
            k = xylem_k_ratio(vc, p_25, f_vis) * _k;
        else
            k = _kh * _k;
        end
        p_end -= _fl / k + _pg;
        r_all += 1 / k;
    end

    k_all = 1 / r_all;

    return p_end, k_all
end

###############################################################################
#
# Calculate xylem end risk factor (k_ratio)
#
###############################################################################
"""
    xylem_risk(hs::LeafHydraulics{FT}, flow::FT) where {FT<:AbstractFloat}

Evaluate the hydraulic risk at the end of leaf xylem, given
- `hs` [`LeafHydraulics`](@ref) type struct
- `flow` Flow rate (per leaf area)
"""
function xylem_risk(
            hs::LeafHydraulics{FT},
            flow::FT
) where {FT<:AbstractFloat}
    @unpack f_st, f_vis, vc = hs;

    p_25 = end_pressure(hs, flow) / hs.f_st;
    T_25 = xylem_k_ratio(vc, p_25, f_vis);

    return T_25
end

max_stars_repo_name: climate-machine/Land

max_forks_repo_name: climate-machine/Land