text:
# knowledge graph notes for mathematics
import KnowledgeGraphNote.Concept

math_notes = [
    Concept("radon measure",
            "https://mathworld.wolfram.com/RadonMeasure.html",
            ["borel measure", "compact set"],
            raw"a radon measure is a Borel measure that is finite on compact sets."
            ),

    Concept("Borel measure",
            "https://mathworld.wolfram.com/BorelMeasure.html",
            ["borel sigma-algebra", "topological space",  "measure", "measurable"],
            raw"if f is the borel sigma-algebra on some topological space, then a measure m:f->r is said to be a borel measure (or borel probability measure). for a borel measure, all continuous functions are measurable."
            ),

    Concept("borel sigma-algebra",
            "https://mathworld.wolfram.com/BorelSigma-Algebra.html",
            ["sigma-algebra", "topological space", "set", "open set", "closed set"],
            raw"a sigma-algebra which is related to the topology of a set. the borel sigma-algebra is defined to be the sigma-algebra generated by the open sets (or equivalently, by the closed sets)."
            ),

    Concept("topological space",
            "https://mathworld.wolfram.com/TopologicalSpace.html",
            ["empty set", "set", "open set"],
            raw"""
A topological space, also called an abstract topological space, is a set X together with a collection of open subsets T that satisfies the four conditions:

1. The empty set emptyset is in T.

2. X is in T.

3. The intersection of a finite number of sets in T is also in T.

4. The union of an arbitrary number of sets in T is also in T.

Alternatively, T may be defined to be the closed sets rather than the open sets, in which case conditions 3 and 4 become:

3. The intersection of an arbitrary number of sets in T is also in T.

4. The union of a finite number of sets in T is also in T.

These axioms are designed so that the traditional definitions of open and closed intervals of the real line continue to be true. For example, the restriction in (3) can be seen to be necessary by considering  $\cap_{n=1}^{\infty}(-1/n,1/n)={0}$, where an infinite intersection of open intervals is a closed set.

In the chapter "Point Sets in General Spaces" Hausdorff (1914) defined his concept of a topological space based on the four Hausdorff axioms (which in modern times are not considered necessary in the definition of a topological space)."""
            ),

    Concept("sigma-algebra",
            "https://mathworld.wolfram.com/Sigma-Algebra.html",
            ["set", "power set"],
            raw"""
Let X be a set. Then a sigma-algebra F is a nonempty collection of subsets of X such that the following hold:

1. X is in F.

2. If A is in F, then so is the complement of A.

3. If $A_n$ is a sequence of elements of F, then the union of the $A_n$s is in F.

If S is any collection of subsets of X, then we can always find a sigma-algebra containing S, namely the power set of X. By taking the intersection of all sigma-algebras containing S, we obtain the smallest such sigma-algebra. We call the smallest sigma-algebra containing S the sigma-algebra generated by S.
"""
            ),
    Concept("probability simplex",
            "https://web.stanford.edu/~boyd/cvxbook/bv_cvxbook.pdf",
            ["simplex"],
            raw"""
The probability simplex is the (n − 1) dimensional simplex determined by the unit vectors e1, . . . , en ∈ Rn. It is the set of vectors that satisfy
x >= 0, 1^T x <= 1.
Vectors in the probability simplex correspond to probability distributions on a set with n elements, with xi interpreted as the probability of the ith element.
    """
            ),
    Concept("simplex",
            "https://web.stanford.edu/~boyd/cvxbook/bv_cvxbook.pdf",
            ["affinely independent", "linearly independent"],
            raw"""
Simplexes are another important family of polyhedra. Suppose the k + 1 points v0, . . . , vk ∈ Rn are affinely independent, which means v1 − v0, . . . , vk − v0 are linearly independent. The simplex determined by them is given by
C = conv{v0, . . . , vk} = {θ0v0 + · · · + θkvk | θ >= 0, 1^Tθ = 1}
where 1 denotes the vector with all entries one. The affine dimension of this simplex is k, so it is sometimes referred to as a k-dimensional simplex in Rn
                    """
            ),
    Concept("linearly independent",
            "https://mathworld.wolfram.com/LinearlyDependentVectors.html",
            [],
            raw"""
n vectors $X_1$, $X_2$, ..., $X_n$ are linearly dependent iff there exist scalars $c_1$, $c_2$, ..., $c_n$, not all zero, such that
 $sum_(i=1)^nc_iX_i=0$.
If no such scalars exist, then the vectors are said to be linearly independent.
                    """
            ),
    Concept("set",
            "https://mathworld.wolfram.com/Set.html",
            [],
            raw"A set is a finite or infinite collection of objects in which order has no significance, and multiplicity is generally also ignored (unlike a list or multiset). Members of a set are often referred to as elements and the notation a in A is used to denote that a is an element of a set A. The study of sets and their properties is the object of set theory."
            ),
    Concept("open set",
            "https://mathworld.wolfram.com/OpenSet.html",
            ["metric space", "neighborhood", ],
            raw"""
Let S be a subset of a metric space. Then the set S is open if every point in S has a neighborhood lying in the set. An open set of radius r and center $x_0$ is the set of all points x such that $|x-x_0|<r$, and is denoted $D_r(x_0)$. In one-space, the open set is an open interval. In two-space, the open set is a disk. In three-space, the open set is a ball.
"""
            ),
    Concept("Metric Space",
            "https://mathworld.wolfram.com/MetricSpace.html",
            ["metric"],
            raw"""
A metric space is a set S with a global distance function (the metric g) that, for every two points x,y in S, gives the distance between them as a nonnegative real number g(x,y). A metric space must also satisfy

1. g(x,y)=0 iff x=y,

2. g(x,y)=g(y,x),

3. The triangle inequality g(x,y)+g(y,z)>=g(x,z).
"""
            ),
    Concept("Metric",
            "https://mathworld.wolfram.com/Metric.html",
            [],
            raw"""
A nonnegative function g(x,y) describing the "distance" between neighboring points for a given set. A metric satisfies the triangle inequality

 g(x,y)+g(y,z)>=g(x,z) 	(1)
and is symmetric, so

 g(x,y)=g(y,x). 	(2)
A metric also satisfies

 g(x,x)=0, 	(3)
as well as the condition that g(x,y)=0 implies x=y. If this latter condition is dropped, then g(x,y) is called a pseudometric instead of a metric.
"""
            ),
    Concept("Tensor",
            "https://mathworld.wolfram.com/Tensor.html",
            [],
            raw"""
An nth-rank tensor in m-dimensional space is a mathematical object that has n indices and m^n components and obeys certain transformation rules. Each index of a tensor ranges over the number of dimensions of space. However, the dimension of the space is largely irrelevant in most tensor equations (with the notable exception of the contracted Kronecker delta). Tensors are generalizations of scalars (that have no indices), vectors (that have exactly one index), and matrices (that have exactly two indices) to an arbitrary number of indices.

Tensors provide a natural and concise mathematical framework for formulating and solving problems in areas of physics such as elasticity, fluid mechanics, and general relativity.

The notation for a tensor is similar to that of a matrix (i.e., $A=(a_(ij))$), except that a tensor $a_(ijk...)$, $a^(ijk...)$, $a_i^(jk)$..., etc., may have an arbitrary number of indices. In addition, a tensor with rank r+s may be of mixed type (r,s), consisting of r so-called "contravariant" (upper) indices and s "covariant" (lower) indices. Note that the positions of the slots in which contravariant and covariant indices are placed are significant so, for example, $a_(munu)^lambda$ is distinct from $a_mu^(nulambda)$.
"""
            ),
    Concept("Closed set",
            "https://mathworld.wolfram.com/ClosedSet.html",
            ["set", "metric space", "set closure", "neighborhood", "Limit Point"],
            raw"""
There are several equivalent definitions of a closed set. Let S be a subset of a metric space. A set S is closed if

1. The complement of S is an open set,

2. S is its own set closure,

3. Sequences/nets/filters in S that converge do so within S,

4. Every point outside S has a neighborhood disjoint from S.

The point-set topological definition of a closed set is a set which contains all of its limit points. Therefore, a closed set C is one for which, whatever point x is picked outside of C, x can always be isolated in some open set which doesn't touch C.
"""
            ),
    Concept("measurable",
            "https://mathworld.wolfram.com/MeasurableSet.html",
            ["sigma-algebra"],
            raw"""
If F is a sigma-algebra and A is a subset of X, then A is called measurable if A is a member of F. X need not have, a priori, a topological structure. Even if it does, there may be no connection between the open sets in the topology and the given sigma-algebra.
"""
            )
]

# the library will ignore case differences in names and dependency names.
duplicate_notes = [
    Concept("a", "empty link", [], "empty text"),
    Concept("b", "empty link", [], "empty text"),
    Concept("A", "empty link", [], "empty text"),
]

# circular definition: https://en.wikipedia.org/wiki/Circular_definition
# avoid them
cyclic_notes = [
    Concept("a", "empty link", ["B"], "empty text"),
    Concept("b", "empty link", ["c"], "empty text"),
    Concept("c", "empty link", ["a"], "empty text"),
]

#using the example from https://eli.thegreenplace.net/2015/directed-graph-traversal-orderings-and-applications-to-data-flow-analysis/
dfs_postordering_test1 = [
    Concept("A", "", ["T", "B", "C"], ""),
    Concept("B", "", ["D"], ""),
    Concept("C", "", ["B", "E"], ""),
    Concept("D", "", [], ""),
    Concept("E", "", ["D"], ""),
    Concept("T", "", ["B"], ""),
]

#using the example from https://en.wikipedia.org/wiki/Depth-first_search#Vertex_orderings
dfs_postordering_test2 = [
    Concept("A", "", ["B", "C"], ""),
    Concept("B", "", ["D"], ""),
    Concept("C", "", ["B", "E"], ""),
    Concept("D", "", [], ""),
]

dfs_postordering_test3 = [
    Concept("A", "", ["B", "C"], ""),
    Concept("B", "", ["E"], ""),
    Concept("C", "", ["D"], ""),
    Concept("D", "", ["F"], ""),
    Concept("E", "", ["F"], ""),
    Concept("F", "", [], ""),
]

max_stars_repo_name: bluepine/KnowledgeGraphNote

max_issues_repo_name: bluepine/KnowledgeGraphNote

max_forks_repo_name: bluepine/KnowledgeGraphNote


text:
#=

  https://edu.swi-prolog.org/mod/assign/view.php?id=249
  """
  Medical diagnosis

  Develop an expert system for medical diagnosis.

  Consider three diseases: flu, gastroenteritis and bronchitis.

  A priori, flu has probability 0.3, gastroenteritis 0 0.2, and bronchitis 0.25.

  If you have the flu, you may have the following symptons, associated with their
  probabilities (symptoms are not mutually exclusive):

    fever, 0.8
    cough 0.6
    sore throat 0.5
    headaches 0.4
    aches 0.7

  If you have gastroenteritis, you may have the following symptons, associated with their
  probabilities (symptoms are not mutually exclusive):

    diarrhea 0.8
    aches 0.7
    nausea 0.4
    fatigue 0.3

  If you have bronchitis, you may have the following symptons, associated with their probabilities
  (symptoms are not mutually exclusive):

    cough 0.8
    fatigue 0.7
    fever 0.3

  Compute the probability of each disease given that the patient has the symptoms fever and aches.

  Do the same supposing the patient adds that he also experiences fatigue.
  """

  Cf ~/cplint/course_medical_diagnosis.pl
     ~/blog/medical_diagnosis.blog
     ~/psi/medical_diagnosis.psi
     ~/webppl/medical_diagnosis.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function medical_diagnosis()
    #=
      Consider three diseases: flu, gastroenteritis and bronchitis.

      A priori, flu has probability 0.3, gastroenteritis 0.2, and bronchitis 0.25.
    =#

    flu ~ flip(0.3)
    gastroenteritis ~ flip(0.2)
    bronchitis ~ flip(0.25)

    #=
      If you have the flu, you may have the following symptons, associated with their
      probabilities (symptoms are not mutually exclusive):

      fever, 0.8
      cough 0.6
      sore throat 0.5
      headaches 0.4
    aches 0.7
    =#

    fever ~ flu ? flip(0.8) : ( bronchitis ? flip(0.3) : flip(0.0))
    cough ~ flu ? flip(0.6) : ( bronchitis ? flip(0.8) : flip(0.0))
    sore_throat ~ flu ? flip(0.5) : flip(0.0)
    headaches ~ flu ? flip(0.4) : flip(0.0)
    aches ~ flu ? flip(0.7) : (gastroenteritis ? flip(0.7) : flip(0.0))

    #=
      If you have gastroenteritis, you may have the following symptons, associated with their
      probabilities (symptoms are not mutually exclusive):

      diarrhea 0.8
      aches 0.7
      nausea 0.4
      fatigue 0.3
    =#

    diarrhea ~ gastroenteritis ? flip(0.8) : flip(0.0)
    nausea  ~ gastroenteritis ? flip(0.4) : flip(0.0)
    fatigue ~ gastroenteritis ? flip(0.3) : ( bronchitis ? flip(0.7) : flip(0.0))

    #=
      If you have bronchitis, you may have the following symptons, associated with their probabilities
      (symptoms are not mutually exclusive):

      cough 0.8
      fatigue 0.7
      fever 0.3
    =#
    cough ~ bronchitis ? flip(0.8) : flip(0.0)
    fatigue ~ bronchitis ? flip(0.7) : flip(0.0)
    fever ~ bronchitis ? flip(0.3) : flip(0.0)

    # The observations
    true ~ Dirac(fever)
    true ~ Dirac(aches)

end

model = medical_diagnosis()
num_chns = 4

