text:
using CloudCovErr
import ImageFiltering
import Distributions
import StatsBase
using Random
using LinearAlgebra

export gen_mask_staticPSF! #
export gen_mask_staticPSF2! #
export prelim_infill! #
export add_sky_noise! #
export add_noise! #
export findmaxpsf #
export kstar_circle_mask #
export im_subrng #

"""
    kstar_circle_mask(Np;rlim=256) -> circmask

Generates a Bool mask for pixels beyond a given (squared) radius of the center of an image.

# Arguments:
- `Np`: size of image stamp

# Keywords:
- `rlim`: squared radius (in pixels^2) beyond which pixels should be masked (default 256)

# Outputs:
- `circmask`: static Bool mask used for assigning pixels beyond some radius of the stellar center as "ignored"
"""
function kstar_circle_mask(Np;rlim=256)
    halfNp = (Np-1) ÷ 2
    x = (-halfNp:halfNp)' .* ones(Int,Np)
    y = ones(Int,Np)' .* (-halfNp:halfNp)
    R = x.^2 .+ y.^2
    return R.>rlim
end

"""
    findmaxpsf(psfstamp1;thr=20) -> flim

Computes the flux a star must have so that the PSF-based masking using `thr`
would require a larger stamp area. Used for computational savings.

# Arguments:
- `psfstamp1`: a small image of a representative PSF

# Keywords:
- `thr`: threshold used to determine which pixels are bright enough to be "hidden"

# Outputs:
- `flim`: maximum flux that can be masked by `thr` without exceeding the PSF stamp footprint
"""
function findmaxpsf(psfstamp1;thr=20)
    thr/maximum([psfstamp1[:,1]...,psfstamp1[1,:]...,psfstamp1[:,end]...,psfstamp1[end,:]...])
end

"""
    im_subrng(jx,jy,cx,cy,sx,sy,px0,py0,stepx,stepy,padx,pady,tilex,tiley) -> xrng, yrng, star_ind

Computes the flux a star must have so that the PSF-based masking using `thr`
would require a larger stamp area. Used for computational savings.

# Arguments:
- `jx`: tile index along x
- `jy`: tile index along y
- `cx`: list of stellar x-coordinates
- `cy`: list of stellar y-coordinates
- `sx`: size of original image in x
- `sy`: size of original image in y
- `px0`: maximal padding in x to account for stars outside image
- `py0`: maximal padding in y to account for stars outside image
- `stepx`: tiling step size in x
- `stepy`: tiling step size in y
- `padx`: tile padding in x required to account for local stamp size, sample size, and pixels outside the image
- `pady`: tile padding in x required to account for local stamp size, sample size, and pixels outside the image
- `tilex`: total number of tile divisions along x
- `tiley`: total number of tile divisions along y

# Outputs:
- `xrng`: slicing range of the tile in x
- `yrng`: slicing range of the tile in y
- `star_ind`: Bool mask of all stars falling within the tile (subimage)
"""
function im_subrng(jx,jy,cx,cy,sx,sy,px0,py0,stepx,stepy,padx,pady,tilex,tiley)
    lowbx = (1 + (jx-1)*stepx)
    uppbx = (1 + jx*stepx-1)
    lowby = (1 + (jy-1)*stepy)
    uppby = (1 + jy*stepy-1)
    # pad right and top by 2 to fix ragged
    if uppbx > sx
        uppbx = sx
    end
    if uppby > sy
        uppby = sy
    end
    xrng = (lowbx-padx):(uppbx+padx)
    yrng = (lowby-pady):(uppby+pady)

    if jx == 1
        lowbx-=px0
    end
    if jx == tilex
        uppbx+=px0
    end
    if jy == 1
        lowby-=py0
    end
    if jy == tiley
        uppby+=py0
    end

    star_ind = findall((lowbx-0.5 .< cx .<= uppbx+0.5) .& (lowby-0.5 .< cy .<= uppby+0.5))

    return xrng, yrng, star_ind
end

"""
    gen_mask_staticPSF!(maskd, psfstamp, x_stars, y_stars, flux_stars, thr=20)

Generate a mask for an input image (which is usually an image of model residuals)
that excludes the cores of stars (which are often mismodeled). In this function,
we use a fixed PSF `psfstamp` for all sources, and adjust the masking fraction based on the
stellar flux and a threshold `thr`. A more general position dependent PSF model could be
used with a slight generalization of this function, but is likely overkill for the problem
of making a mask.

# Arguments:
- `maskd`: bool image to which mask will be added (bitwise or)
- `psfstamp`: simple 2D array of a single PSF to be used for the whole image
- `x_stars`: list of source x positions
- `y_stars`: list of source y positions
- `flux_stars`: list of source fluxes

# Keywords:
- `thr`: threshold used for flux-dependent masking (default 20)
"""
function gen_mask_staticPSF!(bmaskd, psfstamp, x_stars, y_stars, flux_stars; thr=20)
    (sx, sy) = size(bmaskd)
    (psx, psy) = size(psfstamp)
    Δx = (psx-1)÷2
    Δy = (psy-1)÷2
    Nstar = size(x_stars)[1]
    # assumes x/y_star is one indexed
    for i=1:Nstar
        fluxt=abs(flux_stars[i])
        x_star = round(Int64, x_stars[i])
        y_star = round(Int64, y_stars[i])
        @views mskt = (psfstamp .> thr/fluxt)[maximum([1,2+Δx-x_star]):minimum([1+sx-x_star+Δx,psx]),maximum([1,2+Δy-y_star]):minimum([1+sy-y_star+Δy,psy])]
        @views bmaskd[maximum([1,x_star-Δx]):minimum([x_star+Δx,sx]),maximum([1,y_star-Δy]):minimum([y_star+Δy,sy])] .|= mskt
    end
end

"""
    gen_mask_staticPSF2!(maskd, psfstamp, psfstamp1, x_stars, y_stars, flux_stars, thr=20)

Generate a mask for an input image (which is usually an image of model residuals)
that excludes the cores of stars (which are often mismodeled). In this function,
we use a small fixed PSF `psfstamp1` for all faint sources, and adjust the masking
fraction based on the stellar flux and a threshold `thr`. Only for source bright
enough to need a larger PSF stamp do we use `psfstamp`, which saves some computational
cost.

# Arguments:
- `maskd`: bool image to which mask will be added (bitwise or)
- `psfstamp`: simple 2D array of a single PSF to be used for bright stars in the whole image
- `psfstamp1`: simple 2D array of a single PSF to be used for faint stars in the whole image
- `x_stars`: list of source x positions
- `y_stars`: list of source y positions
- `flux_stars`: list of source fluxes

# Keywords:
- `thr`: threshold used for flux-dependent masking (default 20)
"""
function gen_mask_staticPSF2!(bmaskd, psfstamp, psfstamp1, x_stars, y_stars, flux_stars; thr=20)
    (sx, sy) = size(bmaskd)
    (psx, psy) = size(psfstamp)
    Δx = (psx-1)÷2
    Δy = (psy-1)÷2
    (psx1, psy1) = size(psfstamp1)
    Δx1 = (psx1-1)÷2
    Δy1 = (psy1-1)÷2
    Nstar = size(x_stars)[1]
    flim = findmaxpsf(psfstamp1;thr=thr)
    # assumes x/y_star is one indexed
    for i=1:Nstar
        fluxt = abs(flux_stars[i])
        x_star = round(Int64, x_stars[i])
        y_star = round(Int64, y_stars[i])
        if fluxt > flim
            pxrange = maximum([1,2+Δx-x_star]):minimum([1+sx-x_star+Δx,psx])
            pyrange = maximum([1,2+Δy-y_star]):minimum([1+sy-y_star+Δy,psy])
            bxrange = maximum([1,x_star-Δx]):minimum([x_star+Δx,sx])
            byrange = maximum([1,y_star-Δy]):minimum([y_star+Δy,sy])
            @views mskt = (psfstamp[pxrange,pyrange] .> thr/fluxt)
        else
            pxrange = maximum([1,2+Δx1-x_star]):minimum([1+sx-x_star+Δx1,psx1])
            pyrange = maximum([1,2+Δy1-y_star]):minimum([1+sy-y_star+Δy1,psy1])
            bxrange = maximum([1,x_star-Δx1]):minimum([x_star+Δx1,sx])
            byrange = maximum([1,y_star-Δy1]):minimum([y_star+Δy1,sy])
            @views mskt = (psfstamp1[pxrange,pyrange] .> thr/fluxt)
        end
        #println((size(mskt),size(bmaskd[bxrange,byrange])))
        @views bmaskd[bxrange,byrange] .|= mskt
        # FIX ME: worth triple checking these relative indexings (remove inbounds for testing when you do that!!)
    end
end

"""
    prelim_infill!(testim,bmaskim,bimage,bimageI,testim2,bmaskim2,goodpix,ccd;widx=19,widy=19,ftype::Int=32,widmult=1.4)

This intial infill replaces masked pixels with a guess based on a smoothed
boxcar. For large masked regions, the smoothing scale is increased. If this
iteration takes too long/requires too strong of masking, the masked pixels
are replaced with the median of the image.

We use 3 copies of the input image and mask image. The first
is a view (with reflective boundary condition padding) with the pixels to be infilled
replaced with zeros, the second is allocated to hold various smoothings of the image,
and the third holds the output image which contains our best infill guess. A final Bool
array of size corresponding to the image is used to keep track of pixels that have safe
infill values.

# Arguments:
- `testim`: input image which requires infilling
- `bmaskim`: input mask indicating which pixels require infilling
- `bimage`: preallocated array for smoothed version of input image
- `bimageI`: preallocated array for smoothed mask counting the samples for each estimate
- `testim2`: inplace modified ouptut array for infilled version of image input
- `bmaskim2`: inplace modified mask to keep track of which pixels still need infilling
- `goodpix`: preallocated array for Bool indexing pixels with good infill
- `ccd`: string name of FITS extension for verbose cmdline printing

# Keywords:
- `widx`: initial size of boxcar smoothing window in x (default 19)
- `widy`: initial size of boxcar smoothing window in y (default 19)
- `ftype::Int`: determine the Float precision, 32 is Float32, otherwise Float64
- `widmult`: multiplicative factor for increasing the smoothing scale at each iteration step
"""
function prelim_infill!(testim,bmaskim,bimage,bimageI,testim2,bmaskim2,goodpix,ccd;widx=19,widy=19,ftype::Int=32,widmult=1.4)
    if ftype == 32
        T = Float32
    else
        T = Float64
    end
    Δx = (widx-1)÷2
    Δy = (widy-1)÷2
    (sx, sy) = size(testim)

    widxMax = round(Int,(widmult^10)*(widx-1)/2)*2+1
    widyMax = round(Int,(widmult^10)*(widy-1)/2)*2+1
    ΔxMax = (widxMax-1)÷2
    ΔyMax = (widyMax-1)÷2

    #the masked entries in testim must be set to 0 so they drop out of the mean
    testim[bmaskim] .= 0;
    bmaskim2 .= copy(bmaskim)
    testim2 .= copy(testim)

    #hopefully replace with the reflected indexedx arrays
    in_image = ImageFiltering.padarray(testim,ImageFiltering.Pad(:reflect,(ΔxMax,ΔxMax)));
    in_mask = ImageFiltering.padarray(.!bmaskim,ImageFiltering.Pad(:reflect,(ΔyMax,ΔyMax)));

    #loop to try masking at larger and larger smoothing to infill large holes
    cnt=0
    while any(bmaskim2) .& (cnt .< 10)
        # this double deep view should be only 1 deep ideally... need the internal unwrap
        @views in_image1 = in_image[(1-Δx):(sx+Δx),(1-Δy):(sy+Δy)]
        @views in_mask1 = in_mask[(1-Δx):(sx+Δx),(1-Δy):(sy+Δy)]
        (sx1, sy1) = size(in_image1)
        tot = zeros(T,sx1)
        totI = zeros(Int,sx1)
        boxsmooth!(bimage,in_image1,tot,widx,widy)
        boxsmooth!(bimageI,in_mask1,totI,widx,widy)

        goodpix .= (bimageI .> 10)

        testim2[bmaskim2 .& goodpix] .= (bimage./bimageI)[bmaskim2 .& goodpix]
        bmaskim2[goodpix] .= false

        # update loop params
        cnt+=1
        widx*=1.4
        widy*=1.4
        widx = round(Int,(widx-1)/2)*2+1
        widy = round(Int,(widy-1)/2)*2+1
        Δx = (widx-1)÷2
        Δy = (widy-1)÷2
    end
    println("Infilling $ccd completed after $cnt rounds with final width (widx,widy) = ($widx,$widy)")
    flush(stdout)

    #catastrophic failure fallback
    if cnt == 10
        testim2[bmaskim2] .= StatsBase.median(testim)
        println("Infilling Failed Badly")
        flush(stdout)
    end
    return
end

"""
    add_sky_noise!(testim2,maskim,skyim,gain;seed=2021)

Adds noise to the infill that matches the Poisson noise of a rough estimate for
the sky background. A random seed to set a local random generator is provided for
reproducible unit testing.

# Arguments:
- `testim2`: input image which had infilling
- `maskim`: mask of pixels which were infilled
- `skyim`: rough estimate of sky background counts
- `gain`: gain of detector to convert from photon count noise to detector noise

# Keywords:
- `seed`: random seed for random generator
"""
function add_sky_noise!(testim2,maskim,skyim,gain;seed=2021)
    rng = MersenneTwister(seed)
    for i in eachindex(testim2)
        if maskim[i]
            intermed = -(rand(rng, Distributions.Poisson(convert(Float64,gain*(abs.(skyim[i]-testim2[i])))))/gain.-skyim[i])
            testim2[i] = intermed
        end
    end
