text:
module HeatCapacity

# This implements different methods to specify heat capacity
#
# If you want to add a new method here, feel free to do so.
# Remember to also export the function name in GeoParams.jl (in addition to here)

using Parameters, LaTeXStrings, Unitful
using ..Units
using GeoParams: AbstractMaterialParam, AbstractMaterialParamsStruct
import Base.show, GeoParams.param_info
using ..MaterialParameters: MaterialParamsInfo

abstract type AbstractHeatCapacity{T} <: AbstractMaterialParam end

export  compute_heatcapacity,               # calculation routines
        compute_heatcapacity!,              # in-place routine
        ConstantHeatCapacity,               # constant
        T_HeatCapacity_Whittington,          # T-dependent heat capacity
        param_info

include("../Utils.jl")
include("../Computations.jl")
# Constant Heat Capacity -------------------------------------------------------
"""
    ConstantHeatCapacity(cp=1050J/mol/kg)

Set a constant heat capacity:
```math
    cp  = cst
```
where ``cp`` is the thermal heat capacity [``J/kg/K``].
"""
@with_kw_noshow struct ConstantHeatCapacity{T,U} <: AbstractHeatCapacity{T}
    cp::GeoUnit{T,U}        =   1050J/kg/K                # heat capacity
end
ConstantHeatCapacity(args...) = ConstantHeatCapacity(convert.(GeoUnit,args)...)

function param_info(s::ConstantHeatCapacity) # info about the struct
    return MaterialParamsInfo(Equation =   L"c_p = cst")
end

# Calculation routine
function (s::ConstantHeatCapacity)(; kwargs...)
    @unpack_val cp   = s
    return cp
end

compute_heatcapacity(s::ConstantHeatCapacity{_T}; kwargs...) where _T = s()

# Calculation routine
function compute_heatcapacity!(cp_array::AbstractArray{_T}, s::ConstantHeatCapacity{_T}; kwargs...) where _T
    Threads.@threads for i in eachindex(cp_array)
        @inbounds cp_array[i] = s()
    end
    return nothing
end

# compute_heatcapacity!(cp_array::AbstractArray{_T}, s::ConstantHeatCapacity{_T}, args) where _T = compute_heatcapacity!(cp_array, s; args...)

# Print info
function show(io::IO, g::ConstantHeatCapacity)
    print(io, "Constant heat capacity: cp=$(UnitValue(g.cp))")
end
#-------------------------------------------------------------------------

# Temperature dependent heat capacity -------------------------------
"""
    T_HeatCapacity_Whittington()

Sets a temperature-dependent heat capacity following the parameterization of Whittington et al. (2009), Nature:
```math
    Cp = (a + b T - c/T^2)/m
```

where ``Cp`` is the heat capacity [``J/kg/K``], and ``a,b,c`` are parameters that dependent on the temperature `T`:
- a = 199.50 J/mol/K    if T<= 846 K
- a = 199.50 J/mol/K    if T> 846 K
- b = 0.0857J/mol/K^2   if T<= 846 K
- b = 0.0323J/mol/K^2   if T> 846 K
- c = 5e6J/mol*K        if T<= 846 K
- c = 47.9e-6J/mol*K    if T> 846 K
- molmass =   0.22178kg/mol

Note that this is slightly different than the equation in the manuscript, as Cp is in J/kg/K (rather than ``J/mol/K`` as in eq.3/4 of the paper)
"""
@with_kw_noshow struct T_HeatCapacity_Whittington{T,U1,U2,U3,U4,U5} <: AbstractHeatCapacity{T}
    # Note: the resulting curve was visually compared with Fig. 2 of the paper
    a0::GeoUnit{T,U1}          =   199.5J/mol/K                # prefactor for low T       (T<= 846 K)
    a1::GeoUnit{T,U1}          =   229.32J/mol/K               # prefactor for high T      (T>  846 K)
    b0::GeoUnit{T,U2}          =   0.0857J/mol/K^2             # linear term for low T     (T<= 846 K)
    b1::GeoUnit{T,U2}          =   0.0323J/mol/K^2             # linear term for high T    (T>  846 K)
    c0::GeoUnit{T,U3}          =   5e6J/mol*K                  # quadratic term for low T  (T<= 846 K)
    c1::GeoUnit{T,U3}          =   47.9e-6J/mol*K              # quadratic term for high T (T>  846 K)
    molmass::GeoUnit{T,U4}     =   0.22178kg/mol               # average molar mass
    Tcutoff::GeoUnit{T,U5}     =   846K                        # cutoff temperature
end
T_HeatCapacity_Whittington(args...) = T_HeatCapacity_Whittington(convert.(GeoUnit,args)...)

function param_info(s::T_HeatCapacity_Whittington) # info about the struct
    return MaterialParamsInfo(Equation =   L"c_p = (a + b*T - c/T^2)/m")
end

# Calculation routine
function (s::T_HeatCapacity_Whittington{_T})(; T::_T=zero(_T), kwargs...) where _T
    @unpack_val a0,a1,b0,b1,c0,c1, molmass, Tcutoff   = s

    cp = a0/molmass

    if T <= Tcutoff
        a,b,c = a0,b0,c0
    else
        a,b,c = a1,b1,c1
    end

    cp = (a + b*T - c/T^2)/molmass

    return cp
end

"""
    compute_heatcapacity!(cp_array::AbstractArray{_T, N},s::T_HeatCapacity_Whittington{_T}, T::_T=zero(_T), P::_T=zero(_T)) where {_T,N}

Computes T-dependent heat capacity in-place
"""
function compute_heatcapacity!(cp_array::AbstractArray{_T, N},s::T_HeatCapacity_Whittington{_T}; T::AbstractArray{_T, N}, kwargs...) where {_T,N}
    @inbounds for i in eachindex(T)
        cp_array[i] = s(T=T[i])
    end
    return nothing
end

# add methods programatically
for myType in (:ConstantHeatCapacity, :T_HeatCapacity_Whittington)
    @eval begin
        (s::$(myType))(args)= s(; args...)
        compute_heatcapacity(s::$(myType), args) = s(args)
        compute_heatcapacity!(k::AbstractArray{_T,N}, s::$(myType){_T}, args) where {_T,N} = compute_heatcapacity!(k, s; args...)
    end
end

"""
    compute_heatcapacity!(cp_array::AbstractArray{_T,N},s::T_HeatCapacity_Whittington{_T}, T::AbstractArray{_T,N},P::AbstractArray{_T,N}) where {_T,N}

Computes T-dependent heat capacity in-place
"""
#=
function compute_heatcapacity!(cp_array::AbstractArray{_T,N},s::T_HeatCapacity_Whittington{_T}, T::AbstractArray{_T,N},P::AbstractArray{_T,N}) where {_T,N}

    compute_heatcapacity(s, cp_array,T)

    return nothing
end
=#

# Print info
function show(io::IO, g::T_HeatCapacity_Whittington)
#    print(io, "T-dependent heat capacity: cp/$(UnitValue(g.molmass))=$(UnitValue(g.a0)) + $(UnitValue(g.b0))*T - $(UnitValue(g.c0))/T^2 (for T<=$(UnitValue(g.Tcutoff))); ");
#    print(io, " cp/$(UnitValue(g.molmass))=$(UnitValue(g.a1)) + $(UnitValue(g.b1))*T - $(UnitValue(g.c1))/T^2 (for T>$(UnitValue(g.Tcutoff))) \n");
    print(io, "T-dependent heat capacity following Whittington et al. (2009) for average crust). \n");
end
#-------------------------------------------------------------------------

#-------------------------------------------------------------------------
# Heat capacity from phase diagram

# to be implemented - see density implementation

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

# Help info for the calculation routines
"""
    Cp = compute_heatcapacity(s:<AbstractHeatCapacity, P, T)

Returns the heat capacity `Cp` at any temperature `T` and pressure `P` using any of the heat capacity laws implemented.

Currently available:
- ConstantHeatCapacity
- T\\_HeatCapacity_Whittington

# Example
Using dimensional units
```julia
julia> T  = (250:100:1250)*K;
julia> cp = T_HeatCapacity_Whittington()
julia> Cp = ComputeHeatCapacity(0,T,cp)
11-element Vector{Unitful.Quantity{Float64, 𝐋² 𝚯⁻¹ 𝐓⁻², Unitful.FreeUnits{(kg⁻¹, J, K⁻¹), 𝐋² 𝚯⁻¹ 𝐓⁻², nothing}}}:
  635.4269997294616 J kg⁻¹ K⁻¹
  850.7470171764261 J kg⁻¹ K⁻¹
  962.0959598489883 J kg⁻¹ K⁻¹
 1037.5420433770641 J kg⁻¹ K⁻¹
 1097.3517921966488 J kg⁻¹ K⁻¹
 1149.2745563671706 J kg⁻¹ K⁻¹
 1157.7915050948404 J kg⁻¹ K⁻¹
 1172.3554874197264 J kg⁻¹ K⁻¹
 1186.9194697445964 J kg⁻¹ K⁻¹
  1201.483452069455 J kg⁻¹ K⁻¹
 1216.0474343943067 J kg⁻¹ K⁻¹
```

"""
compute_heatcapacity()

"""
    Cp = ComputeHeatCapacity(T::Any, s::AbstractHeatCapacity)

Computes heat capacity if only temperature (and not pressure) is specified
"""
# compute_heatcapacity(s::AbstractHeatCapacity, T::AbstractArray{_T}) where _T =  compute_heatcapacity(s,similar(T), T)
# compute_heatcapacity!(cp_array::AbstractArray{_T}, s::AbstractHeatCapacity, T::AbstractArray{_T}) where _T =  compute_heatcapacity!(cp_array,s,similar(T), T)

"""
    Cp = ComputeHeatCapacity(s::ConstantHeatCapacity)

Returns heat capacity if we are sure that we will only employ constant heat capacity in the simulation
"""

# Computational routines needed for computations with the MaterialParams structure
function compute_heatcapacity(s::AbstractMaterialParamsStruct, args)
    return s.HeatCapacity[1](args)
end

"""
    compute_heatcapacity!(Cp::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct}, Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat})

In-place computation of heat capacity `Cp` for the whole domain and all phases, in case a vector with phase properties `MatParam` is provided, along with `P` and `T` arrays.
This assumes that the `Phase` of every point is specified as an Integer in the `Phases` array.

______________________________________________________________________________________________

compute_heatcapacity!(Cp::AbstractArray{_T, N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct, 1}, PhaseRatios::AbstractArray{_T, M}, P::AbstractArray{_T, N},T::AbstractArray{_T, N})

In-place computation of heat capacity `Cp` for the whole domain and all phases, in case a vector with phase properties `MatParam` is provided, along with `P` and `T` arrays.
This assumes that the `PhaseRatio` of every point is specified as an Integer in the `PhaseRatios` array, which has one dimension more than the data arrays (and has a phase fraction between 0-1)

"""
compute_heatcapacity(args...) = compute_param(compute_heatcapacity, args...)
compute_heatcapacity!(args...) = compute_param!(compute_heatcapacity, args...)

# In case just temperature is provided
function compute_heatcapacity!(Cp::AbstractArray{_T, ndim}, MatParam::NTuple{N,AbstractMaterialParamsStruct}, Phases::AbstractArray{_I, ndim}, T::AbstractArray{_T, ndim}) where {_T,ndim,N,_I<:Integer}
    compute_param!(compute_heatcapacity,Cp,MatParam,Phases,nothing,T)
end

end


text:
using Interpolations, Colors

type ImageStitcher
    image
    tform
    xlim
    ylim
end

ImageStitcher(image, tform) = ImageStitcher(image, tform, Float64[], Float64[])

# in a module we need to export the names, to be imported via the `using` command
export ImageStitcher

"""
 Syntax:       Is, alpha = stitchImages(It);
               Is, alpha = stitchImages(...,'dim',dim,...);
               Is, alpha = stitchImages(...,'b0',b0,...);
               Is, alpha = stitchImages(...,'view',view,...);
               Is, alpha = stitchImages(...,'order',order,...);
               Is, alpha = stitchImages(...,'mode',mode,...);

 Inputs:       It is an n x 1 array of structs, where

                   It(i).image is an Nyi x Nxi x Nc uint8 image

                   It(i).tform is the 3 x 3 transformation matrix that
                   maps It(i).image into the coordinates specified by the
                   mode argument below

               [OPTIONAL] dim = [h, w] are the desired height and width,
               respectively, in pixels, of the stitched image. By default
               dim is set to preserve the input resolutions

               [OPTIONAL] b0 is a Nc x 1 vector of double values in
               [0, 255] specifying the background (i.e., missing data)
               color. The default value is b0 = 255 * ones(Nc,1)

               [OPTIONAL] view = {'left','center','right','default'}
               specifies what view to construct the stitched image from.
               The default is view = 'default' (world coordinates)

               [OPTIONAL] order = {'natural','reverse'} specifies the
               order in which to overlay the stitched images. In natural
               ordering, the closest image is displayed on top. The
               default is order = 'natural'

               [OPTIONAL] mode = {'absolute','relative'} specifies the
               type of transform used. When mode = 'absolute', It(i).tform
               should map to "world" coordinates. When mode = 'relative',
               It(i).tform should map It(i).image into the coordinates of
               It(i-1).image. Here, It(1).tform can be eye(3) or another
               transformation that maps to the desired world coordinates.
               The default value is mode = 'relative'

 Outputs:      Is is a h x w x Nc uint8 matrix containing the stitched
               image

               alpha is a h x w uint8 matrix containing the alpha channel
               (transparency) values for generating a nice transparent png
               of the stitched image

 Author:       Brian Moore
               brimoor@umich.edu

 Date:         September 30, 2015
"""
function stitchImages(It;
        view = "default",
        order = "natural",
        mode = "relative",
        dim = [],
        b0 = oneunit(eltype(It[1].image))
    )

    # Convert to absolute coordinates, if necessary
    if mode == "relative"
        toAbsoluteCoordinates(It)
    end

    # Apply view
    applyView(It, view)

    # Apply order
    applyOrder(It, order)

    # Compute stitched image limits
    xlim, ylim = computeStitchedLimits(It)

    # Get sample points
    x, y, w, h = getSamplePoints(xlim, ylim, dim)

    # Stitch images
    Is = fill(RGBA(NaN, NaN, NaN, NaN), h, w)
    for i = length(It):-1:1
        overlayImage(Is, It[i], x, y)
    end

    # Fill background
    fillBackground(Is, b0)

    Is