# HH has problem with this!
# chns = sample(model, MH(), MCMCThreads(), 10_000, num_chns)
# chns = sample(model, MH(), 40_000)

# chns = sample(model, PG(20), MCMCThreads(), 1000, num_chns)
chns = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chns)

# Note: IS don't generate chns the same way as MH, PG, and SMC!
# chns = sample(model, IS(), MCMCThreads(), 1000, num_chns)

display(chns)

show_var_dist_pct(chns,:flu)
show_var_dist_pct(chns,:gastroenteritis)
show_var_dist_pct(chns,:bronchitis)

ext: jl

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
export DQNLearner

mutable struct DQNLearner{Tq,Tt,Tf,R<:AbstractRNG} <: Any
    approximator::Tq
    target_approximator::Tt
    loss_func::Tf
    min_replay_history::Int
    update_freq::Int
    update_step::Int
    target_update_freq::Int
    sampler::NStepBatchSampler
    rng::R
    # for logging
    loss::Float32
    is_enable_double_DQN::Bool
end

"""
    DQNLearner(;kwargs...)

See paper: [Human-level control through deep reinforcement learning](https://www.nature.com/articles/nature14236)

# Keywords

- `approximator`::[`AbstractApproximator`](@ref): used to get Q-values of a state.
- `target_approximator`::[`AbstractApproximator`](@ref): similar to `approximator`, but used to estimate the target (the next state).
- `loss_func`: the loss function.
- `γ::Float32=0.99f0`: discount rate.
- `batch_size::Int=32`
- `update_horizon::Int=1`: length of update ('n' in n-step update).
- `min_replay_history::Int=32`: number of transitions that should be experienced before updating the `approximator`.
- `update_freq::Int=4`: the frequency of updating the `approximator`.
- `target_update_freq::Int=100`: the frequency of syncing `target_approximator`.
- `stack_size::Union{Int, Nothing}=4`: use the recent `stack_size` frames to form a stacked state.
- `traces = SARTS`: set to `SLARTSL` if you are to apply to an environment of `FULL_ACTION_SET`.
- `rng = Random.GLOBAL_RNG`
- `is_enable_double_DQN = Bool`: enable double dqn, enabled by default.
"""
function DQNLearner(;
    approximator::Tq,
    target_approximator::Tt,
    loss_func::Tf,
    stack_size::Union{Int,Nothing} = nothing,
    γ::Float32 = 0.99f0,
    batch_size::Int = 32,
    update_horizon::Int = 1,
    min_replay_history::Int = 32,
    update_freq::Int = 1,
    target_update_freq::Int = 100,
    traces = SARTS,
    update_step = 0,
    rng = Random.GLOBAL_RNG,
    is_enable_double_DQN::Bool = true,
) where {Tq,Tt,Tf}
    copyto!(approximator, target_approximator)
    sampler = NStepBatchSampler{traces}(;
        γ = γ,
        n = update_horizon,
        stack_size = stack_size,
        batch_size = batch_size,
    )
    DQNLearner(
        approximator,
        target_approximator,
        loss_func,
        min_replay_history,
        update_freq,
        update_step,
        target_update_freq,
        sampler,
        rng,
        0.0f0,
        is_enable_double_DQN,
    )
end

Functors.functor(x::DQNLearner) = (Q = x.approximator, Qₜ = x.target_approximator),
y -> begin
    x = @set x.approximator = y.Q
    x = @set x.target_approximator = y.Qₜ
    x
end

function (learner::DQNLearner)(env)
    env |>
    state |>
    x ->
        Flux.unsqueeze(x, ndims(x) + 1) |>
        x ->
            send_to_device(device(learner), x) |>
            learner.approximator |>
            vec |>
            send_to_host
end

function RLBase.update!(learner::DQNLearner, batch::NamedTuple)
    Q = learner.approximator
    Qₜ = learner.target_approximator
    γ = learner.sampler.γ
    loss_func = learner.loss_func
    n = learner.sampler.n
    batch_size = learner.sampler.batch_size
    is_enable_double_DQN = learner.is_enable_double_DQN
    D = device(Q)

    s, a, r, t, s′ = (send_to_device(D, batch[x]) for x in SARTS)
    a = CartesianIndex.(a, 1:batch_size)

    if is_enable_double_DQN
        q_values = Q(s′)
    else
        q_values = Qₜ(s′)
    end

    if haskey(batch, :next_legal_actions_mask)
        l′ = send_to_device(D, batch[:next_legal_actions_mask])
        q_values .+= ifelse.(l′, 0.0f0, typemin(Float32))
    end

    if is_enable_double_DQN
        selected_actions = dropdims(argmax(q_values, dims = 1), dims = 1)
        q′ = Qₜ(s′)[selected_actions]
    else
        q′ = dropdims(maximum(q_values; dims = 1), dims = 1)
    end

    G = r .+ γ^n .* (1 .- t) .* q′

    gs = gradient(params(Q)) do
        q = Q(s)[a]
        loss = loss_func(G, q)
        ignore() do
            learner.loss = loss
        end
        loss
    end

    update!(Q, gs)
end

max_stars_repo_name: jbrea/ReinforcementLearning

max_issues_repo_name: jbrea/ReinforcementLearning

max_forks_repo_name: jbrea/ReinforcementLearning


text:
# Lebesgue measure

export Lebesgue
struct Lebesgue{X} <: AbstractMeasure end

Lebesgue(X) = Lebesgue{X}()

basemeasure(μ::Lebesgue{X}) where {X} = μ

isprimitive(::Lebesgue) = true

sampletype(::Lebesgue{X}) where{X} = X


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _hadec2altaz(ha::T, dec::T, lat::T, ws::Bool) where {T<:AbstractFloat}
    sh, ch = sincos(deg2rad(ha))
    sd, cd = sincos(deg2rad(dec))
    sl, cl = sincos(deg2rad(lat))

    x = -ch*cd*sl + sd*cl
    y = -sh*cd
    z = ch*cd*cl + sd*sl
    r = hypot(x, y)

    # Now get altitude, azimuth
    az  = rad2deg(mod2pi(atan(y, x)))
    alt = rad2deg(atan(z, r))
    # Convert azimuth to West from South, if desired
    if ws
        az = mod(az + 180, 360)
    end
    return alt, az
end

"""
    hadec2altaz(ha, dec, lat[, ws=true]) -> alt, az

### Purpose ###

Convert Hour Angle and Declination to Horizon (Alt-Az) coordinates.

### Explanation ###

Can deal with the NCP singularity.  Intended mainly to be used by program
`eq2hor`.

### Arguments ###

Input coordinates may be either a scalar or an array, of the same dimension.

* `ha`: the local apparent hour angle, in degrees.  The hour angle is the time
  that right ascension of 0 hours crosses the local meridian.  It is
  unambiguously defined.
* `dec`: the local apparent declination, in degrees.
* `lat`: the local geodetic latitude, in degrees, scalar or array.
* `ws` (optional boolean keyword): if true, the output azimuth is measured West
  from South.  The default is to measure azimuth East from North.

`ha` and `dec` can be given as a 2-tuple `(ha, dec)`.

### Output ###

2-tuple `(alt, az)`

* `alt`: local apparent altitude, in degrees.
* `az`: the local apparent azimuth, in degrees.

The output coordinates are always floating points and have the same type (scalar
or array) as the input coordinates.

### Example ###

Arcturus is observed at an apparent hour angle of 336.6829 and a declination of
19.1825 while at the latitude of +43° 4' 42''.  What are the local altitude and
azimuth of this object?

### Notes ###

`altaz2hadec` converts Horizon (Alt-Az) coordinates to Hour Angle and
Declination.

Code of this function is based on IDL Astronomy User's Library.
"""
hadec2altaz(ha::Real, dec::Real, lat::Real; ws::Bool=false) =
    _hadec2altaz(promote(float(ha), float(dec), float(lat))..., ws)

hadec2altaz(hadec::Tuple{Real, Real}, lat::Real; ws::Bool=false) =
    hadec2altaz(hadec..., lat, ws=ws)

function hadec2altaz(ha::AbstractArray{R}, dec::AbstractArray{<:Real},
                     lat::AbstractArray{<:Real}; ws::Bool=false) where {R<:Real}
    @assert length(ha) == length(dec) == length(lat)
    typeha = float(R)
    alt = similar(ha, typeha)
    az  = similar(ha, typeha)
    for i in eachindex(ha)
        alt[i], az[i] = hadec2altaz(ha[i], dec[i], lat[i], ws=ws)
    end
    return alt, az
end


text:
using TensorBoardLogger, ValueHistories, DelimitedFiles

function load_ho_values(path)
    (data, header) = readdlm(path, ','; header=true)
end

function loaddata(path; kwargs...)
    hist = MVHistory()
    TensorBoardLogger.map_summaries(path; kwargs...) do tag, step, val
        push!(hist, Symbol(tag), step, val)
    end
    hist
end

function getvalues(data, name; smoothing_window=100, sample_range)
    values = data[name]
    smoothed = smooth(values.values, smoothing_window)
    return smoothed[sample_range]
end

function gettime(data)
    data["training/reward"].iterations / (24 * 3600)
end

function smooth(x, window)
    y = similar(x)
    total = 0
    for i in eachindex(x)
        total += x[i]
        if i > window
            total -= x[i-window]
            y[i] = total / window
        else
            y[i] = total / i
        end
    end
    y
end

function read_all(exp_path)
    dict = Dict{String,String}()
    for (root, dirs, files) in walkdir(exp_path)
        for d in dirs
            if d == "tb_log"
                dict[splitpath(root)[end-1]] = joinpath(root, d)
            end
        end
    end
    dict
end

max_stars_repo_name: albheim/ServiceMeshControl

max_issues_repo_name: albheim/ServiceMeshControl

max_forks_repo_name: albheim/ServiceMeshControl


text:
#InspectDR: Top-level interface for Cairo layer
#-------------------------------------------------------------------------------

#Render base annotation (titles, axis labels, ...)
function render_baseannotation(ctx::CairoContext, rplot::RPlot2D, plot::Plot2D)
	lyt = plot.layout.values #WANTCONST
Cairo.save(ctx)
	drawrectangle(ctx, rplot.bb, lyt.frame_canvas)
Cairo.restore(ctx)
	render_baseannotation(ctx, rplot, lyt, plot.annotation)
end

#Render axes & grid lines (graph area)
function render_axes(ctx::CairoContext, rplot::RPlot2D, plot::Plot2D)
	lyt = plot.layout.values #WANTCONST

	#Render titles, axis labels, ...
	render_baseannotation(ctx, rplot, plot.layout.values, plot.annotation)

	nstrips = length(rplot.strips)
	for rstrip in rplot.strips
		render_xticklabels = (rstrip.istrip == nstrips) #Only bottom-most graph
	Cairo.save(ctx)
		drawrectangle(ctx, rstrip.bb, lyt.frame_data)
	Cairo.restore(ctx)
		render_grid(ctx, rstrip, lyt, rstrip.grid)
		render_axes(ctx, rstrip, lyt, render_xticklabels)
	end
	return
end

#Render actual plot data
#-------------------------------------------------------------------------------
function render_data(ctx::CairoContext, rstrip::RStrip2D, plot::Plot2D)
Cairo.save(ctx)
	setclip(ctx, rstrip.bb)

	render(ctx, rstrip, plot.display_data_heat)
	render(ctx, rstrip, plot.display_data)
	if plot.displayNaN
		rendernans(ctx, rstrip, plot.data)
	end
Cairo.restore(ctx)
end
function render_data(ctx::CairoContext, rplot::RPlot2D, plot::Plot2D)
	for rstrip in rplot.strips
		render_data(ctx, rstrip, plot)
	end
end

#Render secondary plot annotation & redraw frame:
#-------------------------------------------------------------------------------
function render_userannotation(ctx::CairoContext, rstrip::RStrip2D, plot::Plot2D)
Cairo.save(ctx)
	setclip(ctx, rstrip.bb)
	render(ctx, rstrip, plot.userannot)
	render(ctx, rstrip, plot.parentannot)
Cairo.restore(ctx)

	render_graphframe(ctx, rstrip, plot.layout.values.frame_data) #Redraw frame
	#TODO: should we be more careful about setclip() - instead of redrawing frame?
end
function render_userannotation(ctx::CairoContext, rplot::RPlot2D, plot::Plot2D)
	for rstrip in rplot.strips
		render_userannotation(ctx, rstrip, plot)
	end
end

#Render entire plot using provided context (static images):
#-------------------------------------------------------------------------------
#TODO: Should base API be written such that user provides databb, and have
#      InspectDR calculate outwards instead?
function render(ctx::CairoContext, plot::Plot2D, bb::BoundingBox)
	refresh!(plot.layout) #Copy in values from defaults
	update_ddata(plot) #Also computes new extents
	rplot = RPlot2D(plot, bb)

	_reset(ctx) #Ensure known state.
	render_baseannotation(ctx, rplot, plot)
	render_axes(ctx, rplot, plot)
	render_legends(ctx, rplot, plot)
	render_data(ctx, rplot, plot)
	render_userannotation(ctx, rplot, plot)

	return rplot #So parent/caller can position objects/format data
end