end

"""
    add_noise!(testim2,gain;seed=2021)

Adds noise to an image that matches the Poisson noise of the pixel counts.
A random seed to set a local random generator is provided for reproducible unit testing.

# Arguments:
- `testim2`: input image which had infilling
- `gain`: gain of detector to convert from photon count noise to detector noise

# Keywords:
- `seed`: random seed for random generator
"""
function add_noise!(testim2,gain;seed=2021)
    rng = MersenneTwister(seed)
    for i in eachindex(testim2)
        @inbounds testim2[i] = rand(rng, Distributions.Poisson(convert(Float64,gain*testim2[i])))/gain
    end
end


text:
###here is Pody layer
const tanh_a = (0, 1, 0, 1/3, 0, 2/15)
Aaf(n) = tanh_a[1 : n]

pody(x) = max(-1, min(1, evalpoly(x, (0, 1, 0, 1/3))))

struct PD{F, M<:AbstractMatrix, B}
    W::M
    b::B
    σ::F
    function PD(W, b)
        return new{typeof(pody), typeof(W), typeof(b)}(W, b, pody)
    end
end

function (layer::PD)(x::AbstractVecOrMat)
    W, b, σ = layer.W, layer.b, layer.σ
    σ.(W * x .+ b)
end

function Base.show(io::IO, a::PD)
    o, i = size(a.W)
    print("PD(IO:$i => $o, pody:[0, 1, 0, 1/3] * x)")
end

Flux.trainable(a::PD) = a.a, a.W
Flux.@functor PD
###end of Pody layer

###here is Data Keeper
struct DataKeeper
    data::Array{Tuple, 1}
    function DataKeeper()
        new(Array{Tuple}(undef, 0))
    end
end

function Base.push!(DK::DataKeeper, x::Tuple)
    if length(DK.data) == 100
        popfirst!(DK.data)
    end
    push!(DK.data, x)
end
###end of Data Keeper

ext: jl


text:
using FiniteDiff, LinearAlgebra, SparseArrays, Test, LinearAlgebra,
      BlockBandedMatrices, ArrayInterfaceCore, BandedMatrices,
      ArrayInterfaceBlockBandedMatrices

fcalls = 0
function f(dx,x)
  global fcalls += 1
  for i in 2:length(x)-1
    dx[i] = x[i-1] - 2x[i] + x[i+1]
  end
  dx[1] = -2x[1] + x[2]
  dx[end] = x[end-1] - 2x[end]
  nothing
end
function oopf(x)
  global fcalls += 1
  vcat([-2x[1]+x[2]],[x[i-1]-2x[i]+x[i+1] for i in 2:length(x)-1],[x[end-1]-2x[end]])
end

function second_derivative_stencil(N)
  A = zeros(N,N)
  for i in 1:N, j in 1:N
      (j-i==-1 || j-i==1) && (A[i,j]=1)
      j-i==0 && (A[i,j]=-2)
  end
  A
end

J = FiniteDiff.finite_difference_jacobian(oopf,rand(30))
@test J ≈ second_derivative_stencil(30)
_J = sparse(J)
@test fcalls == 31

_J2 = similar(_J)
fcalls = 0
FiniteDiff.finite_difference_jacobian!(_J2,f,rand(30),colorvec=repeat(1:3,10))
@test fcalls == 4
@test _J2 ≈ _J

_J2 = similar(_J)
fcalls = 0
FiniteDiff.finite_difference_jacobian!(_J2,f,rand(30),Val{:central},colorvec=repeat(1:3,10))
@test fcalls == 6
@test _J2 ≈ _J

_J2 = similar(_J)
fcalls = 0
FiniteDiff.finite_difference_jacobian!(_J2,f,rand(30),Val{:complex},colorvec=repeat(1:3,10))
@test fcalls == 3
@test _J2 ≈ _J

_J2 = similar(_J)
_denseJ2 = collect(_J2)
fcalls = 0
FiniteDiff.finite_difference_jacobian!(_denseJ2,f,rand(30),colorvec=repeat(1:3,10),sparsity=_J2)
@test fcalls == 4
@test sparse(_denseJ2) ≈ _J

_J2 = similar(_J)
_denseJ2 = collect(_J2)
fcalls = 0
FiniteDiff.finite_difference_jacobian!(_denseJ2,f,rand(30),Val{:central},colorvec=repeat(1:3,10),sparsity=_J2)
@test fcalls == 6
@test sparse(_denseJ2) ≈ _J

_J2 = similar(_J)
_denseJ2 = collect(_J2)
fcalls = 0
FiniteDiff.finite_difference_jacobian!(_denseJ2,f,rand(30),Val{:complex},colorvec=repeat(1:3,10),sparsity=_J2)
@test fcalls == 3
@test sparse(_denseJ2) ≈ _J

_J2 = similar(Tridiagonal(_J))
fcalls = 0
FiniteDiff.finite_difference_jacobian!(_J2,f,rand(30),colorvec=repeat(1:3,10))
@test fcalls == 4
@test _J2 ≈ _J

_J2 = similar(_J2)
fcalls = 0
FiniteDiff.finite_difference_jacobian!(_J2,f,rand(30),Val{:central},colorvec=repeat(1:3,10))
@test fcalls == 6
@test _J2 ≈ _J

_J2 = similar(_J2)
fcalls = 0
FiniteDiff.finite_difference_jacobian!(_J2,f,rand(30),Val{:complex},colorvec=repeat(1:3,10))
@test fcalls == 3
@test _J2 ≈ _J

_Jb = BandedMatrices.BandedMatrix(similar(_J2),(1,1))
FiniteDiff.finite_difference_jacobian!(_Jb, f, rand(30), colorvec=repeat(1:3,10))
@test _Jb ≈ _J

_Jtri = Tridiagonal(similar(_J2))
FiniteDiff.finite_difference_jacobian!(_Jtri, f, rand(30), colorvec=repeat(1:3,10))
@test _Jtri ≈ _J

#https://github.com/JuliaDiff/FiniteDiff.jl/issues/67#issuecomment-516871956
function f(out, x)
	x = reshape(x, 100, 100)
	out = reshape(out, 100, 100)
	for i in 1:100
		for j in 1:100
			out[i, j] = x[i, j] + x[max(i -1, 1), j] + x[min(i+1, size(x, 1)), j] +  x[i, max(j-1, 1)]  + x[i, min(j+1, size(x, 2))]
		end
	end
	return vec(out)
end
x = rand(10000)
Jbbb = BandedBlockBandedMatrix(Ones(10000, 10000), fill(100, 100), fill(100, 100), (1, 1), (1, 1))
Jsparse = sparse(Jbbb)
colorsbbb = ArrayInterfaceCore.matrix_colors(Jbbb)
FiniteDiff.finite_difference_jacobian!(Jbbb, f, x, colorvec=colorsbbb)
FiniteDiff.finite_difference_jacobian!(Jsparse, f, x, colorvec=colorsbbb)
@test Jbbb ≈ Jsparse
Jbb = BlockBandedMatrix(similar(Jsparse),fill(100, 100), fill(100, 100),(1,1));
colorsbb = ArrayInterfaceCore.matrix_colors(Jbb)
FiniteDiff.finite_difference_jacobian!(Jbb, f, x, colorvec=colorsbb)
@test Jbb ≈ Jsparse

# Non-square Jacobian test.
# The Jacobian of f_nonsquare! has size (n, 2*n).
function f_nonsquare!(y, x)
    global fcalls += 1
    @assert length(x) == 2*length(y)
    n = length(x) ÷ 2
    x1 = @view x[1:n]
    x2 = @view x[n+1:end]

    @. y = (x1 .- 3).^2 .+ x1.*x2 .+ (x2 .+ 4).^2 .- 3
    return nothing
end

n = 4
x0 = vcat(ones(n).*(1:n) .+ 0.5, ones(n).*(1:n) .+ 1.5)
y0 = zeros(n)
rows = vcat([i for i in 1:n], [i for i in 1:n])
cols = vcat([i for i in 1:n], [i+n for i in 1:n])
sparsity = sparse(rows, cols, ones(length(rows)))
colorvec = vcat(fill(1, n), fill(2, n))

J_nonsquare1 = zeros(size(sparsity))
FiniteDiff.finite_difference_jacobian!(J_nonsquare1, f_nonsquare!, x0)

J_nonsquare2 = similar(sparsity)
for method in [Val(:forward), Val(:central), Val(:complex)]
    cache = FiniteDiff.JacobianCache(copy(x0), copy(y0), copy(y0), method; sparsity, colorvec)
    global fcalls = 0
    FiniteDiff.finite_difference_jacobian!(J_nonsquare2, f_nonsquare!, x0, cache)
    if method == Val(:central)
        @test fcalls == 2*maximum(colorvec)
    elseif method == Val(:complex)
        @test fcalls == maximum(colorvec)
    else
        @test fcalls == maximum(colorvec) + 1
    end
    @test isapprox(J_nonsquare2, J_nonsquare1; rtol=1e-6)
end


text:
# https://hevodata.com/blog/how-to-set-up-kafka-on-ubuntu-16-04/
# Using Python
# pip install kafka-python
using PyCall
using FileIO

function messageCallback(message)
    @show message
end

py"""
from kafka import KafkaConsumer
consumer = KafkaConsumer(bootstrap_servers='localhost:9092', auto_offset_reset='earliest', consumer_timeout_ms=1000)
"""
consumer = py"consumer"
consumer[:subscribe](["rawImageStream"])

py"""
from kafka import KafkaConsumer
def processMessages(consumer, juliaCallback):
    #consumer = KafkaConsumer(bootstrap_servers='localhost:9092',
    #                             auto_offset_reset='earliest',
    #                             consumer_timeout_ms=1000)
    consumer.subscribe(['rawImageStream'])
    for message in consumer:
        print(message)
        juliaCallback(message)
    #consumer.close()
"""
processMessages = py"processMessages"
processMessages(consumer, messageCallback)
consumer[:close]()


text:
# This file is a part of NeighborhoodApproximationIndex.jl

module NeighborhoodApproximationIndex
using SimilaritySearch

include("refs.jl")
include("knr.jl")
include("search.jl")
include("opt.jl")
end

max_stars_repo_name: sadit/NeighborhoodApproximationIndex

max_issues_repo_name: sadit/NeighborhoodApproximationIndex

max_forks_repo_name: sadit/NeighborhoodApproximationIndex


text:
let
    epc = Epoch(2018,1,1,12,0,0)

    a   = R_EARTH + 500e3
    oe  = [a, 0, 0, 0, 0, sqrt(GM_EARTH/a)]
    eci = sOSCtoCART(oe, use_degrees=true)

    r_rtn   = rRTNtoECI(eci)
    r_rtn_t = rECItoRTN(eci)

    r = r_rtn * r_rtn_t

    tol = 1e-8
    @test isapprox(r[1, 1], 1.0, atol=tol)
    @test isapprox(r[1, 2], 0.0, atol=tol)
    @test isapprox(r[1, 3], 0.0, atol=tol)

    @test isapprox(r[2, 1], 0.0, atol=tol)
    @test isapprox(r[2, 2], 1.0, atol=tol)
    @test isapprox(r[2, 3], 0.0, atol=tol)

    @test isapprox(r[3, 1], 0.0, atol=tol)
    @test isapprox(r[3, 2], 0.0, atol=tol)
    @test isapprox(r[3, 3], 1.0, atol=tol)
end

let
    epc = Epoch(2018,1,1,12,0,0)

    oe  = [R_EARTH + 500e3, 0, 0, 0, 0, 0]
    eci = sOSCtoCART(oe, use_degrees=true)

    xt = deepcopy(eci) + [100, 0, 0, 0, 0, 0]

    x_rtn = sECItoRTN(eci, xt)

    tol = 1e-8
    @test isapprox(x_rtn[1], 100.0, atol=tol)
    @test isapprox(x_rtn[2], 0.0, atol=tol)
    @test isapprox(x_rtn[3], 0.0, atol=tol)
    @test isapprox(x_rtn[4], 0.0, atol=tol)
    @test isapprox(x_rtn[5], 0.0, atol=0.5)
    @test isapprox(x_rtn[6], 0.0, atol=tol)

    xt2 = sRTNtoECI(eci, x_rtn)

    @test isapprox(xt[1], xt2[1], atol=tol)
    @test isapprox(xt[2], xt2[2], atol=tol)
    @test isapprox(xt[3], xt2[3], atol=tol)
    @test isapprox(xt[4], xt2[4], atol=tol)
    @test isapprox(xt[5], xt2[5], atol=tol)
    @test isapprox(xt[6], xt2[6], atol=tol)
end


text:
module Comms

using LCMCore

mutable struct CommsT <: LCMType
    utime::Int64
    format::String
    format_version_major::Int32
    format_version_minor::Int32
    num_bytes::Int32
    data::Vector{UInt8}
end

@lcmtypesetup(CommsT,
    data => (num_bytes, )
)

function CommsT(utime::Integer, format::String, format_version_major::Integer, format_version_minor::Integer, data::Vector{UInt8})
    CommsT(utime, format, format_version_major, format_version_minor, length(data), data)
end

end


text:
using HDF5, FileIO, JLD2, Test, Random, DSGE, DelimitedFiles
path = dirname(@__FILE__)
writing_output = false

m = AnSchorfheide()

save = normpath(joinpath(dirname(@__FILE__),"save"))
m <= Setting(:saveroot, save)