end

export stitchImages

#--------------------------------------------------------------------------
function toAbsoluteCoordinates(It)
    # Map all transforms to coordinates of It(1)
    for i = 2:length(It)
        It[i].tform = It[i].tform * It[i - 1].tform;
    end
end

#--------------------------------------------------------------------------
function applyView(It, view)
    # Get transformed limits
    Ns = length(It)
    xc = zeros(Ns)
    for i = 1:Ns
        xlimi, _ = getOutputLimits(It[i].image, It[i].tform)
        xc[i] = mean(xlimi)
    end

    # Get ordering
    view = lowercase(view)
    if "left" == view
        # Left view
        idx = sort(xc)
    elseif "center" == view
        # Center view
        idx = sort(abs(xc - mean(xc)));
    elseif "right" == view
        # Right view
        idx = sort(xc, rev = true)
    elseif "default" == view
        # Use input ordering
        idx = 1:Ns
    else
        # Unsupported view
        error("Unsupported view $view")
    end

    # Apply ordering
    It[1:end] = It[idx]
    if view != "default"
        H1 = It[1].tform
        for i = 1:Ns
            It[i].tform = It[i].tform / H1
        end
    end
end
#--------------------------------------------------------------------------
function applyOrder(It, order)
    # Apply order
    order = lowercase(order)
    if order == "natural"
        # Natural ordering
        # Empty
    elseif order == "reverse"
        # Reverse ordering
        reverse!(It)
    else
        # Unsupported order
        error("Unsupported order: $order")
    end
end
#--------------------------------------------------------------------------
function computeStitchedLimits(It)
    # Compute limits
    minx = Inf;
    maxx = -Inf;
    miny = Inf;
    maxy = -Inf;
    for i = 1:length(It)
        xlimi, ylimi = getOutputLimits(It[i].image, It[i].tform)
        It[i].xlim = xlimi
        It[i].ylim = ylimi
        minx = min(minx, xlimi[1])
        maxx = max(maxx, xlimi[2])
        miny = min(miny, ylimi[1])
        maxy = max(maxy, ylimi[2])
    end
    xlim = [floor(Int, minx), ceil(Int, maxx)]
    ylim = [floor(Int, miny), ceil(Int, maxy)]
    xlim, ylim
end
#--------------------------------------------------------------------------
function getOutputLimits(I, H)
    # Compute limits of transformed image
    Ny, Nx = size(I)
    X = [1 Nx Nx 1]
    Y = [1 1 Ny Ny]
    Xt, Yt = applyTransform(X, Y, H)
    xlim = [minimum(Xt), maximum(Xt)]
    ylim = [minimum(Yt), maximum(Yt)]
    xlim, ylim
end

#--------------------------------------------------------------------------
function applyTransform(X, Y, H)
    # Apply transformation
    sz = size(X)
    n = length(X)
    tmp = [X[:] Y[:] ones(n)] * H
    Xt = reshape(tmp[:, 1] ./ tmp[:, 3], sz)
    Yt = reshape(tmp[:, 2] ./ tmp[:, 3], sz)
    Xt, Yt
end

#--------------------------------------------------------------------------
function applyInverseTransform(Xt, Yt, H)
    # Apply inverse transformation
    sz = size(Xt)
    n = length(Xt)
    tmp = [Xt[:] Yt[:] ones(n)] / H
    X = reshape(tmp[:, 1] ./ tmp[:, 3], sz)
    Y = reshape(tmp[:, 2] ./ tmp[:, 3], sz)
    X, Y
end
#--------------------------------------------------------------------------
function getSamplePoints(xlim, ylim, dim)
    # Get sample dimensions
    if isempty(dim)
        w = (diff(xlim) + 1)[1]
        h = (diff(ylim) + 1)[1]
    else
        h, w = dim
    end

    # Limit resolution to a reasonable value, if necessary
    MAX_PIXELS = 2000 * 2000
    w, h = limitRes(w, h, MAX_PIXELS)

    # Compute sample points
    x = linspace(xlim[1], xlim[2], w)
    y = linspace(ylim[1], ylim[2], h)
    x, y, w, h
end

#--------------------------------------------------------------------------
function limitRes(w, h, lim)
    if w * h <= lim
        # No rescaling needed
        return w, h
    end

    # Rescale to meet limit
    kappa = w / h
    w = round(Int, sqrt(lim * kappa))
    h = round(Int, sqrt(lim / kappa))
    warn("Output resolution too large, rescaling to $h x $w")
    w, h
end

#--------------------------------------------------------------------------
function overlayImage(Is, It, x, y)
    # Overlay image
    If = fillImage(It, x, y)
    mask = (.~)(any(isnan.(If), 3))
    Isj = Is[:, :]
    Ifj = If[:, :]
    Isj[mask] = Ifj[mask]
    Is[:, :] = Isj
end

#--------------------------------------------------------------------------
function fillImage(It, x, y)
    # Parse inputs
    w = length(x)
    h = length(y)

    # Get active coordinates
    xIdx1 = findlast( x-> x <= It.xlim[1], x)
    xIdx2 = findfirst(x-> x >= It.xlim[2], x)

    yIdx1 = findlast( y-> y <= It.ylim[1], y)
    yIdx2 = findfirst(y-> y >= It.ylim[2], y)
    wa = xIdx2 + 1 - xIdx1
    ha = yIdx2 + 1 - yIdx1

    # Compute inverse transformed coordinates
    # part overlapping with It.image
    Xta = [x[i] for j in yIdx1:yIdx2, i in xIdx1:xIdx2]
    Yta = [y[j] for j in yIdx1:yIdx2, i in xIdx1:xIdx2]
    # coordinates of target image transformed to It.image
    Xa, Ya = applyInverseTransform(Xta, Yta, It.tform)

    # Compute active image
    Ia = fill(RGBA(NaN, NaN, NaN, NaN), ha, wa)
    imgip = interpolate(It.image, BSpline(Linear()), OnGrid())
    for xi = 1:wa, yi = 1:ha
        if xi <= size(Xa, 2) && yi <= size(Xa, 1)
            xt, yt = Xa[yi, xi], Ya[yi, xi]
            if xt <= size(imgip, 2) && yt <= size(imgip, 1) && yt >= 1 && xt >= 1
                Ia[yi, xi] = imgip[yt, xt]
            end
        end
    end

    # Embed into full image
    If = fill(RGBA(NaN, NaN, NaN, NaN), h, w)
    If[yIdx1:yIdx2, xIdx1:xIdx2] = Ia
    If
end
#--------------------------------------------------------------------------
function fillBackground(Is, b0)
    # Fill background
    mask = any(isnan.(Is), 3)
    Isj = Is[:, :]
    Isj[mask] = b0
    Is[:, :] = Isj

    # alpha = zeros(N0f8, size(mask))
    # alpha[~mask] = 1.0
end

ext: jl


text:
# GV post filters
# Reference:
# Silén, Hanna, et al. "Ways to Implement Global Variance in Statistical Speech
# Synthesis." INTERSPEECH. 2012.

immutable VarianceScaling
    σ²::Vector{Float64}
end

function fvpostf!(vs::VarianceScaling, src::AbstractMatrix)
    D = size(src, 1)
    μ = mean(src, 2)
    src[:,:] = sqrt(vs.σ² ./ var(src[1:D,:], 2)) .* (src .- μ) .+ μ
    src
end

function fvpostf(vs::VarianceScaling, src::AbstractMatrix)
    filtered = copy(src)
    fvpostf!(vs, filtered)
    filtered
end


text:
using Plots
using Random
Random.seed!(1)

# soft rect
soft_rect(z) = log(1.0 + exp(z))
zz = range(-5.0, stop = 5.0, length = 100)
plot(zz, soft_rect.(zz))
plot(zz, 5.0 * tanh.(zz))

# ddp
include_ddp()

# Model
include_model("double_integrator")

model = DoubleIntegratorContinuous{RK3, FixedTime}(2, 1, 1)
n = model.n
m = model.m

struct MultipleModel{I, T} <: Model{I, T}
	n::Vector{Int}
	m::Vector{Int}
	d::Vector{Int}

	model::Model

	N::Int
	p::Int
end

function multiple_model(model, T, N; p = 0)
	n = model.n
	m = model.m
	d = model.d

	n = [N * n + (t != 1 ? p + N * m : 0) for t = 1:T]
	m = [N * m + (t == 1 ? p : 0) for t = 1:T-1]
	d = [N * d for t = 1:T-1]

	MultipleModel{typeof(model).parameters...}(n, m, d, model, N, p)
end

# Policy
dp = 1
p_policy = dp * model.n * model.n + dp * model.n + dp * model.n * dp * model.n + dp * model.n + model.m * dp * model.n + model.m
# p_policy = dp * model.n * model.n + dp * model.n + model.m * dp * model.n + model.m
# p_policy = model.m * model.n + model.m

function policy(θ, x, t, n, m)
	# 2 layer
	p_policy = dp * model.n * model.n + dp * model.n + dp * model.n * dp * model.n + dp * model.n + model.m * dp * model.n + model.m

	K1 = reshape(view(θ, 1:(dp * n) * n), dp * n, n)
	k1 = view(θ, (dp * n) * n .+ (1:(dp * n)))

	K2 = reshape(view(θ, dp * n * n + dp * n .+ (1:(dp * n * dp * n))), dp * n, dp * n)
	k2 = view(θ, dp * n * n + dp * n + dp * n * dp * n .+ (1:(dp * n)))

	Ko = reshape(view(θ, dp * n * n + dp * n + dp * n * dp * n + dp * n .+ (1:(m * dp * n))), m, dp * n)
	ko = view(θ, dp * n * n + dp * n + dp * n * dp * n + dp * n + m * dp * n .+ (1:m))

	# z1 = soft_rect.(K1 * x + k1)
	# z2 = soft_rect.(K2 * z1 + k2)
	z1 = tanh.(K1 * x + k1)
	z2 = tanh.(K2 * z1 + k2)
	zo = Ko * z2 + ko

	## 1 layer
	# p_policy = dp * model.n * model.n + dp * model.n + model.m * dp * model.n + model.m
	#
	# K1 = reshape(view(θ, 1:(dp * n) * n), dp * n, n)
	# k1 = view(θ, (dp * n) * n .+ (1:(dp * n)))
	#
	# z1 = soft_rect.(K1 * x + k1)
	# Ko = reshape(view(θ, dp * n * n + dp * n .+ (1:m * (dp * n))), m, dp * n)
	# ko = view(θ, dp * n * n + dp * n + m * dp * n .+ (1:m))
	#
	# zo = Ko * z1 + ko

	## affine
	# p_policy = model.m * model.n + model.m
	#
	# Ko = reshape(view(θ, 1:m * n), m, n)
	# ko = view(θ, m * n .+ (1:m))
	#
	# zo = Ko * x + ko

	return zo
end

function f(model::DoubleIntegratorContinuous, x, u, w)
    [x[2];
	 (1.0 + w[1]) * (u[1]
	 	+ policy(view(u, 1 .+ (1:p_policy)), x, nothing, model.n, model.m)[1])]
end

function fd(models::MultipleModel, x, u, w, h, t)
	N = models.N

	n = models.n[t]
	m = models.m[t]
	d = models.d[t]

	ni = models.model.n
	mi = models.model.m
	di = models.model.d

	p = models.p

	x⁺ = []
	u_policy = []

	if t == 1
		θ = view(u, N * mi .+ (1:p))
	else
		θ = view(x, N * ni + N * mi .+ (1:p))
	end

	for i = 1:N
		xi = view(x, (i - 1) * ni .+ (1:ni))
		ui = [u[(i - 1) * mi .+ (1:mi)]; θ]
		wi = view(w, (i - 1) * di .+ (1:di))
		push!(x⁺, fd(models.model, xi, ui, wi, h, t))
		push!(u_policy, policy(θ, xi, t, ni, mi))
	end

	return vcat(x⁺..., u_policy..., θ)
end

# Time
T = 11
h = 0.1
tf = h * (T - 1)
N = 5
models = multiple_model(model, T, N, p = p_policy)

x_ref = [[0.0; 0.0] for t = 1:T]
u_ref = [zeros(model.m) for t = 1:T-1]
_xT = [0.0; 0.0]
xT = [vcat([_xT for i = 1:N]..., zeros(t == 1 ? 0 : N * model.m), zeros(t == 1 ? 0 : models.p)) for t = 1:T]

# Initial conditions, controls, disturbances
x1 = zeros(models.n[1])
x1_add = [1.0, -1.0, 2.0, -2.0, 0.0]#3.0, -3.0, 4.0, -4.0, 0.0]
for i = 1:N
	x1[(i - 1) * (model.n) + 1] = 0.0 + x1_add[i]
end

ū = [1.0 * randn(models.m[t]) for t = 1:T-1]
wi = [0.0, 0.0, 0.0, 0.0, 0.0]#, 0.0, 0.0, 0.0, 0.0]

@assert length(wi) == N
w = [vcat(wi...) for t = 1:T-1]

# Rollout
x̄ = rollout(models, x1, ū, w, h, T)

# Objective
_R = 1.0e-1 * ones(N * model.m)
_Ru = 10000.0 * ones(N * model.m)

Q = [(t < T ?
	 Diagonal(vcat([h * [1.0e-1; 1.0e-1] for i = 1:N]..., (t == 1 ? zeros(0) : _R)..., 1.0 * ones(t == 1 ? 0 : models.p)))
	: Diagonal(vcat([[100.0; 100.0] for i = 1:N]..., (t == 1 ? zeros(0) : _R)..., 1.0 * ones(t == 1 ? 0 : models.p)))) for t = 1:T]
q = [-2.0 * Q[t] * xT[t] for t = 1:T]

R = [Diagonal(t == 1 ? [h * _Ru; 1.0 * ones(models.p)] : _Ru) for t = 1:T-1]
r = [zeros(models.m[t]) for t = 1:T-1]