#Render entire plot plot using CairoBufferedPlot (GUI optimization):
#-------------------------------------------------------------------------------
#TODO: split in parts so that GUI can be refreshed before refreshing data
#TODO: Or maybe send in final context so it can be done here.
function render(bplot::CairoBufferedPlot, plot::Plot2D, bb::BoundingBox, refreshdata::Bool)
	refresh!(plot.layout) #Copy in values from defaults
	refreshed = false
	if refreshdata
		update_ddata(plot) #Also computes new extents
		#TODO: don't refresh data if !refresh && bounding box did not change??
	end
	rplot = RPlot2D(plot, bb)

	ctx = CairoContext(bplot.surf)
	_reset(ctx); clear(ctx)
	render_baseannotation(ctx, rplot, plot)
	render_axes(ctx, rplot, plot)
	render_legends(ctx, rplot, plot)

	if refreshdata #Update data surface:
		ctxdata = CairoContext(bplot.data)
		_reset(ctxdata); clear(ctxdata)
		render_data(ctxdata, rplot, plot)
		Cairo.destroy(ctxdata)
	end

	#Now render data portion onto main context:
	Cairo.set_source_surface(ctx, bplot.data, 0, 0)
	Cairo.paint(ctx) #Applies contents of bplot.data

	render_userannotation(ctx, rplot, plot)
	Cairo.destroy(ctx)

	return rplot #So parent/caller can position objects/format data
end

#Last line


text:
module ValidateRCall
using Compat: @warn

# Environment variable to avoid boring R package builds
mustCrossvalidate = haskey(ENV, "JULIA_MUST_CROSSVALIDATE") && ENV["JULIA_MUST_CROSSVALIDATE"] == "1"

# Only run R on unix or when R is installed because JULIA_MUST_CROSSVALIDATE is set to 1
global skipR = !mustCrossvalidate &&
    !((VERSION < v"0.7.0-" && is_unix()) ||
      (VERSION >= v"0.7.0-" && Sys.isunix()))
try
    skipR && error("Skipping R testing...")
    using RCall
    global skipR = false
catch
    global skipR = true
    if mustCrossvalidate && VERSION < v"0.7.0-"
        error("R not installed, but JULIA_MUST_CROSSVALIDATE is set")
    else
        @warn "R or appropriate Phylo package not installed, skipping R cross-validation."
    end
end

!skipR && include("run_rcall.jl")

end


text:
module TestMeshGraphics
using Test

using CairoMakie
using CombinatorialSpaces
using CombinatorialSpaces.MeshGraphics
using StaticArrays: SVector

const Point3D = SVector{3,Float64}

CairoMakie.AbstractPlotting.inline!(true)

s = EmbeddedDeltaSet2D(joinpath(@__DIR__, "assets", "square.obj"))
sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s)
subdivide_duals!(sd, Barycenter())

# Test Graphs
###############

for ds in [s, sd]
  fig, ax, ob = wireframe(ds)
  @test fig isa Figure
  wireframe!(ds)
  @test fig isa Figure
  fig, ax, ob = mesh(ds)
  @test fig isa Figure
  mesh!(ds)
  @test fig isa Figure
  fig, ax, ob = scatter(ds)
  @test fig isa Figure
  scatter!(ds)
  @test fig isa Figure
end

end


text:
# time: 2022-03-26 12:42:15 UTC
# mode: julia
	soiefjoisejfseoijf
# time: 2022-03-26 12:42:18 UTC
# mode: julia
	1+1


text:
using RICoverAlg

# Queueing problems
m = 2;  n = 40 ; lmd = 1.0
g = ones(n); g /=sum(g)
U = zeros(m, n)
for i in (1:n)
   U[2,i] = 28 - i
end
# shiftterm = (((maximum(U)+minimum(U))/2))   ## reduce extremes of U for numerical stability
Z = exp.((U.-0)/lmd)

@test convertInputsToFloats(U,g,1) == (U,g,1.)
@test convertInputsToFloats([2 1],g, 1//10) == ([2. 1.],g,.1)
#@test convertInputsToFloats([2; 1],1, 1//10) == ([2.; 1.],1.,.1) # how to convert a number to Array if necessary
#@test cancelZeroStates([4 5; 6 7], [0 1]) == ([5; 7], [1], [2]) # why is this test not working?
@test checkForTriviality([2. 1.],[.1; 0.9],1) == (ones(1), log(2.)*.1+log(1.)*.9, 0, "Trivial: only one action is possible")
@test checkForTriviality([2. ; 1.],[1],2) == ([1. ; 0.], log(2.0), 0, "Trivial: only one state is possible")

@test risolve(U, g, 10)[1] == [0.0, 1.0]
@test risolve(reverse(U,dims=1), g, 10)[3][1:47] == "Found boundary solution (two alternatives only)"

# Matejka & McKay 2015, RBT
U = [0 1 0 1;
    0 0 1 1;
    1/2 1/2 1/2 1/2]
lmd = 0.4
g = [1+1; 1-1; 1-1; 1+1]./4
@test convertInputsToFloats(U,g,lmd) == (U./1.,g,lmd)
@test cancelZeroStates(U,g ) == (U[:,[1,4]],[.5, .5 ],[1, 4])
@test transform_utilities(U[:,[1,4]],lmd) == ( exp.((U[:,[1,4]].-.0)/lmd)   ,  .0   )


text:
# test to make sure our geodeics have constant hamiltonians over time
include("../src/flow_ode.jl")
include("../src/rfuncs.jl")
include("../src/targets.jl")

tmpx = [rand(50), randn(70)/10 + .8]
tmpy = tmpx + rand(size(tmpx)) .* 1.1
N = length(tmpx)

kappa = Array{Float64,1}[[(tmpx[i]- minimum(tmpx))/maximum(tmpx), (tmpy[i]-minimum(tmpy))/maximum(tmpy)] for i=1:N]
eta_coeff = Array{Float64,1}[ randn(size(kappa[i]))/100.0 for i in 1:length(kappa)]

sigma = 0.1

function henergy(kappa, eta_coeff, sigma)
	sum = 0.0
	N = length(kappa)
	for i=1:N, j=1:N
		sum += dot(eta_coeff[i], eta_coeff[j]) * R(kappa[i], kappa[j], sigma)
	end
	sum
end
function forward_test!(ham, kappa, eta_coeff, sigma)
	# this function records the time varying hamiltonian and pushes the values to ham
	stepsODE = 1000
	epsilon = 1.0/stepsODE
	Nkappa = length(kappa)
	# initialize the delta s
	dEtaDt   = Array{Float64,1}[zero(kappa[1]) for i in 1:Nkappa]
	dKappaDt = Array{Float64,1}[zero(kappa[1]) for i in 1:Nkappa]
	for counter = 1:stepsODE
		d_eta_dt!(dEtaDt,  eta_coeff, kappa, sigma, epsilon)
		d_kappa_dt!(dKappaDt, eta_coeff, kappa, sigma, epsilon)
		for k = 1:Nkappa
			eta_coeff[k] +=  dEtaDt[k]
			kappa[k]    +=  dKappaDt[k]
		end
		push!(ham, henergy(kappa, eta_coeff, sigma))
	end
end

ham = Float64[henergy(kappa, eta_coeff, sigma)]
forward_test!(ham, kappa, eta_coeff, sigma)

using PyCall
@pyimport matplotlib.pyplot as plt
plt.plot(push!(ham,0.0))
plt.show()

max_stars_repo_name: EthanAnderes/MleDensityWarp

max_issues_repo_name: EthanAnderes/MleDensityWarp

max_forks_repo_name: EthanAnderes/MleDensityWarp


text:
## Coordinate relaxation  in general

""" Set x[j] s.t. we solve for row i in the
linear system Ax = b. This updates x
in place"""
function rowsolve!(A,b,x,i,j)
  rhs = b[i]
  for k=1:size(A,2)
    if k!=j
      rhs -= A[i,k]*x[k]
    end
  end
  x[j] = rhs/A[i,j]
end
""" Solve Ax=b with coordinate relaxation sweeps """
function solve_linsys_coord_relax(A,b,nsweeps)
  x = copy(b)
  for iter=1:nsweeps
    for j=1:size(A,1)
      rowsolve!(A,b,x,j,j)
    end
  end
  return x
end
A = [2 1; 1 2]
b = [1,1.5]
x = solve_linsys_coord_relax(A,b,5)
##
xtrue = A\b
norm(x-xtrue)
##
for nsteps = 1:10
  x = solve_linsys_coord_relax(A,b,nsteps)
  @show x
end


text:
@inline function ExaTronKernels.dscal(n::Int,da::Float64,dx,incx::Int,tx)

    # Ignore incx for now.
    if tx <= n
        @inbounds dx[tx] = da*dx[tx]
    end
    @synchronize

    return
end


text:
module AoB

using Cropbox
using ..Radish
using TimeZones

ND = @config (
    :Calendar => (
        :init => ZonedDateTime(2007, 9, 1, tz"Asia/Seoul"),
        :last => ZonedDateTime(2008, 8, 31, tz"Asia/Seoul"),
    ),
    :Weather => (:store => Radish.loadwea(Radish.datapath("2007.wea"), tz"Asia/Seoul")),
    :Phenology => (:planting_date => ZonedDateTime(2007, 11, 1, tz"Asia/Seoul")),
)

tz = tz"America/Los_Angeles"

KMSP = @config (
# # CV PHYL ILN GLN LL LER SG SD LTAR LTARa LIR Topt Tceil critPPD
# KM1 134 4 10 65.0 4.70 1.84 122 0 0.4421 0.1003 22.28 34.23 12
    :Phenology => (
        optimal_temperature = 22.28, # Topt
        ceiling_temperature = 34.23, # Tceil
        critical_photoperiod = 12, # critPPD
        #initial_leaves_at_harvest = , # ILN
        maximum_leaf_initiation_rate = 0.1003, # LIR
        # storage_days = 100, # SD
        storage_temperature = 5,
        maximum_phyllochron_asymptote = 0.4421, # LTARa
        leaves_generic = 10, # GLN
    ),
    :Leaf => (
        maximum_elongation_rate = 4.70, # LER
        minimum_length_of_longest_leaf = 65.0, # LL
        # stay_green = , # SG
    ),
    :Carbon => (
# # Rm Yg
# 0.012 0.8
        maintenance_respiration_coefficient = 0.012, # Rm
        synthesis_efficiency = 0.8, # Yg
        partitioning_table = [
        # root leaf sheath scape bulb
          0.00 0.00   0.00  0.00 0.00 ; # seed Radish before germination
          0.35 0.30   0.25  0.00 0.10 ; # vegetative stage between germination and scape initiation
          0.15 0.15   0.10  0.25 0.35 ; # period between scape initiation and scape appearance
          0.05 0.10   0.00  0.35 0.50 ; # period after scape appearance before removal (scape stays intact)
          0.05 0.00   0.00  0.00 0.95 ; # period after scape removal (scape appeared and subsequently removed)
          0.00 0.00   0.00  0.00 0.00 ; # dead
        ],
    ),
)
KM = @config (KMSP, (
    :Phenology => (initial_leaves_at_harvest = 4,), # ILN
    :Leaf => (stay_green = 1.84,), # SG
    :Meta => (cultivar = :KM,),
))
SP = @config (KMSP, (
    :Phenology => (initial_leaves_at_harvest = 6,), # ILN
    :Leaf => (stay_green = 1.47,), # SG
    :Meta => (cultivar = :SP,),
))

CUH = @config (
# # LAT LONG ALT
# 47.66 122.29 20.0
    :Location => (
        latitude = 47.66, # LAT
        longitude = 122.29, #LONG
        altitude = 20.0, # ALT
    ),
# # CO2 timestep
# 390 60
    :Weather => (
        CO2 = 390, # CO2
    ),
    :Plant => (initial_planting_density = 55,), # PD0
)
CUH_2013 = @config (CUH, (
    :Weather => (
        store = Radish.loadwea(Radish.datapath("CUH/2013.wea"), tz),
    ),
    :Calendar => (
        init = ZonedDateTime(2013, 10, 30, tz), # Y1 bgn
        last = ZonedDateTime(2014, 7, 28, tz), #Y2 end
    ),
    :Meta => (
        year = 2013,
    )
))
CUH_2014 = @config (CUH, (
    :Weather => (
        store = Radish.loadwea(Radish.datapath("CUH/2014.wea"), tz),
    ),
    :Calendar => (
        # 2014.wea starts from 2014-09-01 01:00, not 00:00
        init = ZonedDateTime(2014, 9, 1, 1, tz), # Y1 bgn
        last = ZonedDateTime(2015, 7, 7, tz), #Y2 end
    ),
    :Meta => (
        year = 2014,
    )
))

P1 = @config :Meta => (; planting_group = 1)
P2 = @config :Meta => (; planting_group = 2)

CUH_2013_P1 = @config (CUH_2013, P1, (
    :Phenology => (
        storage_days = 122, # SD
        planting_date = ZonedDateTime(2013, 10, 30, tz), # Y1 sow
    ),
))
CUH_2013_P2 = @config (CUH_2013, P2, (
    :Phenology => (
        storage_days = 170, # SD
        planting_date = ZonedDateTime(2013, 12, 17, tz), # Y1 sow
    ),
))
CUH_2014_P1 = @config (CUH_2014, P1, (
    :Phenology => (
        storage_days = 93, # SD
        planting_date = ZonedDateTime(2014, 10, 1, tz), # Y1 sow
    ),
))
CUH_2014_P2 = @config (CUH_2014, P2, (
    :Phenology => (
        storage_days = 143, # SD
        planting_date = ZonedDateTime(2014, 11, 20, tz), # Y1 sow
    ),
))

KM_2013_P1_SR0 = @config (KM, CUH_2013_P1, (
    :Phenology => (
        emergence_date = ZonedDateTime(2013, 12, 29, tz), # Y1 emg
        scape_removal_date = nothing, # Y2 SR
    ),
))
KM_2013_P2_SR0 = @config (KM, CUH_2013_P2, (
    :Phenology => (
        emergence_date = ZonedDateTime(2014, 1, 26, tz), # Y1 emg
        scape_removal_date = nothing, # Y2 SR
    ),
))
KM_2014_P1_SR0 = @config (KM, CUH_2014_P1, (
    :Phenology => (
        emergence_date = ZonedDateTime(2014, 10, 26, tz), # Y1 emg
        scape_removal_date = nothing, # Y2 SR
    ),
))
KM_2014_P2_SR0 = @config (KM, CUH_2014_P2, (
    :Phenology => (
        emergence_date = ZonedDateTime(2014, 12, 30, tz), # Y1 emg
        scape_removal_date = nothing, # Y2 SR
    ),
))

SP_2013_P1_SR0 = @config (SP, CUH_2013_P1, (
    :Phenology => (
        emergence_date = ZonedDateTime(2013, 11, 14, tz), # Y1 emg
        scape_removal_date = nothing, # Y2 SR
    ),
))
SP_2013_P2_SR0 = @config (SP, CUH_2013_P2, (
    :Phenology => (
        emergence_date = ZonedDateTime(2014, 1, 6, tz), # Y1 emg
        scape_removal_date = nothing, # Y2 SR
    ),
))
SP_2014_P1_SR0 = @config (SP, CUH_2014_P1, (
    :Phenology => (
        emergence_date = ZonedDateTime(2014, 10, 6, tz), # Y1 emg
        scape_removal_date = nothing, # Y2 SR
    ),
))
SP_2014_P2_SR0 = @config (SP, CUH_2014_P2, (
    :Phenology => (
        emergence_date = ZonedDateTime(2014, 11, 30, tz), # Y1 emg
        scape_removal_date = nothing, # Y2 SR
    ),
))

using CSV
using DataFrames
using DataFramesMeta
import Dates
loaddata(cv, y, p) = begin
    ps = CSV.File(Radish.datapath("CUH/PhenologySummary.csv")) |> DataFrame |> unitfy
    @chain ps begin
        @subset(:cultivar .== cv, :year .== y, :planting_group .== p)
        @select(:DAP, :leaves)
    end
end

calibrate_LTAR(cv, y, p) = begin
    c = eval(Symbol("$(cv)_$(y)_P$(p)_SR0")) #HACK for now
    obs = loaddata(cv, y, p)
    f(s) = s.DAP' in obs.DAP && Dates.hour(s.calendar.time') == 12
    #r = simulate(Radish.Model;
    calibrate(Radish.Model, obs;
        config=c,
        stop="calendar.count",
        index=:DAP,
        target=:Leaves => :leaves_appeared,
        parameters=:Phenology => :LTARa_max => (0.01, 0.80),
        snap=f,
        optim=(:MaxSteps => 100,),
    )