data = h5read("$path/../../reference/smc.h5", "data")

m <= Setting(:n_particles, 400)
m <= Setting(:n_Φ, 100)
m <= Setting(:λ, 2.0)
m <= Setting(:n_smc_blocks, 1)
m <= Setting(:use_parallel_workers, true)
m <= Setting(:step_size_smc, 0.5)
m <= Setting(:n_mh_steps_smc, 1)
m <= Setting(:resampler_smc, :polyalgo)
m <= Setting(:target_accept, 0.25)

m <= Setting(:mixture_proportion, .9)
m <= Setting(:adaptive_tempering_target_smc, false)
m <= Setting(:resampling_threshold, .5)
m <= Setting(:smc_iteration, 0)
m <= Setting(:use_chand_recursion, true)

@everywhere Random.seed!(42)

println("Estimating AnSchorfheide Model... (approx. 2 minutes)")
DSGE.smc2(m, data, verbose = :none) # us.txt gives equiv to periods 95:174 in our current dataset
println("Estimation done!")

test_file = load(rawpath(m, "estimate", "smc_cloud.jld2"))
test_cloud  = test_file["cloud"]
test_w      = test_file["w"]
test_W      = test_file["W"]
# test_z = test_file["z"]

if writing_output
    jldopen("$path/../../reference/smc_cloud_fix=true.jld2", true, true, true, IOStream) do file
        write(file, "cloud", test_cloud)
        write(file, "w", test_w)
        write(file, "W", test_W)
        write(file, "z", test_z)
    end
end

saved_file  = load("$path/../../reference/smc_cloud_fix=true.jld2")
saved_cloud = old_to_new_cloud(DSGE.Cloud(saved_file["cloud"]))
saved_w     = saved_file["w"]
saved_W     = saved_file["W"]

####################################################################
cloud_fields = fieldnames(typeof(test_cloud))
@testset "ParticleCloud Fields: AnSchorf" begin
    @test @test_matrix_approx_eq SMC.get_vals(test_cloud) SMC.get_vals(saved_cloud)
    @test @test_matrix_approx_eq SMC.get_loglh(test_cloud) SMC.get_loglh(saved_cloud)
    @test length(test_cloud.particles) == length(saved_cloud.particles)
    @test test_cloud.tempering_schedule == saved_cloud.tempering_schedule
    @test test_cloud.ESS ≈ saved_cloud.ESS
    @test test_cloud.stage_index == saved_cloud.stage_index
    @test test_cloud.n_Φ == saved_cloud.n_Φ
    @test test_cloud.resamples == saved_cloud.resamples
    @test test_cloud.c == saved_cloud.c
    @test test_cloud.accept == saved_cloud.accept
end

test_particle  = test_cloud.particles[1,:]
saved_particle = saved_cloud.particles[1,:]
N = length(test_particle)
@testset "Individual Particle Fields Post-SMC: AnSchorf" begin
    @test test_particle[1:SMC.ind_para_end(N)] ≈ saved_particle[1:SMC.ind_para_end(N)]
    @test test_particle[SMC.ind_loglh(N)]      ≈ saved_particle[SMC.ind_loglh(N)]
    @test test_particle[SMC.ind_logpost(N)]    ≈ saved_particle[SMC.ind_logpost(N)]
    @test test_particle[SMC.ind_logprior(N)]   ≈ saved_particle[SMC.ind_logprior(N)]
    @test test_particle[SMC.ind_old_loglh(N)] == saved_particle[SMC.ind_old_loglh(N)]
    @test test_particle[SMC.ind_accept(N)]    == saved_particle[SMC.ind_accept(N)]
    @test test_particle[SMC.ind_weight(N)]     ≈ saved_particle[SMC.ind_weight(N)]
end

@testset "Weight Matrices: AnSchorf" begin
    @test @test_matrix_approx_eq test_w saved_w
    @test @test_matrix_approx_eq test_W saved_W
end

######################################################################
# Test: SmetsWoutersOrig Model, no blocks
######################################################################
m = SmetsWoutersOrig()

save = normpath(joinpath(dirname(@__FILE__),"save"))
m <= Setting(:saveroot, save)