obj = StageCosts([QuadraticCost(Q[t], q[t],
	t < T ? R[t] : nothing, t < T ? r[t] : nothing) for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		Q = obj.cost[t].Q
		q = obj.cost[t].q
	    R = obj.cost[t].R
		r = obj.cost[t].r
        return (x' * Q * x + q' * x + u' * R * u + r' * u) ./ N
    elseif t == T
		Q = obj.cost[T].Q
		q = obj.cost[T].q
        return (x' * Q * x + q' * x) ./ N
    else
        return 0.0
    end
end

# Constraints
ns = models.N * models.model.n
ms = models.N * models.model.m
p_con = [t == T ? 0 * ns : 0 * ms for t = 1:T]
ul = [-Inf]
uu = [Inf]
info_t = Dict()
info_T = Dict()#:xT => _xT)

con_set = [StageConstraint(p_con[t], t < T ? info_t : info_T) for t = 1:T]

function c!(c, cons::StageConstraints, x, u, t)
	T = cons.T
	N = models.N
	n = models.model.n
	m = models.model.m
	p = models.p
	np = n + p
	ns = N * n
	ms = N * m

	# if t < T
	# 	c[1:ms] = view(u, 1:ms)
	# end

	# if t == T
	# 	for i = 1:N
	# 		c[(i - 1) * n .+ (1:n)] .= view(x, (i - 1) * n .+ (1:n)) - cons.con[T].info[:xT]
	# 	end
	# end
end

prob = problem_data(models, obj, con_set, copy(x̄), copy(ū), w, h, T,
	n = models.n, m = models.m, d = models.d)

# Solve
@time constrained_ddp_solve!(prob,
    max_iter = 1000, max_al_iter = 10,
	ρ_init = 1.0, ρ_scale = 10.0,
	con_tol = 1.0e-5)

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)
x̄i = [x̄[t][1:model.n] for t = 1:T]

# individual trajectories
x_idx = [(i - 1) * model.n .+ (1:model.n) for i = 1:N]
u_idx = [N * model.n + (i - 1) * model.m .+ (1:model.m) for i = 1:N]

# Visualize
x_idxs = vcat(x_idx...)
u_idxs = vcat(u_idx...)

# state
plot(hcat([xT[t][x_idxs] for t = 1:T]...)',
    width = 2.0, color = :black, label = "")
plot!(hcat([x[t][x_idxs] for t = 1:T]...)', color = :magenta, label = "")

# control
plot(hcat([x[t][u_idxs] for t = 2:T]..., x[T][u_idxs], )',
    width = 2.0, color = :black, label = "", linetype = :steppost)

# verify solution
uθ = u[1][models.N * model.m .+ (1:models.p)]
xθ = [x[t][models.N * model.n + models.N * model.m .+ (1:models.p)] for t = 2:T]

policy_err = []
for t = 2:T
	push!(policy_err, norm(xθ[t-1] - uθ, Inf))
end
@show maximum(policy_err)

slack_err = []
for t = 1:T-1
	if t > 1
		push!(slack_err, norm(ū[t], Inf))
	else
		push!(slack_err, norm(ū[t][1:models.N * model.m], Inf))
	end
end
@show maximum(slack_err)

# Simulate policy
include(joinpath(@__DIR__, "simulate.jl"))

# Policy
θ = u[1][models.N * model.m .+ (1:models.p)]

# Model
model_sim = DoubleIntegratorContinuous{RK3, FixedTime}(model.n, model.m, model.d)
x1_sim = copy(x1[1:model.n])
T_sim = 10 * T

# Time
tf = h * (T - 1)
t = range(0, stop = tf, length = T)
t_sim = range(0, stop = tf, length = T_sim)
dt_sim = tf / (T_sim - 1)

# Simulate
N_sim = 10
x_sim = []
u_sim = []
J_sim = []
Random.seed!(1)
for k = 1:N_sim
	wi_sim = rand(range(-3.0, stop = 3.0, length = 1000))
	println("w: $wi_sim")
	# w_sim = [wi_sim for t = 1:T-1]
	w_sim = [wi_sim for t = 1:T-1]

	x_nn, u_nn, J_nn, Jx_nn, Ju_nn = simulate_policy(
		model_sim,
		θ,
		x_ref, u_ref,
		[_Q[1:model.n, 1:model.n] for _Q in Q], [_R[1:model.m, 1:model.m] for _R in R],
		T_sim, h,
		[wi_sim; 0.0],
		[zeros(model.d) for t = 1:T-1],
		ul = ul,
		uu = uu)

	push!(x_sim, x_nn)
	push!(u_sim, u_nn)
	push!(J_sim, J_nn)
end

# Visualize
idx = (1:2)
plt = plot(t, hcat(x_ref...)[idx, :]',
	width = 2.0, color = :black, label = "",
	xlabel = "time (s)", ylabel = "state",
	title = "double integrator (J_avg = $(round(mean(J_sim), digits = 3)), N_sim = $N_sim)")

for xs in x_sim
	plt = plot!(t_sim, hcat(xs...)[idx, :]',
	    width = 1.0, color = :magenta, label = "")
end
display(plt)

plt = plot(
	label = "",
	xlabel = "time (s)", ylabel = "control",
	title = "double integrator (J_avg = $(round(mean(J_sim), digits = 3)), N_sim = $N_sim)")
for us in u_sim
	plt = plot!(t_sim, hcat(us..., us[end])',
		width = 1.0, color = :magenta, label = "",
		linetype = :steppost)
end
display(plt)


text:
### This particular simulation allows particles to live indefenitively
# without erasing when arriving to the boundary.

### Domain [mm]
x_max = 1

### Models ###
# The considered PSF for the considered imaging tecnique
	# sigma also defines the resolution of the imaging system.
	sigma = 0.04
	filter = (x, y) -> exp(-(x^2 + y^2)/2/sigma^2)
	filter_dx = (x, y) -> -x/sigma^2*filter(x, y)
	filter_dy = (x, y) -> -y/sigma^2*filter(x, y)
# The number of time samples for reconstruction: 2K+1
	K = 2
# the maximum allowed speed for any considered particle
	v_max = 15 #[mm/s]
# The sampling rate
	tau = 1/500

### Noise ###
sigma_noise = 0.01

### Medium ###
# total number of the experimets measurements (gives the length of the experience)
	n_im = 1000
# Vessels separation
	dx = 0.02

### Inverse problem ###
# Considered norm to minimize
	lp_norm = 1
# Threshold on when to consider that there was a difference on the number of particles
# at two consecutive time steps.
	jump_threshold = 0.1

### Sequence ###
# Probability of a particle to activate or deactivate
	activation_probability=0.02


text:
######### Stan program example  ###########

using StanSample
using StatsPlots

bernoullimodel = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

observed_data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

# Default for tmpdir is to create a new tmpdir location
# To prevent recompilation of a Stan progam, choose a fixed location,
tmpdir=joinpath(@__DIR__, "tmp")

sm = SampleModel("bernoulli", bernoullimodel,
  method=StanSample.Sample(save_warmup=true, num_warmup=1000,
  num_samples=1000, thin=1, adapt=StanSample.Adapt(delta=0.85)),
  tmpdir=tmpdir);

(sample_file, log_file) = stan_sample(sm, data=observed_data);

if !isnothing(sample_file)
  chns = read_samples(sm)

  # Describe the results
  println()
  show(chns)
  println()

  # Optionally, convert to a DataFrame
  DataFrame(chns, showall=true, sorted=true, append_chains=true) |> display
  println()

  # Look at effective sample saize
  ess(chns) |> display
  println()

  # Check if StatsPlots is available and show basic MCMCChains plots
  if isdefined(Main, :StatsPlots)
    cd(@__DIR__) do
      p1 = plot(chns)
      savefig(p1, "traceplot.pdf")
      p2 = pooleddensity(chns)
      savefig(p2, "pooleddensity.pdf")
    end
  end

  # Ceate a ChainDataFrame
  summary_df = read_summary(sm)
  summary_df[:theta, [:mean, :ess]]
end


text:
# build GAP binary for Julia 1.7 (and possibly beyond)
include("../common.jl")

platforms = configure(v"1.7", v"1.7.0")

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("GMP_jll", compat="6.2.0"),
    Dependency("Readline_jll", compat="8.1.1"),
    Dependency("Zlib_jll"),
    BuildDependency(PackageSpec(name="libjulia_jll", version=v"1.7.0")),
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"7", init_block=init_block, julia_compat="1.7")

max_stars_repo_name: singularitti/Yggdrasil

max_issues_repo_name: singularitti/Yggdrasil

max_forks_repo_name: singularitti/Yggdrasil


text:
export calculate_psd

function calculate_psd(A::AbstractArray, sampling_rate::Float64)
    N = size(A, 1)
    psd = Array{Float64}(undef, N)
    delta = 1. / sampling_rate

    for i in 1:N
        psd[i] = 2 * delta / N * abs2(A[i])
    end

    return psd
end


text:
using LAR

include("test/merge.jl")
include("test/spacePartition.jl")

using PyCall
@pyimport triangle

function chainCoords(EW,rows)
	edges = hcat([EW[:,e] for e in rows]...)
	for k=1:length(edges)
		v1,v2 = find(edges .== edges[k])
		if (v1 + v2) % 2 == 0
			# swap column elements, for column that contains v2
			col = v2 ÷ 2
			if v2 % 2 != 0
				col +=1
			end
			edges[2,col],edges[1,col] = edges[1,col],edges[2,col]
		end
	end
	out = Int[]
	for k=1:size(edges,2)
		if edges[1,k] < edges[2,k]
			push!(out,1)
		else
			push!(out,-1)
		end
	end
	out
end

function boundary_2(EW,FW)
	I,J,V = Int[],Int[],Int[]
	FE = crossRelation(FW,EW)
	for (f,rows) in enumerate(FE)
		append!(I,rows)
		append!(J,[f for k=1:length(rows)])
		append!(V,chainCoords(EW,rows))
	end
	return sparse(I,J,V)
end

function edgeSlopeOrdering(VE,V,EV)
 	VE_sorted = Array{Array{Int64,1},1}()
 	forwardBackward = Tuple{Array{Int64,1},Dict{Int64,Int64},Dict{Int64,Int64}}[]
 	for (v,ve) in enumerate(VE)
 		ve_angle = []
 		if ve != []
 			for edge in ve
 				w = pop!(setdiff(Set(EV[:,edge]), Set(v)))
 				(x,y) = V[:,w] - V[:,v]
 				angle = atan2(y,x)
 				push!(ve_angle, 180*angle/pi)
 			end
 		end
 		pairs = sort(collect(zip(ve_angle,ve)))
 		sortedEdges = [pair[2] for pair in pairs]
 		push!(VE_sorted, sortedEdges)
 		theNexts,thePrevs = nextPrev(sortedEdges)
 		next = Dict(zip(sortedEdges,theNexts))
 		prev = Dict(zip(sortedEdges,thePrevs))
 		push!(forwardBackward, (sortedEdges,next,prev))
 	end
 	forwardBackward
end

function larTriangulate(W::Array{Float64,2},FW::Array{Int64,2},EW::Array{Int,2})
	FV = [face for face in FW]
	return larTriangulate(W,FV,EW)
end

function larTriangulate(W::Array{Float64,2},FW::Array{Array{Int64,1},1},EW::Array{Int,2})
	FE = crossRelation(FW,EW)
	TW = Array{Array{Int64,2},1}()
	for (f,face) in enumerate(FW)
		@show (f,face)," in larTriangulate"
		# submanifold mapping and 2D projection
		vs = hcat([W[:,v] for v in face]...)
		fv = [[v for v=1:length(face)]]
		ev = Array{Array{Int64,1},1}([])
		for e in FE[f]
			edge = Int[]
			push!( edge, find(face .== EW[1,e])[1] )
			push!( edge, find(face .== EW[2,e])[1] )
			push!(ev,edge)
		end
		ev = hcat(ev ...)
		M = submanifoldMapping(vs,fv,1)
		U = vcat(vs,ones((1,size(vs,2))))
		Y = M * U
		Z = Y[ 1:2, : ]
		# prepare data for triangle.triangulation
		polygon = Dict(
			#"holes" => Array{Float64,2}(),
			"segments" => Array{Int32,2}(ev-1)',
			#"vertex_attributes" => hcat([[1.0] for k=1:size(Z,2)]...)',
			"vertices" => Array{Float64,2}(Z'),
		)
		tria = triangle.triangulate(PyObject(polygon), PyObject("p"))
		triangles = [face[v+1] for v in tria["triangles"]]
		push!(TW,triangles)
	end
	return TW
end

function boundary_3(W,EW,FW)
	BND_2 = boundary_2(EW,FW)
	m = size(BND_2,1)
	n = size(BND_2,2)
	EF = crossRelation(EW,FW)

	#forwardBackward = edgeSlopeOrdering(VE,V,EV)

end

v,(vv,ev,fv,cv) = p.larCuboids((3,3,1),true)
V = hcat([Array{Float64,1}(v[k,:]) for k=1:size(v,1)]...)
FV = hcat([Array{Int64,1}(fv[k,:]+1) for k=1:size(fv,1)]...)
EV = hcat([Array{Int64,1}(ev[k,:]+1) for k=1:size(ev,1)]...)
model1 = Any[V,FV,EV]

W = hcat([V[:,k] + [.5;0.5;0.5] for k=1:size(V,2)]...)
#W = rotate([0,0,π/6],W)
FW = copy(FV)
EW = copy(EV)
model2 = Any[W,FW,EW]

models = [model1,model2]
X,FX,EX = cmerge(models)

viewexploded(X,FX)
viewexploded(X,EX)

V,FV,EV = deepcopy((X,FX,EX))
W,FW,EW = spacePartition(V,FV,EV)

viewexploded(W,EW)
viewexploded(W,FW)

#full(boundary_2(EW,FW))
#viewLarIndices(W,EW,FW,2)

TW = larTriangulate(W,FW,EW)
viewexploded(W,vcat(TW...)')


text:
# Patterns consists of 2 or more channels. This represents just one of the channels.
mutable struct PatternChannel <: AbstractBitStream
    base::BaseData
    pattern::AbstractBitStream # Could be a RegularPatternStream

    function PatternChannel(pattern::AbstractBitStream, id::Int64)
        o = new()
        o.base = BaseData{UInt8}()
        o.base.id = id
        o.pattern = pattern
        o
    end
end

function output(stream::PatternChannel)
    output(stream.pattern, stream.base.id)
end


text:
@testset "step(r)" begin
    for (r,b) in ((OneToMRange(10), OneTo(10)),
                  (OneToSRange(UInt(10)), OneTo(UInt(10))),
                  (UnitMRange(1, 10), UnitRange(1, 10)),
                  (UnitSRange(1.,10.), UnitRange(1.0, 10.0)),
                  (StepMRange(1, 2, 10), StepRange(1, 2, 10)),
                  (StepSRange(UInt32(1), UInt32(2), UInt32(10)), StepRange(UInt32(1), UInt32(2), UInt32(10))),
                  (mrange(1.0, step=2.0, stop=10.0), 1.0:2.0:10.0),
                  (srange(UInt32(1), step=UInt32(2), stop=UInt32(10)), UInt32(1):UInt32(2):UInt32(10)),
                 )
        @test @inferred(step(r)) === step(b)
        @test @inferred(step_hp(r)) === step_hp(b)
        if b isa StepRangeLen
            @test @inferred(stephi(r)) == stephi(b)
            @test @inferred(steplo(r)) == steplo(b)
        end
    end
end


text:
using NeuralQuantum, QuantumOpticsBase, ProgressMeter

N = 7
g = 0.4
V = 2.0

hilb = HomogeneousSpin(N,1//2)
hilb = HomogeneousFock(N,2)

ops = []
H = LocalOperator(hilb)

Sx = LocalOperator(hilb)
Sy = LocalOperator(hilb)
Sz = LocalOperator(hilb)

for i=1:N
    global H += g/2.0 * sigmax(hilb, i)
    global H += V/4.0 * sigmaz(hilb, i) * sigmaz(hilb, mod(i, N)+1)

    global Sx += sigmax(hilb, i)/N
    global Sy += sigmay(hilb, i)/N
    global Sz += sigmaz(hilb, i)/N

    push!(ops, sigmam(hilb, i))
end

liouv = liouvillian(H, ops)

sampl = MetropolisSampler(LocalRule(), 125, N, burn=100)
#sampl = ExactSampler(5000)
algo  = SR(ϵ=(0.001), algorithm=sr_cholesky)
#algo  = Gradient()

net  = NDM(Float64, N, 1, 1, af_softplus)
is = BatchedSampler(net, sampl, liouv, algo; batch_sz=16)
add_observable!(is, "Sx", Sx)
add_observable!(is, "Sy", Sy)
add_observable!(is, "Sz", Sz)

optimizer = Optimisers.Descent(0.01)

Evalues = Float64[];
Eerr = Float64[];
for i=1:200
    ldata, prec = sample!(is)
    ob = compute_observables(is)

    println("$i - $ldata")

    push!(Evalues, ldata.mean)
    push!(Eerr, ldata.error)
    grad = precondition!(prec, algo, i)
    Optimisers.update!(optimizer, net, grad)
end

plot(Evalues, yerr=Eerr, yscale=:log10)


text:
using DistributedFactorGraphs
struct TestInferenceVariable1 <: InferenceVariable end
a = DFGVariable(:a, TestInferenceVariable1())
b = copy(a)

d1 = getSolverData(a)
pop!(solverDataDict(a), :default)
d2 = getSolverData(b)

d1.dims = 6
d2.dims


text:
############################################################
# joCoreBlockConstructors ##################################
############################################################

# joBlockDiag
include("joCoreBlockConstructors/joBlockDiag.jl")

# joDict
include("joCoreBlockConstructors/joDict.jl")

# joStack
include("joCoreBlockConstructors/joStack.jl")

# joBlock
include("joCoreBlockConstructors/joBlock.jl")


text:
# ============================================================
# I / O Methods
# ============================================================

export read_cfg,
       read_dict,
       write_jld,
       write_vtk,
       plot_line,
       plot_contour

"""
    read_cfg(filename::T) where {T<:AbstractString}

Read configuration into dictionary

* @arg filename: configuration text file
* @return vars: dictionary with values of variables

"""
function read_cfg(filename::T) where {T<:AbstractString}
    D = read_dict(filename)

    if haskey(D, :boundary)
        D[:boundary] = begin
            if parse(Int, D[:space][1]) == 1
                [D[:boundary], D[:boundary]]
            elseif parse(Int, D[:space][1]) == 2
                [D[:boundary], D[:boundary], D[:boundary], D[:boundary]]
            end
        end
    elseif haskey(D, :boundary4)
        D[:boundary] = [D[:boundary1], D[:boundary2], D[:boundary3], D[:boundary4]]
    elseif haskey(D, :boundary2)
        D[:boundary] = begin
            if parse(Int, D[:space][1]) == 1
                [D[:boundary1], D[:boundary2]]
            elseif parse(Int, D[:space][1]) == 2
                [D[:boundary1], D[:boundary1], D[:boundary2], D[:boundary2]]
            end
        end
    end

    return D
end

"""
    read_dict(filename::T, allowed) where {T<:AbstractString}
    read_dict(filename::T) where {T<:AbstractString}

Read text into dictionary

* @args filename: configuration text file
* @args allowed: keywords in range
* @return vars: dictionary with values of variables

"""
function read_dict(filename::T, allowed) where {T<:AbstractString}

    @info "reading config from $filename"
    println("--------------------------------------------------------------")
    f = open(filename)
    vars = Dict{String,Any}()

    for line in eachline(f)
        # skip comments
        if length(line) == 0 || line[1] == '#'
            continue
        end

        #print("\t");
        #println(line)
        var, val = split(line, "=")
        stripped = strip(var)
        if stripped in allowed
            println(line)

            #vars[stripped] = parse(Float64, val)
            #vars[stripped] = strip(val)
            tmp = tryparse(Float64, val)
            if isa(tmp, Nothing)
                vars[stripped] = strip(val)
            else
                vars[stripped] = isinteger(tmp) ? Int(tmp) : tmp
            end
        end
    end
    println("--------------------------------------------------------------")
    println("")
    return vars

end

function read_dict(filename::T) where {T<:AbstractString}

    @info "reading config from $filename"
    println("--------------------------------------------------------------")
    f = open(filename)
    vars = Dict{Symbol,Any}()

    for line in eachline(f)
        if length(line) == 0 || line[1] == '#'
            continue
        end
        println(line)

        var, val = split(line, "=")
        stripped = strip(var)
        stripped = Symbol(stripped)
        val = split(val, "#")[1]

        tmp = tryparse(Float64, val)
        if isa(tmp, Nothing)
            vars[stripped] = strip(val)
        else
            vars[stripped] = isinteger(tmp) ? Int(tmp) : tmp
        end
    end
    println("--------------------------------------------------------------")
    println("")
    return vars

end

"""
    write_jld(KS, ctr, t)

Write data into JLD2
"""
function write_jld(KS::X, ctr::Y, t = 0) where {X<:AbstractSolverSet,Y}

    strIter = string(t)
    fileOut = KS.outputFolder * "data/t=" * strIter * ".jld2"

    save(fileOut, Dict("set" => KS, "ctr" => ctr, "t" => t))

end

"""
    write_vtk(points, cells, cdata, pdata = zeros(axes(points, 1)))

Write data into VTK
"""
function write_vtk(
    points::T,
    cells,
    cdata,
    pdata = zeros(axes(points, 1)),
) where {T<:AbstractMatrix}
    mcells = [MeshCell(VTKCellTypes.VTK_TRIANGLE, cells[i, :]) for i in axes(cells, 1)]
    vtkfile = vtk_grid("sol", permutedims(points), mcells)

    len = size(cdata, 2)
    cname = Array{String}(undef, len)
    for i in eachindex(cname)
        cname[i] = "c" * string(i)
        vtkfile[cname[i], VTKCellData()] = cdata[:, i]
    end

    len = size(pdata, 2)
    pname = Array{String}(undef, len)
    for i in eachindex(pname)
        pname[i] = "p" * string(i)
        vtkfile[pname[i], VTKPointData()] = pdata[:, i]
    end

    outfiles = vtk_save(vtkfile)

    return nothing
end

"""
    plot_line(KS, ctr; backend = :plots)

Plot solution profiles
"""
function plot_line(
    KS::X,
    ctr::Y,
) where {X<:AbstractSolverSet,Y<:AbstractArray{<:AbstractControlVolume,1}}
    plot(KS, ctr)
end

@recipe function plot_line(
    KS::X,
    ctr::Y,
) where {X<:AbstractSolverSet,Y<:AbstractArray{<:AbstractControlVolume,1}}

    # solution
    pltx = KS.pSpace.x[1:KS.pSpace.nx]
    plty = zeros(KS.pSpace.nx, 3)
    for i in eachindex(pltx)
        for j = 1:2
            plty[i, j] = ctr[i].prim[j]
        end

        plty[i, 3] = 1.0 / ctr[i].prim[end]
    end

    # attributes
    xguide --> "x"
    :linewidth --> 1.5

    @series begin
        label := "ρ"
        pltx, plty[:, 1]
    end

    if ctr[1].w isa AbstractArray
        @series begin
            label := "u"
            pltx, plty[:, 2]
        end
        @series begin
            label := "T"
            pltx, plty[:, 3]
        end
    end

    # user-defined
    c = get(plotattributes, :linewidth, :auto)

    nothing

end

"""
    plot_contour(KS, ctr; backend = :plots)

Plot solution contours
"""
function plot_contour(
    KS::X,
    ctr::Y;
    color = :inferno,
    backend = :plots::Symbol,
) where {X<:AbstractSolverSet,Y<:AbstractArray{<:AbstractControlVolume,2}}

    sol = zeros(size(ctr[1].w, 1), KS.pSpace.nx, KS.pSpace.ny)
    for i in axes(sol, 2)
        for j in axes(sol, 3)
            for k = 1:size(sol, 1)-1
                sol[k, i, j] = ctr[i, j].prim[k]
            end
            sol[end, i, j] = 1.0 / ctr[i, j].prim[end]
        end
    end

    p1 = contourf(
        KS.pSpace.x[1:KS.pSpace.nx, 1],
        KS.pSpace.y[1, 1:KS.pSpace.ny],
        sol[1, :, :]',
        color = color,
    )
    p2 = contourf(
        KS.pSpace.x[1:KS.pSpace.nx, 1],
        KS.pSpace.y[1, 1:KS.pSpace.ny],
        sol[2, :, :]',
        color = color,
    )
    p3 = contourf(
        KS.pSpace.x[1:KS.pSpace.nx, 1],
        KS.pSpace.y[1, 1:KS.pSpace.ny],
        sol[3, :, :]',
        color = color,
    )
    p4 = contourf(
        KS.pSpace.x[1:KS.pSpace.nx, 1],
        KS.pSpace.y[1, 1:KS.pSpace.ny],
        sol[4, :, :]',
        color = color,
    )

    plot(p1, p2, p3, p4, layout = (2, 2))

end


text:
function test_poly_cvx(s)
    global solver = s

    test_dir = dirname(@__FILE__)

    @testset "Convex Polynomial Models" begin end
end


text:
# using Revise
using Test
using JDF, WeakRefStrings, DataFrames

@testset "JDF.jl WeakRefStrings.StringArrays" begin
    a = StringVector(["a", "b", "a", missing, "c"])
    io = open("io.jdf", "w")
    metadata = compress_then_write(a, io)
    close(io)

    io = open("io.jdf", "r")
    aa = column_loader(StringVector{String}, io, metadata)
    close(io)

    df = DataFrame(a = a)
    savejdf("pls_del.jdf", df)
    df2 = loadjdf("pls_del.jdf")

    rm("io.jdf", force = true)
    rm("pls_del.jdf", force = true, recursive = true)
end

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 4ea0ccfa-9622-11eb-1cf0-e9ae2f927dd2
html"""

Section 2.7

 Discrete and Continuous

Lecture Video

"""

# ╔═╡ d155ea12-9628-11eb-347f-7754a33fd403
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
        Pkg.PackageSpec(name="HypertextLiteral", version="0.6"),
        Pkg.PackageSpec(name="LightGraphs", version="1"),
        Pkg.PackageSpec(name="GraphPlot", version="0.4"),
		Pkg.PackageSpec(name="SpecialFunctions", version="1")
    ])
    using Plots, PlutoUI, HypertextLiteral, LightGraphs, GraphPlot, Printf, SpecialFunctions
end

# ╔═╡ 01506de2-918a-11eb-2a4d-c554a6e54631
TableOfContents(title="📚 Table of Contents", aside=true)

# ╔═╡ 877deb2c-702b-457b-a54b-f27c277928d4
md"""
## Julia concepts
- printing using css fancy stuff (not really julia)
## pedagogical concepts
- curiosity based learning
- bridging what is often two different communities
"""

# ╔═╡ ee349b52-9189-11eb-2b86-b5dc15ebe432
md"""
###  Discrete and Continuous
"""

# ╔═╡ 43e39a6c-918a-11eb-2408-93563b4fb8c1
md"""
An exact technical definition of discrete and continuous can be difficult,
Nonetheless the idea of discrete mathematics is associated with  finite or countably many values that are isolated.  The set {1,...,n} and the integers are both discrete sets.
"""

# ╔═╡ 719a4c8c-9615-11eb-3dd7-7fb786f7fa17
md"""
DISCRETE MATH OBJECTS (examples):

1. Finite Sets: [1,2,...,100]

2. Infinite Discrete Sets: ``\mathbb{Z}`` = integers = ``\{\ldots,-2,-1,0,1,2,\ldots\}``

3. Graphs:
"""

# ╔═╡ 45ecee7e-970e-11eb-22fd-01f56876684e
gplot( barabasi_albert(150, 2) )

# ╔═╡ 52fa7f18-757a-4bf5-b851-32a1fca9c378
Pkg.add("GraphPlot")

# ╔═╡ 61ffe0f2-9615-11eb-37d5-f9e30a31c111
md"""
By contrast, entire intervals or the whole real line is associated with continuous mathematics.

In fact, mathematicians have worked painstakingly to define these terms, inventing fields such as point set topology that among other things can recognize the discrete, and analysis to rigorously work with the continuous.
"""

# ╔═╡ 627f6db6-9617-11eb-0453-a1f9e341ecfe
md"""
**Continuous real line:**

![continuous](data:image/png;base64,[encoded data, 9660 characters])
"""

# ╔═╡ 091a8a44-918c-11eb-2ee3-9be84a311afd
md"""
### Heard in the hallways: I only like discrete math.  I only like continuous math.
"""

# ╔═╡ 173b44ea-918c-11eb-116b-0bbaeffc3fe2
md"""
It is not unusual for students (and professors) to gravitate towards the discrete or the continuous.  We wish to point out, that the discrete and the continuous are so closely related, that it is worthwhile to be comfortable with both.  Up until fairly recently, much of computer science was often associated with discrete mathemtics, while computational science and engineering was associated with physical systems, hence continuous mathematics.

$(blue("That is blurring these days:"))  The popularity of machine learning has brought continuous optimization ideas such as gradient descent into the world of computer science and the impact of the physical world on us all (e.g. climate change, pandemics) is motivating applications in computer science.  The newfound popularity of Data science and statistics is also mixing the discrete with the continuous.

"""

# ╔═╡ a3f005a8-9617-11eb-1503-75c31ec54f70
md"""
$(blue("Continuous math often lets you replace complicated large systems
		with lots of details with a simpler abstraction that is easier to work with."))
"""

# ╔═╡ 870cdf5f-f896-4060-9548-5d9c1749d100
md"""
$(blue("The combination of continuous and discrete is often more useful than either one alone."))
"""

# ╔═╡ d9dfe7c5-9211-4707-bb33-a3ff258e10f4
md"""
$(blue("Machine Learning, Epidemics, climate change, etc. show how critical continuous math is these days."))
"""

# ╔═╡ 5c536430-9188-11eb-229c-e7feba62d257
md"""
### Indexing and Function Evaluation

Analogy: ``v_i`` (ith element of v) vs. $f(x)$  (evaluate f at x):

These are different, right?
In the one you are extracting an element, and in the other you are applying what in some high schools might be called a "function machine."

However, a moment's thought tells you that a vector is really a discrete function, in that the argument can take on the values i = 1,2,...,n and the evaluation is $v_i$. That's a function.

In fact, think of a range object such as 2:2:20.  You could think of this as just a shorthand for the vector [2,4,...,20] but in fact when you index into this "vector" like thing, you are indeed explicitly evaluating a function i.e. i->2i.
"""

# ╔═╡ 1e8ea849-40b7-41fd-b17f-cd2d991d5c24
[2:2:20;]  # this expands the "iterator" into an ordinary vector

# ╔═╡ 679a39ee-99a5-4211-9adc-8296d499e37e
[2:2:20;][7] # Extracts an element from Memory (of course there is an address calculation)

# ╔═╡ 2c64f98d-dc84-4fa5-81ce-25b319ff9583
(2:2:20)[7] # Compute 2*4 (more or less)

# ╔═╡ 0a379cae-386d-4daa-ab6f-9d0424c1cdc1
begin
	f(x)=2x
	f(7)     # Compute 2*7
end

# ╔═╡ 890c0fa2-c247-4f14-84f6-2bed69d0f0c5
md"""
Any which way $v$ is a function "machine" whose input is $\{1,2,3,4,5,6,7,8,9,10\}$
"""

# ╔═╡ 68b60d09-acee-48d8-8bb1-7ab4faa6b785
gr()

# ╔═╡ 40095ad2-961f-11eb-1f23-83d1a381ace7
md"""
### Area
"""

# ╔═╡ ed71b026-9565-11eb-1058-d77efe114562
md"""
Area of a circle using regular polygons:
"""

# ╔═╡ 3b84bb0a-9566-11eb-1c1f-e30ca7330c09
md"""
n = $(@bind sides Slider(3:100, show_value=true, default=6))
"""

# ╔═╡ 02784976-9566-11eb-125c-a7f1f1bafd6b
begin
	θ = (0:.01:1)*2π
		plot( cos.(θ),sin.(θ), ratio=1, axis=false, legend=false, ticks=false, lw=4, color=:black, fill=false)
	plot!( cos.(θ),sin.(θ), ratio=1, axis=false, legend=false, ticks=false, lw=4, color=:white, fill=true, alpha=.6)

	ϕ = (0:sides)*2π/sides
	for i=1:sides
	   plot!( Shape( [0,cos(ϕ[i]),cos(ϕ[i+1])],[0,sin(ϕ[i]),sin(ϕ[i+1])]), fill=true,lw=0)
	end
	title!("Area = ($sides/2)sin(2π/$sides) ≈  $(area(sides)/π )  π")
end

# ╔═╡ f20da096-9712-11eb-2a67-cd33f6ab8750
area(s) = (s/2) * sin(2π/s)

# ╔═╡ 6fd93018-c33b-4682-91c3-7a20a41d9b03
area0 = area.( 2 .^ (2:10) )

# ╔═╡ a306559f-e095-4f6d-94e8-b0be160e77fa
π

# ╔═╡ ea29e286-4b4a-4291-a093-cd942ba46e49
md"""
A carefully chosen convolution: [-1/3,4/3]
"""

# ╔═╡ 103c93ae-8175-4996-ab8f-5d537691defc
area1 = [ 4/3 * area0[i+1] .-  1/3 * area0[i] for i = 1:length(area0)-1 ]

# ╔═╡ 686904c9-1cc4-4476-860b-159e56471e38
function colorgoodbad(should_be, given)
	indexofmistake = something(
		findfirst(collect(should_be) .!== collect(given)),
		length(given)+1,
	)
	@htl("""
		$(given[1:indexofmistake-1])$(given[indexofmistake:end])
		""")
end

# ╔═╡ bcfd1585-8161-43a2-8b19-ed654df2e0e1
colorgoodbad(string(float(π)) , string(22/7))

# ╔═╡ d2d1366b-9b6d-4e54-a0c4-7087f5f063c4
pyramid( [area0,area1], horizontal = true)

# ╔═╡ a76ac67b-27b9-4e2b-9fca-61480dca5264
area2 = [16/15 * area1[i+1] .-  1/15 * area1[i] for i = 1:length(area1)-1 ]

# ╔═╡ c742742a-765b-4eb5-bd65-dc0cd6328255
md"""
Another carefully chosen convolution: [-1/15,16/15], do you see the pattern?
"""

# ╔═╡ 6577e546-8f0b-413a-a8bb-b9c12803199d
pyramid([area0,area1,area2], horizontal = true)

# ╔═╡ 5273fe09-fe38-4c88-b84a-51af17cff906
big(π)

# ╔═╡ 4dd03325-2498-4fe7-9212-f964081a0300
area3 = [64/63 * area2[i+1] .-  1/63 * area2[i] for i = 1:length(area2)-1 ]

# ╔═╡ 626242ea-544c-49fc-9884-c70dd6800902
area4 = [128/127 * area3[i+1] .-  1/127 * area3[i] for i = 1:length(area3)-1 ]

# ╔═╡ 893a56b0-f5d0-4f8d-ba15-1048180a7e53
pyramid([pp.(area0), pp.(area1), pp.(area2), pp.(area3), pp.(area4)], horizontal = true)

# ╔═╡ dbccc2d5-c2af-48c4-8726-a95c09da78ae
md"""
Why does this work?
"""

# ╔═╡ 82a407b6-0ecb-4011-a0f6-bc9e1f51393f
md"""
Area(s) = ``(s/2) \sin (2\pi/2) = \pi- \frac{2\pi^3}{3} s^{-2}  + \frac{2\pi^5}{15}s^{-4} -  \frac{4\pi^7}{315} s^{-6} + \ldots``

as `` s \rightarrow \infty ``.
"""

# ╔═╡ 5947dc80-9714-11eb-389d-1510a1137a50
md"""
Area(s) = `` \pi  -  c_1 /s^2  + c_2 / s^4  - \ldots``

Area(2s) = `` \pi  -  c_1 /(4s^2)  + c_2 / (16s^4)  - \ldots``
"""

# ╔═╡ db8121ec-8546-4f1e-8153-cff5b4df39df
md"""
Think about taking (4/3) Area(2s) - (1/3) Area(s).

Now we have
``\pi + c s^{-4}`` as the leading term so doubling the s approximately reduces the area error by 16, when before it was only 4. etc.
"""

# ╔═╡ d4f83a20-62cf-47f1-a622-d5c4c34e4813
areab(s) = (s/2) * sin(big(2)*big(π)/s)

# ╔═╡ 01631c38-9713-11eb-30bf-3d23cf0d3dc8
begin
	area0b = areab.(big.([2,4,8,16,32,62,128,256,512,1024,2048,4096,8192,16384,32768,65536]))
	area0b[end]
end

# ╔═╡ 553bdb0a-9714-11eb-1646-413a969d6884
begin
	area1b = [ 4//3 * area0b[i+1] .-  1//3 * area0b[i] for i = 1:length(area0b)-1 ]
	area1b[end]
end

# ╔═╡ 453f2585-157d-490a-9d1c-0b02939d0a11
begin
	area2b = [16//15 * area1b[i+1] .-  1//15 * area1b[i] for i = 1:length(area1b)-1 ]
	area2b[end]
end

# ╔═╡ bc1efddd-c959-4407-9a86-ba73a64508a8
begin
	area3b = [64//63 * area2b[i+1] .-  1//63 * area2b[i] for i = 1:length(area2b)-1 ]
	area3b[end]
end

# ╔═╡ 516b69d8-5d94-4b4d-9596-2db0dfbf4038
begin
	area4b = [256//255 * area3b[i+1] .-  1//255 * area3b[i] for i = 1:length(area3b)-1 ]
	area4b[end]
end

# ╔═╡ cec13915-8adb-4627-b220-591377239997
begin
	area5b = [1024//1023 * area4b[i+1] .-  1//1023 * area4b[i] for i = 1:length(area4b)-1 ]
	area5b[end]
end

# ╔═╡ 6d954628-6290-4867-8144-dd486551545d
begin
	area6b = [4096//4095 * area5b[i+1] .-  1//4095 * area5b[i] for i = 1:length(area5b)-1 ]
	area6b[end]
end

# ╔═╡ 00478b2c-5dcc-44fc-a7be-a3dadf6300e7
begin
	area7b = [16336//16335 * area6b[i+1] .-  1//16335 * area6b[i] for i = 1:length(area6b)-1 ]
	area7b[end]
end

# ╔═╡ 37fc6e56-9714-11eb-1427-b75613800366
big(π)

# ╔═╡ 4a072870-961f-11eb-1215-17efa0013873
md"""
Area using inscribed squares
"""

# ╔═╡ de9066e2-d5eb-49e3-be71-edda8e8e31dd
@bind s Slider(2:40, show_value=true)

# ╔═╡ 4d4705d0-9568-11eb-085c-0fc556c4cfe7
let

    plot()
	for i=-s:s
		plot!([i/s,i/s],[-1,1],color=RGB(0,1,0),lw=1)
		plot!([-1,1],[i/s,i/s],color=RGB(0,1,0),lw=1)
	end
		P = plot!( cos.(θ),sin.(θ), ratio=1, axis=false, legend=false, ticks=false, lw=3, color=:black)
	plot!(P)

	h = 1/s
	a = 0

	xx=  floor(√2/2h)
	x = xx*h
	y=x
	plot!( Shape([-x, -x, x ,x],[-y, y ,y, -y]), color=RGB(1,0,0),alpha=.7)

	a = a+Int(2*xx)^2

	 for i=-s:(-xx-1), j=-s:(-1)
	   x = i*h
	    y = j*h
	   if (x^2+y^2≤1) & ( (x+h)^2+(y+h)^2 ≤1) & (x^2+(y+h)^2 ≤1) & ((x+h)^2+y^2 ≤1)
	 	 plot!( Shape([x, x, x+h ,x+h],[y, y+h ,y+h, y]), color=:blue)
		 plot!( Shape([-x-h, -x-h, -x ,-x],[y, y+h ,y+h, y]), color=:blue)
	     plot!( Shape([x, x, x+h ,x+h],[-y-h, -y ,-y, -y-h]), color=:blue)
		 plot!( Shape([-x-h, -x-h, -x ,-x],[-y-h, -y ,-y, -y-h]), color=:blue)
		 plot!( Shape([y, y+h ,y+h, y],[x, x, x+h ,x+h]), color=:blue)
		 plot!( Shape([-y-h, -y ,-y, -y-h],[x, x, x+h ,x+h]), color=:blue)
		 plot!( Shape([y, y+h ,y+h, y],[-x-h, -x-h, -x ,-x]), color=:blue)
		 plot!( Shape([-y-h, -y ,-y, -y-h],[-x-h, -x-h, -x ,-x]), color=:blue)
	 		a += 8
	 	end
	 end

	xlabel!("s  =  $s")

	title!( "$(a//s^2) =  $(a*h^2/π) π")
	plot!()

end

# ╔═╡ 155241b0-9646-11eb-180e-89c8651536c6
@bind j Slider(1:9, show_value=true, default=6)

# ╔═╡ 4f845436-9646-11eb-2445-c12746a9e556
begin
	N  = 1024
	h =  1/N
	v = randn(N)
end

# ╔═╡ 31d56008-9646-11eb-1985-2b68af354773
J = N ÷ 2^j

# ╔═╡ 1761187e-9645-11eb-3778-b132f856696d
begin
	plot()
	c = [0;cumsum(v)] .* √h
	plot!((0:N)./N,c)
	plot!( (0:J:N)./N,   c[1:J:end],legend=false,m=:o,ms=1, color=:red, lw=5)
	plot!(ylims=(-2,2))
end

# ╔═╡ 1e18f95c-cd53-4ede-8d93-572c81f872da
md"""
A random walk is a discrete random function. It is defined at grid points.  Brownian motion is a continuous random function.  It is defined on an entire interval.
If one has an instance of a Brownian motion, you can say its exact value at, say, .7.
If one looks at the random variable that represents evaluation at .7 it is a normal distribution.
"""

# ╔═╡ c32e0f9c-918e-11eb-1cf9-a340786db24a
md"""
Alan's essay:

In what sense does the continuous even exist?  The fact of the matter is that there are limits that give the same answer no matter how you get there, and these limits
are important to us. For example, no matter how you cover an area, by little rectangles, the sum always converges to what we intuitively call area.
The normal distribution is interesting in that no matter which starting finite distribution we might take, if add n independent copies and normalize to variance 1 we get the same limit.  Again, there are so many ways to start, and yet we always end up with the same thing.  Continuous mathematics is full of so many examples, where discrete objects end up behaving the same.

Indeed what happens as discrete objects get larger and larger, their complexity gets out of control if one wants to keep track of every detail, but they get simpler in their aggregate behavior.
"""

# ╔═╡ 632eea46-9710-11eb-1abe-85da8d9c30a9
f(x,t) =  exp(-x^2/t)/√(π*t)

# ╔═╡ 9c519eca-9710-11eb-20dc-3f76801545d1
@bind t Slider(.01:.01:8, show_value=true)

# ╔═╡ 7c4b82c8-9710-11eb-101e-53616e278289
begin
	x = -3 : .01 : 3
	plot( x, f.(x,t), ylims=(0,1), legend=false)
end

# ╔═╡ 021d7e9a-9711-11eb-063b-11441afa2e69
begin
	surface(-2:.05:2, .2:.01:1, f, alpha=.4, c=:Reds, legend=false)
	for t = .2 : .1 :1
		plot!(-2:.05:2, fill(t,length(-2:.05:2)),  f.(-2:.05:2,t), c=:black)
	end
	xlabel!("x")
	ylabel!("t")
	plot!()
end

# ╔═╡ 653a1cd6-9711-11eb-2284-770f00bebeec
plotly()

# ╔═╡ bb8dc4fe-918d-11eb-2bde-bb00c47a1c27
md"""
### Sum and a (Definite) Integral
"""

# ╔═╡ c4a3bf6c-918d-11eb-1d50-911f83b6df81
md"""
### Cumsum and an Indefinite Integral
"""

# ╔═╡ d99dc494-918d-11eb-2733-29ce93ba584e
md"""
### Discrete (Finite Differencing) Filters and Derivatives/Gradients
"""

# ╔═╡ 0fb84ff2-918e-11eb-150f-8dad121c87bc
md"""
### Discrete and continuous convolutions e.g. probability densities
"""

# ╔═╡ a7c5ef96-918d-11eb-0632-f94386eb64f2
md"""
### Discrete Random Walks and Brownian Motion
"""

# ╔═╡ 75672be6-918d-11eb-1e10-07fbcc72abbd
md"""
### Binomial Distribution and the Normal Distribution
"""

# ╔═╡ 906758c6-918d-11eb-08ae-b3c4f7870b4e
md"""
### Discrete and the Continuous Fourier Transform
"""

# ╔═╡ b972b218-970c-11eb-1949-535830e20990
P = [ binomial(n,k) for n=0:5,k=0:5]

# ╔═╡ 80682786-970d-11eb-223b-b5f762b19c24
P*P'

# ╔═╡ ed6d7404-970c-11eb-13ee-5f5a454d2222
begin
	A = (1 ./beta.( (1:6)', 0:5) ) ./ [1;1:5;]
	A[1,:] .= 1
	round.(Int,A)
end

# ╔═╡ 0e6cab25-70f8-46ab-a5ab-8542e232274e
function blue(s::String)
	HTML(" $(s)  ")
end

# ╔═╡ c03d45f8-9188-11eb-2e11-0fafa39f253d
function pyramid(rows::Vector{<:Vector};
		horizontal=false,
		padding_x=8,
		padding_y=2,
	)

	style = (; padding = "$(padding_y)px $(padding_x)px")
	render_row(xs) = @htl("""
<padder></padder>$([@htl("$(x)") for x in xs])<padder></padder>
""")

	@htl("""

		$(render_row.(rows))

		""")
end

# ╔═╡ 43d20d56-d56a-47a8-893e-f726c1a99651
pp(x) =  colorgoodbad( string(float(π)) , (@sprintf "%.15f" x) )

# ╔═╡ Cell order:
# ╟─4ea0ccfa-9622-11eb-1cf0-e9ae2f927dd2
# ╠═d155ea12-9628-11eb-347f-7754a33fd403
# ╠═01506de2-918a-11eb-2a4d-c554a6e54631
# ╟─877deb2c-702b-457b-a54b-f27c277928d4
# ╟─ee349b52-9189-11eb-2b86-b5dc15ebe432
# ╟─43e39a6c-918a-11eb-2408-93563b4fb8c1
# ╟─719a4c8c-9615-11eb-3dd7-7fb786f7fa17
# ╠═45ecee7e-970e-11eb-22fd-01f56876684e
# ╠═52fa7f18-757a-4bf5-b851-32a1fca9c378
# ╟─61ffe0f2-9615-11eb-37d5-f9e30a31c111
# ╟─627f6db6-9617-11eb-0453-a1f9e341ecfe
# ╟─091a8a44-918c-11eb-2ee3-9be84a311afd
# ╟─173b44ea-918c-11eb-116b-0bbaeffc3fe2
# ╟─a3f005a8-9617-11eb-1503-75c31ec54f70
# ╟─870cdf5f-f896-4060-9548-5d9c1749d100
# ╟─d9dfe7c5-9211-4707-bb33-a3ff258e10f4
# ╟─5c536430-9188-11eb-229c-e7feba62d257
# ╠═1e8ea849-40b7-41fd-b17f-cd2d991d5c24
# ╠═679a39ee-99a5-4211-9adc-8296d499e37e
# ╠═2c64f98d-dc84-4fa5-81ce-25b319ff9583
# ╠═0a379cae-386d-4daa-ab6f-9d0424c1cdc1
# ╟─890c0fa2-c247-4f14-84f6-2bed69d0f0c5
# ╠═68b60d09-acee-48d8-8bb1-7ab4faa6b785
# ╟─40095ad2-961f-11eb-1f23-83d1a381ace7
# ╟─ed71b026-9565-11eb-1058-d77efe114562
# ╟─3b84bb0a-9566-11eb-1c1f-e30ca7330c09
# ╟─02784976-9566-11eb-125c-a7f1f1bafd6b
# ╠═f20da096-9712-11eb-2a67-cd33f6ab8750
# ╠═6fd93018-c33b-4682-91c3-7a20a41d9b03
# ╠═a306559f-e095-4f6d-94e8-b0be160e77fa
# ╟─ea29e286-4b4a-4291-a093-cd942ba46e49
# ╠═103c93ae-8175-4996-ab8f-5d537691defc
# ╠═686904c9-1cc4-4476-860b-159e56471e38
# ╠═bcfd1585-8161-43a2-8b19-ed654df2e0e1
# ╠═d2d1366b-9b6d-4e54-a0c4-7087f5f063c4
# ╠═a76ac67b-27b9-4e2b-9fca-61480dca5264
# ╟─c742742a-765b-4eb5-bd65-dc0cd6328255
# ╠═6577e546-8f0b-413a-a8bb-b9c12803199d
# ╠═5273fe09-fe38-4c88-b84a-51af17cff906
# ╠═4dd03325-2498-4fe7-9212-f964081a0300
# ╠═626242ea-544c-49fc-9884-c70dd6800902
# ╠═893a56b0-f5d0-4f8d-ba15-1048180a7e53
# ╟─dbccc2d5-c2af-48c4-8726-a95c09da78ae
# ╟─82a407b6-0ecb-4011-a0f6-bc9e1f51393f
# ╟─5947dc80-9714-11eb-389d-1510a1137a50
# ╟─db8121ec-8546-4f1e-8153-cff5b4df39df
# ╠═d4f83a20-62cf-47f1-a622-d5c4c34e4813
# ╠═01631c38-9713-11eb-30bf-3d23cf0d3dc8
# ╠═553bdb0a-9714-11eb-1646-413a969d6884
# ╠═453f2585-157d-490a-9d1c-0b02939d0a11
# ╠═bc1efddd-c959-4407-9a86-ba73a64508a8
# ╠═516b69d8-5d94-4b4d-9596-2db0dfbf4038
# ╠═cec13915-8adb-4627-b220-591377239997
# ╠═6d954628-6290-4867-8144-dd486551545d
# ╠═00478b2c-5dcc-44fc-a7be-a3dadf6300e7
# ╠═37fc6e56-9714-11eb-1427-b75613800366
# ╟─4a072870-961f-11eb-1215-17efa0013873
# ╠═de9066e2-d5eb-49e3-be71-edda8e8e31dd
# ╟─4d4705d0-9568-11eb-085c-0fc556c4cfe7
# ╠═155241b0-9646-11eb-180e-89c8651536c6
# ╠═4f845436-9646-11eb-2445-c12746a9e556
# ╠═31d56008-9646-11eb-1985-2b68af354773
# ╠═1761187e-9645-11eb-3778-b132f856696d
# ╟─1e18f95c-cd53-4ede-8d93-572c81f872da
# ╟─c32e0f9c-918e-11eb-1cf9-a340786db24a
# ╠═632eea46-9710-11eb-1abe-85da8d9c30a9
# ╠═9c519eca-9710-11eb-20dc-3f76801545d1
# ╠═7c4b82c8-9710-11eb-101e-53616e278289
# ╠═021d7e9a-9711-11eb-063b-11441afa2e69
# ╠═653a1cd6-9711-11eb-2284-770f00bebeec
# ╟─bb8dc4fe-918d-11eb-2bde-bb00c47a1c27
# ╟─c4a3bf6c-918d-11eb-1d50-911f83b6df81
# ╟─d99dc494-918d-11eb-2733-29ce93ba584e
# ╟─0fb84ff2-918e-11eb-150f-8dad121c87bc
# ╟─a7c5ef96-918d-11eb-0632-f94386eb64f2
# ╟─75672be6-918d-11eb-1e10-07fbcc72abbd
# ╟─906758c6-918d-11eb-08ae-b3c4f7870b4e
# ╠═b972b218-970c-11eb-1949-535830e20990
# ╠═80682786-970d-11eb-223b-b5f762b19c24
# ╠═ed6d7404-970c-11eb-13ee-5f5a454d2222
# ╠═0e6cab25-70f8-46ab-a5ab-8542e232274e
# ╟─c03d45f8-9188-11eb-2e11-0fafa39f253d
# ╟─43d20d56-d56a-47a8-893e-f726c1a99651


text:
module MortalityTables
import XMLDict
import DataStructures

include("tableProcessor.jl")

"""
    qx(table::MortalityTable, issueAge::Int, duration::Int)

Given a mortality table, an issue age, and a duration, returns the appropriate select or ultimate qx.

If the table is a select and ultimate table, it will return the select mortality during the select period.
If the table does not have select rates, will just return the rate for the age at the given duration.
"""

function qx(table::MortalityTable, issueAge::Int, duration::Int)
    try
        table.select[issueAge+1,duration] # +1 because julia index starts with 1
    catch err
        if isa(err,BoundsError)
            missing
        else
            throw(err)
        end
    end
end

"""
    qx(table::XTbMLTable, age)

Given a mortality table and an age returns the qx.

If the table is a select and ultimate table, it will return the select mortality during the ulitmate period because duration was not given.

If the table does not have select rates, will just return the rate for the given age.

"""
function qx(table::MortalityTable, age)
    try
        table.ultimate[age+1] # +1 because julia index starts with 1
    catch err
        if isa(err,BoundsError)
            missing
        else
            throw(err)
        end
    end
end

end # module

export Tables,loadXTbMLTable, qx


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   General structures related to orbits.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export Orbit, KeplerianElements, OrbitStateVector

"""
    Orbit

Abstract type of an orbit representation.

"""
abstract type Orbit{Tepoch, T} end

"""
    KeplerianElements{Tepoch, T}

This structure defines the orbit in terms of the Keplerian elements.

# Fields

- `t::Tepoch`: Epoch.
- `a::T`: Semi-major axis [m].
- `e::T`: Eccentricity [ ].
- `i::T`: Inclination [rad].
- `Ω::T`: Right ascension of the ascending node [rad].
- `ω::T`: Argument of perigee [rad].
- `f::T`: True anomaly [rad].
"""
@with_kw_noshow struct KeplerianElements{Tepoch, T} <: Orbit{Tepoch, T}
    t::Tepoch
    a::T
    e::T
    i::T
    Ω::T
    ω::T
    f::T
end

function KeplerianElements(
    t::Number,
    a::T1,
    e::T2,
    i::T3,
    Ω::T4,
    ω::T5,
    f::T6
) where {T1<:Number, T2<:Number, T3<:Number, T4<:Number, T5<:Number, T6<:Number}
    T = promote_type(T1, T2, T3, T4, T5, T6)
    return KeplerianElements{typeof(t), T}(t, a, e, i, Ω, ω, f)
end

"""
    OrbitStateVector{Tepoch, T}

Store the state vector representation of an orbit.

# Fields

- `t::Tepoch`: Epoch [Julian Day].
- `r::SVector{3, T}`: Position vector [m].
- `v::SVector{3, T}`: Velocity vector [m/s].
- `a::SVector{3, T}`: Acceleration vector [m/s²].
"""
@with_kw_noshow struct OrbitStateVector{Tepoch, T} <: Orbit{Tepoch, T}
    t::Tepoch        = 0
    r::SVector{3, T} = SVector{3, T}(0, 0, 0)
    v::SVector{3, T} = SVector{3, T}(0, 0, 0)
    a::SVector{3, T} = SVector{3, T}(0, 0, 0)
end


text:
module Rules
    using TermInterface
    using ..Util
    using ..EMatchCompiler
    import ..closure_generator

    include("rule_types.jl")
    include("rule_cache.jl")
    include("rule_dsl.jl")
    include("exports.jl")
end


text:
export devec,demat,mat

doc"""
`ArraySpace` used to represent array-valued expansions in `space`.  The
coefficients are of each entry are interlaced.
"""
immutable ArraySpace{S,n,T,DD,dim} <: Space{T,DD,dim}
     space::S
     dimensions::NTuple{n,Int}
#      # for AnyDomain() usage
    ArraySpace(sp::S,dims)=new(sp,dims)
    ArraySpace(d::Domain,dims)=new(S(d),dims)
end

InterlaceIterator(sp::ArraySpace) = InterlaceIterator(fill(dimension(sp.space),length(sp)))
interlacer(sp::ArraySpace) = InterlaceIterator(sp)

typealias VectorSpace{S,T,DD,dim} ArraySpace{S,1,T,DD,dim}
typealias MatrixSpace{S,T,DD,dim} ArraySpace{S,2,T,DD,dim}

ArraySpace(S::Space,n::Tuple{Vararg{Int}}) =
    ArraySpace{typeof(S),length(n),basistype(S),
               domaintype(S),domaindimension(S)}(S,n)
ArraySpace(S::Space,n::Integer) = ArraySpace(S,(n,))
ArraySpace(S::Space,n,m) =
    ArraySpace{typeof(S),2,basistype(S),
               domaintype(S),domaindimension(S)}(S,(n,m))
ArraySpace(d::Domain,n...) = ArraySpace(Space(d),n...)

Base.length{SS}(AS::ArraySpace{SS,1}) = AS.dimensions[1]
Base.length(AS::ArraySpace) = *(AS.dimensions...)

Base.length{AS<:ArraySpace}(f::Fun{AS}) = length(space(f))

Base.size(AS::ArraySpace)=AS.dimensions
Base.size(AS::ArraySpace,k)=AS.dimensions[k]

Base.size{AS<:ArraySpace}(f::Fun{AS},k...) = size(space(f),k...)

function Base.reshape(AS::VectorSpace,k,j)
    @assert length(AS)==k*j
    ArraySpace(AS.space,(k,j))
end

dimension(AS::ArraySpace) = dimension(AS.space)*length(AS)

domain(AS::ArraySpace)=domain(AS.space)

## transforms

transform{SS,V}(AS::ArraySpace{SS,1},vals::Vector{Vector{V}})=transform(AS,transpose(hcat(vals...)))

function transform{SS,T,V<:Number}(AS::ArraySpace{SS,1,T},M::Array{V,2})
    n=length(AS)

    @assert size(M,2)==n
    plan = plan_transform(AS.space,M[:,1])
    cfs=Vector{V}[transform(AS.space,M[:,k],plan)  for k=1:size(M,2)]

    C=zeros(coefficient_type(T,V),mapreduce(length,max,cfs)*size(M,2))

    for k=1:size(M,2)
        C[k:n:end]=cfs[k]
    end
    C
end

# transform of array is same order as vectorizing and then transforming
transform{SS,n,V}(AS::ArraySpace{SS,n},vals::Vector{Array{V,n}})=transform(vec(AS),map(vec,vals))

Base.vec(AS::ArraySpace)=ArraySpace(AS.space,length(AS))
function Base.vec{S<:Space,V,T,DD,d}(f::Fun{VectorSpace{S,V,DD,d},T})
    n=length(space(f))
    Fun{S,T}[Fun(f.coefficients[j:n:end],space(f).space) for j=1:n]
end
Base.vec{AS<:ArraySpace,T}(f::Fun{AS,T})=vec(Fun(f.coefficients,vec(space(f))))

mat{AS<:ArraySpace,T}(f::Fun{AS,T})=reshape(vec(f),size(space(f))...)

# mat(f,1) vectorizes columnwise
function mat{S,V,T,DD,d}(f::Fun{MatrixSpace{S,V,DD,d},T},j::Integer)
    @assert j==1
    m=mat(f)
    r=Array(Fun{VectorSpace{S,V,DD,d},T},1,size(m,2))
    for k=1:size(m,2)
        r[1,k]=devec(m[:,k])
    end
    r
end

spaces(A::ArraySpace) = fill(A.space,A.dimensions)

TupleSpace{SS}(A::ArraySpace{SS,1}) = TupleSpace(spaces(A))

Base.getindex{S,V,DD,d}(f::Fun{VectorSpace{S,V,DD,d}},k...)=vec(f)[k...]
Base.getindex{S,V,DD,d}(f::Fun{MatrixSpace{S,V,DD,d}},k...)=mat(f)[k...]

Base.getindex(S::ArraySpace,k...)=S.space

Base.start(S::ArraySpace) = 1
Base.next(S::ArraySpace,k) = S.space,k+1
Base.done(S::ArraySpace,k) = k>length(S)
Base.endof(S::ArraySpace) = length(S)

#support tuple set
for OP in (:(Base.start),:(Base.done),:(Base.endof))
    @eval $OP{SS<:ArraySpace}(f::Fun{SS},k...)=$OP(space(f),k...)
end

Base.next{SS<:ArraySpace}(f::Fun{SS},k)=f[k],k+1

function devec{F<:Fun}(v::Vector{F})
    sps=map(space,v)
    if spacesequal(sps)
        Fun(vec(coefficientmatrix(v).'),ArraySpace(first(sps),length(v)))
    else
        Fun(vec(coefficientmatrix(v).'),TupleSpace(sps))
    end
end

devec(v::Vector{Any})=devec([v...])

function devec{S<:Space}(spl::Vector{S})
    #TODO: Redesign
    @assert spacescompatible(spl)
    ArraySpace(first(spl),length(spl))
end

function demat{S,T}(v::Array{Fun{S,T}})
    ff=devec(vec(v))  # A vectorized version
    Fun(coefficients(ff),ArraySpace(space(ff).space,size(v)...))
end

demat(v::Vector{Any})=devec(v)

function demat{S,T,V,DD,d}(A::Array{Fun{VectorSpace{S,T,DD,d},V},2})
    @assert size(A,1)==1

    M=Array(Fun{S,V},length(space(A[1])),size(A,2))
    for k=1:size(A,2)
        M[:,k]=vec(A[k])
    end
    demat(M)
end

Fun{F<:Fun}(V::AbstractVector{F})=devec(V)
Fun{F<:Fun}(V::AbstractMatrix{F})=demat(V)

Fun(v::Vector{Any},sp::ArraySpace) = devec(map(f->Fun(f,sp.space),v))

function union_rule{S,n,T,DD,dim,S2,T2,DD2}(a::ArraySpace{S,n,T,DD,dim},b::ArraySpace{S2,n,T2,DD2,dim})
    if a.dimensions==b.dimensions
        sp=union(a.space,b.space)
        if !isa(sp,NoSpace)
            return ArraySpace(sp,a.dimensions)
        end
    end

    NoSpace()
end

## routines

spacescompatible(AS::ArraySpace,BS::ArraySpace)=size(AS)==size(BS) && spacescompatible(AS.space,BS.space)
canonicalspace(AS::ArraySpace)=ArraySpace(canonicalspace(AS.space),size(AS))
evaluate(f::AbstractVector,S::ArraySpace,x)=map(g->evaluate(g,x),mat(Fun(f,S)))

for OP in (:(Base.transpose),)
    @eval $OP{AS<:ArraySpace,T}(f::Fun{AS,T})=demat($OP(mat(f)))
end

Base.reshape{AS<:ArraySpace}(f::Fun{AS},k...)=Fun(f.coefficients,reshape(space(f),k...))

Base.diff{AS<:ArraySpace,T}(f::Fun{AS,T},n...)=demat(diff(mat(f),n...))

## conversion

function coefficients(f::Vector,a::VectorSpace,b::VectorSpace)
    n=length(a)
    @assert n==length(b)
    A=a.space;B=b.space
    ret=Array(eltype(f),length(f))
    for k=1:n
        ret[k:n:end]=coefficients(f[k:n:end],A,B)
    end
    ret
end

coefficients{F<:Fun}(Q::Vector{F},rs::VectorSpace) = vec(coefficientmatrix(Q,rs.space).')

Fun{FF<:Fun}(f::Vector{FF},d::VectorSpace) = Fun(coefficients(f,d),d)
Fun{FF<:Fun}(f::Matrix{FF},d::MatrixSpace) = Fun(coefficients(f,d),d)

## constructor

# change to ArraySpace
Fun{AS<:ArraySpace}(f::Fun{AS},d::ArraySpace)=Fun(coefficients(f,d),d)
Fun{AS<:ArraySpace}(f::Fun{AS},d::Space)=Fun(f,ArraySpace(d,space(f).dimensions))

# columns are coefficients
Fun{T<:Number}(M::Array{T,2},sp::Space)=devec([Fun(M[:,k],sp) for k=1:size(M,2)])

# Automatically change to ArraySpace
# A is interpreted as coefficients
function Fun{T<:Number}(A::Array{T,2},sp::VectorSpace)
    n=length(sp)
    m=size(A,2)

    cfs=Array(T,m*n*(div(size(A,1),n)+1))
    for k=1:size(A,1),j=1:m
        cfs[m*n*div(k-1,n)+mod(k-1,n)+(j-1)*n+1]=A[k,j]
    end
    # pad with zeros
    for k=size(A,1)+1:n*(div(size(A,1),n)+1),j=1:m
        cfs[m*n*div(k-1,n)+mod(k-1,n)+(j-1)*n+1]=zero(T)
    end
    Fun(cfs,ArraySpace(sp.space,n,size(A,2)))
end

Base.ones{T<:Number}(::Type{T},A::ArraySpace)=demat(fill(ones(T,A.space),A.dimensions...))
Base.ones(A::ArraySpace)=demat(fill(ones(A.space),A.dimensions...))

## calculus

for op in (:differentiate,:integrate,:(Base.cumsum))
    @eval $op{V<:ArraySpace}(f::Fun{V})=demat(map($op,mat(f)))
end

function Base.det{A<:ArraySpace,V}(f::Fun{A,V})
    @assert size(space(f))==(2,2)
    m=mat(f)
    m[1,1]*m[2,2]-m[1,2]*m[2,1]
end

function Base.inv{A<:ArraySpace,T}(V::Fun{A,T})
    @assert size(space(V),1)==size(space(V),2)
    M=Multiplication(V,ArraySpace(space(V).space,size(space(V),1)))
    M\eye(size(space(V),2))
end

## Algebra

for OP in (:*,:.*,:+,:-)
    @eval begin
        $OP{T<:Number,AS<:ArraySpace,V}(A::Array{T},f::Fun{AS,V})=demat($OP(A,mat(f)))
        $OP{T<:Number,AS<:ArraySpace,V}(f::Fun{AS,V},A::Array{T})=demat($OP(mat(f),A))
        $OP{T,S,AS<:ArraySpace,V}(A::Vector{Fun{S,T}},f::Fun{AS,V})=demat($OP(A,mat(f)))
        $OP{T,S,AS<:ArraySpace,V}(f::Fun{AS,V},A::Vector{Fun{S,T}})=demat($OP(mat(f),A))
        $OP{T,S,AS<:ArraySpace,V}(A::Array{Fun{S,T}},f::Fun{AS,V})=demat($OP(A,mat(f)))
        $OP{T,S,AS<:ArraySpace,V}(f::Fun{AS,V},A::Array{Fun{S,T}})=demat($OP(mat(f),A))
        $OP{AS<:ArraySpace,V}(A::UniformScaling,f::Fun{AS,V})=demat($OP(A,mat(f)))
        $OP{AS<:ArraySpace,V}(f::Fun{AS,V},A::UniformScaling)=demat($OP(mat(f),A))
    end
end

for OP in (:*,:.*)
    @eval $OP{BS<:ArraySpace,T,AS<:ArraySpace,V}(A::Fun{BS,T},f::Fun{AS,V})=demat($OP(mat(A),mat(f)))
end

## ConstantVectorSpace

typealias ConstantVectorSpace VectorSpace{ConstantSpace{AnyDomain},RealBasis,AnyDomain,1}

function Base.vec{V,TT,DD,d,T}(f::Fun{SumSpace{Tuple{ConstantVectorSpace,V},TT,DD,d},T},k)
    m=length(space(f)[1])
    if k≤m
        Fun(f.coefficients[k],ConstantSpace())
    else
        Fun(f.coefficients[m+1:end],space(f)[2])
    end
end

Base.vec{V,TT,DD,d,T}(f::Fun{SumSpace{Tuple{ConstantVectorSpace,V},TT,DD,d},T})=Any[vec(f,k) for k=1:length(space(f)[1])+1]

## linsolve

linsolve{S,T,DD,dim}(A::QROperator,b::Fun{MatrixSpace{S,T,DD,dim}};kwds...) =
    linsolve(A,mat(b);kwds...)

function linsolve{S,T,DD,dim}(A::Operator,b::Fun{MatrixSpace{S,T,DD,dim}};kwds...)
    if dimension(domain(A)) > 1
        pdesolve(A,b;kwds...)
    elseif isambiguous(domainspace(A))
        A=choosespaces(A,Fun(b[:,1]))  # use only first column
        if isambiguous(domainspace(A))
            error("Cannot infer spaces")
        end
        linsolve(A,b;kwds...)
    else
        linsolve(qrfact(A),b;kwds...)
    end
end


text:
module RickTracy

using DataStructures, JLD

__init__() = begin
    global _num_trace_locations = 0
    global location_counts = DefaultDict{String, Int}(1) #number of times tracepoint at each location has been hit (but not necessarily logged)
    global watched_exprs = OrderedDict{String, Bool}()
    global happysnaps = Vector{TraceItem}()
    global autowatch = true
end

include("types.jl")
include("snapdat.jl")
include("snapall.jl")
include("traceviews.jl")
include("io.jl")
include("replay.jl")
include("clarence.jl")
include("utils.jl")

#=
1) http://julialang.org/blog/2016/10/StructuredQueries https://davidagold.github.io/StructuredQueries.jl/latest/man/guide.html
2) http://www.david-anthoff.com/Query.jl/stable/
3) https://github.com/FugroRoames/TypedTables.jl
4) https://github.com/JuliaComputing/IndexedTables.jl
=#

end

max_stars_repo_name: JobJob/RickTracy


text:
# using CytofResearch
# using Distributions
# using Flux, Flux.Tracker
include(joinpath(@__DIR__, "util.jl"))
include(joinpath(@__DIR__, "logprob_lam.jl"))

function sample_lam(state::CytofResearch.VB.StateMP, y::Vector{Matrix{Float64}},
                    c::CytofResearch.VB.Constants)
  _, s, ys, _ = CytofResearch.VB.rsample(state, y, c)
  lp_lam= logprob_lam(data(s), Tracker.data.(ys), c)
  lam = [[CytofResearch.MCMC.wsample_logprob(lp_lam[i][n, :])
          for n in 1:c.N[i]] for i in 1:c.I]

  # Relabel noisy class as 0
  for i in 1:c.I
    for n in 1:c.N[i]
      if lam[i][n] > c.K
        lam[i][n] = 0
      end
    end
  end

  return Vector{Int8}.(lam)
end

VV = Vector{Vector{I}} where I
VVV = Vector{VV{I}}  where I

function lam_f(lam_samps::VVV{T}, i::Integer, n::Integer, f::Function) where {T <: Integer}
  B = length(lam_samps)
  return f(lam_samps[b][i][n] for b in 1:B)
end

function lam_f(lam_samps::VVV{T}, f::Function) where {T <: Integer}
  I = length(lam_samps[1])
  N = length.(lam_samps[1])
  B = length(lam_samps)

  return [[lam_f(lam_samps, i, n, f) for n in 1:N[i]] for i in 1:I]
end

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
"""
    pyplotshow([fig]; close=true, [format])

Show the matplotlib/pyplot/seaborn/etc figure `fig`, or all open figures if not given, using Julia's display mechanism.

If `close` is true, the figure is also closed.

The `format` specifies the file format of the generated image.
By default this is `pyplot.rcParams["savefig.format"]`.
It can be one of `"png"`, `"jpg"`, `"jpeg"`, `"tif"`, `"tiff"`, `"svg"` or `"pdf"`.
"""
function pyplotshow(fig; close::Bool = true, format::String = "")
    @py ```
    import matplotlib.pyplot as plt, io
    fig = $fig
    if not isinstance(fig, plt.Figure):
        fig = plt.figure(fig)
    buf = io.BytesIO()
    format = $format
    if not format:
        format = plt.rcParams.get("savefig.format", "png")
    if format not in ["png", "jpg", "jpeg", "tif", "tiff", "svg", "pdf"]:
        raise ValueError("Unsupported format: {}".format(format))
    fig.savefig(buf, format=format)
    $(data::Vector{UInt8}) = buf.getvalue()
    if $close:
        plt.close(fig)
    $(format::String) = format
    ```
    if format == "png"
        display(MIME"image/png"(), data)
    elseif format in ("jpg", "jpeg")
        display(MIME"image/jpeg"(), data)
    elseif format in ("tif", "tiff")
        display(MIME"image/tiff"(), data)
    elseif format == "svg"
        display(MIME"image/svg+xml"(), String(data))
    elseif format == "pdf"
        display(MIME"application/pdf"(), data)
    else
        error()
    end
    nothing
end
function pyplotshow(; opts...)
    @py ```
    import sys
    plt = sys.modules.get("matplotlib.pyplot", None)
    $(fignums::Vector{Int}) = [] if plt is None else plt.get_fignums()
    ```
    for fig in fignums
        pyplotshow(fig; opts...)
    end
end
export pyplotshow


text:
function  poisson(theta, y, x)
	lambda = exp.(x*theta)
    n = size(y,1)
    c = zeros(n)
    for i = 1:n
        c[i] = logabsgamma(y[i])[1]
    end
    logdensity = y .* (x*theta) .- lambda .- c
end

max_stars_repo_name: Vkamalvanshi/Econometrics

max_issues_repo_name: Vkamalvanshi/Econometrics

max_forks_repo_name: Vkamalvanshi/Econometrics


text:
### MCSamplerState represents the Monte Carlo samplers' internal state ("local variables")

abstract type MCSamplerState{F<:VariateForm} end

abstract type MHSamplerState{F<:VariateForm} <: MCSamplerState{F} end

abstract type HMCSamplerState{F<:VariateForm} <: MCSamplerState{F} end

abstract type LMCSamplerState{F<:VariateForm} <: MCSamplerState{F} end

### Abstract Monte Carlo sampler typesystem

## Root Monte Carlo sampler

abstract type MCSampler end

## Family of samplers based on Metropolis-Hastings

abstract type MHSampler <: MCSampler end

## Family of Hamiltonian Monte Carlo samplers

abstract type HMCSampler <: MCSampler end

## Family of Langevin Monte Carlo samplers

abstract type LMCSampler <: MCSampler end

tuner_state(parameter::Parameter, sampler::MCSampler, tuner::VanillaMCTuner) = BasicMCTune(NaN, 0, 0, tuner.period)

tuner_state(parameter::Parameter, sampler::MHSampler, tuner::VanillaMCTuner) = BasicMCTune(1., 0, 0, tuner.period)

tuner_state(parameter::Parameter, sampler::HMCSampler, tuner::VanillaMCTuner) =
  BasicMCTune(sampler.leapstep, 0, 0, tuner.period)

tuner_state(parameter::Parameter, sampler::LMCSampler, tuner::VanillaMCTuner) =
  BasicMCTune(sampler.driftstep, 0, 0, tuner.period)

tuner_state(parameter::Parameter, sampler::MHSampler, tuner::AcceptanceRateMCTuner) = BasicMCTune(1., 0, 0, tuner.period)

tuner_state(parameter::Parameter, sampler::HMCSampler, tuner::AcceptanceRateMCTuner) =
  BasicMCTune(sampler.leapstep, 0, 0, tuner.period)

tuner_state(parameter::Parameter, sampler::LMCSampler, tuner::AcceptanceRateMCTuner) =
  BasicMCTune(sampler.driftstep, 0, 0, tuner.period)

function set_diagnosticindices!(sstate::MCSamplerState, diagnostickeys::Vector{Symbol}, samplerkeys::Vector{Symbol})
  if !isempty(diagnostickeys)
    dindices = map(k -> findfirst(diagnostickeys, k), samplerkeys)

    for (k, i) in zip(samplerkeys, dindices)
      if i != 0
        sstate.diagnosticindices[k] = i
      end
    end
  end
end

function reset!(
  pstate::ParameterState{Continuous, Univariate},
  x::Real,
  parameter::Parameter{Continuous, Univariate},
  sampler::HMCSampler
)
  pstate.value = x
  parameter.uptogradlogtarget!(pstate)
end

function reset!(
  pstate::ParameterState{Continuous, Multivariate},
  x::RealVector,
  parameter::Parameter{Continuous, Multivariate},
  sampler::HMCSampler
)
  pstate.value = copy(x)
  parameter.uptogradlogtarget!(pstate)
end

reset!(tune::BasicMCTune, sampler::MCSampler, tuner::VanillaMCTuner) =
  ((tune.accepted, tune.proposed, tune.totproposed, tune.rate) = (0, 0, tuner.period, NaN))

function reset!(tune::BasicMCTune, sampler::HMCSampler, tuner::AcceptanceRateMCTuner)
  tune.step = sampler.leapstep
  (tune.accepted, tune.proposed, tune.totproposed, tune.rate) = (0, 0, tuner.period, NaN)
end

function reset!(tune::BasicMCTune, sampler::LMCSampler, tuner::AcceptanceRateMCTuner)
  tune.step = sampler.driftstep
  (tune.accepted, tune.proposed, tune.totproposed, tune.rate) = (0, 0, tuner.period, NaN)
end

reset!(
  sstate::MCSamplerState{F},
  pstate::ParameterState{S, F},
  parameter::Parameter{S, F},
  sampler::MCSampler,
  tuner::MCTuner
) where {S<:ValueSupport, F<:VariateForm} =
  reset!(sstate.tune, sampler, tuner)

hamiltonian(logtarget::Real, momentum::Real) = logtarget-0.5*abs2(momentum)

hamiltonian(logtarget::Real, momentum::RealVector) = logtarget-0.5*dot(momentum, momentum)

function leapfrog!(
  pstate::ParameterState{Continuous, Univariate},
  pstate0::ParameterState{Continuous, Univariate},
  momentum0::Real,
  step::Real,
  gradlogtarget!::Function
)
  local momentum::Real

  momentum = momentum0+0.5*step*pstate0.gradlogtarget
  pstate.value = pstate0.value+step*momentum
  gradlogtarget!(pstate)
  momentum = momentum+0.5*step*pstate.gradlogtarget

  return momentum
end

function leapfrog!(
  pstate::ParameterState{Continuous, Multivariate},
  momentum::RealVector,
  pstate0::ParameterState{Continuous, Multivariate},
  momentum0::RealVector,
  step::Real,
  gradlogtarget!::Function
)
  momentum[:] = momentum0+0.5*step*pstate0.gradlogtarget
  pstate.value[:] = pstate0.value+step*momentum
  gradlogtarget!(pstate)
  momentum[:] = momentum+0.5*step*pstate.gradlogtarget
end

function initialize_step!(
  pstate::ParameterState{Continuous, Univariate},
  pstate0::ParameterState{Continuous, Univariate},
  momentum0::Real,
  step0::Real,
  gradlogtarget!::Function,
  ::Type{DualAveragingMCTuner}
)
  local momentum::Real
  local oldhamiltonian::Real
  local newhamiltonian::Real
  local ratio::Real
  local a::Real
  local step::Real = step0

  oldhamiltonian = hamiltonian(pstate0.logtarget, momentum0)

  momentum = leapfrog!(pstate, pstate0, momentum0, step, gradlogtarget!)
  newhamiltonian = hamiltonian(pstate.logtarget, momentum)

  ratio = newhamiltonian-oldhamiltonian
  a = 2*(exp(ratio) > 0.5)-1

  while exp(ratio)^a > 2^(-a)
    step = (2^a)*step
    momentum = leapfrog!(pstate, pstate, momentum, step, gradlogtarget!)
    newhamiltonian = hamiltonian(pstate.logtarget, momentum)

    ratio = newhamiltonian-oldhamiltonian
  end

  return step
end

function initialize_step!(
  pstate::ParameterState{Continuous, Multivariate},
  momentum::RealVector,
  pstate0::ParameterState{Continuous, Multivariate},
  momentum0::RealVector,
  step0::Real,
  gradlogtarget!::Function,
  ::Type{DualAveragingMCTuner}
)
  local oldhamiltonian::Real
  local newhamiltonian::Real
  local ratio::Real
  local a::Real
  local step::Real = step0

  oldhamiltonian = hamiltonian(pstate0.logtarget, momentum0)

  leapfrog!(pstate, momentum, pstate0, momentum0, step, gradlogtarget!)
  newhamiltonian = hamiltonian(pstate.logtarget, momentum)

  ratio = newhamiltonian-oldhamiltonian
  a = 2*(exp(ratio) > 0.5)-1

  while exp(ratio)^a > 2^(-a)
    step = (2^a)*step
    leapfrog!(pstate, momentum, pstate, moment, step, gradlogtarget!)
    newhamiltonian = hamiltonian(pstate.logtarget, momentum)

    ratio = newhamiltonian-oldhamiltonian
  end

  return step
end


text:
"""
    SARSOPSolver

Base solver type for SARSOP. Contains an options dictionary with the following entries:

* 'fast': use fast (but very picky) alternate parser for .pomdp files
* 'randomization': turn on randomization for the sampling algorithm
* 'precision': run ends when target precision is reached
* 'timeout':  [sec] If running time exceeds the specified value, pomdpsol writes out a policy and terminates
* 'memory': [MB] If memory usage exceeds the specified value, pomdpsol writes out a policy and terminates
* 'trial-improvement-factor': terminates when the gap between bounds reaches this value
* 'policy-interval':  the time interval between two consecutive write-out of policy files
* `verbose`: [true/false] whether to print output from the solver. Default: true
"""
@with_kw mutable struct SARSOPSolver <: Solver
    fast::Bool = false
    randomization::Bool = false
    precision::Float64 = 1e-3
    timeout::Union{Nothing, Float64} = nothing
    memory::Union{Nothing, Float64} = nothing
    trial_improvement_factor::Float64 = 0.5
    policy_interval::Union{Nothing, Float64} = nothing
    verbose::Bool = true
    policy_filename::String = "policy.out"
    pomdp_filename::String = "model.pomdpx"
end

function get_solver_options(solver::SARSOPSolver)
    options_list = AbstractString[]
    if solver.fast
        push!(options_list, "--fast")
    end
    if solver.randomization
        push!(options_list, "--randomization")
    end
    push!(options_list, "--precision")
    push!(options_list, string(solver.precision))
    if solver.timeout != nothing
        push!(options_list, "--timeout")
        push!(options_list, string(solver.timeout))
    end
    if solver.memory != nothing
        push!(options_list, "--memory")
        push!(options_list, string(solver.memory))
    end
    push!(options_list, "--trial-improvement-factor")
    push!(options_list, string(solver.trial_improvement_factor))
    if solver.policy_interval != nothing
        push!(options_list, "--policy-interval")
        push!(options_list, string(solver.policy_interval))
    end
    return options_list
end

"""
    solve(solver, pomdp, policy)

Runs pomdpsol using the options in 'solver' on 'pomdp',
and writes out a .policy xml file specified by 'policy'.
"""
function POMDPs.solve(solver::SARSOPSolver, pomdp::POMDP; kwargs...)
    if !isempty(kwargs)
        @warn("Keyword args for solve(::SARSOPSolver, ::POMDP) are no longer supported. Use the options in the SARSOPSolver")
    end
    pomdp_file = POMDPFile(pomdp, solver.pomdp_filename, verbose=solver.verbose)
    options = get_solver_options(solver)
    pomdpsol() do pomdpsol_path
        if solver.verbose
            run(`$pomdpsol_path $(pomdp_file.filename) --output $(solver.policy_filename) $options`)
        else
            success(`$pomdpsol_path $(pomdp_file.filename) --output $(solver.policy_filename) $options`)
        end
    end
    alphas = POMDPAlphas(solver.policy_filename)
    action_map = broadcast(x -> getindex(ordered_actions(pomdp), x), alphas.alpha_actions .+ 1)
    return AlphaVectorPolicy(pomdp, alphas.alpha_vectors, action_map)
end

POMDPs.solve(solver::SARSOPSolver, mdp::MDP) = mdp_error()

# utilities

"""
    load_policy(pomdp::POMDP, file_name::AbstractString)

Load a policy from an xml file output by SARSOP.
"""
function load_policy(pomdp::POMDP, file_name::AbstractString)
    alphas = nothing
    if isfile(file_name)
        alphas = POMDPAlphas(file_name)
    else
        error("Policy file ", file_name, " does not exist")
    end
    action_map = broadcast(x -> getindex(ordered_actions(pomdp), x), alphas.alpha_actions .+ 1)
    return AlphaVectorPolicy(pomdp, alphas.alpha_vectors, action_map)
end

"""
    to_pomdpx(pomdp::POMDPFile)

Convert a .pomdp file to a .pomdpx file.
"""
function to_pomdpx(pomdp::POMDPFile)
    @assert(splitext(pomdp.filename)[2] == ".pomdp")
    pomdpconvert() do pomdpconvert_path
        run(`$pomdpconvert_path $(pomdp.filename)`)
    end
end

mdp_error() = error("SARSOP is designed to solve POMDPs and is not set up to solve MDPs; consider using DiscreteValueIteration.jl to solve MDPs.")

POMDPLinter.@POMDP_require solve(solver::SARSOPSolver, pomdp::POMDP) begin
    P = typeof(pomdp)
    S = statetype(P)
    A = actiontype(P)
    O = obstype(P)
    @req discount(::P)
    @subreq ordered_states(pomdp)
    @subreq ordered_actions(pomdp)
    @subreq ordered_observations(pomdp)
    @req transition(::P,::S,::A)
    @req reward(::P,::S,::A)
    @req stateindex(::P,::S)
    @req obsindex(::P,::O)
    @req observations(::P)
    @req observation(::P,::A,::S)
    @req actionindex(::P, ::A)
    @req actions(::P, ::S)
    as = actions(pomdp)
    ss = states(pomdp)
    a = first(as)
    s = first(ss)
    dist = transition(pomdp, s, a)
    D = typeof(dist)
    @req support(::D)
    @req pdf(::D,::S)
end