end
#calibrate_LTAR(:KM, 2014, 2)

end


text:
using FlightMechanics
using FlightMechanics.Models
using Dierckx

##----------------------------------------------------------------------------------------------------
## imports

# Aerodynamics
import FlightMechanics.Models:
    calculate_aerodynamics

# Propulsion
import FlightMechanics.Models:
    get_pfm, get_cj, get_power, get_efficiency, get_tanks,
    get_engine_position, get_engine_orientation, get_engine_gyro_effects,
    calculate_engine

# Aircraft Model
import FlightMechanics.Models:
    get_name, get_wing_area, get_wing_span, get_chord, get_arp,
    get_empty_mass_props, get_payload_mass_props,
    set_stick_lon!, set_stick_lat!, set_pedals!, set_thtl!,
    get_controls_ranges_trimmer

##----------------------------------------------------------------------------------------------------
## exports

# Aerodynamics
export C310Aerodynamics,
    calculate_aerodynamics

# Propulsion
export C310Engine, C310EngineLeft, C310EngineRight,
    get_pfm, get_cj, get_power, get_efficiency, get_tanks,
    get_engine_position, get_engine_orientation,
    calculate_engine

# FCS
export C310FCS,
    set_thtl1!, set_thtl2!

# Aircraft Model
export C310,
    get_name,
    get_payload_mass_props,
    get_empty_mass_props,
    get_wing_area,
    get_wing_span,
    get_chord,
    get_arp

##----------------------------------------------------------------------------------------------------
## structs

# Aerodynamics
struct C310Aerodynamics<:Aerodynamics       # TODO: CDα, CLα, CYβ... could be included here
    wind_pfm::PointForcesMoments
    wind_coeff_pfm::PointForcesMoments
    body_pfm::PointForcesMoments
    body_coeff_pfm::PointForcesMoments
end

# Propulsion
abstract type C310Engine<:Engine end

struct C310EngineLeft<:C310Engine
    pfm::PointForcesMoments
    cj::Number
    power::Number
    efficiency::Number
    tanks::Array{RigidSolid, 1}
    h::Array{T, 1} where T<:Number  # angular momentum [kg·m²/s]
end

struct C310EngineRight<:C310Engine
    pfm::PointForcesMoments
    cj::Number
    power::Number
    efficiency::Number
    tanks::Array{RigidSolid, 1}
    h::Array{T, 1} where T<:Number  # angular momentum [kg·m²/s]
end

# FCS
mutable struct C310FCS<:FCS
    # Cabin controls
    stick_longitudinal::RangeControl
    stick_lateral::RangeControl
    pedals::RangeControl
    thtl::RangeControl

    # Aerodynamic surfaces
    de::RangeControl
    da::RangeControl
    dr::RangeControl
    # Engine
    t1::RangeControl
    t2::RangeControl
    # Configuration
    allow_out_of_range::Bool
    throw_error_on_out_of_range::Bool

end

# Aircraft Model
struct C310<:Aircraft
    mass_props::RigidSolid
    pfm::PointForcesMoments

    aerodynamics::C310Aerodynamics
    propulsion::Propulsion
    fcs::FCS
end

##----------------------------------------------------------------------------------------------------
## Aerodynamics

# Drag
#             α (rad)   cDα
tablecDα = [-3.14      0.01  ;
            -1.51      1.    ;
            -0.785     0.707 ;
            -0.0873    0.0041;
            -0.0698    0.0013;
            -0.0524    0.0001;
            -0.0349    0.0003;
            -0.0175    0.0020;
             0.0000    0.0052;
             0.0175    0.0099;
             0.0349    0.0162;
             0.0524    0.0240;
             0.0698    0.0334;
             0.0873    0.0442;
             0.1047    0.0566;
             0.1222    0.0706;
             0.1396    0.0860;
             0.1571    0.0962;
             0.1745    0.1069;
             0.1920    0.1180;
             0.2094    0.1298;
             0.2269    0.1424;
             0.2443    0.1565;
             0.2618    0.1727;
             0.2793    0.1782;
             0.2967    0.1716;
             0.3142    0.1618;
             0.3316    0.1475;
             0.3491    0.1097;
             0.785     0.707 ;
             1.57      1.00  ;
             3.14      0.01  ]
splcDα = Spline1D(tablecDα[:, 1], tablecDα[:, 2], k=1, bc="nearest")

cDα(aero::C310Aerodynamics, α) = splcDα(α)
cD0(aero::C310Aerodynamics) = 0.0280
cDde(aero::C310Aerodynamics, de) = 0.0
cDβ(aero::C310Aerodynamics, β) = 0.1400 * β

# Side
#              β (rad)    cYβ
tablecYβ = [-0.3490     0.2120;
             0.0000     0.0000;
             0.3490    -0.2120]

splcYβ = Spline1D(tablecYβ[:, 1], tablecYβ[:, 2], k=1, bc="nearest")

cYp(aero::C310Aerodynamics, p) = 0.1410 * p
cYr(aero::C310Aerodynamics, r) = 0.3550 * r
cYdr(aero::C310Aerodynamics, dr) = 0.2300 * dr
cYβ(aero::C310Aerodynamics, β) = splcYβ(β)

# Lift
#             α (rad)   cLα
tablecLα = [-3.1416     0.0   ;
            -2.356      0.97  ;
            -1.57       0.0   ;
            -0.785     -0.970 ;
            -0.3000    -1.0100;
            -0.2820    -1.1720;
            -0.2470    -1.0870;
            -0.2120    -0.9480;
            -0.1780    -0.8100;
             0.0000     0.0000;
             0.1040     0.4760;
             0.1390     0.6160;
             0.1740     0.7540;
             0.2090     0.8570;
             0.2440     0.9470;
             0.2790     0.8990;
             0.3000     0.8060;
             0.3140     0.5690;
             0.3490     0.2880;
             0.6        0.904 ;
             0.785      0.97  ;
             1.57       0.0   ;
             2.356     -0.97  ;
             3.14159	 0     ]
splcLα = Spline1D(tablecLα[:, 1], tablecLα[:, 2], k=1, bc="nearest")

#              α (rad)   cLαde
tablecLde = [-1.57       0.0   ;
              0.0000    -0.8100;
              0.0873    -0.9000;
              0.1152    -0.7700;
              1.57       0.0   ]
splcLde = Spline1D(tablecLde[:, 1], tablecLde[:, 2], k=1, bc="nearest")

#          α (rad)   cLααdot
tablecLαdot = [0.0000    5.3000;
          0.0873    4.5000;
          0.1152    4.1000]
splcLαdot = Spline1D(tablecLαdot[:, 1], tablecLαdot[:, 2], k=1, bc="nearest")

#            α (rad)   cLqα
tablecLq = [0.0000    9.7000;
            0.0873    8.8000;
            0.1152    8.4000]
splcLq = Spline1D(tablecLq[:, 1], tablecLq[:, 2], k=1, bc="nearest")

cL0(aero::C310Aerodynamics) = 0.280
cLα(aero::C310Aerodynamics, α) = splcLα(α)
cLde(aero::C310Aerodynamics, de, α) = splcLde(α) * de
cLαdot(aero::C310Aerodynamics, αdot, α) = splcLαdot(α) * αdot
cLq(aero::C310Aerodynamics, q, α) = splcLq(α) * q

# Roll
#             β (rad)    cLβ
tablecLβ = [-0.3490     0.0382;
             0.0000     0.0000;
             0.3490    -0.0382]
# TODO: not a good idea that Spline is created once and again
splcLβ = Spline1D(tablecLβ[:, 1], tablecLβ[:, 2], k=1, bc="nearest")

clp(aero::C310Aerodynamics, p) = -0.7500 * p
clr(aero::C310Aerodynamics, r) = 0.0729 * r
clda(aero::C310Aerodynamics, da) = 0.1720 * da
cldr(aero::C310Aerodynamics, dr) = 0.0192 * dr
clβ(aero::C310Aerodynamics, β) = splcLβ(β)

# Pitch
cm0(aero::C310Aerodynamics) = 0.070
cmα(aero::C310Aerodynamics, α) = -0.989 * α
cmαdot(aero::C310Aerodynamics, αdot) = -12.7000 * αdot
cmq(aero::C310Aerodynamics, q) =  -80.0000 * q
cmde(aero::C310Aerodynamics, de) = -2.2600 * de

# Yaw
cnβ(aero::C310Aerodynamics, β) = 0.1000 * β
cnp(aero::C310Aerodynamics, p) = -0.0257 * p
cnr(aero::C310Aerodynamics, r) = -0.3000 * r
cnda(aero::C310Aerodynamics, da) = -0.0168 * da
cndr(aero::C310Aerodynamics, dr) = -0.1152 * dr

# FORCES
function calculate_aerodynamics(
    ac::Aircraft, aero::C310Aerodynamics, aerostate::AeroState, state::State
    )

    ARP = get_arp(ac)

    qinf = get_qinf(aerostate)
    Sw = get_wing_area(ac)
    b = get_wing_span(ac)
    c = get_chord(ac)

    fcs = get_fcs(ac)

    de = get_value(fcs.de)  # rad
    da = get_value(fcs.da)  # rad
    dr = get_value(fcs.dr)  # rad

    α = get_alpha(aerostate)  # rad
    β = get_beta(aerostate)    # rad
    tas = get_tas(aerostate)  # m/s
    αdot = get_alpha_dot(aerostate)  # rad/s

    c2v = c / (2*tas)
    b2v = b / (2*tas)

    p, q, r = get_body_ang_velocity(state)  # rad/s

    cD = cD0(aero) + cDα(aero, α) + cDde(aero, de) + cDβ(aero, β)
    cY = cYβ(aero, β) + cYdr(aero, dr) +
         b2v * (cYp(aero, p)  + cYr(aero, r))
    cL = cL0(aero) + cLα(aero, α) + cLde(aero, de, α) +
         c2v * (cLαdot(aero, αdot, α) + cLq(aero, q, α))

    cl = clβ(aero, β) + clda(aero, da) + cldr(aero, dr) +
         b2v * (clp(aero, p) + clr(aero, r))
    cm = cm0(aero) + cmα(aero, α) + cmde(aero, de) +
         c2v * (cmαdot(aero, αdot) + cmq(aero, q))
    cn = cnβ(aero, β) + cnda(aero, da) + cndr(aero, dr) +
         b2v * (cnp(aero, p) + cnr(aero, r))

    pfm_wind = PointForcesMoments(ARP,
                                  qinf*Sw*[-cD, cY, -cL],
                                  qinf*Sw*[b*cl, c*cm, b*cn])

    aerodynamics_from_wind_total(pfm_wind, α, β, qinf, Sw, b, c, aero)
end

##----------------------------------------------------------------------------------------------------
## Propulsion

function C310EngineLeft()
    C310EngineLeft(
        PointForcesMoments(zeros(3), zeros(3), zeros(3)),
        0,
        0,
        0,
        [PointMass(225 * LB2KG, [35, -209.8, 28.3] .* IN2M),
         PointMass(95 * LB2KG, [35, -41.6, 11.7] .* IN2M)],
        [0.0, 0.0, 0.0]
        )
end

function C310EngineRight()
    C310EngineRight(
        PointForcesMoments(zeros(3), zeros(3), zeros(3)),
        0,
        0,
        0,
        [PointMass(225 * LB2KG, [35, 209.8, 28.3] .* IN2M),
         PointMass(95 * LB2KG, [35, 41.6, 11.7] .* IN2M)],
        [0.0, 0.0, 0.0]
        )