data = readdlm("reference/YY.txt")
data = Matrix{Float64}(data')

m <= Setting(:n_particles, 12000)
m <= Setting(:n_Φ, 500)
m <= Setting(:λ, 2.0)
m <= Setting(:n_smc_blocks, 1)
m <= Setting(:use_parallel_workers, true)
m <= Setting(:step_size_smc, 0.5)
m <= Setting(:n_mh_steps_smc, 1)
m <= Setting(:resampler_smc, :polyalgo)
m <= Setting(:target_accept, 0.25)

m <= Setting(:mixture_proportion, .9)
m <= Setting(:adaptive_tempering_target_smc, false)
m <= Setting(:resampling_threshold, .5)
m <= Setting(:smc_iteration, 0)
m <= Setting(:use_chand_recursion, true)

@everywhere Random.seed!(42)
smc(m, data, verbose = :none, run_test = true)

test_file = load(rawpath(m, "estimate", "smc_cloud.jld2"))
test_cloud  = test_file["cloud"]
test_w      = test_file["w"]
test_W      = test_file["W"]
test_z      = test_file["z"]

if writing_output
    jldopen("reference/smc_sw_cloud_fix=true_blocks=1.jld2", true, true, true, IOStream) do file
        write(file, "cloud", test_cloud)
        write(file, "w", test_w)
        write(file, "W", test_W)
        write(file, "z", test_z)
    end
end

saved_file   = load("reference/smc_sw_cloud_fix=true_blocks=1.jld2")
saved_cloud  = saved_file["cloud"]
saved_w      = saved_file["w"]
saved_W      = saved_file["W"]

#####################################################################
cloud_fields = fieldnames(typeof(test_cloud))
@testset "ParticleCloud Fields: SW, No Blocking" begin
    @test @test_matrix_approx_eq DSGE.get_vals(test_cloud) DSGE.get_vals(saved_cloud)
    @test @test_matrix_approx_eq DSGE.get_loglh(test_cloud) DSGE.get_loglh(saved_cloud)
    @test length(test_cloud.particles) == length(saved_cloud.particles)
    @test test_cloud.tempering_schedule == saved_cloud.tempering_schedule
    @test test_cloud.ESS == saved_cloud.ESS
    @test test_cloud.stage_index == saved_cloud.stage_index
    @test test_cloud.n_Φ == saved_cloud.n_Φ
    @test test_cloud.resamples == saved_cloud.resamples
    @test test_cloud.c == saved_cloud.c
    @test test_cloud.accept == saved_cloud.accept
end

test_particle = test_cloud.particles[1]
saved_particle = saved_cloud.particles[1]
particle_fields = fieldnames(typeof(test_particle))
@testset "Individual Particle Fields Post-SMC: SW, No Blocking" begin
    @test test_particle.weight == saved_particle.weight
    @test test_particle.keys == saved_particle.keys
    @test test_particle.value == saved_particle.value
    @test test_particle.loglh == saved_particle.loglh
    @test test_particle.logpost == saved_particle.logpost
    @test test_particle.old_loglh == saved_particle.old_loglh
    @test test_particle.accept == saved_particle.accept
end

@testset "Weight Matrices: SW, No Blocking" begin
    @test @test_matrix_approx_eq test_w saved_w
    @test @test_matrix_approx_eq test_W saved_W
end

######################################################################
# Test: SmetsWoutersOrig Model, blocks = 3
######################################################################
m = SmetsWoutersOrig()

save = normpath(joinpath(dirname(@__FILE__),"save"))
m <= Setting(:saveroot, save)

data = readdlm("reference/YY.txt")
data = Matrix{Float64}(data')

m <= Setting(:n_particles, 12000)
m <= Setting(:n_Φ, 500)
m <= Setting(:λ, 2.0)
m <= Setting(:n_smc_blocks, 3)
m <= Setting(:use_parallel_workers, true)
m <= Setting(:step_size_smc, 0.5)
m <= Setting(:n_mh_steps_smc, 1)
m <= Setting(:resampler_smc, :polyalgo)
m <= Setting(:target_accept, 0.25)

m <= Setting(:mixture_proportion, .9)
m <= Setting(:adaptive_tempering_target_smc, false)
m <= Setting(:resampling_threshold, .5)
m <= Setting(:smc_iteration, 0)
m <= Setting(:use_chand_recursion, true)

@everywhere Random.seed!(42)
smc(m, data, verbose = :none, run_test = true)

test_file = load(rawpath(m, "estimate", "smc_cloud.jld2"))
test_cloud  = test_file["cloud"]
test_w      = test_file["w"]
test_W      = test_file["W"]
test_z = test_file["z"]

if writing_output
    jldopen("reference/smc_sw_cloud_fix=true_blocks=3.jld2", true, true, true, IOStream) do file
        write(file, "cloud", test_cloud)
        write(file, "w", test_w)
        write(file, "W", test_W)
        write(file, "z", test_z)
    end
end

saved_file   = load("reference/smc_sw_cloud_fix=true_blocks=3.jld2")
saved_cloud  = saved_file["cloud"]
saved_w      = saved_file["w"]
saved_W      = saved_file["W"]
####################################################################
cloud_fields = fieldnames(typeof(test_cloud))
@testset "ParticleCloud Fields: SW, Blocks=3" begin
    @test @test_matrix_approx_eq DSGE.get_vals(test_cloud) DSGE.get_vals(saved_cloud)
    @test @test_matrix_approx_eq DSGE.get_loglh(test_cloud) DSGE.get_loglh(saved_cloud)
    @test length(test_cloud.particles) == length(saved_cloud.particles)
    @test test_cloud.tempering_schedule == saved_cloud.tempering_schedule
    @test test_cloud.ESS == saved_cloud.ESS
    @test test_cloud.stage_index == saved_cloud.stage_index
    @test test_cloud.n_Φ == saved_cloud.n_Φ
    @test test_cloud.resamples == saved_cloud.resamples
    @test test_cloud.c == saved_cloud.c
    @test test_cloud.accept == saved_cloud.accept
end

test_particle = test_cloud.particles[1]
saved_particle = saved_cloud.particles[1]
particle_fields = fieldnames(typeof(test_particle))
@testset "Individual Particle Fields Post-SMC: SW, Blocks=3" begin
    @test test_particle.weight == saved_particle.weight
    @test test_particle.keys == saved_particle.keys
    @test test_particle.value == saved_particle.value
    @test test_particle.loglh == saved_particle.loglh
    @test test_particle.logpost == saved_particle.logpost
    @test test_particle.old_loglh == saved_particle.old_loglh
    @test test_particle.accept == saved_particle.accept
end

@testset "Weight Matrices: SW, Blocks=3" begin
    @test @test_matrix_approx_eq test_w saved_w
    @test @test_matrix_approx_eq test_W saved_W
end

rmprocs(procs())


text:
using Test
using Jin

@Jin.register plus_one(x)=x+1
@Jin.register plus_two(x)=x+2

@Jin.configurable function f(;x=5)
    return x+1
end

@Jin.configurable function g(value;name="default")
    if name=="default"
        return value
    elseif name=="from_config"
        return -value
    else
        throw("Unrecognized name")
    end
end

@Jin.configurable function outer_fn(x;kernel=x->x)
    return -kernel(x)
end

@test f()==6
@test f(x=10)==11
@test g(5)==5
@test outer_fn(3)== -3
@test outer_fn(3, kernel=plus_one) == -4
@test outer_fn(3, kernel=plus_two) == -5

config_filename = joinpath(@__DIR__, "test_config.jin")
Jin.load_config(config_filename)
println(Jin.jind)

@test f()==16
@test f(x=10)==11
@test g(5)==-5
@test outer_fn(3) == -5

Jin.reset()
@test f()==6
@test g(5)==5


text:
# IRI_ass_kfold_include.jl

include("../tcm_code/tcm_transactionsToCounts.jl")

include("../tcm_code/tcm_estimateMNL.jl")
include("../tcm_code/tcm_predictMNL.jl")

include("../tcm_code/tcm_estimateMMNL.jl")
include("../tcm_code/tcm_estimateMMNL_kcv.jl")
include("../tcm_code/tcm_predictMMNL.jl")

include("../tcm_code/tcm_estimateHALOMNL.jl")
include("../tcm_code/tcm_predictHALOMNL.jl")

include("../tcm_code/tcm_estimateRanking.jl")
include("../tcm_code/tcm_predictRanking.jl")

include("../tcm_code/tcm_estimateForest.jl")
include("../tcm_code/tcm_estimateForest_kcv.jl")

include("../tcm_code/tcm_predictForest.jl")

include("../tcm_code/tcm_evaluateKL.jl")

include("IRI_ass_kfold_helperfn.jl")

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
module InteractiveDynamics

using Makie, Observables

include("colors.jl")
include("utils.jl")
using Requires

function __init__()
    @require DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" begin
        # include("numericdata/plot_dataset.jl")
        # include("numericdata/trajectory_highlighter.jl")
        include("chaos/orbitdiagram.jl")
        include("chaos/poincareclick.jl")
        include("chaos/trajanim.jl")
        include("chaos/cobweb.jl")
        include("chaos/brainscan.jl")
    end
    @require DynamicalBilliards = "4986ee89-4ee5-5cef-b6b8-e49ba721d7a5" begin
        include("billiards/defs_plotting.jl")
        include("billiards/defs_animating.jl")
        include("billiards/interactive_billiard.jl")
    end
    @require Agents = "46ada45e-f475-11e8-01d0-f70cc89e6671" begin
        include("agents/stepping.jl")
        include("agents/plots_videos.jl")
        include("agents/interactive_parameters.jl")
    end
end

end


text:
#Load some packages and add workers
# using Pkg
# Pkg.add("Images")
# Pkg.add("ImageIO")
# Pkg.add("QuartzImageIO")
# Pkg.add("ImageMagick")

using Images
using Distributed

using LinearAlgebra
using Statistics

#Load the package on all workers
# addprocs(4)
# @everywhere

# niw mixturemodels
using DPMMSubClusters

function segment_rgb(img_path,out_path,rgb_prior_multiplier,xy_prior_multiplier,nu)
    img = load(img_path)
    #Change to channel view
    img_channels = channelview(img)
    z,x,y = size(img_channels)
    #Create input
    input_arr = zeros(5,x*y)
    for i=1:x
        for j=1:y
            input_arr[:,(i-1)*y+j]=[img_channels[1,i,j],img_channels[2,i,j],img_channels[3,i,j],i,j]
        end
    end

    #Create HyperParams
    #The rgb,xy multiplier allows us to both play with the weight of the xy/rgb
    data_cov = cov(input_arr')
    data_cov[4:5,1:3] .= 0
    data_cov[1:3,4:5] .= 0

    data_cov[1:3,1:3] .*= rgb_prior_multiplier
    data_cov[4:5,4:5] .*= xy_prior_multiplier

    data_mean = mean(input_arr,dims = 2)[:]

    hyper_params = DPMMSubClusters.niw_hyperparams(0.5,
            data_mean,
            nu, #nu
            data_cov)
    #Run the model
    labels,clusters,weights = DPMMSubClusters.fit(input_arr,hyper_params,50000.0,iters = 200, verbose = true)

    #Get the cluster color means
    color_means = [x.μ[1:3] for x in clusters]

    segemnated_image = zeros(3,x,y)
    for i=1:x
        for j=1:y
            segemnated_image[:,i,j] = color_means[labels[(i-1)*y+j]]
        end
    end
    segemnated_image = colorview(RGB,segemnated_image)
    save(out_path, segemnated_image)
end

function segment_rgbd(img_path,depth_path,out_path,rgb_prior_multiplier,xyz_prior_multiplier,nu)
    img = load(img_path)
    depth = load(depth_path)
    #Change to channel view
    img_channels = channelview(img)
    rgb,x,y = size(img_channels)
    #Create input
    input_arr = zeros(6,x*y)
    for i=1:x
        for j=1:y
            input_arr[:,(i-1)*y+j]=[img_channels[1,i,j],img_channels[2,i,j],img_channels[3,i,j],i,j,depth[i,j]]
        end
    end
    #Create HyperParams
    #The rgb,xy multiplier allows us to both play with the weight of the xy/rgb
    # rgb_prior_multiplier = 30
    # xyz_prior_multiplier = 1
    data_cov = cov(input_arr')
    data_cov[4:6,1:3] .= 0
    data_cov[1:3,4:6] .= 0
    data_cov[1:3,1:3] .*= rgb_prior_multiplier
    data_cov[4:6,4:6] .*= xyz_prior_multiplier
    data_mean = mean(input_arr,dims = 2)[:]
    hyper_params = DPMMSubClusters.niw_hyperparams(0.5,
            data_mean,
            nu, #nu
            data_cov)
    #Run the model
    labels,clusters,weights = DPMMSubClusters.fit(input_arr,hyper_params,50000.0,iters = 200, verbose = true)
    #Get the cluster color means
    color_means = [x.μ[1:3] for x in clusters]
    segemnated_image = zeros(3,x,y)
    for i=1:x
        for j=1:y
            segemnated_image[:,i,j] = color_means[labels[(i-1)*y+j]]
        end
    end
    segemnated_image = colorview(RGB,segemnated_image)
    save(out_path, segemnated_image)
end

img_path = "desk_1_1.png"
depth_path = "desk_1_1_depth.png"
rgb_prior_multiplier = 10
xy_prior_multiplier =1
nu = 128
# segment_rgb(img_path,"rgb_result.png",rgb_prior_multiplier,xy_prior_multiplier, nu)
segment_rgbd(img_path,depth_path,"rgbd_result.png",rgb_prior_multiplier,xy_prior_multiplier, nu)


text:
@testset "Transient Responses" begin

    fname = joinpath(dirname(@__FILE__), "..", "..", "data", "test_Hz19.5-testing.bdf")
    fname2 = joinpath(dirname(@__FILE__), "..", "..", "data", "tmp", "test_Hz19.5-copy.bdf")
    fname_out = joinpath(dirname(@__FILE__), "..", "..", "data", "tmp", "testwrite.bdf")
    cp(fname, fname2, force = true)  # So doesnt use .mat file

    @testset "Read file" begin

        s = read_TR(fname)

        @test isa(s, NeuroimagingMeasurement)
        @test isa(s, EEG)
        @test isa(s, TR)

        s = read_TR(fname2)

        @test samplingrate(s) == 8192.0u"Hz"
        @test samplingrate(Int, s) == 8192
        @test isa(samplingrate(Float64, s), AbstractFloat) == true
        @test isa(samplingrate(Int, s), Int) == true

        @test isapprox(maximum(s.data[:, 2]), 54.5939; atol = 0.1)
        @test isapprox(minimum(s.data[:, 4]), -175.2503; atol = 0.1)

        s = read_TR(fname, valid_triggers = [8])
        s = read_TR(fname, min_epoch_length = 8388)
        s = read_TR(fname, max_epoch_length = 8389)

        s = read_TR(fname)
        original_length = length(s.triggers["Index"])
        s = read_TR(fname, remove_first = 10)
        @test original_length - 10 == length(s.triggers["Index"])

        s = read_TR(fname, max_epochs = 12)
        @test length(s.triggers["Index"]) == 12

        s = read_TR(fname)
        s.triggers = extra_triggers(s.triggers, 1, 7, 0.7, samplingrate(Float64, s))
        @test length(s.triggers["Index"]) == 56
    end

    @testset "Show" begin

        s = read_TR(fname)
        show(s)

    end

    @testset "Plot" begin
        s = read_TR(fname)
        Neuroimaging.plot(s)
        Neuroimaging.plot(s, "Cz")
        Neuroimaging.plot(s, ["Cz"])
        Neuroimaging.plot(s, ["Cz", "10Hz_SWN_70dB_R"])
    end

end


text:
# https://github.com/gnuradio/SigMF/blob/7996ac17d06f690ab12d02d95968d3974f712b49/extensions/signal.sigmf-ext.md

Base.@kwdef struct SignalDetailExtension
    type::Union{String, Nothing} = nothing
    class::Union{String, Nothing} = nothing
    standard::Union{String, Nothing} = nothing
    carrier_variant::Union{String, Nothing} = nothing
    symbol_variant::Union{String, Nothing} = nothing
    order::Union{UInt, Nothing} = nothing
    duplexing::Union{String, Nothing} = nothing
    multiplexing::Union{String, Nothing} = nothing
    multiple_access::Union{String, Nothing} = nothing
    spreading::Union{String, Nothing} = nothing
    channel_bw::Union{Float64, Nothing} = nothing
    channel::Union{UInt, Nothing} = nothing
    class_variant::Union{String, Nothing} = nothing
end

Base.@kwdef struct SignalEmitterExtension
    seid::Union{UInt, Nothing} = nothing
    manufacturer::Union{String, Nothing} = nothing
    power_tx::Union{Float64, Nothing} = nothing
    power_eirp::Union{Float64, Nothing} = nothing
    geolocation::Union{String, Nothing} = nothing
end

StructTypes.StructType(::Type{SignalDetailExtension}) = StructTypes.Struct()
StructTypes.StructType(::Type{SignalEmitterExtension}) = StructTypes.Struct()


text:
# Define convolutional layer:
struct Conv; w; stride; padding; f;end #Define convolutional layer
(c::Conv)(x)= c.f.(conv4(c.w,x; stride = c.stride, padding = c.padding))
Conv(w1::Int,w2::Int,cx::Int,cy::Int,st,pd,f) = Conv(param(w1,w2,cx,cy),st,pd,f)

# Define dense layer:
struct Dense; w; b; f; end
(d::Dense)(x) = d.f.(d.w * mat(x) .+ d.b)
Dense(i::Int,o::Int,f=relu) = Dense(param(o,i), param0(o), f)

# Define pooling layer:
struct Pool; size; stride; pad; end # Define pool layer
(p::Pool)(x) = pool(x; window = p.size, stride = p.stride, padding=p.pad) #pool function

# Define a chain of layers and a loss function:
struct Chain
    layers
    Chain(layers...) = new(layers)
end
(c::Chain)(x) = (for l in c.layers; x = l(x); end; x) #chain function
(c::Chain)(x,y) = nll(c(x),y) # loss function


text:
using HSA

if !isdefined(:get_testagent)

function get_testagent(agents = HSA.all_agents())
	spectre_idx = findfirst(a -> begin
		name = HSA.agent_info_name(a)
		# Creating queues for the CPU fails for some reason (Generic Error)
		return startswith(name, "Spectre")
	end, agents)
	return agents[spectre_idx]
end

macro with_agents(args...)
	quote
    if size(HSA.all_agents(),1) > 0
		$(esc(args...))
	else
        # no agents
		$(esc(:(@pending)),:(x => y))
    end
    end
end

macro with_codegen(context)
    if context.head != :call
        error("use with_codegen only on facts(...) or context(...) calls")
    end

    title = context.args[3]
	quote
    if HSA.has_hsa_codegen()
		$context
	else
        # no codegen
        title = $title
        warn("codegen is unavailable: skipping '$title'")
    end
    end
end

function get_softqueue(agent, s = HSA.Signal(value = 0), r_out = Ref{Nullable{HSA.Region}}())
    regs = HSA.regions(agent)

	gl_idx = findfirst(x -> begin
        seg = HSA.region_info_segment(x)
		if seg == HSA.RegionSegmentGlobal
			flags = HSA.region_info_global_flags(x)

			if (flags & HSA.RegionGlobalFlagKernarg) != 0
				return true
			end
		end
		return false
	end, regs)

	gl_reg = regs[gl_idx]

	r_out.x = gl_reg

	return HSA.Queue(gl_reg, 0x04, s)
end
end


text:
include("../src/datastats.jl")

# normal (variables on columns)
X = rand(40, 10)
ds = DataStats(10)
updatestats(ds, X[1:20, :])
updatestats(ds, X[21:end, :])

@assert all((cov(X) - cov(ds)) .< 0.0001)

# transposed (variables on rows)
X = rand(40, 10)
ds = DataStats(10, 2)
updatestats(ds, X')

@assert all((cov(X) - cov(ds)) .< 0.0001)

println("All OK")


text:
using Kinetic, ProgressMeter, Plots, LinearAlgebra, JLD2

###
# initialize kinetic solver
###

cd(@__DIR__)
D = Dict{Symbol,Any}()
begin
    D[:matter] = "gas"
    D[:case] = "sod"
    D[:space] = "1d2f1v"
    D[:flux] = "kfvs"
    D[:collision] = "bgk"
    D[:nSpecies] = 1
    D[:interpOrder] = 2
    D[:limiter] = "vanleer"
    D[:boundary] = "fix"
    D[:cfl] = 0.8
    D[:maxTime] = 0.2

    D[:x0] = 0.0
    D[:x1] = 1.0
    D[:nx] = 100
    D[:pMeshType] = "uniform"
    D[:nxg] = 1

    D[:umin] = -5.0
    D[:umax] = 5.0
    D[:nu] = 100
    D[:vMeshType] = "rectangle"
    D[:nug] = 0

    D[:knudsen] = 1e-4
    D[:mach] = 0.0
    D[:prandtl] = 1.0
    D[:inK] = 2.0
    D[:omega] = 0.81
    D[:alphaRef] = 1.0
    D[:omegaRef] = 0.5
end

ks = SolverSet(D)
ctr, face = init_fvm(ks, ks.ps, :dynamic_array; structarray = true)
for i in eachindex(face)
    face[i].fw .= 0.0
    #face[i].ff .= 0.0
end

function flux_tmp!(
    fw::X,
    wL::Y,
    wR::Y,
    γ::Real,
    inK::Real,
    μᵣ::Real,
    ω::Real,
    dt::Real,
    dxL::Real,
    dxR::Real,
    swL::Y,
    swR::Y,
) where {X<:AbstractArray{<:AbstractFloat,1},Y<:AbstractArray{<:AbstractFloat,1}}

    primL = conserve_prim(wL .+ swL .* dxL, γ)
    primR = conserve_prim(wR .- swR .* dxR, γ)

    Mu1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Mu2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)

    w =
        primL[1] .* moments_conserve(MuL1, Mxi1, 0, 0) .+
        primR[1] .* moments_conserve(MuR2, Mxi2, 0, 0)
    prim = conserve_prim(w, γ)
    tau =
        vhs_collision_time(prim, μᵣ, ω)

    faL = pdf_slope(primL, swL, inK)
    sw = -primL[1] .* moments_conserve_slope(faL, Mu1, Mxi1, 1)
    faTL = pdf_slope(primL, sw, inK)

    faR = pdf_slope(primR, swR, inK)
    sw = -primR[1] .* moments_conserve_slope(faR, Mu2, Mxi2, 1)
    faTR = pdf_slope(primR, sw, inK)

    Mu, Mxi, MuL, MuR = gauss_moments(prim, inK)
    sw0L = (w .- wL) ./ dxL
    sw0R = (wR .- w) ./ dxR
    gaL = pdf_slope(prim, sw0L, inK)
    gaR = pdf_slope(prim, sw0R, inK)
    sw =
        -prim[1] .* (
            moments_conserve_slope(gaL, MuL, Mxi, 1) .+
            moments_conserve_slope(gaR, MuR, Mxi, 1)
        )
    # ga = pdf_slope(prim, sw, inK)
    # sw = -prim[1] .* moments_conserve_slope(ga, Mu, Mxi, 1)
    gaT = pdf_slope(prim, sw, inK)

    # time-integration constants
    Mt = zeros(5)
    Mt[4] = tau * (1.0 - exp(-dt / tau))
    Mt[5] = -tau * dt * exp(-dt / tau) + tau * Mt[4]
    Mt[1] = dt - Mt[4]
    Mt[2] = -tau * Mt[1] + Mt[5]
    Mt[3] = 0.5 * dt^2 - tau * Mt[1]

    # flux related to central distribution
    Muv = moments_conserve(Mu, Mxi, 1, 0)
    MauL = moments_conserve_slope(gaL, MuL, Mxi, 2)
    MauR = moments_conserve_slope(gaR, MuR, Mxi, 2)
    # Mau = moments_conserve_slope(ga, MuR, Mxi, 2)
    MauT = moments_conserve_slope(gaT, Mu, Mxi, 1)

    fw .=
        Mt[1] .* prim[1] .* Muv .+ Mt[2] .* prim[1] .* (MauL .+ MauR) .+
        Mt[3] .* prim[1] .* MauT
    # fw .= Mt[1] .* prim[1] .* Muv .+ Mt[2] .* prim[1] .* Mau .+ Mt[3] .* prim[1] .* MauT

    # flux related to upwind distribution
    MuvL = moments_conserve(MuL1, Mxi1, 1, 0)
    MauL = moments_conserve_slope(faL, MuL1, Mxi1, 2)
    MauLT = moments_conserve_slope(faTL, MuL1, Mxi1, 1)

    MuvR = moments_conserve(MuR2, Mxi2, 1, 0)
    MauR = moments_conserve_slope(faR, MuR2, Mxi2, 2)
    MauRT = moments_conserve_slope(faTR, MuR2, Mxi2, 1)

    @. fw +=
        Mt[4] * primL[1] * MuvL - (Mt[5] + tau * Mt[4]) * primL[1] * MauL -
        tau * Mt[4] * primL[1] * MauLT + Mt[4] * primR[1] * MuvR -
        (Mt[5] + tau * Mt[4]) * primR[1] * MauR - tau * Mt[4] * primR[1] * MauRT
    # @. fw += Mt[4] * primL[1] * MuvL + Mt[4] * primR[1] * MuvR

    return nothing

end

t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(ks.set.maxTime ÷ dt) + 1
res = zero(ks.ib.wL)
for iter = 1:nt
    reconstruct!(ks, ctr)

    for i in eachindex(face)
        flux_gks!(
            face[i].fw,
            ctr[i-1].w .+ ctr[i-1].sw .* ks.ps.dx[i-1] / 2,
            ctr[i].w .- ctr[i].sw .* ks.ps.dx[i] / 2,
            ks.gas.γ,
            ks.gas.K,
            ks.gas.μᵣ,
            ks.gas.ω,
            dt,
            ks.ps.dx[i-1] / 2,
            ks.ps.dx[i] / 2,
            ctr[i-1].sw,
            ctr[i].sw,
        )
        #=flux_hll!(
            face[i].fw,
            ctr[i-1].w,
            ctr[i].w,
            ks.gas.γ,
            dt,
        )=#
    end

    for i = 1:ks.ps.nx
        @. ctr[i].w += (face[i].fw - face[i+1].fw) / ks.ps.dx[i]
        ctr[i].prim .= conserve_prim(ctr[i].w, ks.gas.γ)
    end

    t += dt
    #if t > ks.set.maxTime || maximum(res) < 5.e-7
    #    break
    #end
end

plot_line(ks, ctr)

sol = zeros(ks.ps.nx, 3)
for i = 1:ks.ps.nx
    sol[i, :] .= ctr.prim[i]
    sol[i, 3] = 1 / sol[i, 3]
end

plot(ks.ps.x[1:ks.ps.nx], sol0[:, 3])
plot!(ks.ps.x[1:ks.ps.nx], sol[:, 3], line=:dash)

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
"""
    HypertextLiteral

The `HypertextLiteral` module exports the `@htl` macro which implements
interpolation aware of hypertext escape context. It also provides for
escaping of JavaScript within the `<script>` tag.

```jldoctest
julia> v = "<1 Brown \\\"M&M's\\\"!";

julia> @htl "<span>\$v</span>"
<span><1 Brown "M&M's"!</span>

julia> @htl "<script>console.log(\$v)</script>"
<script>console.log("<1 Brown \\\"M&M's\\\"!")</script>
```

This escaping of Julia values to JavaScript values is done with `js`
function, which is not exported by default.

```jldoctest
julia> v = "<1 Brown \\\"M&M's\\\"!";

julia> @htl "<div onclick='alert(\$(HypertextLiteral.js(v)))'>"
<div onclick='alert("<1 Brown \\"M&M's\\"!")'>
```

There is also a non-standard string literal, `@htl_str` that is not
exported. It can be used with dynamically constructed templates.

See also: [`@htl`](@ref), [`HypertextLiteral.@htl_str`](@ref)
"""
module HypertextLiteral

export @htl, @htl_str

include("primitives.jl") # Wrap, Unwrap, EscapeProxy
include("macro.jl")      # @htl macro and `Result` object
include("convert.jl")    # runtime conversion of objects
include("style.jl")      # printing of content within a style tag
include("script.jl")     # printing of content within a script tag
include("lexer.jl")      # interpolate string to macro expression
include("rewrite.jl")    # macro optimizations called by interpolate

end


text:
using Statistics
using FFTW
using Plots
using BenchmarkTools
using Profile
using LinearAlgebra
using Measures
using ProgressMeter
using Plots
theme(:dark)
using Revise
push!(LOAD_PATH, pwd()*"/main")
using eqws

filter_hash = fink_filter_hash(1,8)
a = rand(256,256)
d = eqws_compute_convmap(a,filter_hash)

heatmap(d[57])

d2 = sum.(d)

d1 = eqws_compute(a,filter_hash)

count(d1 .≈ d2)

maximum(abs.(d1.-d2)) < 1e-15

rod_test = rod_image(1,1,30,35,8)

d1 = eqws_compute(rod_test,filter_hash)
d = eqws_compute_convmap(rod_test,filter_hash)
d2 = sum.(d)
count(d1 .≈ d2)
test2 = (d1 .≈ d2)
argmin(test2)

d1[57]
d2[57]

d1[57]
d2[57]

maximum(abs.(d1.-d2)) < 1e-16

scatter(d1 .≈ d2)

b = [im_rd_0_1[:,:,i]  for i in 1:size(im_rd_0_1)[3]]

b[1]

filter_hash["S1_iso_mat"]*b

filter_hash["S1_iso_mat"]

filter_hash["S1_iso_mat"]*b

function test_compute(i)
    a = rand(256,256)
    return eqws_compute_convmap(a,filter_hash)
end

using Distributed

test = pmap(test_compute,1:10)

DHC_out = hcat(test...)

DHC_out[3,:]

49*69000

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
@testset "dffunctions_test.jl" begin
    # test reading out a csv file into a DataFrame, then interpolate it to produce a function
    spG = LaserLab.CsvG(';',',')  # spanish: decimals represented with ',' delimited with ';'
    enG = LaserLab.CsvG(',','.')  # english
    df = LaserLab.load_df_from_csv("/Users/jj/JuliaProjects/LaserLab/data/objectives", "LMU-40X-NUV.csv", enG)
    wl=200.0:20:800.0
    fdf =LaserLab.dftof(wl, df, "Transmission(%)")
    @test all([fdf(df[1,"Wavelength(µm)"]) == df[1,"Transmission(%)"] for i in 1:31])

    green_area = LaserLab.qpdf(fdf, 400.0, 600.0)
    red_area = LaserLab.qpdf(fdf, 600.0, 800.0)
    @test isapprox(green_area/red_area, 1.88, rtol=0.01)

    total_area = LaserLab.qpdf(fdf, 200.0, 800.0)
    N, dfpdf = LaserLab.ftopdf(wl, fdf)
    @test total_area ≈ N
end

max_stars_repo_name: jjgomezcadenas/laserlab

max_issues_repo_name: jjgomezcadenas/laserlab

max_forks_repo_name: jjgomezcadenas/laserlab


text:
"""
  Prokaryotic auto-regulation network
"""
function autoreg(;
  gene = 10, P2_gene = 0, RNA = 0, P = 0, P2 = 0,
  k1  :: Float64 = 0.01,
  k2  :: Float64 = 10.0,
  k3  :: Float64 = 1.0,
  k3r :: Float64 = 1.0,
  k4  :: Float64 = 1.0,
  k4r :: Float64 = 10.0,
  k5  :: Float64 = 0.1,
  k6  :: Float64 = 0.01)
  m = Network("auto-regulation")

  m <= Species("gene",    gene)
  m <= Species("P2_gene", P2_gene)
  m <= Species("RNA",     RNA)
  m <= Species("P",       P)
  m <= Species("P2",      P2)

  m <= Reaction("transcription",              k1,  "gene --> gene + RNA")
  m <= Reaction("translation",                k2,  "RNA --> RNA + P")
  m <= Reaction("dimerization",               k3,  "P + P --> P2")
  m <= Reaction("dissociation",               k3r, "P2 --> P + P")
  m <= Reaction("repression binding",         k4,  "gene + P2 --> P2_gene")
  m <= Reaction("reverse repression binding", k4r, "P2_gene --> gene + P2")
  m <= Reaction("RNA degradation",            k5,  "RNA --> 0")
  m <= Reaction("protein degradation",        k6,  "P --> 0")

  return m
end

model = autoreg()


text:
using DataStructures
using GRIB

export FileIndex

"""
A mutable store for indices of a GRIB file

TODO: Should probably change this to a immutable struct
"""
mutable struct FileIndex
    "Version number used when saving/hashing index files, should change if the indexing structure changes breaking backwards-compatibility"
    allowed_protocol_version::VersionNumber

    "Path to the file the index belongs to"
    grib_path::String
    "Path to the index cache file"
    index_path::String

    "Array containing all of the index keys"
    index_keys::Array{String,1}
    "Array containing pairs of `offsets[HeaderTuple(header_values)] => offset_field`"
    offsets::Array{Pair{T where T <: NamedTuple,Int},1}
    "Array containing the length of each message in the GRIB file"
    message_lengths::Array{Int,1}
    "Dictionary of all of the loaded header values in the GRIB file"
    header_values::OrderedDict{String,Array}

    "Filters used when creating the file index"
    filter_by_keys::Dict

    FileIndex() = new(v"0.0.0")
end

function FileIndex(
    grib_path::String,
    index_keys::Array{String,1};
    filter_by_keys::Dict=Dict(),
)::FileIndex
    fileindex = FileIndex()
    fileindex.grib_path = grib_path
    fileindex.index_keys = index_keys

    if isfile(index_path(fileindex))
        from_indexfile!(fileindex)
    else
        from_gribfile!(fileindex)
        get_header_values!(fileindex)
    end

    if length(filter_by_keys) > 0
        fileindex = filter(fileindex, filter_by_keys...)
    end

    return fileindex
end

Base.getindex(obj::FileIndex, key) = obj.header_values[key]

function index_path(index::FileIndex)
    index_keys_hash = hash(
        join([index.index_keys..., index.allowed_protocol_version])
    )
    index_keys_hash = string(index_keys_hash, base=16)

    return "$(index.grib_path).$index_keys_hash.idx"
end

index_path!(index::FileIndex) = setfield!(index, :index_path, index_path(index))

function save_indexfile(index::FileIndex)
    throw("unimplemented")
end

function from_indexfile!(index::FileIndex)
    throw("unimplemented")
end

function from_gribfile!(index::FileIndex)
    offsets = OrderedDict{T where T <: NamedTuple,Int}()
    count_offsets = Dict{Int,Int}()

    index_keys = index.index_keys
    index_key_count = length(index_keys)
    index_key_symbols = Tuple(Symbol.(index_keys))
    HeaderTuple = NamedTuple{index_key_symbols}

    #  TODO: Time function to see if it is worth optimising
    #  based on gribfile.nmessages w/ known-length arrays
    #  more, or if I/O overhead too large
    GribFile(index.grib_path) do f
        message_lengths = Array{Int,1}(undef, f.nmessages)
        for (nmessage, message) in enumerate(f)
            header_values = Array{Any}(undef, index_key_count)
            for (i, key) in enumerate(index_keys)
                value = read_message(message, key)

                # #  Useful for debugging
                # value_naive = haskey(message, key) ? message[key] : missing
                # value_naive = value_naive isa Array ? Tuple(value_naive) : value_naive
                # value_naive = key == "time" ? from_grib_date_time(message) : value_naive
                # if any(ismissing.([value_naive, value]))
                #     if !all(ismissing.([value_naive, value]))
                #         println("$key: $value vs $value_naive")
                #     end
                # elseif value_naive != value
                #     println("$key: $value vs $value_naive")
                # end

                header_values[i] = value
            end

            offset = Int(message["offset"])
            if offset in keys(count_offsets)
                count_offsets[offset] += 1
                offset_field = (offset, count_offsets[offset])
            else
                count_offsets[offset] = 0
                offset_field = offset
            end

            message_lengths[nmessage] = Int(message["totalLength"])
            offsets[HeaderTuple(header_values)] = offset_field
        end
        index.message_lengths = message_lengths
    end

    index.offsets = collect(pairs(offsets))
end

function get_header_values!(index::FileIndex)
    header_values = OrderedDict{String,Array}()
    for key in index.index_keys
        header_values[key] = unique([
            offset[1][Symbol(key)]
            for offset
            in index.offsets
        ])
    end

    index.header_values = header_values
end

function getone(index::FileIndex, item)
    values = index[item]

    if length(values) != 1
        throw("Expected 1 value for $(item), found $(length(values)) instead")
    end

    return values[1]
end

function first(index::FileIndex)
    GribFile(index.grib_path) do file
        first_offset = index.offsets[1][2][1]
        #  There is a discrepancy between how offsets are defined and used
        #  in cfgrib with the GRIB file seek method and in the Julia GRIB
        #  package, in Julia seek seeks through the messages themselves not
        #  the actual offset values. Here we use the cumulative sum of the
        #  message lengths to work out which message an offset value is in.
        message_length_cumsum = cumsum(index.message_lengths)
        offset_message_index = findfirst(message_length_cumsum .> first_offset) - 1
        seek(file, offset_message_index)
        return Message(file)
    end
end

function filter_offsets(index::FileIndex; query...)
    filtered_offsets = Array{Pair{Any,Any},1}()

    for (header_values, offset_values) in index.offsets
        for (k, v) in query
            if ismissing(header_values[k]) && ismissing(v)
                append!(filtered_offsets, [Pair(header_values, offset_values)])
            elseif ismissing(header_values[k]) || header_values[k] != v
                break
            else
                append!(filtered_offsets, [Pair(header_values, offset_values)])
            end
        end
    end

    return filtered_offsets
end

function filter(index::FileIndex; query...)::FileIndex
    filtered_offsets = filter_offsets(index; query...)

    filtered_index = deepcopy(index)
    filtered_index.offsets = filtered_offsets
    filtered_index.filter_by_keys = query

    get_header_values!(filtered_index)

    return filtered_index
end

function filter!(index::FileIndex; query...)
    filtered_offsets = filter_offsets(index; query...)

    index.offsets = filtered_offsets
    index.filter_by_keys = query

    get_header_values!(index)
end


text:
function check_dimension(::Type{MathUtils}, arg0::jint, arg1::jint)
    return jcall(MathUtils, "checkDimension", void, (jint, jint), arg0, arg1)
end

function check_finite(::Type{MathUtils}, arg0::Vector{jdouble})
    return jcall(MathUtils, "checkFinite", void, (Vector{jdouble},), arg0)
end

function check_finite(::Type{MathUtils}, arg0::jdouble)
    return jcall(MathUtils, "checkFinite", void, (jdouble,), arg0)
end

function check_not_null(::Type{MathUtils}, arg0::Object)
    return jcall(MathUtils, "checkNotNull", void, (Object,), arg0)
end

function check_not_null(::Type{MathUtils}, arg0::Object, arg1::Localizable, arg2::Vector{Object})
    return jcall(MathUtils, "checkNotNull", void, (Object, Localizable, Vector{Object}), arg0, arg1, arg2)
end

function check_range_inclusive(::Type{MathUtils}, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(MathUtils, "checkRangeInclusive", void, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function check_range_inclusive(::Type{MathUtils}, arg0::jlong, arg1::jlong, arg2::jlong)
    return jcall(MathUtils, "checkRangeInclusive", void, (jlong, jlong, jlong), arg0, arg1, arg2)
end

function copy_sign(::Type{MathUtils}, arg0::jbyte, arg1::jbyte)
    return jcall(MathUtils, "copySign", jbyte, (jbyte, jbyte), arg0, arg1)
end

function copy_sign(::Type{MathUtils}, arg0::jint, arg1::jint)
    return jcall(MathUtils, "copySign", jint, (jint, jint), arg0, arg1)
end

function copy_sign(::Type{MathUtils}, arg0::jlong, arg1::jlong)
    return jcall(MathUtils, "copySign", jlong, (jlong, jlong), arg0, arg1)
end

function copy_sign(::Type{MathUtils}, arg0::jshort, arg1::jshort)
    return jcall(MathUtils, "copySign", jshort, (jshort, jshort), arg0, arg1)
end

function equals(::Type{MathUtils}, arg0::jdouble, arg1::jdouble)
    return jcall(MathUtils, "equals", jboolean, (jdouble, jdouble), arg0, arg1)
end

function hash(::Type{MathUtils}, arg0::Vector{jdouble})
    return jcall(MathUtils, "hash", jint, (Vector{jdouble},), arg0)
end

function hash(::Type{MathUtils}, arg0::jdouble)
    return jcall(MathUtils, "hash", jint, (jdouble,), arg0)
end

function max(::Type{MathUtils}, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(MathUtils, "max", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function min(::Type{MathUtils}, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(MathUtils, "min", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function normalize_angle(::Type{MathUtils}, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(MathUtils, "normalizeAngle", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function normalize_angle(::Type{MathUtils}, arg0::jdouble, arg1::jdouble)
    return jcall(MathUtils, "normalizeAngle", jdouble, (jdouble, jdouble), arg0, arg1)
end

function reduce(::Type{MathUtils}, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(MathUtils, "reduce", jdouble, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end


text:
### A Pluto.jl notebook ###
# v0.16.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

# ╔═╡ 1d6260d4-f663-11ea-03da-efe9ed63f9bd
# Explicit use of own environment instead of a local one for each notebook
begin
	using Pkg
	cd(joinpath(dirname(@__FILE__),".."))
    Pkg.activate(pwd())
    using NativeSVG
	using Plots
end

# ╔═╡ e6ff0e98-f662-11ea-03a7-e3d09e6272a6
md"""# Physical Modelling

Port of [Think Complexity chapter 7](http://greenteapress.com/complexity2/html/index.html) by Allen Downey."""

# ╔═╡ 2abf49e0-f663-11ea-25f3-2f9229de732e
md"""## Diffusion

In 1952 Alan Turing published a paper called “The chemical basis of morphogenesis”, which describes the behavior of systems involving two chemicals that diffuse in space and react with each other. He showed that these systems produce a wide range of patterns, depending on the diffusion and reaction rates, and conjectured that systems like this might be an important mechanism in biological growth processes, particularly the development of animal coloration patterns.

Turing’s model is based on differential equations, but it can be implemented using a cellular automaton.

Before we get to Turing’s model, we’ll start with something simpler: a diffusion system with just one chemical. We’ll use a 2-D CA where the state of each cell is a continuous quantity (usually between 0 and 1) that represents the concentration of the chemical.

We’ll model the diffusion process by comparing each cell with the average of its neighbors. If the concentration of the center cell exceeds the neighborhood average, the chemical flows from the center to the neighbors. If the concentration of the center cell is lower, the chemical flows the other way.

We’ll use a diffusion constant, `r`, that relates the difference in concentration to the rate of flow:"""

# ╔═╡ 40373f12-f663-11ea-256c-abd1458a8e85
function applydiffusion(array::Array{Float64, 2}, r::Float64=0.1)
    nr_y, nr_x = size(array)
    out = deepcopy(array)
    for y in 2:nr_y-1
        for x in 2:nr_x-1
            c = array[y-1, x] + array[y, x-1] + array[y, x+1] + array[y+1, x] - 4*array[y, x]
            out[y, x] += r*c
        end
    end
    out
end

# ╔═╡ 52f4e280-f663-11ea-38a4-c52a7b06b564
md"""visualisation:"""

# ╔═╡ 61437b46-f663-11ea-12c6-2ff992812dca
function visualizearray(array::Array{Float64, 2}, dim)
    (nr_y, nr_x) = size(array)
    width = dim * (nr_x - 1)
    height = dim * (nr_y - 1)
    Drawing(width=width, height=height) do
		for (j, y) in enumerate(2:nr_y-1)
			for (i, x) in enumerate(2:nr_x-1)
				gray = 80*(1-array[y, x])+10
				fill = "rgb($gray%,$gray%,$gray%"
				rect(x=i*dim, y=j*dim, width=dim, height=dim, fill=fill)
			end
		end
	end
end

# ╔═╡ 358ac3e2-f666-11ea-2058-11eedeabff5f
mutable struct Diffusion
	array :: Array{Float64, 2}
	function Diffusion()
		diffusion = new(zeros(Float64, 11, 11))
		diffusion.array[5:7, 5:7] = ones(Float64, 3, 3)
		diffusion
	end
end

# ╔═╡ 863a3fa4-f664-11ea-2597-8f18f2a862ac
diffusion = Diffusion();

# ╔═╡ 966fc752-f666-11ea-1794-195433f4cce5
visualizearray(diffusion.array, 30)

# ╔═╡ b0148fb4-f664-11ea-2d64-976658b08661
@bind togglediffusion html"<input type=button value='Next'>"

# ╔═╡ b25969c4-f665-11ea-222b-df9786a710d2
if togglediffusion === "Next"
	diffusion.array = applydiffusion(diffusion.array)
	visualizearray(diffusion.array, 30)
end

# ╔═╡ 95dc78e8-f667-11ea-1c07-8f5cc11b011b
let
	diffusion2 = Diffusion();
	for _ in 1:1000
    	diffusion2.array = applydiffusion(diffusion2.array)
	end
	visualizearray(diffusion2.array, 30)
end

# ╔═╡ c03905f2-f667-11ea-3484-b111f7c14f60
md"""## Reaction-Diffusion

Now let’s add a second chemical."""

# ╔═╡ cf7c1a2c-f667-11ea-358f-df431ec27476
function applyreactiondiffusion(
        a::Array{Float64, 2},
        b::Array{Float64, 2},
        ra::Float64=0.5, rb::Float64=0.25,
        f::Float64=0.055, k::Float64=0.062)
    nr_y, nr_x = size(a)
    a_out = deepcopy(a)
    b_out = deepcopy(b)
    for y in 2:nr_y-1
        for x in 2:nr_x-1
            reaction = a[y, x] * b[y, x]^2
            ca = 0.25*(a[y-1, x] + a[y, x-1] + a[y, x+1] + a[y+1, x]) - a[y, x]
            cb = 0.25*(b[y-1, x] + b[y, x-1] + b[y, x+1] + b[y+1, x]) - b[y, x]
            a_out[y, x] += ra*ca - reaction + f * (1 - a[y, x])
            b_out[y, x] += rb*cb + reaction - (f+k) * b[y, x]
        end
    end
    a_out, b_out
end

# ╔═╡ 0a815d80-f668-11ea-0dd3-9155b36a4134
md"""- `ra`: The diffusion rate of A (analogous to `r` in the previous section).
- `rb`: The diffusion rate of B. In most versions of this model, `rb` is about half of `ra`.
- `f`: The “feed” rate, which controls how quickly A is added to the system.
- `k`: The “kill” rate, which controls how quickly B is removed from the system.

`ca` and `ca` are the result of applying a diffusion kernel to A and B. Multiplying by `ra` and `rb` yields the rate of diffusion into or out of each cell.

The term `a*b^2` represents the rate that A and B react with each other. Assuming that the reaction consumes A and produces B, we subtract this term in the first equation and add it in the second.

The term `f * (1-a)` determines the rate that A is added to the system. Where A is near 0, the maximum feed rate is `f`. Where A approaches 1, the feed rate drops off to zero.

Finally, the term `(f+k) * b` determines the rate that B is removed from the system. As B approaches 0, this rate goes to zero.

As long as the rate parameters are not too high, the values of A and B usually stay between 0 and 1."""

# ╔═╡ 23f98030-f668-11ea-389b-ff911adaadfd
mutable struct ReactionDiffusion
	a::Array{Float64, 2}
    b::Array{Float64, 2}
	function ReactionDiffusion()
		reactiondiffusion = new(ones(Float64, 258, 258), rand(Float64, 258, 258)*0.1)
		reactiondiffusion.b[129-12:129+12, 129-12:129+12] += ones(Float64, 25, 25)*0.1
		reactiondiffusion
	end
end

# ╔═╡ 3c706f38-f669-11ea-2266-f743535073b5
reactiondiffusion = ReactionDiffusion();

# ╔═╡ 1b551fc2-f669-11ea-0d2f-ad0e089592d8
visualizearray(reactiondiffusion.b, 2)

# ╔═╡ e1988846-f668-11ea-1b26-8f7dca1bb4e7
@bind togglereactiondiffusion html"<input type=button value='Next'>"

# ╔═╡ dda066c0-f668-11ea-3096-bdea16f40db7
# f = 0.035, 0.055, 0.039 k = 0.057, 0.062, 0.065
if togglereactiondiffusion === "Next"
	for _ in 1:500
		reactiondiffusion.a, reactiondiffusion.b = applyreactiondiffusion(reactiondiffusion.a, reactiondiffusion.b, 0.5, 0.25, 0.035, 0.057)
	end
	visualizearray(reactiondiffusion.b, 2)
end

# ╔═╡ 1734b35e-f66a-11ea-10b4-317d12008809
md"""Since 1952, observations and experiments have provided some support for Turing’s conjecture. At this point it seems likely, but not yet proven, that many animal patterns are actually formed by reaction-diffusion processes of some kind."""

# ╔═╡ 40f84b2e-f66a-11ea-1dc7-3ded9e11af06
md"""## Percolation

Percolation is a process in which a fluid flows through a semi-porous material. Examples include oil in rock formations, water in paper, and hydrogen gas in micropores. Percolation models are also used to study systems that are not literally percolation, including epidemics and networks of electrical resistors.

We’ll explore a 2-D CA that simulates percolation.

- Initially, each cell is either “porous” with probability `q` or “non-porous” with probability `1-q`.
- When the simulation begins, all cells are considered “dry” except the top row, which is “wet”.
- During each time step, if a porous cell has at least one wet neighbor, it becomes wet. Non-porous cells stay dry.
- The simulation runs until it reaches a “fixed point” where no more cells change state.

If there is a path of wet cells from the top to the bottom row, we say that the CA has a “percolating cluster”.

Two questions of interest regarding percolation are (1) the probability that a random array contains a percolating cluster, and (2) how that probability depends on `q`."""

# ╔═╡ 77dcd828-f66a-11ea-0942-c5d613d828b7
function applypercolation(array::Array{Float64, 2})
    nr_y, nr_x = size(array)
    out = deepcopy(array)
    for y in 2:nr_y-1
        for x in 2:nr_x-1
            if out[y, x] > 0.0
                c = array[y-1, x] + array[y, x-1] + array[y, x+1] + array[y+1, x]
                if c ≥ 0.5
                    out[y, x] = 0.5
                end
            end
        end
    end
    out
end

# ╔═╡ 818e5534-f66a-11ea-07b2-7722084e1e07
md"visualisation:"

# ╔═╡ 44c8b17a-f66b-11ea-39af-f554593e33eb
mutable struct Wall
	array::Array{Float64, 2}
	function Wall(n, q)
    	array = zeros(Float64, n+2, n+2)
		array[2, 2:n+1] = ones(Float64, n)*0.5
		array[3:n+1, 2:n+1] = rand(Float64, n-1, n)
		for y in 3:n+1
			for x in 2:n+1
				if array[y, x] < q
					array[y, x] = 0.1
				else
					array[y, x] = 0.0
				end
			end
		end
		new(array)
	end
end

# ╔═╡ d9ccf8f8-f66b-11ea-2c81-0d771c2e900e
wall = Wall(100, 0.62);

# ╔═╡ e492630e-f66b-11ea-3aed-f37a400243be
visualizearray(wall.array, 8)

# ╔═╡ e048231c-f66b-11ea-191d-8517a6c65bc5
@bind togglepercolation html"<input type=button value='Next'>"

# ╔═╡ 11847488-f66c-11ea-1aca-4b1613dbfb8e
if togglepercolation === "Next"
	wall.array = applypercolation(wall.array)
	visualizearray(wall.array, 8)
end

# ╔═╡ 3563b208-f66c-11ea-3c1a-77f799805d06
md"""## Phase Change

Now let’s test whether a random array contains a percolating cluster."""

# ╔═╡ 2b9c8a86-f66c-11ea-0447-1fb7ba83d3a1
function testpercolation(array::Array{Float64, 2}, vis=false)
    numberwet = count(x->x==0.5, array[3:101, 2:101])
    while true
        array = applypercolation(array)
        if count(x->x==0.5, array[101, 2:101]) > 0
			if vis
				return true, visualizearray(array, 8)
			else
            	return true
			end
        end
        newnumberwet = count(x->x==0.5, array[3:101, 2:101])
        if numberwet == newnumberwet
            if vis
				return false, visualizearray(array, 8)
			else
            	return false
			end
        end
        numberwet = newnumberwet
    end
end

# ╔═╡ c71f0f24-f66c-11ea-0ba7-d765f3998ed9
let
	wall = Wall(100, 0.5)
	result, drawing = testpercolation(wall.array, true)
	drawing
end

# ╔═╡ e51771ec-f66c-11ea-1606-fd159b8d216f
md"""To estimate the probability of a percolating cluster, we generate many random arrays and test them."""

# ╔═╡ 1254b21a-f66e-11ea-2b36-7985687edd02
function estimateprob(;n=100, q=0.5, iters=100)
    t = Bool[]
    for _ in 1:iters
        wall = Wall(n, q)
        push!(t, testpercolation(wall.array))
    end
    count(x->x, t) / iters
end

# ╔═╡ dda2cd12-f66c-11ea-3f31-837c9a749b5c
estimateprob(q = 0.50)

# ╔═╡ 2b8081ae-f66e-11ea-28dd-d160aa6c80ef
md"""We can estimate the critical value more precisely using a random walk. Starting from an initial value of `q`, we construct a wall and check whether it has a percolating cluster. If so, `q` is probably too high, so we decrease it. If not, `q` is probably too low, so we increase it."""

# ╔═╡ 3d494330-f66e-11ea-118c-7f36892d863e
function findcritical(;n=100, q=0.5, iters=100)
    qs = [q]
    for _ in 1:iters
        wall = Wall(n, q)
        if testpercolation(wall.array)
            q -= 0.004
        else
            q += 0.004
        end
        push!(qs, q)
    end
    qs
end

# ╔═╡ 501878b4-f66e-11ea-270f-d7542715acbb
findcritical()

# ╔═╡ 59540966-f66e-11ea-2cac-6fac6c62910b
md"""With `n=100` the mean of `qs` is about 0.59; this value does not seem to depend on `n`.

The rapid change in behavior near the critical value is called a phase change by analogy with phase changes in physical systems, like the way water changes from liquid to solid at its freezing point.

A wide variety of systems display a common set of behaviors and characteristics when they are at or near a critical point. These behaviors are known collectively as critical phenomena."""

# ╔═╡ 87c7816a-f66e-11ea-02fe-1752fee2c7eb
md"""## Fractals

To understand fractals, we have to start with dimensions.

For simple geometric objects, dimension is defined in terms of scaling behavior. For example, if the side of a square has length `l`, its area is `l^2`. The exponent, 2, indicates that a square is two-dimensional. Similarly, if the side of a cube has length `l`, its volume is `l^3`, which indicates that a cube is three-dimensional.

More generally, we can estimate the dimension of an object by measuring some kind of size (like area or volume) as a function of some kind of linear measure (like the length of a side).

As an example, I’ll estimate the dimension of a 1-D cellular automaton by measuring its area (total number of “on” cells) as a function of the number of rows."""

# ╔═╡ af69c480-f66e-11ea-3e4a-1db5d58bdd7b
function inttorule1dim(val::UInt8)
    digs = BitArray(digits(val, base=2))
    for i in length(digs):7
        push!(digs, false)
    end
    digs
end

# ╔═╡ d6cab994-f66e-11ea-1734-dd3ab2d22e92
function applyrule1dim(rule::BitArray{1}, bits::BitArray{1})
    val = 1 + bits[3] + 2*bits[2] + 4*bits[1]
    rule[val]
end

# ╔═╡ daa15a5a-f66e-11ea-1983-51732624404b
function step1dim(x₀::BitArray{1}, rule::BitArray{1}, steps::Int64)
    xs = [x₀]
    len = length(x₀)
    for i in 1:steps
        x = copy(x₀)
        for j in 2:len-1
            x[j] = applyrule1dim(rule, xs[end][j-1:j+1])
        end
        push!(xs, x)
    end
    xs
end

# ╔═╡ e047f8b8-f66e-11ea-05af-398e5d22ccf8
function visualize1dim(res, dim)
    width = dim * (length(res[1]) + 1)
    height = dim * (length(res) + 1)
    Drawing(width = width, height = height) do
        for (i, arr) in enumerate(res)
            for (j, val) in enumerate(arr)
                fill = if val "grey" else "lightgrey" end
                rect(x = j*dim, y = i*dim, width = dim, height = dim, fill = fill)
            end
        end
    end
end

# ╔═╡ ea25224a-f66e-11ea-1c90-b30dba111177
begin
	x₀ = falses(65)
	x₀[33] = true
	res = step1dim(x₀, inttorule1dim(UInt8(18)), 31) # 20, 50, 18
	visualize1dim(res, 10)
end

# ╔═╡ 070b57aa-f66f-11ea-1a12-bd1529268d4f
md"""I’ll estimate the dimension of these CAs with the following function, which counts the number of on cells after each time step."""

# ╔═╡ 041baf7c-f66f-11ea-270b-2179bd949a67
function countcells(rule, n=501)
    x₀ = falses(2*n+3)
    x₀[n+2] = true
    res = step1dim(x₀, inttorule1dim(UInt8(rule)), n)
    cells = [1]
    for i in 2:n
        push!(cells, cells[end]+sum(line->count(cell->cell, line), res[i]))
    end
    cells
end

# ╔═╡ 73fd7b5e-f66f-11ea-12cd-dd36181cf956
begin
	n = 501;
	rule = 20 # 20, 50, 18
	plot(1:n, 1:n, xaxis=:log, yaxis=:log, label="d = 1")
	plot!(1:n, (1:n).^2, xaxis=:log, yaxis=:log, label="d = 2")
	plot!(1:n, countcells(rule, n), xaxis=:log, yaxis=:log, label="rule $rule")
end

# ╔═╡ 81657f1a-f66f-11ea-183b-0ba3f1c3a38e
md"""Rule 20 (left) produces 3 cells every 2 time steps, so the total number of cells after $i$ steps is $y = 1.5i$. Taking the log of both sides, we have $\log y = \log 1.5 + \log i$, so on a log-log scale, we expect a line with slope 1. In fact, the estimated slope of the line is 1.01.

Rule 50 (center) produces i+1 new cells during the ith time step, so the total number of cells after $i$ steps is $y = i^2 + i$. If we ignore the second term and take the log of both sides, we have $\log y \approx 2 \log i$, so as $i$ gets large, we expect to see a line with slope 2. In fact, the estimated slope is 1.97.

Finally, for Rule 18 (right), the estimated slope is about 1.57, which is clearly not 1, 2, or any other integer. This suggests that the pattern generated by Rule 18 has a “fractional dimension”; that is, it is a fractal.

This way of estimating a fractal dimension is called box-counting."""

# ╔═╡ feefb6e4-f66f-11ea-2998-a30311bee9f4
md"""## Fractals and Percolation

Now let’s get back to percolation models.

To estimate their fractal dimension, we can run CAs with a range of sizes, count the number of wet cells in each percolating cluster, and then see how the cell counts scale as we increase the size of the array."""

# ╔═╡ 21d918ee-f670-11ea-0359-cbb1809bab9b
function percolationwet(array::Array{Float64, 2})
    numberwet = count(x->x==0.5, array[3:end-1, 2:end-1])
    while true
        array = applypercolation(array)
        if count(x->x==0.5, array[end-1, 2:end-1]) > 0
            break
        end
        newnumberwet = count(x->x==0.5, array[3:end-1, 2:end-1])
        if numberwet == newnumberwet
            break
        end
        numberwet = newnumberwet
    end
    numberwet
end

# ╔═╡ 3206914e-f670-11ea-3072-d162f584fbf0
let
	res = Float64[]
	sizes = 10:10:200
	q = 0.4 # 0.4, 0.8, 0.596
	for size in sizes
		wall = Wall(size, q)
		push!(res, percolationwet(wall.array))
	end
	plot(sizes, sizes, xaxis=:log, yaxis=:log, label="d = 1")
	plot!(sizes, (sizes).^2, xaxis=:log, yaxis=:log, label="d = 2")
	plot!(sizes, res, xaxis=:log, yaxis=:log, seriestype=:scatter, label="q = $q")
end

# ╔═╡ 4da9c3ba-f670-11ea-1627-d7c5e6b36742
md"""The dots show the number of cells in each percolating cluster. The slope of a line fitted to these dots is often near 1.85, which suggests that the percolating cluster is, in fact, fractal when `q` is near the critical value.

When `q` is larger than the critical value, nearly every porous cell gets filled, so the number of wet cells is close to `q * size^2`, which has dimension 2.

When `q` is substantially smaller than the critical value, the number of wet cells is proportional to the linear size of the array, so it has dimension 1."""

# ╔═╡ Cell order:
# ╟─e6ff0e98-f662-11ea-03a7-e3d09e6272a6
# ╟─1d6260d4-f663-11ea-03da-efe9ed63f9bd
# ╟─2abf49e0-f663-11ea-25f3-2f9229de732e
# ╠═40373f12-f663-11ea-256c-abd1458a8e85
# ╟─52f4e280-f663-11ea-38a4-c52a7b06b564
# ╠═61437b46-f663-11ea-12c6-2ff992812dca
# ╠═358ac3e2-f666-11ea-2058-11eedeabff5f
# ╠═863a3fa4-f664-11ea-2597-8f18f2a862ac
# ╠═966fc752-f666-11ea-1794-195433f4cce5
# ╠═b0148fb4-f664-11ea-2d64-976658b08661
# ╠═b25969c4-f665-11ea-222b-df9786a710d2
# ╠═95dc78e8-f667-11ea-1c07-8f5cc11b011b
# ╟─c03905f2-f667-11ea-3484-b111f7c14f60
# ╠═cf7c1a2c-f667-11ea-358f-df431ec27476
# ╟─0a815d80-f668-11ea-0dd3-9155b36a4134
# ╠═23f98030-f668-11ea-389b-ff911adaadfd
# ╠═3c706f38-f669-11ea-2266-f743535073b5
# ╠═1b551fc2-f669-11ea-0d2f-ad0e089592d8
# ╠═e1988846-f668-11ea-1b26-8f7dca1bb4e7
# ╠═dda066c0-f668-11ea-3096-bdea16f40db7
# ╟─1734b35e-f66a-11ea-10b4-317d12008809
# ╟─40f84b2e-f66a-11ea-1dc7-3ded9e11af06
# ╠═77dcd828-f66a-11ea-0942-c5d613d828b7
# ╟─818e5534-f66a-11ea-07b2-7722084e1e07
# ╠═44c8b17a-f66b-11ea-39af-f554593e33eb
# ╠═d9ccf8f8-f66b-11ea-2c81-0d771c2e900e
# ╠═e492630e-f66b-11ea-3aed-f37a400243be
# ╠═e048231c-f66b-11ea-191d-8517a6c65bc5
# ╠═11847488-f66c-11ea-1aca-4b1613dbfb8e
# ╟─3563b208-f66c-11ea-3c1a-77f799805d06
# ╠═2b9c8a86-f66c-11ea-0447-1fb7ba83d3a1
# ╠═c71f0f24-f66c-11ea-0ba7-d765f3998ed9
# ╟─e51771ec-f66c-11ea-1606-fd159b8d216f
# ╠═1254b21a-f66e-11ea-2b36-7985687edd02
# ╠═dda2cd12-f66c-11ea-3f31-837c9a749b5c
# ╟─2b8081ae-f66e-11ea-28dd-d160aa6c80ef
# ╠═3d494330-f66e-11ea-118c-7f36892d863e
# ╠═501878b4-f66e-11ea-270f-d7542715acbb
# ╟─59540966-f66e-11ea-2cac-6fac6c62910b
# ╟─87c7816a-f66e-11ea-02fe-1752fee2c7eb
# ╠═af69c480-f66e-11ea-3e4a-1db5d58bdd7b
# ╠═d6cab994-f66e-11ea-1734-dd3ab2d22e92
# ╠═daa15a5a-f66e-11ea-1983-51732624404b
# ╠═e047f8b8-f66e-11ea-05af-398e5d22ccf8
# ╠═ea25224a-f66e-11ea-1c90-b30dba111177
# ╟─070b57aa-f66f-11ea-1a12-bd1529268d4f
# ╠═041baf7c-f66f-11ea-270b-2179bd949a67
# ╠═73fd7b5e-f66f-11ea-12cd-dd36181cf956
# ╟─81657f1a-f66f-11ea-183b-0ba3f1c3a38e
# ╟─feefb6e4-f66f-11ea-2998-a30311bee9f4
# ╠═21d918ee-f670-11ea-0359-cbb1809bab9b
# ╠═3206914e-f670-11ea-3072-d162f584fbf0
# ╟─4da9c3ba-f670-11ea-1627-d7c5e6b36742


text:
export AbstractStage,
    PreExperimentStage,
    PostExperimentStage,
    PreEpisodeStage,
    PostEpisodeStage,
    PreActStage,
    PostActStage

abstract type AbstractStage end

struct PreExperimentStage <: AbstractStage end
struct PostExperimentStage <: AbstractStage end
struct PreEpisodeStage <: AbstractStage end
struct PostEpisodeStage <: AbstractStage end
struct PreActStage <: AbstractStage end
struct PostActStage <: AbstractStage end

(p::AbstractPolicy)(::AbstractStage, ::AbstractEnv) = nothing

RLBase.optimise!(::AbstractPolicy) = nothing


text:
"""
    render(m::Union{MDP,POMDP}, step::NamedTuple)

Return a renderable representation of the step in problem `m`.

The renderable representation may be anything that has `show(io, mime, x)`
methods. It could be a plot, svg, Compose.jl context, Cairo context, or image.

# Arguments
`step` is a `NamedTuple` that contains the states, action, etc. corresponding
to one transition in a simulation. It may have the following fields:
- `t`: the time step index
- `s`: the state at the beginning of the step
- `a`: the action
- `sp`: the state at the end of the step (s')
- `r`: the reward for the step
- `o`: the observation
- `b`: the belief at the
- `bp`: the belief at the end of the step
- `i`: info from the model when the state transition was calculated
- `ai`: info from the policy decision
- `ui`: info from the belief update

Keyword arguments are reserved for the problem implementer and can be used to control appearance, etc.

# Important Notes
- `step` may not contain all of the elements listed above, so `render` should check for them and render only what is available
- `o` typically corresponds to `sp`, so it is often clearer for POMDPs to render `sp` rather than `s`.
"""
function render(m::Union{MDP,POMDP}, step)
    @warn("No implementation of POMDPModelTools.render(m::$m, step) found. Falling back to text default.", maxlog=1)
    io = IOBuffer()
    ioc = IOContext(io, :short=>true)
    try
        for (k, v) in pairs(step)
            print(ioc, k)
            print(ioc, ": ")
            show(ioc, v)
            println(ioc)
        end
    finally
        println(ioc, """

            Please implement POMDPModelTools.render(m::$(typeof(m)), step) to enable visualization.
            """)
    end
    return String(take!(io))
end

render(m::Union{MDP, POMDP}) = render(m, NamedTuple())


text:
let
    solver = VanillaAStar()
    set_iteration_limit!(solver,100)
    set_verbosity!(solver,0)
    for f in mapf_test_problems()
        mapf = f()
        s0 = get_start(mapf,1)
        sF = get_goal(mapf,1)
        env = build_env(mapf,initialize_root_node(mapf),1)
        path, cost = CRCBS.a_star!(solver,env,s0)
        @test get_distance(env,s0,sF) == cost
        reset_solver!(solver)
    end
end


text:
"""
```
transition(m::PoolModel{T}) where {T<:AbstractFloat}
```

Assign transition equation

```
x_t = (1-ρ) μ + ρ x_{t-1} + sqrt{1 - ρ^2} σ ϵ_t
ϵ_t ∼ iid N(0,1), x_0 ∼ N(μ,σ^2)
λ_t = Φ(x_t)
```
where Φ(⋅) is the cdf of a N(0,1) random variable, F_ϵ is the distribution of ϵ_t,
and F_λ is the distribution of λ(x_0).
"""
function transition(m::PoolModel{T}) where {T<:AbstractFloat}
    @inline Φ(x::Vector{Float64}, ϵ::Vector{Float64}) = abs.([0;1] .-
                                                             (cdf.(Normal(), (1 - m[:ρ].value) .*
                                                                   m[:μ].value .+ m[:ρ].value .*
                                                                   quantile(Normal(),x[1]) .+
                                                                   sqrt(1 - m[:ρ].value^2) .*
                                                                   m[:σ].value .* ϵ)))
    F_ϵ = Normal(0.,1.)
    F_λ = Uniform(0.,1.)
    return Φ, F_ϵ, F_λ
end


text:
println("Compile libraries on main processor..."); flush(stdout)
include("imports.jl")
println("Finished loading libraries."); flush(stdout)

using Distributed
rmprocs(workers())
addprocs(32)
println("Compile libraries on workers ..."); flush(stdout)
@everywhere include("imports.jl")
println("Finished loading libraries."); flush(stdout)

Ks = collect(2:9)
markers = [:CD3z, :EOMES, :Perforin, :Granzyme_A, :Siglec7, :LAG3, :CD56, :CD57]
skewtmix = [true, false]
sims = dict_list(Dict(:K => Ks, :marker => markers, :skewtmix => skewtmix))

# Run simulations in parallel.
@time res = pmap(run, sims, on_error=identity)

# Print results status.
foreach(z -> println("$(z[1]) => $(z[2])"), zip(sims, res))

# Post process
@time pp_res = pmap(postprocess, sims)

# Combine results, plot DIC.
for marker in markers
  imdir = mkpath(joinpath(Info.resultsdir_datastudy, simname, "img"))
  plot(size=plotsize)
  xlabel!("K")
  ylabel!("DIC")
  for stm in skewtmix
    dics = Float64[]
    for K in Ks
      info_path = joinpath(Info.resultsdir_datastudy, simname,
                           savename(Dict(:marker => marker, :K => K, :skewtmix => stm)),
                           "info.txt")
      model_info = open(f->read(f, String), info_path)
      _dic = parse_dic(model_info)
      append!(dics, _dic)
    end
    modelname = stm ? "Skew-t mixture" : "Normal mixture"

    # FIXME: This was hardcoded...
    pos = marker in (:LAG3, ) ? :bottomleft : :topright

    plot!(Ks, dics, marker=:square, ms=8, label=modelname, legend=pos)
  end
  savefig(joinpath(imdir, "dic_marker=$(marker).pdf"))
  closeall()
end

# Send results to S3
Util.s3sync(from="$(Info.resultsdir_datastudy)/$(simname)",
            to="$(Info.awsbucket_datastudy)/$(simname)",
            tags=`--exclude '*.nfs'`)


text:
using Pkg

Pkg.activate("../..")

using ColonyRuntime
using ColonyRuntime.Crypto

function fib(n)
  if n == 0
    return 0
  elseif n == 1
    return 1
  else
    return fib(n-2) + fib(n-1)
  end
end

function main(args)
  colony_prvkey = args[1]
  colonyid = Crypto.id(colony_prvkey)
  server = ColonyRuntime.ColoniesServer("localhost", 8080)

  # register a runtime
  println("- registering a new runtime to colony " * colonyid)
  runtime_prvkey = Crypto.prvkey()
  runtime = ColonyRuntime.Runtime(Crypto.id(runtime_prvkey), "fibonacci_solver", "fibonacci_solver", colonyid, "AMD Ryzen 9 5950X (32) @ 3.400GHz", 32, 80326, "NVIDIA GeForce RTX 2080 Ti Rev. A", 1)
  runtime = ColonyRuntime.addruntime(server, runtime, colony_prvkey)

  # and approve it so that can use the api
  println("- approving runtime " * runtime.runtimeid)
  ColonyRuntime.approveruntime(server, runtime.runtimeid, colony_prvkey)

  # request a waiting process
  try
    println("- assign process")
    assigned_process = ColonyRuntime.assignprocess(server, colonyid, runtime_prvkey)
    fibonacci_num = parse(Int64, assigned_process.attributes[1].value)
    println("  fibonacci_num: ", fibonacci_num)
    res = fib(fibonacci_num)
    println("  result: ", res)

    # add an attribute to the process
    println("- add result attribute")
    attribute = ColonyRuntime.Attribute(assigned_process.processid, "result", string(res))
    ColonyRuntime.addattribute(server, attribute, runtime_prvkey)

    # close the process
    println("- close process")
    ColonyRuntime.closeprocess(server, assigned_process.processid, true, runtime_prvkey)
  catch err
    println(err)
    println("No waiting process found")
    return
  end

end

main(ARGS)


text:
# rewrite of [TensorFlow Probability Example](https://medium.com/tensorflow/an-introduction-to-probabilistic-programming-now-available-in-tensorflow-probability-6dcc003ca29e) in `ADCME`
using Statistics
using ADCME
using PyCall
using PyPlot
using DelimitedFiles

if !("challenger_data.csv" in readdir("."))
    download("https://raw.githubusercontent.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/master/Chapter2_MorePyMC/data/challenger_data.csv", "challenger_data.csv")
end

F = readdlm("challenger_data.csv",',', header=true)[1][:,2:3]
I = ones(Bool, size(F,1))
for i = 1:size(F,1)
    if F[i,2]=="NA" || F[i,2] == "Challenger Accident"
        I[i] = false
    end
end
F = Array{Float32}(F[I,:]);

sess = Session()

temperature = constant(F[:,1], dtype=Float32)
D = constant(F[:,2], dtype=Float32)

function challenger_joint_log_prob(D, temperature, α, β)
    rv_alpha = Normal(loc=0., scale=1000.)
    rv_beta = Normal(loc=0., scale=1000.)
    logistic_p = 1.0/(1+exp(β*temperature+α ))
    rv_observed = Bernoulli(probs = logistic_p)

    return log_prob(rv_alpha, α) + log_prob(rv_beta, β) +
        sum(log_prob(rv_observed, D))
end

number_of_steps = 40000
burnin = 38000
initial_chain_state = [
    constant(0.0, dtype=Float32),
    constant(0.0, dtype=Float32)
]

function unnormalized_posterior_log_prob(args...)
    challenger_joint_log_prob(D, temperature,args...)
end

unconstraining_bijectors = [Identity(), Identity()]

variable_scope("mcmc", reuse=AUTO_REUSE) do
global step_size = get_variable(
        "step_size",
        initializer=constant(0.5, dtype=Float32),
        trainable=false
    )
end

ham=TransformedTransitionKernel(HamiltonianMonteCarlo(
        target_log_prob_fn=unnormalized_posterior_log_prob,
        num_leapfrog_steps=2,
        step_size_update_fn=make_simple_step_size_update_policy(),
        state_gradients_are_stopped=true,
        step_size=step_size),
        bijector=unconstraining_bijectors)

states,is_accepted_, kernel_results = sample_chain(
    num_results = number_of_steps,
    num_burnin_steps = burnin,
    current_state=initial_chain_state,
    kernel=ham
    )

init(sess)
posterior_α, posterior_β, is_accepted =  run(sess, [states[1], states[2], is_accepted_])

println("Accepted rate=", sum(is_accepted)/length(is_accepted))

function logistic(x, beta, alpha)
    return 1.0 ./ (1.0 .+ exp.(beta*x .+ alpha) )
end

figure()
xval = LinRange(minimum(F[:,1])-1.0, maximum(F[:,1])+1.0, 150)
y1 = logistic(xval, posterior_β[end-2000], posterior_α[end-2000])
y2 = logistic(xval, posterior_β[end-8], posterior_α[end-8])
y3 = logistic(xval, mean(posterior_β), mean(posterior_α))
scatter(F[:,1], F[:,2])
plot(xval, y1, label="posterior1")
plot(xval, y2, label="posterior2")
plot(xval, y3, label="mean")
ylim(-0.1,1.1)
legend()

# challenger
figure()
plt.hist(logistic(31, posterior_β[end-2000:end], posterior_α[end-2000:end]), 50)