end

#Getters
get_pfm(eng::C310Engine) = eng.pfm
get_cj(eng::C310Engine) = eng.cj
get_power(eng::C310Engine) = eng.power
get_efficiency(eng::C310Engine) = eng.efficiency
get_tanks(eng::C310Engine) = eng.tanks

get_engine_position(prop::C310EngineLeft) = [-27.5, -70, 15.5] .* IN2M
get_engine_position(prop::C310EngineRight) = [-27.5, 70, 15.5] .* IN2M
get_engine_orientation(prop::C310Engine) = [0, 0, 0] .* DEG2RAD

function calculate_engine_power(eng::C310Engine, fcs::FCS, aerostate::AeroState, state::State)
   # Max power at sea level
   Pmax0 = 260 * HP2WAT
   # Power assumed proportional to thrust lever
   P0 = Pmax0 * get_thrust_setting(eng, fcs)
   # Assume supercharged engine, no variation of power with altitude until crtic altitude
   Pm = P0 * 1.0
   # Calculate fuel consumption also
   cj = 0.0
   return Pm, cj
end

function calculate_propeller_thrust(
    eng::C310Engine, fcs::FCS, aerostate::AeroState, state::State, Pm::Number
    )
   # Assume constant propulsive efficiency
   ηp = 0.75
   # Calculate thrust
   thrust = ηp * Pm / get_tas(aerostate)
   return thrust, ηp
end

function calculate_engine(
    eng::C310Engine, fcs::FCS, aerostate::AeroState, state::State; consume_fuel=false
    )
    Pm, cj = calculate_engine_power(eng, fcs, aerostate, state)
    thrust, ηp = calculate_propeller_thrust(eng, fcs, aerostate, state, Pm)

    if consume_fuel==true
        error("fuel consumption not implemented")
    else
        tanks = get_tanks(eng)
    end

    # TODO: Torque calcualtion is missing
    pfm = PointForcesMoments(get_engine_position(eng),
                             [thrust, 0, 0],
                             [0, 0, 0]
                             )

    return typeof(eng)(pfm, cj, Pm, ηp, tanks, get_engine_gyro_effects(eng))
end

get_thrust_setting(eng::C310EngineLeft, fcs::FCS) = get_value(fcs.t1)
get_thrust_setting(eng::C310EngineRight, fcs::FCS) = get_value(fcs.t2)

##----------------------------------------------------------------------------------------------------
## FCS

C310FCS(stick_lon, stick_lat, pedals, thtl, de, da, dr, t1, t2) = C310FCS(
    stick_lon, stick_lat, pedals, thtl, de, da, dr, t1, t2, false, false
    )

C310FCS() = C310FCS(
    # Cabin Inputs
    RangeControl(0.0, [0, 1]),  # stick_longitudinal
    RangeControl(0.0, [0, 1]),  # stick_lateral
    RangeControl(0.0, [0, 1]),  # pedals
    RangeControl(0.0, [0, 1]),  # thtl
    # Controls
    RangeControl(0.0, [-25.0, 35.0] .* DEG2RAD),  # elevator
    RangeControl(0.0, [-18.0, 14.0] .* DEG2RAD),  # ailerons
    RangeControl(0.0, [-27.0, 27.0] .* DEG2RAD),  # rudder
    RangeControl(0.0, [0.0, 1.0]),                # t1
    RangeControl(0.0, [0.0, 1.0])                 # t1
    )

function set_stick_lon!(fcs::C310FCS, value)
    set_value!(fcs.stick_longitudinal, value)
    min, max = get_value_range(fcs.de)
    range = max - min
    allow_out_of_range = get_allow_out_of_range_inputs(fcs)
    throw_error = get_throw_error_on_out_of_range_inputs(fcs)
    set_value!(fcs.de, min + range * value, allow_out_of_range, throw_error)
end

function set_stick_lat!(fcs::C310FCS, value)
    set_value!(fcs.stick_lateral, value)
    min, max = get_value_range(fcs.da)
    range = max - min
    allow_out_of_range = get_allow_out_of_range_inputs(fcs)
    throw_error = get_throw_error_on_out_of_range_inputs(fcs)
    set_value!(fcs.da, min + range * value, allow_out_of_range, throw_error)
end

function set_pedals!(fcs::C310FCS, value)
    set_value!(fcs.pedals, value)
    min, max = get_value_range(fcs.dr)
    range = max - min
    allow_out_of_range = get_allow_out_of_range_inputs(fcs)
    throw_error = get_throw_error_on_out_of_range_inputs(fcs)
    set_value!(fcs.dr, min + range * value, allow_out_of_range, throw_error)
end

function set_thtl1!(fcs::C310FCS, value)
    set_value!(fcs.thtl, value)
    min, max = get_value_range(fcs.t1)
    range = max - min
    allow_out_of_range = get_allow_out_of_range_inputs(fcs)
    throw_error = get_throw_error_on_out_of_range_inputs(fcs)
    set_value!(fcs.t1, min + range * value, allow_out_of_range, throw_error)
end

function set_thtl2!(fcs::C310FCS, value)
    set_value!(fcs.thtl, value)
    min, max = get_value_range(fcs.t2)
    range = max - min
    allow_out_of_range = get_allow_out_of_range_inputs(fcs)
    throw_error = get_throw_error_on_out_of_range_inputs(fcs)
    set_value!(fcs.t2, min + range * value, allow_out_of_range, throw_error)
end

function set_thtl!(fcs::C310FCS, value)
    set_thtl1!(fcs, value)
    set_thtl2!(fcs, value)
end

function get_controls_ranges_trimmer(fcs::C310FCS)
    [get_value_range(fcs.stick_longitudinal),
     get_value_range(fcs.stick_lateral),
     get_value_range(fcs.pedals),
     get_value_range(fcs.thtl)]
end

##----------------------------------------------------------------------------------------------------
## Aircraft Model

function C310()
    pfm0 = PointForcesMoments(zeros(3), zeros(3), zeros(3))
    aero0 = C310Aerodynamics(pfm0, pfm0, pfm0, pfm0)
    engine_right = C310EngineRight()
    engine_left = C310EngineLeft()
    propulsion0 = Propulsion(
            PointForcesMoments(zeros(3), zeros(3), zeros(3)),
            0, 0, 0,
            [engine_right, engine_left],
            [0., 0., 0.]
            )

    fcs = C310FCS()
    # mass properties cannot be retrieved until ac is created... so:
    ac = C310(RigidSolid(0, zeros(3), zeros(3, 3)),
              pfm0,
              aero0,
              propulsion0,
              fcs)
    mass = get_fuel_mass_props(get_propulsion(ac)) +
           get_empty_mass_props(ac) +
           get_payload_mass_props(ac)

    C310(mass, pfm0, aero0, propulsion0, fcs)
end

# Name
function get_name(ac::C310)
    return "Cessna 310"
end

# GEOMETRIC PROPERTIES
get_wing_area(ac::C310) = 175.0 * FT2M^2
get_wing_span(ac::C310) = 36.5 * FT2M
get_chord(ac::C310) = 4.9 * FT2M
# Aerodynamic Reference Point
get_arp(ac::C310) = [46, 0, 8.6] .* IN2M

# MASS PROPERTIES
function get_empty_mass_props(ac::C310)
    RigidSolid(
        2950 * LB2KG,                            # Empty mass
        [46, 0, 8.6] .* IN2M,                    # Empty CG
        [8884      0      0;                     # Empty inertia
            0   1939      0;
            0      0  11001]  .* SLUGFT2_2_KGM2
    )
end

# Crew and cargo
get_pilot_mass_props(ac::C310) = PointMass(180 * LB2KG, [37, -14, 24] .* IN2M)
get_copilot_mass_props(ac::C310) = PointMass(180 * LB2KG, [37, 14, 24] .* IN2M)
get_lugage_mass_props(ac::C310) = PointMass(100 * LB2KG, [90, -0, 24] .* IN2M)

get_payload_mass_props(ac::C310) = (get_pilot_mass_props(ac) +
                                    get_copilot_mass_props(ac) +
                                    get_lugage_mass_props(ac)
                                    )

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


text:
using GLM, CSV, DataFrames, PyCall, PyPlot
PyDict(matplotlib."rcParams")["axes.formatter.use_locale"] = true
ihd17r = CSV.File("data/IHME-GBD_2017_DATA-f3f91654-1.csv") |> DataFrame
ihd17paf = CSV.File("data/IHME-GBD_2017_DATA-cba77f73-1.csv") |> DataFrame
gbd17iso3166_ctries = CSV.File("data/gbd17iso3166_ctries.csv") |> DataFrame
ihd17r[:lrate] = log.(ihd17r[:val] ./ 10^5)
ihd17paf[:lnegpaf] = log.(1 .- ihd17paf[:val])
sexes = ["män", "kvinnor"]

function ihd17sex(s)
	r_sex = ihd17r[ihd17r[:sex].==s, :]
	paf_sex = ihd17paf[ihd17paf[:sex].==s, :]
	r_sex_sub = DataFrame(location = r_sex[:location], lrate = r_sex[:lrate])
	paf_sex_sub = DataFrame(location = paf_sex[:location], lnegpaf = paf_sex[:lnegpaf])
    rpaf = join(r_sex_sub, paf_sex_sub, on = :location)
	Dict(:sex => s, :df => join(rpaf, gbd17iso3166_ctries, on = :location))
end

function pafrate_ols(sdict)
	ols = lm(@formula(lrate~lnegpaf), sdict[:df])
	Dict(:ols => ols, :sdict => sdict)
end

function poplot(podict)
	close("all")
	df = podict[:sdict][:df]
	ols = podict[:ols]
	co1lab = round(coef(ols)[1], digits=3)
	co2lab = round(coef(ols)[2]*-1, digits=3)
	reglab = latexstring(replace("$co1lab-$co2lab x", "."=>"{,}"))
	plot(df[:lnegpaf], df[:lrate], ".")
	plot(df[:lnegpaf], predict(ols), "-", label = reglab)
	for(i, iso3166) in enumerate(df[:iso3166])
		annotate(iso3166, (df[:lnegpaf][i], df[:lrate][i]))
	end
	grid(1)
	legend(loc=3)
	title("Ischemisk hjärtsjukdom 2017 50\u201369 $(sexes[podict[:sdict][:sex]])")
	xlabel("log(1\u2212PAF)")
	ylabel("log(dödstal)")
end

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
function test_grid(g::Grid)
    for c in g.cells
        @test c.obj <= 1.0
        @test c.obj >= 0.0
        @test c.color <= 1.0
        @test c.color >= 0.0
    end
    for s in g.starts
        @test s >= 1
        @test s <= length(g.cells)
    end
    for e in g.exits
        @test e >= 1
        @test e <= length(g.cells)
    end
end

@testset "Empty grid" begin
    g = Grid("maps/default.yaml")
    test_grid(g)
end

@testset "Default grid" begin
    g = Grid("maps/default.yaml")
    test_grid(g)
end

@testset "T-maze grid" begin
    g = Grid("maps/tmaze.yaml")
    test_grid(g)
    nfood = 0
    for c in g.cells
        if RoboGrid.object(c) == "food"
            nfood += 1
        end
    end
    @test nfood == 2
end

@testset "Cross grid" begin
    g = Grid("maps/cross.yaml")
    test_grid(g)
end


text:
@testset "affine.jl" begin

    map = NLS_Fit.Map_Affine(1=>2,3=>5)

    @test parameter_size(map) == 2

    θ = Float32[1, 1]
    X_hat = Int[1,2,3]
    X = eval_map(map,X_hat,θ)

    @test eltype(X) == Float32
    @test X ≈ Float32[2,3.5,5]

end


text:
export Disk,𝔻

# T is Real or Complex
# TT is (Real,Real) or Complex
immutable Disk{T,TT} <: BivariateDomain{T}
    center::TT
    radius::Float64
end

Disk(c::Complex,r)=Disk{typeof(c),typeof(c)}(c,r)
Disk{T<:Real}(c::Tuple{T,T},r)=Disk{T,typeof(c)}(c,r)

Disk()=Disk((0.,0.),1.)
Disk(::AnyDomain)=Disk(NaN,(NaN,NaN))

const 𝔻=Disk()

isambiguous(d::Disk)=isnan(d.radius) && all(isnan,d.center)

#canonical is rectangle [r,0]x[-π,π]
# we assume radius and centre are zero for now
fromcanonical{T<:Real}(D::Disk{T},x,t)=D.radius*x*cos(t)+D.center[1],D.radius*x*sin(t)+D.center[2]
tocanonical{T<:Real}(D::Disk{T},x,y)=sqrt((x-D.center[1])^2+(y-D.center[2])^2)/D.radius,atan2(y-D.center[2],x-D.center[1])
checkpoints(d::Disk)=[fromcanonical(d,(.1,.2243));fromcanonical(d,(-.212423,-.3))]

# function points(d::Disk,n,m,k)
#     ptsx=0.5*(1-gaussjacobi(n,1.,0.)[1])
#     ptst=points(PeriodicInterval(),m)
#
#     Float64[fromcanonical(d,x,t)[k] for x in ptsx, t in ptst]
# end

∂(d::Disk)=Circle(Complex(d.center...),d.radius)


text:
using GenericInstruments
GI = GenericInstruments
####################################
# Instantiate obj
niDmm1 = GI.DMM.INSTR(:NI4071,"PXI1Slot5")
# Init
GI.DMM.connect!(niDmm1)
GI.DMM.get_idn(niDmm1)
GI.DMM.disconnect!(niDmm1)

# pending...


text:
using CmdStan, DiffEqBayes, OrdinaryDiffEq, ParameterizedFunctions,
      ModelingToolkit, RecursiveArrayTools, Distributions, Random, Test

# Uncomment for local testing only, make sure MCMCChains and StatsPlots are available
using MCMCChains, StatsPlots

#Random.seed!(123)
cd(@__DIR__)
isdir("tmp") && rm("tmp", recursive=true)

include("../../stan_inference.jl")

println("\nOne parameter case 1\n")
f1 = @ode_def begin
  dx = a*x - x*y
  dy = -3y + x*y
end a
u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5]
prob1 = ODEProblem(f1,u0,tspan,p)
sol = solve(prob1,Tsit5())
t = collect(range(1,stop=10,length=10))
randomized = VectorOfArray([(sol(t[i]) + .5randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)
priors = [truncated(Normal(0.7,1),0.1,2)]

bayesian_result = stan_inference(prob1,t,data,priors;num_samples=2000, nchains=4,
  num_warmup=1000,likelihood=Normal)

sdf  = CmdStan.read_summary(bayesian_result.model)
@test sdf[sdf.parameters .== :theta1, :mean][1] ≈ 1.5 atol=3e-1

# Uncomment for local chain inspection
chn = CmdStan.convert_a3d(bayesian_result.chains, bayesian_result.cnames, Val(:mcmcchains))
plot(chn)
!isdir("tmp") && mkdir("tmp")
savefig("$(@__DIR__)/one_par_case_01.png")

v1 = [sdf[sdf.parameters .== Symbol("u_hat[$i,1]"), :mean][1] for i in 1:10]
v2 = [sdf[sdf.parameters .== Symbol("u_hat[$i,2]"), :mean][1] for i in 1:10]
qa_prey = zeros(10, 3);
qa_pred = zeros(10, 3);
achn = Array(chn);
for i in 1:10
  qa_prey[i, :] = quantile(achn[:, 4+i], [0.055, 0.5, 0.945])
  qa_pred[i, :] = quantile(achn[:, 14+i], [0.055, 0.5, 0.945])
end
p1 = plot(1:10, v1; ribbon=(qa_prey[:, 1], qa_prey[:, 3]), color=:lightgrey, leg=false)
title!("Prey u_hat 89% quantiles")
plot!(v1, lab="prey", xlab="time", ylab="prey", color=:darkred)

p2 = plot(1:10, v2; ribbon=(qa_pred[:, 1], qa_pred[:, 3]), color=:lightgrey, leg=false)
title!("Preditors u_hat 89% quantiles")
plot!(v2, lab="pred", xlab="time", ylab="preditors", color=:darkblue)

plot(p1, p2, layout=(2,1))
savefig("$(@__DIR__)/one_par_pred_prey_01.png")

println("\nOne parameter case 2\n")
priors = [
  truncated(Normal(1.,0.5), 0.1, 3),
  truncated(Normal(1.,0.5), 0.1, 3),
  truncated(Normal(1.5,0.5), 0.1, 4)
]
bayesian_result = stan_inference(prob1,t,data,priors;num_samples=2000, nchains=4,
  num_warmup=1000,likelihood=Normal,sample_u0=true)

sdf  = CmdStan.read_summary(bayesian_result.model)
@test sdf[sdf.parameters .== :theta1, :mean][1] ≈ 1. atol=3e-1
@test sdf[sdf.parameters .== :theta2, :mean][1] ≈ 1. atol=3e-1
@test sdf[sdf.parameters .== :theta3, :mean][1] ≈ 1.5 atol=3e-1

# Uncomment for local chain inspection
chn = CmdStan.convert_a3d(bayesian_result.chains, bayesian_result.cnames, Val(:mcmcchains))
plot(chn)
!isdir("tmp") && mkdir("tmp")
savefig("$(@__DIR__)/one_par_case_02.png")

v1 = [sdf[sdf.parameters .== Symbol("u_hat[$i,1]"), :mean][1] for i in 1:10]
v2 = [sdf[sdf.parameters .== Symbol("u_hat[$i,2]"), :mean][1] for i in 1:10]
qa_prey = zeros(10, 3);
qa_pred = zeros(10, 3);
achn = Array(chn);
for i in 1:10
  qa_prey[i, :] = quantile(achn[:, 8+i], [0.055, 0.5, 0.945])
  qa_pred[i, :] = quantile(achn[:, 18+i], [0.055, 0.5, 0.945])
end
p1 = plot(1:10, v1; ribbon=(qa_prey[:, 1], qa_prey[:, 3]), color=:lightgrey, leg=false)
title!("Prey u_hat 89% quantiles")
plot!(v1, lab="prey", xlab="time", ylab="prey", color=:darkred)

p2 = plot(1:10, v2; ribbon=(qa_pred[:, 1], qa_pred[:, 3]), color=:lightgrey, leg=false)
title!("Preditors u_hat 89% quantiles")
plot!(v2, lab="pred", xlab="time", ylab="preditors", color=:darkblue)

plot(p1, p2, layout=(2,1))
savefig("$(@__DIR__)/one_par_pred_prey_02.png")

println("\nOne parameter case 3\n")

sol = solve(prob1,Tsit5(),save_idxs=[1])
randomized = VectorOfArray([(sol(t[i]) + .01 * randn(1)) for i in 1:length(t)])
data = convert(Array,randomized)
priors = [truncated(Normal(1.5,0.5), 0.1,2)]
bayesian_result = stan_inference(prob1,t,data,priors;num_samples=2000, nchains=4,
  num_warmup=1000,likelihood=Normal,save_idxs=[1])

sdf  = CmdStan.read_summary(bayesian_result.model)
@test sdf[sdf.parameters .== :theta1, :mean][1] ≈ 1.5 atol=3e-1

# Uncomment for local chain inspection
chn = CmdStan.convert_a3d(bayesian_result.chains, bayesian_result.cnames, Val(:mcmcchains))
plot(chn)
savefig("$(@__DIR__)/one_par_case_03.png")

v1 = [sdf[sdf.parameters .== Symbol("u_hat[$i,1]"), :mean][1] for i in 1:10]
v2 = [sdf[sdf.parameters .== Symbol("u_hat[$i,2]"), :mean][1] for i in 1:10]
qa_prey = zeros(10, 3);
qa_pred = zeros(10, 3);
achn = Array(chn);
for i in 1:10
  qa_prey[i, :] = quantile(achn[:, 3+i], [0.055, 0.5, 0.945])
  qa_pred[i, :] = quantile(achn[:, 13+i], [0.055, 0.5, 0.945])
end
p1 = plot(1:10, v1; ribbon=(qa_prey[:, 1], qa_prey[:, 3]), color=:lightgrey, leg=false)
title!("Prey u_hat 89% quantiles")
plot!(v1, lab="prey", xlab="time", ylab="prey", color=:darkred)

p2 = plot(1:10, v2; ribbon=(qa_pred[:, 1], qa_pred[:, 3]), color=:lightgrey, leg=false)
title!("Preditors u_hat 89% quantiles")
plot!(v2, lab="pred", xlab="time", ylab="preditors", color=:darkblue)

plot(p1, p2, layout=(2,1))
savefig("$(@__DIR__)/one_par_pred_prey_03.png")

println("\nOne parameter case 4\n")

priors = [truncated(Normal(1.,0.5), 0.1, 3),
  truncated(Normal(1.5,0.5), 0.1, 4)]
bayesian_result = stan_inference(prob1,t,data,priors;num_samples=2000,
  nchains=4, num_warmup=1000,likelihood=Normal,save_idxs=[1],sample_u0=true)

sdf  = CmdStan.read_summary(bayesian_result.model)
@test sdf[sdf.parameters .== :theta1, :mean][1] ≈ 1. atol=3e-1
@test sdf[sdf.parameters .== :theta2, :mean][1] ≈ 1.5 atol=3e-1

# Uncomment for local chain inspection
chn = CmdStan.convert_a3d(bayesian_result.chains, bayesian_result.cnames, Val(:mcmcchains))
plot(chn)
savefig("$(@__DIR__)/one_par_case_04.png")

v1 = [sdf[sdf.parameters .== Symbol("u_hat[$i,1]"), :mean][1] for i in 1:10]
v2 = [sdf[sdf.parameters .== Symbol("u_hat[$i,2]"), :mean][1] for i in 1:10]
qa_prey = zeros(10, 3);
qa_pred = zeros(10, 3);
achn = Array(chn);
for i in 1:10
  qa_prey[i, :] = quantile(achn[:, 4+i], [0.055, 0.5, 0.945])
  qa_pred[i, :] = quantile(achn[:, 14+i], [0.055, 0.5, 0.945])
end
p1 = plot(1:10, v1; ribbon=(qa_prey[:, 1], qa_prey[:, 3]), color=:lightgrey, leg=false)
title!("Prey u_hat 89% quantiles")
plot!(v1, lab="prey", xlab="time", ylab="prey", color=:darkred)

p2 = plot(1:10, v2; ribbon=(qa_pred[:, 1], qa_pred[:, 3]), color=:lightgrey, leg=false)
title!("Preditors u_hat 89% quantiles")
plot!(v2, lab="pred", xlab="time", ylab="preditors", color=:darkblue)

plot(p1, p2, layout=(2,1))
savefig("$(@__DIR__)/one_par_pred_prey_04.png")


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

## Project Euler
#
#  problems: http://projecteuler.net/problems
#  solutions: https://code.google.com/p/projecteuler-solutions/wiki/ProjectEulerSolutions

#1: 233168
@test sum(filter(n->(n%3==0)|(n%5==0),1:999)) == 233168

#2: 4613732
function euler2(n)
    t, i, j = 0, 1, 2
    while j <= n
          t += j
          i, j = j, i+j
          i, j = j, i+j
          i, j = j, i+j
    end
    return t
end
@test euler2(4000000) == 4613732

#4: 906609
function euler4(n)
    m = 1
    for a = 10^n-1:-1:10^(n-1)
        for b = 10^n-1:-1:max(a,-fld(-m,a))
            p = a*b
            d = digits(p)
            if d == reverse(d) && p > m
                m = p
                b < -fld(-m,a) && break
            end
        end
    end
    return m
end
@test euler4(3) == 906609

#5: 232792560
@test lcm(1:20) == 232792560

#6: 25164150
@test sum(1:100)^2 - sum((1:100).^2) == 25164150

#8: 40824
function euler8(n,m)
    d = digits(n)
    maximum([prod(d[k:k+m-1]) for k=1:length(d)-m+1])
end
let n = 7316717653133062491922511967442657474235534919493496983520312774506326239578318016984801869478851843858615607891129494954595017379583319528532088055111254069874715852386305071569329096329522744304355766896648950445244523161731856403098711121722383113622298934233803081353362766142828064444866452387493035890729629049156044077239071381051585930796086670172427121883998797908792274921901699720888093776657273330010533678812202354218097512545405947522435258490771167055601360483958644670632441572215539753697817977846174064955149290862569321978468622482839722413756570560574902614079729686524145351004748216637048440319989000889524345065854122758866688116427171479924442928230863465674813919123162824586178664583591245665294765456828489128831426076900422421902267105562632111110937054421750694165896040807198403850962455444362981230987879927244284909188845801561660979191338754992005240636899125607176060588611646710940507754100225698315520005593572972571636269561882670428252483600823257530420752963450
    @test euler8(n,5) == 40824
end

#9: 31875000
function euler9(n)
    for a = 1:n, b = 1:n-a, c = n-a-b
        a^2 + b^2 == c^2 && return a*b*c
    end
end
@test euler9(1000) == 31875000

#11: 70600674
function euler11(grid,n)
    m = typemin(eltype(grid))
    for i = n:size(grid,1)-n+1,
        j = n:size(grid,2)-n+1,
        di = -1:1, dj = -1:1
        di == dj == 0 && continue
        idx = sub2ind(size(grid),
                      di==0 ? fill(i,n) : range(i,di,n),
                      dj==0 ? fill(j,n) : range(j,dj,n))
        m = max(m,prod(grid[idx]))
    end
    return m
end
let grid = [
    08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08
    49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00
    81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65
    52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91
    22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80
    24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50
    32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70
    67 26 20 68 02 62 12 20 95 63 94 39 63 08 40 91 66 49 94 21
    24 55 58 05 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72
    21 36 23 09 75 00 76 44 20 45 35 14 00 61 33 97 34 31 33 95
    78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92
    16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57
    86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58
    19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40
    04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66
    88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69
    04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36
    20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16
    20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54
    01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48
]
    @test euler11(grid,4) == 70600674
end

#12: 76576500

#13: 5537376230
let nums = [
    37107287533902102798797998220837590246510135740250
    46376937677490009712648124896970078050417018260538
    74324986199524741059474233309513058123726617309629
    91942213363574161572522430563301811072406154908250
    23067588207539346171171980310421047513778063246676
    89261670696623633820136378418383684178734361726757
    28112879812849979408065481931592621691275889832738
    44274228917432520321923589422876796487670272189318
    47451445736001306439091167216856844588711603153276
    70386486105843025439939619828917593665686757934951
    62176457141856560629502157223196586755079324193331
    64906352462741904929101432445813822663347944758178
    92575867718337217661963751590579239728245598838407
    58203565325359399008402633568948830189458628227828
    80181199384826282014278194139940567587151170094390
    35398664372827112653829987240784473053190104293586
    86515506006295864861532075273371959191420517255829
    71693888707715466499115593487603532921714970056938
    54370070576826684624621495650076471787294438377604
    53282654108756828443191190634694037855217779295145
    36123272525000296071075082563815656710885258350721
    45876576172410976447339110607218265236877223636045
    17423706905851860660448207621209813287860733969412
    81142660418086830619328460811191061556940512689692
    51934325451728388641918047049293215058642563049483
    62467221648435076201727918039944693004732956340691
    15732444386908125794514089057706229429197107928209
    55037687525678773091862540744969844508330393682126
    18336384825330154686196124348767681297534375946515
    80386287592878490201521685554828717201219257766954
    78182833757993103614740356856449095527097864797581
    16726320100436897842553539920931837441497806860984
    48403098129077791799088218795327364475675590848030
    87086987551392711854517078544161852424320693150332
    59959406895756536782107074926966537676326235447210
    69793950679652694742597709739166693763042633987085
    41052684708299085211399427365734116182760315001271
    65378607361501080857009149939512557028198746004375
    35829035317434717326932123578154982629742552737307
    94953759765105305946966067683156574377167401875275
    88902802571733229619176668713819931811048770190271
    25267680276078003013678680992525463401061632866526
    36270218540497705585629946580636237993140746255962
    24074486908231174977792365466257246923322810917141
    91430288197103288597806669760892938638285025333403
    34413065578016127815921815005561868836468420090470
    23053081172816430487623791969842487255036638784583
    11487696932154902810424020138335124462181441773470
    63783299490636259666498587618221225225512486764533
    67720186971698544312419572409913959008952310058822
    95548255300263520781532296796249481641953868218774
    76085327132285723110424803456124867697064507995236
    37774242535411291684276865538926205024910326572967
    23701913275725675285653248258265463092207058596522
    29798860272258331913126375147341994889534765745501
    18495701454879288984856827726077713721403798879715
    38298203783031473527721580348144513491373226651381
    34829543829199918180278916522431027392251122869539
    40957953066405232632538044100059654939159879593635
    29746152185502371307642255121183693803580388584903
    41698116222072977186158236678424689157993532961922
    62467957194401269043877107275048102390895523597457
    23189706772547915061505504953922979530901129967519
    86188088225875314529584099251203829009407770775672
    11306739708304724483816533873502340845647058077308
    82959174767140363198008187129011875491310547126581
    97623331044818386269515456334926366572897563400500
    42846280183517070527831839425882145521227251250327
    55121603546981200581762165212827652751691296897789
    32238195734329339946437501907836945765883352399886
    75506164965184775180738168837861091527357929701337
    62177842752192623401942399639168044983993173312731
    32924185707147349566916674687634660915035914677504
    99518671430235219628894890102423325116913619626622
    73267460800591547471830798392868535206946944540724
    76841822524674417161514036427982273348055556214818
    97142617910342598647204516893989422179826088076852
    87783646182799346313767754307809363333018982642090
    10848802521674670883215120185883543223812876952786
    71329612474782464538636993009049310363619763878039
    62184073572399794223406235393808339651327408011116
    66627891981488087797941876876144230030984490851411
    60661826293682836764744779239180335110989069790714
    85786944089552990653640447425576083659976645795096
    66024396409905389607120198219976047599490197230297
    64913982680032973156037120041377903785566085089252
    16730939319872750275468906903707539413042652315011
    94809377245048795150954100921645863754710598436791
    78639167021187492431995700641917969777599028300699
    15368713711936614952811305876380278410754449733078
    40789923115535562561142322423255033685442488917353
    44889911501440648020369068063960672322193204149535
    41503128880339536053299340368006977710650566631954
    81234880673210146739058568557934581403627822703280
    82616570773948327592232845941706525094512325230608
    22918802058777319719839450180888072429661980811197
    77158542502016545090413245809786882778948721859617
    72107838435069186155435662884062257473692284509516
    20849603980134001723930671666823555245252804609722
    53503534226472524250874054075591789781264330331690
]
    @test sum(digits(sum(nums))[end-9:end].*Int64(10).^(0:9)) == 5537376230
end

#14: 837799
function euler14(m)
    c = zeros(Int,m)
    c[1] = 1
    for n::Int64 = 2:m
        nʹ, d = n, 0
        while nʹ > length(c) || c[nʹ] == 0
            nʹ = iseven(nʹ) ? nʹ>>1 : 3nʹ+1
            d += 1
        end
        d += c[nʹ]
        while n > length(c) || c[n] == 0
            n <= length(c) && (c[n] = d)
            n = iseven(n) ? n>>1 : 3n+1
            d -= 1
        end
    end
    indmax(c)
end
@test euler14(999999) == 837799

#15: 137846528820

#16: 1366
@test sum(digits(big(2)^1000)) == 1366

#17: 21124
#18: 1074
#19: 171

#20: 648
@test sum(digits(factorial(big(100)))) == 648

#21: 31626
#22: 871198282
#23: 4179871
#24: 2783915460
#25: 4782
#26: 983
#27: -59231
#28: 669171001
#29: 9183
#30: 443839
#31: 73682
#32: 45228
#33: 100
#34: 40730
#35: 55
#36: 872187
#37: 748317
#38: 932718654
#39: 840
#40: 210
#41: 7652413
#42: 162
#43: 16695334890
#44: 5482660
#45: 1533776805
#46: 5777
#47: 134043
#48: 9110846700
#49: 296962999629
#50: 997651
#51: 121313
#52: 142857
#53: 4075
#54: 376
#55: 249
#56: 972
#57: 153
#58: 26241
#59: 107359
#60: 26033
#61: 28684
#62: 127035954683
#63: 49
#64: 1322
#65: 272
#66: 661
#67: 7273
#68: 6531031914842725
#69: 510510
#70: 8319823
#71: 428570
#72: 303963552391
#73: 7295372
#74: 402
#75: 161667
#76: 190569291
#77: 71
#78: 55374
#79: 73162890
#80: 40886
#81: 427337
#82: 260324
#83: 425185
#84: 101524
#85: 2772
#86: 1818
#87: 1097343
#88: 7587457
#89: 743
#90: 1217
#91: 14234
#92: 8581146
#93: 1258
#94: 518408346
#95: 14316
#96: 24702
#97: 8739992577
#98: 18769
#99: 709
#100: 756872327473
#101: 37076114526
#102: 228
#103: 20313839404245
#104: 329468
#105: 73702
#106: 21384
#107: 259679
#108: 180180
#109: 38182
#110: 9350130049860600
#111: 612407567715
#112: 1587000
#113: 51161058134250
#114: 16475640049
#115: 168
#116: 20492570929
#117: 100808458960497
#118: 44680
#119: 248155780267521
#120: 333082500
#121: 2269
#122: 1582
#123: 21035
#124: 21417
#125: 2906969179
#126: 18522
#127: 18407904
#128: 14516824220
#129: 1000023
#130: 149253
#131: 173
#132: 843296
#133: 453647705
#134: 18613426663617118
#135: 4989
#136: 2544559
#137: 1120149658760
#138: 1118049290473932
#139: 10057761
#140: 5673835352990
#141: 878454337159
#142: 1006193
#143: 30758397
#144: 354
#145: 608720
#146: 676333270
#147: 846910284
#148: 2129970655314432
#149: 52852124
#150: -271248680
#151: 0.464399
#152: 301
#153: 17971254122360635
#154: 479742450
#155: 3857447
#156: 21295121502550
#157: 53490
#158: 409511334375
#159: 14489159
#160: 16576
#161: 20574308184277971
#162: 3D58725572C62302
#163: 343047
#164: 378158756814587
#165: 2868868
#166: 7130034
#167: 3916160068885
#168: 59206
#169: 178653872807
#170: 9857164023
#171: 142989277
#172: 227485267000992000
#173: 1572729
#174: 209566
#175: 1,13717420,8
#176: 96818198400000
#177: 129325
#178: 126461847755
#179: 986262
#180: 285196020571078987
#181: 83735848679360680
#182: 399788195976
#183: 48861552
#184: 1725323624056
#185: 4640261571849533
#186: 2325629
#187: 17427258
#188: 95962097
#189: 10834893628237824
#190: 371048281
#191: 1918080160
#192: 57060635927998347
#193: 684465067343069
#194: 61190912
#195: 75085391
#196: 322303240771079935
#197: 1.710637717
#198: 52374425
#199: 0.00396087
#200: 229161792008
#201: 115039000
#202: 1209002624
#203: 34029210557338
#204: 2944730
#205: 0.5731441
#206: 1389019170
#207: 44043947822
#208: 331951449665644800
#209: 15964587728784
#210: 1598174770174689458
#211: 1922364685
#212: 328968937309
#213: 330.721154
#214: 1677366278943
#215: 806844323190414
#216: 5437849
#217: 6273134
#218: 0
#219: 64564225042
#220: 139776,963904
#221: 1884161251122450
#222: 1590933
#223: 61614848
#224: 4137330
#225: 2009
#226: 0.11316017
#227: 3780.618622
#228: 86226
#229: 11325263
#230: 850481152593119296
#231: 7526965179680
#232: 0.83648556
#233: 271204031455541309
#234: 1259187438574927161
#235: 1.002322108633
#236: 123/59
#237: 15836928
#238: 9922545104535661
#239: 0.001887854841
#240: 7448717393364181966
#241: 482316491800641154
#242: 997104142249036713
#243: 892371480
#244: 96356848
#245: 288084712410001
#246: 810834388
#247: 782252
#248: 23507044290
#249: 9275262564250418
#250: 1425480602091519
#251: 18946051
#252: 104924.0
#253: 11.492847
#254: 8184523820510
#255: 4.4474011180
#256: 85765680
#257: 139012411
#258: 12747994
#259: 20101196798
#260: 167542057
#261: 238890850232021
#262: 2531.205
#263: 2039506520
#264: 2816417.1055
#265: 209110240768
#266: 1096883702440585
#267: 0.999992836187
#268: 785478606870985
#269: 1311109198529286
#270: 82282080
#271: 4617456485273129588
#272: 8495585919506151122
#273: 2032447591196869022
#274: 1601912348822
#275: 15030564
#276: 5777137137739632912
#277: 1125977393124310
#278: 1228215747273908452
#279: 416577688
#280: 430.088247
#281: 1485776387445623
#282: 1098988351
#283: 28038042525570324
#284: 5a411d7b
#285: 157055.80999
#286: 52.6494571953
#287: 313135496
#288: 605857431263981935
#289: 6567944538
#290: 20444710234716473
#291: 4037526
#292: 3600060866
#293: 2209
#294: 789184709
#295: 4884650818
#296: 1137208419
#297: 2252639041804718029
#298: 1.76882294
#299: 549936643
#300: 8.0540771484375
#301: 2178309
#302: 1170060
#303: 1111981904675169
#304: 283988410192
#305: 18174995535140
#306: 852938
#307: 0.7311720251
#308: 1539669807660924
#309: 210139
#310: 2586528661783
#311: 2466018557
#312: 324681947
#313: 2057774861813004
#314: 132.52756426
#315: 13625242
#316: 542934735751917735
#317: 1856532.8455
#318: 709313889
#319: 268457129
#320: 278157919195482643
#321: 2470433131948040
#322: 999998760323313995
#323: 6.3551758451
#324: 96972774
#325: 54672965
#326: 1966666166408794329
#327: 34315549139516
#328: 260511850222
#329: 199740353/29386561536000
#330: 15955822
#331: 467178235146843549
#332: 2717.751525
#333: 3053105
#334: 150320021261690835
#335: 5032316
#336: CAGBIHEFJDK
#337: 85068035
#338: 15614292
#339: 19823.542204
#340: 291504964
#341: 56098610614277014
#342: 5943040885644
#343: 269533451410884183
#344: 65579304332
#345: 13938
#346: 336108797689259276
#347: 11109800204052
#348: 1004195061
#349: 115384615384614952
#350: 84664213
#351: 11762187201804552
#352: 378563.260589
#353: 1.2759860331
#354: 58065134
#355: 1726545007
#356: 28010159
#357: 1739023853137
#358: 3284144505
#359: 40632119
#360: 878825614395267072
#361: 178476944
#362: 457895958010
#363: 0.0000372091
#364: 44855254
#365: 162619462356610313
#366: 88351299
#367: 48271207
#368: 253.6135092068
#369: 862400558448
#370: 41791929448408
#371: 40.66368097
#372: 301450082318807027
#373: 727227472448913
#374: 334420941
#375: 7435327983715286168
#376: 973059630185670
#377: 732385277
#378: 147534623725724718
#379: 132314136838185
#380: 6.3202e25093
#381: 139602943319822
#382: 697003956
#383: 22173624649806
#384: 3354706415856332783
#385: 3776957309612153700
#386: 528755790
#387: 696067597313468
#388: 831907372805129931
#389: 2406376.3623
#390: 2919133642971
#391: 61029882288
#392: 3.1486734435
#393: 112398351350823112
#394: 3.2370342194
#395: 28.2453753155
#396: 173214653
#397: 141630459461893728
#398: 2010.59096
#399: 1508395636674243,6.5e27330467
#400: 438505383468410633
#401: 281632621
#402: 356019862
#403: 18224771
#404: 1199215615081353
#405: 237696125
#406: 36813.12757207
#407: 39782849136421
#408: 299742733
#409: 253223948
#410:
#411: 9936352
#412: 38788800
#413: 3079418648040719
#414:
#415:
#416:
#417: 446572970925740
#418: 1177163565297340320
#419: 998567458,1046245404,43363922
#420: 145159332
#421:
#422:
#423:
#424:
#425: 46479497324
#426:
#427:
#428:
#429: 98792821
#430: 5000624921.38


text:
# collection formats
const COLL_MULTI = "multi"  # aliased to CSV, as multi is not supported by Requests.jl (https://github.com/JuliaWeb/Requests.jl/issues/140)
const COLL_PIPES = "pipes"
const COLL_SSV = "ssv"
const COLL_TSV = "tsv"
const COLL_CSV = "csv"
const COLL_DLM = Dict{String,String}([COLL_PIPES=>"|", COLL_SSV=>" ", COLL_TSV=>"\t", COLL_CSV=>",", COLL_MULTI=>","])

const DATETIME_FORMATS = (Dates.DateFormat("yyyy-mm-dd HH:MM:SS.sss"), Dates.DateFormat("yyyy-mm-ddTHH:MM:SS.sss"), Dates.DateFormat("yyyy-mm-ddTHH:MM:SSZ"))
const DATE_FORMATS = (Dates.DateFormat("yyyy-mm-dd"),)

const DEFAULT_TIMEOUT_SECS = 5*60

function convert(::Type{DateTime}, str::String)
    # strip off timezone, as Julia DateTime does not parse it
    if '+' in str
        str = split(str, '+')[1]
    end
    # truncate micro/nano seconds to milliseconds, as Julia DateTime does not parse it
    if '.' in str
        uptosec,subsec = split(str, '.')
        if length(subsec) > 3
            str = uptosec * "." * subsec[1:3]
        end
    end
    for fmt in DATETIME_FORMATS
        try
            return DateTime(str, fmt)
        catch
            # try next format
        end
    end
    throw(SwaggerException("Unsupported DateTime format: $str"))
end

function convert(::Type{Date}, str::String)
    for fmt in DATETIME_FORMATS
        try
            return Date(str, fmt)
        catch
            # try next format
        end
    end
    throw(SwaggerException("Unsupported Date format: $str"))
end

struct SwaggerException <: Exception
    reason::String
end

struct ApiException <: Exception
    status::Int
    reason::String
    resp::HTTP.Response

    function ApiException(resp::HTTP.Response; reason::String="")
        isempty(reason) && (reason = get(HTTP.Messages.STATUS_MESSAGES, resp.status, reason))
        new(resp.status, reason, resp)
    end
end

struct Client
    root::String
    headers::Dict{String,String}
    get_return_type::Function   # user provided hook to get return type from response data
    clntoptions::Dict

    function Client(root::String; headers::Dict{String,String}=Dict{String,String}(), get_return_type::Function=(default,data)->default, sslconfig=nothing, require_ssl_verification=true)
        endswith(root, '/') && warn("Root URI ($root) terminates with '/'. Ensure that resource paths do not begin with '/'. This is unconventional.")
        clntoptions = Dict(:sslconfig=>sslconfig, :status_exception=>false, :retries=>0, :require_ssl_verification=>require_ssl_verification)
        new(root, headers, get_return_type, clntoptions)
    end
end

set_user_agent(client::Client, ua::String) = set_header("User-Agent", ua)
set_cookie(client::Client, ck::String) = set_header("Cookie", ck)
set_header(client::Client, name::String, value::String) = (client.headers[name] = value)

struct Ctx
    client::Client
    method::String
    return_type::Type
    resource::String
    auth::Vector{String}

    path::Dict{String,String}
    query::Dict{String,String}
    header::Dict{String,String}
    form::Dict{String,String}
    file::Dict{String,String}
    body::Any
    timeout::Int

    function Ctx(client::Client, method::String, return_type, resource::String, auth, body=nothing; timeout::Int=DEFAULT_TIMEOUT_SECS)
        resource = client.root * resource
        headers = copy(client.headers)
        new(client, method, return_type, resource, auth, Dict{String,String}(), Dict{String,String}(), headers, Dict{String,String}(), Dict{String,String}(), body, timeout)
    end
end

is_json_mime(mime::T) where {T <: AbstractString} = ("*/*" == mime) || occursin(r"(?i)application/json(;.*)?", mime) || occursin(r"(?i)application/(.*)-patch\+json(;.*)?", mime)

function select_header_accept(accepts::Vector{String})
    isempty(accepts) && (return "application/json")
    for accept in accepts
        is_json_mime(accept) && (return accept)
    end
    return join(accepts, ", ")
end

function select_header_content_type(ctypes::Vector{String})
    isempty(ctypes) && (return "application/json")
    for ctype in ctypes
        is_json_mime(ctype) && (return (("*/*" == ctype) ? "application/json" : ctype))
    end
    return ctypes[1]
end

set_header_accept(ctx::Ctx, accepts::Vector{T}) where {T} = set_header_accept(ctx, convert(Vector{String}, accepts))
function set_header_accept(ctx::Ctx, accepts::Vector{String})
    accept = select_header_accept(accepts)
    !isempty(accept) && (ctx.header["Accept"] = accept)
    return nothing
end

set_header_content_type(ctx::Ctx, ctypes::Vector{T}) where {T} = set_header_content_type(ctx, convert(Vector{String}, ctypes))
function set_header_content_type(ctx::Ctx, ctypes::Vector{String})
    ctx.header["Content-Type"] = select_header_content_type(ctypes)
    return nothing
end

set_param(params::Dict{String,String}, name::String, value::Nothing; collection_format=nothing) = nothing

function set_param(params::Dict{String,String}, name::String, value::Union{Nothing,T}; collection_format=nothing) where {T}
    (value === nothing) && return

    if !isa(value, Vector) || isempty(collection_format)
        params[name] = string(value)
    else
        dlm = get(COLL_DLM, collection_format, "")
        isempty(dlm) && throw(SwaggerException("Unsupported collection format $collection_format"))
        params[name] = join(map((x)->string(x), value), dlm)
    end
end

function prep_args(ctx::Ctx)
    kwargs = copy(ctx.client.clntoptions)
    isempty(ctx.file) && (ctx.body === nothing) && isempty(ctx.form) && !("Content-Length" in keys(ctx.header)) && (ctx.header["Content-Length"] = "0")
    isempty(ctx.query) || (kwargs[:query] = ctx.query)
    headers = ctx.header
    body = nothing
    if !isempty(ctx.form)
        headers["Content-Type"] = "application/x-www-form-urlencoded"
        body = HTTP.URIs.escapeuri(ctx.form)
    end
    if !isempty(ctx.file)
        (body === nothing) && (body = Dict())
        idx = 1
        for (_k,_v) in ctx.file
            body["multi$idx"] = HTTP.Multipart(_k, open(_v))
            idx += 1
        end
    end
    if ctx.body !== nothing
        (isempty(ctx.form) && isempty(ctx.file)) || throw(SwaggerException("Can not send both form-encoded data and a request body"))
        if is_json_mime(get(ctx.header, "Content-Type", "application/json"))
            body = to_json(ctx.body)
        elseif ("application/x-www-form-urlencoded" == ctx.header["Content-Type"]) && isa(ctx.body, Dict)
            body = HTTP.URIs.escapeuri(ctx.body)
        elseif isa(ctx.boody, SwaggerModel) && isempty(get(ctx.header, "Content-Type", ""))
            headers["Content-Type"] = "application/json"
            body = to_json(ctx.body)
        else
            body = ctx.body
        end
    end
    # set the timeout
    kwargs[:readtimeout] = ctx.timeout
    return body, kwargs
end

response(::Type{Nothing}, resp::HTTP.Response) = nothing::Nothing
response(::Type{T}, resp::HTTP.Response) where {T <: Real} = response(T, resp.body)::T
response(::Type{T}, resp::HTTP.Response) where {T <: String} = response(T, resp.body)::T
function response(::Type{T}, resp::HTTP.Response) where {T}
    ctype = HTTP.header(resp, "Content-Type", "application/json")
    (length(resp.body) == 0) && return T()
    v = response(T, is_json_mime(ctype) ? JSON.parse(String(resp.body)) : resp.body)
    v::T
end
response(::Type{T}, data::Vector{UInt8}) where {T<:Real} = parse(T, String(data))
response(::Type{T}, data::Vector{UInt8}) where {T<:String} = String(data)::T
response(::Type{T}, data::T) where {T} = data
response(::Type{T}, data) where {T} = convert(T, data)
response(::Type{T}, data::Dict{String,Any}) where {T} = from_json(T, data)::T
response(::Type{T}, data::Dict{String,Any}) where {T<:Dict} = convert(T, data)
response(::Type{Vector{T}}, data::Vector{V}) where {T,V} = [response(T, v) for v in data]

function exec(ctx::Ctx)
    resource_path = replace(ctx.resource, "{format}"=>"json")
    for (k,v) in ctx.path
        resource_path = replace(resource_path, "{$k}"=>v)
    end

    # TODO: use auth_settings for authentication
    body, kwargs = prep_args(ctx)
    if body !== nothing
        resp = HTTP.request(uppercase(ctx.method), HTTP.URIs.URI(resource_path), ctx.header, body; kwargs...)
    else
        resp = HTTP.request(uppercase(ctx.method), HTTP.URIs.URI(resource_path), ctx.header; kwargs...)
    end
    (200 <= resp.status <= 206) || throw(ApiException(resp))

    response(ctx.client.get_return_type(ctx.return_type, resp), resp)
end

name_map(o::T) where {T<:SwaggerModel} = name_map(T)
field_map(o::T) where {T<:SwaggerModel} = field_map(T)

# TODO: will be good to have a comprehensive selector syntax
function get_field(o::T, path...) where {T<:SwaggerModel}
    val = get_field(o, path[1])
    rempath = path[2:end]
    (length(rempath) == 0) && (return val)

    if isa(val, Vector)
        if isa(rempath[1], Integer)
            val = val[rempath[1]]
            rempath = rempath[2:end]
        else
            return [get_field(item, rempath...) for item in val]
        end
    end

    (length(rempath) == 0) && (return val)
    get_field(val, rempath...)
end
get_field(o::T, name::String) where {T<:SwaggerModel} = get_field(o, name_map(o)[name])
get_field(o::T, name::Symbol) where {T<:SwaggerModel} = getfield(o, name)

function isset_field(o::T, path...) where {T<:SwaggerModel}
    ret = isset_field(o, path[1])
    rempath = path[2:end]
    (length(rempath) == 0) && (return ret)

    val = get_field(o, path[1])
    if isa(val, Vector)
        if isa(rempath[1], Integer)
            ret = length(val) >= rempath[1]
            if ret
                val = val[rempath[1]]
                rempath = rempath[2:end]
            end
        else
            return [isset_field(item, rempath...) for item in val]
        end
    end

    (length(rempath) == 0) && (return ret)
    isset_field(val, rempath...)
end
isset_field(o::T, name::String) where {T<:SwaggerModel} = isset_field(o, name_map(o)[name])
isset_field(o::T, name::Symbol) where {T<:SwaggerModel} = (getfield(o, name) !== nothing)

set_field!(o::T, name::String, val) where {T<:SwaggerModel} = set_field!(o, name_map(o)[name], val)
function set_field!(o::T, name::Symbol, val) where {T<:SwaggerModel}
    validate_field(o, name, val)
    FT = fieldtype(T,name)

    if isa(val, FT)
        return setfield!(o, name, val)
    else
        ftval = try
            convert(FT, val)
        catch
            FT(val)
        end
        return setfield!(o, name, ftval)
    end
end

convert(::Type{T}, json::Dict{String,Any}) where {T<:SwaggerModel} = from_json(T, json)
convert(::Type{T}, v::Nothing) where {T<:SwaggerModel} = T()

show(io::IO, model::T) where {T<:SwaggerModel} = print(io, JSON.json(model, 2))
summary(model::T) where {T<:SwaggerModel} = print(io, T)


text:
"""
Acceps a map function that is applied to the `zero_rate` of the curve.

The map function `f` takes the form `f(rate, maturity)` with the following arguments:

    * `rate` is the retult of `zero_rate` applied to underlying curve for `maturity`.

    * `maturity` is the requested zero rate maturity.

The `zero_rate` for a `CurveMap` is implemented as:

```julia
zero_rate(curve::CurveMap, maturity::Date) = curve.f(zero_rate(curve.curve, maturity), maturity)
```

# Example

```julia
vert_x = [11, 15, 19, 23]
vert_y = [0.09, 0.14, 0.19, 0.18] # yield values 9%, 14%, 19%, 18%

# parallel shock of 1%
function map_parallel_1pct(rate, maturity)
    return rate + 0.01
end

dt_curve = Date(2015,08,03)

curve_map = InterestRates.CurveMap("mapped-curve", map_parallel_1pct, InterestRates.IRCurve("dummy-cont-flatforward", InterestRates.Actual360(),
    InterestRates.ContinuousCompounding(), InterestRates.FlatForward(), dt_curve,
    vert_x, vert_y))

# will report zero rate as 10% for maturity 11 days
zero_rate(curve_map, dt_curve + Dates.Day(11)) ≈ 0.1
```
"""
struct CurveMap{M, C<:AbstractIRCurve{M}, F<:Function} <: AbstractIRCurve{M}
    name::String
    f::F
    curve::C
end

for fun in (:curve_get_daycount, :curve_get_compounding, :curve_get_method, :curve_get_date, :curve_get_dtm, :curve_get_zero_rates, :curve_get_model_parameters, :curve_get_spline_fit_on_rates, :curve_get_spline_fit_on_discount_factors)
    @eval begin
        ($fun)(curve::CurveMap) = ($fun)(curve.curve)
    end
end

zero_rate(curve::CurveMap, maturity::Date) = curve.f(zero_rate(curve.curve, maturity), maturity)
curve_get_name(curve::CurveMap) = curve.name


text:
# Practical Machine learning
# Artificial Neural Network
# Chapter 11

function debugInitializeWeights(fan_out, fan_in)
  #   Note that W should be set to a matrix of size(1 + fan_in, fan_out) as
  #   the first row of W handles the "bias" terms

  # Set W to zeros
  W = DEBUGINITIALIZEWEIGHTS(fan_in, fan_out)

  # Initialize W using "sin", this ensures that W is always of the same
  # values and will be useful for debugging
  W = reshape(sin(1:length(W)), size(W)) / 10

  return W
end