text:
const DATA_DIR = joinpath(MODULE_DIR, "..", "data")
using CategoricalArrays

function MMI.matrix(table)
    cols = columns(table)
    a, b = length(cols[1]), length(cols)
    matrix = Matrix{Float64}(undef, a, b)
    @inbounds for i=1:b
        matrix[:, i] = cols[i]
    end
    return matrix
end

"""
load_dataset(fpath)

Load one of standard dataset like Chinatown etc assuming the file is a
comma separated file.

"""
function load_dataset(fname::String)
    fpath = joinpath(DATA_DIR, fname)
    data_raw = load(fpath, header_exists=false)
    Xmatrix = matrix(table(data_raw)) #check if we can use data_raw directly for y
    width = size(Xmatrix)[2]
    #return table(view(Xmatrix, :, 1:(width-1))), CategoricalArray(view(Xmatrix, :, width))   moving away from IndexedTable
    return Xmatrix[:, 1:(width-1)], CategoricalArray(view(Xmatrix, :, width))
end

#NOTE: permutedims(cat(mat, mat2, dims=3),[1,3,2]) For Motions, eg for multi-dimensional dataset.

"""
   `load_ts_file(fpath; return_array=false)`

`load_ts_file` takes the path to `ts` files and returns the table (IndexedTable),
where `fpath` is path to `ts` file located on your computer.
on `return_array=true` reurns an Array.
eg
```julia
   tableT = load_ts_file("/Users/your_user_name/../Adiac_TRAIN.ts")
   arr = load_ts_file("/Users/your_user_name/../Adiac_TRAIN.ts", return_array=true)
```
"""
function load_ts_file(fpath; return_array=false)
#`readuntil` and `readlines` are used in combination so that we read data
# required only for transformation into julia array and table.
    data = open(fpath) do input
          readuntil(input, "@data")
          readlines(input)
    end
    data = data[2:end] # removes  the empty string.
    # split the string and convert each element into Float64.
    arrays = map(i -> parse.(Float64, split(i, r"[:,]")) , data)
    Xmatrix = transpose(hcat(arrays...)) # Creates 2D Array.
    if return_array == true
        return Xmatrix
    else
        return table(eachcol(Xmatrix[:, 1:end-1])...), CategoricalArray(Xmatrix[:, end])
    end
end

"""
   `ts_dataset(dataset; split=false)`
`ts_dataset` takes two inputes, `dataset` is the name of datasets available in the
`data` folder & `split` specifile the train or test part.
Reurns dataset and target variable.
eg
```julia
   X, y = ts_dataset("Adiac")
   X_test, y_test = ts_dataset("Adiac", split="test")
```
"""
function ts_dataset(dataset::String; split=false)
    if split in ["test", "train"]
         fpath = joinpath(DATA_DIR, dataset, string(dataset, "_", uppercase(split), ".ts" ))
         return load_ts_file(fpath)
    elseif split == false
         v = []
         for i in ["test", "train"]
             fpath = joinpath(DATA_DIR, dataset, string(dataset, "_", uppercase(i), ".ts" ))
             push!(v, fpath)
         end
         test, train =  load_ts_file.(v)
         return merge(test[1], train[1]), vcat(test[2], train[2])
    else
        throw(ArgumentError("Invalid `split` value: $split"))
    end
end

"""
   `load_from_tsfile_to_NDArray(fpath)`

`load_from_tsfile_to_NDArray` takes the path to `ts` files and reurns an Array.
eg
```julia
   tableT = load_from_tsfile_to_NDArray("/Users/your_user_name/../BasicMotions_TRAIN.ts")
   arr = load_from_tsfile_to_NDArray("/Users/your_user_name/../BasicMotions_TRAIN.ts")
```
"""
function load_from_tsfile_to_NDArray(fpath)
    data = open(fpath) do input
        readuntil(input, "@data")
        readlines(input)
    end
    data = split.(data[2:end], r"[:]")
    instance = length(data)
    dimensions  = length(data[1]) - 1
    series_dim1 = split(data[1][1], r"[,]")
    serieslength = length(series_dim1)
    NDMatrix = zeros(instance, serieslength, dimensions)
    data = cat(data..., dims=3)
    @inbounds for i=1:instance
        @inbounds for j=1:dimensions
            NDMatrix[i,:,j] = parse.(Float64, split(data[j,1,i], r"[,]"))
        end
    end
    return NDMatrix, categorical(data[dimensions+1,1,:])
end

"""
   `load_NDdataset(dataset::String; split=false)`
`load_NDdataset` takes two inputes, `dataset` is the name of datasets available in the
`data` folder & `split` specifile the train or test part.
Reurns dataset and target variable.
eg
```julia
   X, y = ts_dataset("BasicMotions")
   X_test, y_test = ts_dataset("BasicMotions", split="test")
```
"""
function load_NDdataset(dataset::String; split=false)
    if split in ["test", "train"]
        fpath = joinpath(DATA_DIR, dataset, string(dataset, "_", uppercase(split), ".ts" ))
        X, y = load_from_tsfile_to_NDArray(fpath)
        return X, y
    elseif split == false
        v = []
        for i in ["test", "train"]
            fpath = joinpath(DATA_DIR, dataset, string(dataset, "_", uppercase(i), ".ts" ))
            push!(v, fpath)
        end
        test, train =  load_from_tsfile_to_NDArray.(v)
        return vcat(test[1], train[1]), vcat(test[2], train[2])
    else
        throw(ArgumentError("Invalid `split` value: $split"))
    end
end

univariate_datasets = [
    "ACSF1",
    "Adiac",
    "AllGestureWiimoteX",
    "AllGestureWiimoteY",
    "AllGestureWiimoteZ",
    "ArrowHead",
    "Coffee",
    "Beef",
    "BeetleFly",
    "BirdChicken",
    "BME",
    "Car",
    "CBF",
    "Chinatown",
    "ChlorineConcentration",
    "CinCECGTorso",
    "Coffee",
    "Computers",
    "CricketX",
    "CricketY",
    "CricketZ",
    "Crop",
    "DiatomSizeReduction",
    "DistalPhalanxOutlineCorrect",
    "DistalPhalanxOutlineAgeGroup",
    "DistalPhalanxTW",
    "DodgerLoopDay",
    "DodgerLoopGame",
    "DodgerLoopWeekend",
    "Earthquakes",
    "ECG200",
    "ECG5000",
    "ECGFiveDays",
    "ElectricDevices",
    "EOGHorizontalSignal",
    "EOGVerticalSignal",
    "EthanolLevel",
    "FaceAll",
    "FaceFour",
    "FacesUCR",
    "FiftyWords",
    "Fish",
    "FordA",
    "FordB",
    "FreezerRegularTrain",
    "FreezerSmallTrain",
    "Fungi",
    "GestureMidAirD1",
    "GestureMidAirD2",
    "GestureMidAirD3",
    "GesturePebbleZ1",
    "GesturePebbleZ2",
    "Ham",
    "HandOutlines",
    "Haptics",
    "Herring",
    "InlineSkate",
    "InsectEPGRegularTrain",
    "InsectEPGSmallTrain",
    "InsectWingbeatSound",
    "ItalyPowerDemand",
    "LargeKitchenAppliances",
    "Lightning2",
    "Lightning7",
    "Mallat",
    "Meat",
    "MedicalImages",
    "MelbournePedestrian",
    "MiddlePhalanxOutlineCorrect",
    "MiddlePhalanxOutlineAgeGroup",
    "MiddlePhalanxTW",
    "MixedShapesRegularTrain",
    "MixedShapesSmallTrain",
    "MoteStrain",
    "NonInvasiveFetalECGThorax1",
    "NonInvasiveFetalECGThorax2",
    "OliveOil",
    "OSULeaf",
    "PhalangesOutlinesCorrect",
    "Phoneme",
    "PickupGestureWiimoteZ",
    "PigAirwayPressure",
    "PigArtPressure",
    "PigCVP",
    "PLAID",
    "Plane",
    "PowerCons",
    "ProximalPhalanxOutlineCorrect",
    "ProximalPhalanxOutlineAgeGroup",
    "ProximalPhalanxTW",
    "RefrigerationDevices",
    "Rock",
    "ScreenType",
    "SemgHandGenderCh2",
    "SemgHandMovementCh2",
    "SemgHandSubjectCh2",
    "ShakeGestureWiimoteZ",
    "ShapeletSim",
    "ShapesAll",
    "SmallKitchenAppliances",
    "SmoothSubspace",
    "SonyAIBORobotSurface1",
    "SonyAIBORobotSurface2",
    "StarlightCurves",
    "Strawberry",
    "SwedishLeaf",
    "Symbols",
    "SyntheticControl",
    "ToeSegmentation1",
    "ToeSegmentation2",
    "Trace",
    "TwoLeadECG",
    "TwoPatterns",
    "UMD",
    "UWaveGestureLibraryAll",
    "UWaveGestureLibraryX",
    "UWaveGestureLibraryY",
    "UWaveGestureLibraryZ",
    "Wafer",
    "Wine",
    "WordSynonyms",
    "Worms",
    "WormsTwoClass",
    "Yoga",
]


text:
import Clang.wrap_c
import DataStructures: DefaultDict
import Base.Meta.isexpr
using Match

indexh         = joinpath(JULIA_HOME, "../include/clang-c/Index.h")
clang_includes = [joinpath(JULIA_HOME, "../lib/clang/3.3/include"), joinpath(dirname(indexh), "..")]

av_headers = ["/usr/include/avbin.h"]
for header in av_headers
    println(header)
end

get_header_library(hpath) = :libavbin

get_header_outfile(hpath) = "avbin.jl"

check_use_header(top_h, hpath) = basename(top_h) == "avbin.h"

function rewrite_fn(e::Expr)
    (head, call, body) = @match e Expr(head, [call, body], _) => (head, call, body)

    parms = Any[]
    content = Any[]

    fn_name = call.args[1]
    push!(parms, fn_name)  # function name

    for call_arg in call.args[2:end]
        @match call_arg begin
            Expr(:(::), [sym, Expr(:curly, [:Ptr, :UInt8], _)], _) =>
                begin
                    orig_type = Expr(:curly, :Ptr, :UInt8)
                    _sym = Symbol(string("_", sym))
                    push!(parms,   :($_sym::Union{Ptr,ByteString}))
                    push!(content, :($sym = convert($orig_type, $_sym)))
                end
            Expr(:(::), [sym, Expr(:curly, [:Ptr, target_type], _)], _) =>
                begin
                    orig_type = Expr(:curly, :Ptr, target_type)
                    _sym = Symbol(string("_", sym))
                    push!(parms,   :($_sym::Ptr))
                    push!(content, :($sym = convert($orig_type, $_sym)))
                end
            Expr(:(::), [sym, (:Int32 || :UInt32 || :Int64 || :UInt64 || :Cuint || :Cint || :Csize_t)], _) =>
                begin
                    println(sym)
                    push!(parms, :($sym::Integer))
                end
            _ => push!(parms, call_arg)
        end
    end

    if !isexpr(body, :block)
        push!(content, body)
    else
        append!(content, body.args)
    end

    call = Expr(:call, parms...)
    body = Expr(:block, content...)

    new = Expr(e.head, call, body)

    return new
end

out_path = pwd()

context = wrap_c.init(
    output_file = "avbin.jl",
    common_file = "avbin_h.jl",
    clang_includes = clang_includes,
    header_wrapped = check_use_header,
    header_library = get_header_library,
    header_outputfile = get_header_outfile,
    func_rewriter = rewrite_fn)
wrap_c.wrap_c_headers(context, av_headers)

ext: jl

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
module Radial

export R_vals, set_d1!, set_d2!, filter!

"""
Computes array of compactified radial points

R_vals(
   nx::Int64,
   dr::Float64
   )::Vector{Float64}
"""
function R_vals(nx::Int64, dr::Float64)::Vector{Float64}
   return [dr*i for i=1:nx]
end

"""
Takes first derivative in radial direction.

set_d1!(
   dv,
   v,
   dr::Float64,
   )
"""
function set_d1!(
      dv,
      v,
      dr::Float64,
     )
   inv_dr   = 1.0/dr

   nx, ny = size(dv)

   for j=1:ny
      for i=3:nx-2
         dv[i,j] = (-
                    (1.0/12.0)*v[i+2,j]
                    +
                    (2.0/3.0 )*v[i+1,j]
                    -
                    (2.0/3.0 )*v[i-1,j]
                    +
                    (1.0/12.0)*v[i-2,j]
                   )*inv_dr
      end
      dv[1,j] = (-
                 3.0 *v[5,j]
                 +
                 16.0*v[4,j]
                 -
                 36.0*v[3,j]
                 +
                 48.0*v[2,j]
                 -
                 25.0*v[1,j]
                )*(1.0/12.0)*inv_dr
      dv[2,j] = (+
                 1.0 *v[5,j]
                 -
                 6.0 *v[4,j]
                 +
                 18.0*v[3,j]
                 -
                 10.0*v[2,j]
                 -
                 3.0 *v[1,j]
                )*(1.0/12.0)*inv_dr

      dv[nx,j] = -(-
                   3.0 *v[nx-4,j]
                   +
                   16.0*v[nx-3,j]
                   -
                   36.0*v[nx-2,j]
                   +
                   48.0*v[nx-1,j]
                   -
                   25.0*v[nx,j]
                  )*(1.0/12.0)*inv_dr
      dv[nx-1,j] = -(+
                     1.0 *v[nx-4,j]
                     -
                     6.0 *v[nx-3,j]
                     +
                     18.0*v[nx-2,j]
                     -
                     10.0*v[nx-1,j]
                     -
                     3.0 *v[nx,j]
                    )*(1.0/12.0)*inv_dr
   end
   return nothing
end
"""
Takes second derivative in radial direction

set_d2!(
   dv,
   v,
   dr::Float64
   )
"""
function set_d2!(
      dv,
      v,
      dr::Float64
     )
   inv_dr2 = 1.0/dr^2

   nx, ny = size(dv)

   for j=1:ny
      for i=3:nx-2
         dv[i,j] = (-
                    1.0 *v[i+2,j]
                    +
                    16.0*v[i+1,j]
                    -
                    30.0*v[i  ,j]
                    +
                    16.0*v[i-1,j]
                    -
                    1.0 *v[i-2,j]
                   )*(1.0/12.0)*inv_dr2
      end
      dv[1,j] = (-
                 10.0 *v[6,j]
                 +
                 61.0 *v[5,j]
                 -
                 156.0*v[4,j]
                 +
                 214.0*v[3,j]
                 -
                 154.0*v[2,j]
                 +
                 45.0 *v[1,j]
                )*(1.0/12.0)*inv_dr2
      dv[2,j] = (+
                 1.0 *v[6,j]
                 -
                 6.0 *v[5,j]
                 +
                 14.0*v[4,j]
                 -
                 4.0 *v[3,j]
                 -
                 15.0*v[2,j]
                 +
                 10.0*v[1,j]
                )*(1.0/12.0)*inv_dr2
      dv[nx,j] = (-
                  10.0 *v[nx-5,j]
                  +
                  61.0 *v[nx-4,j]
                  -
                  156.0*v[nx-3,j]
                  +
                  214.0*v[nx-2,j]
                  -
                  154.0*v[nx-1,j]
                  +
                  45.0 *v[nx  ,j]
                 )*(1.0/12.0)*inv_dr2
      dv[nx-1,j] = (+
                    1.0 *v[nx-5,j]
                    -
                    6.0 *v[nx-4,j]
                    +
                    14.0*v[nx-3,j]
                    -
                    4.0 *v[nx-2,j]
                    -
                    15.0*v[nx-1,j]
                    +
                    10.0*v[nx  ,j]
                   )*(1.0/12.0)*inv_dr2
   end
   return nothing
end

"""
6th order Kreiss-Oliger dissipation
"""
function filter!(
      v,
      tmp,
      eps_KO::Float64
     )
   nx, ny = size(v)

   for j=1:ny
      for i=1:nx
         tmp[i,j] = v[i,j]
      end
   end

   for j=1:ny
      for i=4:nx-3
         v[i,j] = (
                         (eps_KO/64.) *tmp[i-3,j]
           +        (-6.0*eps_KO/64.) *tmp[i-2,j]
           +        (15.0*eps_KO/64.) *tmp[i-1,j]
           + (1.0 - (20.0*eps_KO/64.))*tmp[i  ,j]
           +        (15.0*eps_KO/64.) *tmp[i+1,j]
           +        (-6.0*eps_KO/64.) *tmp[i+2,j]
           +             (eps_KO/64.) *tmp[i+3,j]
          )
      end
   end

   return nothing
end

end


text:
module PLCTag
	module libplctag
		using CBinding

		let
			incdir = joinpath(dirname(@__DIR__), "deps/usr/include")
			libdir = joinpath(dirname(@__DIR__), "deps/usr/lib")

			c`-I$(incdir) -fparse-all-comments -L$(libdir) -lplctag`
		end

		const c"int8_t"  = Int8
		const c"int16_t" = Int16
		const c"int32_t" = Int32
		const c"int64_t" = Int64
		const c"uint8_t"  = UInt8
		const c"uint16_t" = UInt16
		const c"uint32_t" = UInt32
		const c"uint64_t" = UInt64

		c"""
			#include <libplctag.h>
		"""ji
	end

	using CBinding
	using .libplctag

	export PLC, PLCRef, PLCBinding

	Base.@kwdef struct PLC
		cpu::String
		protocol::String
		gateway::String
		path::String
		debug::Int = 1
	end

	struct PLCRef{T}
		tag::Ref{Int32}
		count::Int

		function PLCRef{T}(plc::PLC, tagName::String, count::Int = 1; timeout::Int = 1000) where {T}
			isprimitivetype(T) || error("Unable to create a PLCRef for anything other than primitive types")

			path = plcpath(plc, tagName, sizeof(T), count)

			tag = plc_tag_create(path, timeout)
			tag <= 0 && error("Unable to create tag `$(tagName)` of type `$(T)` and count $(count): $(unsafe_string(plc_tag_decode_error(tag)))")

			result = new{T}(Ref{Int32}(tag), count)
			finalizer(t -> plc_tag_destroy(t[]), result.tag)
			return result
		end
	end

	Base.length(ref::PLCRef) = ref.count
	Base.size(ref::PLCRef) = (ref.count,)
	Base.eltype(ref::PLCRef{T}) where {T} = T
	Base.iterate(ref::PLCRef, state = 1) = state > length(ref) ? nothing : (ref[state], state+1)

	function Base.read(ref::PLCRef)
		code = plc_tag_read(ref.tag[], 0)
		code in (c"PLCTAG_STATUS_OK", c"PLCTAG_STATUS_PENDING") || error("Unable to read tag: $(unsafe_string(plc_tag_decode_error(code)))")
	end
	function Base.fetch(ref::PLCRef, ind::Int = 1)
		timedwait(() -> plc_tag_status(ref.tag[]) == c"PLCTAG_STATUS_OK", 1.0, pollint = 0.001) === :ok || error("Failed to complete tag read: $(unsafe_string(plc_tag_decode_error(plc_tag_status(ref.tag[]))))")
		return convert(eltype(ref), plcget(eltype(ref), ref.tag[], ind))
	end

	function Base.write(ref::PLCRef, val, ind::Int = 1)
		plcset(eltype(ref), ref.tag[], ind, convert(eltype(ref), val))
		code = plc_tag_write(ref.tag[], 0)
		code in (c"PLCTAG_STATUS_OK", c"PLCTAG_STATUS_PENDING") || error("Unable to write tag: $(unsafe_string(plc_tag_decode_error(code)))")
	end
	function Base.flush(ref::PLCRef)
		timedwait(() -> plc_tag_status(ref.tag[]) == c"PLCTAG_STATUS_OK", 1.0, pollint = 0.001) === :ok || error("Failed to complete tag write: $(unsafe_string(plc_tag_decode_error(plc_tag_status(ref.tag[]))))")
	end

	function Base.getindex(ref::PLCRef, ind::Int = 1)
		read(ref)
		return fetch(ref, ind)
	end

	function Base.setindex!(ref::PLCRef, val, ind::Int = 1)
		write(ref, val, ind)
		flush(ref)
	end

	plcpath(plc::PLC, tagName::String, size::Int, count::Int) = join(map(((k, v),) -> "$(k)=$(v)", (
		:cpu => plc.cpu,
		:protocol => plc.protocol,
		:gateway => plc.gateway,
		:path => plc.path,
		:debug => plc.debug,
		:elem_size => size,
		:elem_count => count,
		:name => tagName,
	)), '&')

	const PLCPrimitives = Union{
		Bool,
		Int8, UInt8, Int16, UInt16,
		Int32, UInt32, Int64, UInt64,
		Float32, Float64
	}

	for T in (
		:Bool,
		:Int8, :UInt8, :Int16, :UInt16,
		:Int32, :UInt32, :Int64, :UInt64,
		:Float32, :Float64,
	)
		if T === :Bool
			# NOTE: using Int8 (aka SINT on the PLC) to store a Bool
			@eval plcget(::Type{Bool}, tag::Int32, ind::Int)            = plcget(Int8, tag, ind) != 0
			@eval plcset(::Type{Bool}, tag::Int32, ind::Int, val::Bool) = plcset(Int8, tag, ind, convert(Int8, val))
		else
			@eval plcget(::Type{$(T)}, tag::Int32, ind::Int)            = $(Symbol("plc_tag_get_", lowercase(String(T))))(tag, ind-1)
			@eval plcset(::Type{$(T)}, tag::Int32, ind::Int, val::$(T)) = $(Symbol("plc_tag_set_", lowercase(String(T))))(tag, ind-1, val)
		end

		@eval Base.$(T)(plc::PLC, tagName::String) = PLCRef{$(T)}(plc, tagName)[]
	end

	function Base.Vector{T}(plc::PLC, tagName::String, count::Int) where {T<:PLCPrimitives}
		ref = PLCRef{T}(plc, tagName, count)
		return T[ref[i] for i in 1:length(ref)]
	end

	function Base.String(plc::PLC, tagName::String)
		len = PLCRef{Int32}(plc, "$(tagName).LEN")[]
		data = PLCRef{UInt8}(plc, "$(tagName).DATA", len)
		return String(data)
	end
end


text:
module DelayedKalmanFilter

using DelimitedFiles
using DifferentialEquations
using Distributions
using ForwardDiff
using JLD
using LatinHypercubeSampling
using LinearAlgebra

include("kalman_filter_alg.jl")
include("log_likelihood.jl")
include("hes5_ode_steady_state.jl")
include("mcmc_samplers.jl")

export kalman_filter
export calculate_log_likelihood_at_parameter_point, calculate_log_likelihood_and_derivative_at_parameter_point
export log_likelihood_and_derivative_with_prior_and_transformation
export calculate_steady_state_of_ode
export MALA, run_mala_for_dataset

end # module


text:
function model_step!(model::Agents.ABM)
    model.new_trash = sum(model.trash_in_spaces) - model.trash_aggregated
    collect_garbage!(model)
    model.trash_aggregated = sum(model.trash_in_spaces)
    return model
end

function collect_garbage!(model::Agents.ABM)
    model.trash_in_spaces = model.trash_in_spaces .- model.garbage_collection
    model.trash_in_spaces = [i >= 0 ? i : 0 for i in model.trash_in_spaces]
    return model
end


text:
using YaoBlocks, YaoArrayRegister
using Test

@testset "check apply" begin
    r = rand_state(1)
    channel = UnitaryChannel([X, Y, Z], [1, 0, 0])
    @test apply!(copy(r), channel) == apply!(copy(r), X)

    r = rand_state(3)
    @test apply!(copy(r), put(3, 1 => channel)) == apply!(copy(r), put(3, 1 => X))
end

@testset "check mat" begin
    @test mat(UnitaryChannel([X, Y, Z], [1, 0, 0])) == mat(X)
    @test mat(UnitaryChannel([X, Y, Z], [0, 1, 0])) == mat(Y)
    @test mat(UnitaryChannel([X, Y, Z], [0, 0, 1])) == mat(Z)
end

@testset "check compare" begin
    @test UnitaryChannel([X, Y, Z], [1, 0, 0]) == UnitaryChannel([X, Y, Z], [1, 0, 0])
    @test UnitaryChannel([X, Y, Z], [1, 0, 0]) != UnitaryChannel([X, Y, Z], [1, 1, 0])
    @test UnitaryChannel([X, Y, Z], [1, 0, 0]) != UnitaryChannel([X, Y], [1, 1, 0])
end

@testset "check adjoint" begin
    channel = UnitaryChannel([X, Y, Z])
    adj_channel = adjoint(channel)
    @test adjoint.(adj_channel.operators) == channel.operators
end

@testset "check ocuppied locations" begin
    channel = UnitaryChannel(put.(6, [1 => X, 3 => Y, 5 => Z]))
    @test occupied_locs(channel) == [1, 3, 5]
end


text:
# concrete types -
struct VLResult{T}
    value::T
end

max_stars_repo_name: varnerlab/PregnancyCoagDataAnalysisKit

max_issues_repo_name: varnerlab/PregnancyCoagDataAnalysisKit

max_forks_repo_name: varnerlab/PregnancyCoagDataAnalysisKit


text:
module ReachabilityBenchmarks

# ========
# Exports
# ========
export @relpath

# ================================
# Input/Output functionality
# ================================
include("io.jl")

end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


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

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

# For a moment, recall the Navier–Stokes equations for an incompressible fluid, where $\vec{v}$ represents the velocity field:

# $$
# \begin{eqnarray*}
# \nabla \cdot\vec{v} &=& 0 \\
# \frac{\partial \vec{v}}{\partial t}+(\vec{v}\cdot\nabla)\vec{v} &=& -\frac{1}{\rho}\nabla p + \nu \nabla^2\vec{v}
# \end{eqnarray*}
# $$

# The first equation represents mass conservation at constant density. The second equation is the conservation of momentum. But a problem appears: the continuity equation for incompressble flow does not have a dominant variable and there is no obvious way to couple the velocity and the pressure. In the case of compressible flow, in contrast, mass continuity would provide an evolution equation for the density $\rho$, which is coupled with an equation of state relating $\rho$ and $p$.

# In incompressible flow, the continuity equation $\nabla \cdot\vec{v}=0$ provides a *kinematic constraint* that requires the pressure field to evolve so that the rate of expansion $\nabla \cdot\vec{v}$ should vanish everywhere. A way out of this difficulty is to *construct* a pressure field that guarantees continuity is satisfied; such a relation can be obtained by taking the divergence of the momentum equation. In that process, a Poisson equation for the pressure shows up!

# Step 10: 2D Poisson Equation
# ----
# ***

# Poisson's equation is obtained from adding a source term to the right-hand-side of Laplace's equation:

# $$\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = b$$

# So, unlike the Laplace equation, there is some finite value inside the field that affects the solution. Poisson's equation acts to "relax" the initial sources in the field.

# In discretized form, this looks almost the same as [Step 9](./12_Step_9.ipynb), except for the source term:

# $$\frac{p_{i+1,j}^{n}-2p_{i,j}^{n}+p_{i-1,j}^{n}}{\Delta x^2}+\frac{p_{i,j+1}^{n}-2 p_{i,j}^{n}+p_{i,j-1}^{n}}{\Delta y^2}=b_{i,j}^{n}$$

# As before, we rearrange this so that we obtain an equation for $p$ at point $i,j$. Thus, we obtain:

# $$p_{i,j}^{n}=\frac{(p_{i+1,j}^{n}+p_{i-1,j}^{n})\Delta y^2+(p_{i,j+1}^{n}+p_{i,j-1}^{n})\Delta x^2-b_{i,j}^{n}\Delta x^2\Delta y^2}{2(\Delta x^2+\Delta y^2)}$$

# We will solve this equation by assuming an initial state of $p=0$ everywhere, and applying boundary conditions as follows:

# $p=0$ at $x=0, \ 2$ and $y=0, \ 1$

# and the source term consists of two initial spikes inside the domain, as follows:

# $b_{i,j}=100$ at $i=\frac{1}{4}nx, j=\frac{1}{4}ny$

# $b_{i,j}=-100$ at $i=\frac{3}{4}nx, j=\frac{3}{4}ny$

# $b_{i,j}=0$ everywhere else.

# The iterations will advance in pseudo-time to relax the initial spikes. The relaxation under Poisson's equation gets slower and slower as they progress. *Why?*

# Let's look at one possible way to write the code for Poisson's equation. As always, we load our favorite Julia libraries.
# We also want to make some lovely plots in 3D.
# Let's get our parameters defined and the initialization out of the way.
# What do you notice of the approach below?

using Plots

# Parameters
nx = 50
ny = 50
nt  = 100
xmin = 0
xmax = 2
ymin = 0
ymax = 1

dx = (xmax - xmin) / (nx - 1)
dy = (ymax - ymin) / (ny - 1)

# Initialization
p  = zeros((ny, nx))
pd = zeros((ny, nx))
b  = zeros((ny, nx))
x  = LinRange(xmin, xmax, nx)
y  = LinRange(xmin, xmax, ny)

# Source
b[int(ny / 4), int(nx / 4)]  = 100
b[int(3 * ny / 4), int(3 * nx / 4)] = -100

# With that, we are ready to advance the initial guess in pseudo-time. How is the code below different from the function used in [Step 9](./12_Step_9.jl) to solve Laplace's equation?

for it in range(nt)

    pd = copy(p)

    p[2:end,2:end] = (((pd[2:end, 2:] + pd[2:end, :-2]) * dy^2 +
                    (pd[2:, 2:end] + pd[:-2, 2:end]) * dx^2 -
                    b[2:end, 2:end] * dx^2 * dy^2) /
                    (2 * (dx^2 + dy^2)))

    # TODO check indexing
    p[begin, :] = 0
    p[ny-1,  :] = 0
    p[:, begin] = 0
    p[:, nx-1] = 0
end

# TODO PLOT
function plot2D(x, y, p)
    fig = pyplot.figure(figsize=(11, 7), dpi=100)
    ax = fig.gca(projection='3d')
    X, Y = numpy.meshgrid(x, y)
    surf = ax.plot_surface(X, Y, p[:], rstride=1, cstride=1, cmap=cm.viridis,
            linewidth=0, antialiased=False)
    ax.view_init(30, 225)
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')

    fig
end
# TODO Plots
plot2D(x, y, p)

# Ah! The wonders of code reuse! Now, you probably think: "Well, if I've written this neat little function that does something so useful, I want to use it over and over again. How can I do this without copying and pasting it each time? —If you are very curious about this, you'll have to learn about *packaging*. But this goes beyond the scope of our CFD lessons. You'll just have to Google it if you really want to know.

# ***

## Learn More

# To learn more about the role of the Poisson equation in CFD, watch **Video Lesson 11** on You Tube:

#src="https://www.youtube.com/embed/ZjfxA3qq2Lg"

max_stars_repo_name: miguelraz/CFDPython

max_issues_repo_name: miguelraz/CFDJulia

max_forks_repo_name: miguelraz/CFDJulia


text:
using Pkg
Pkg.activate(".")

using IJulia
notebook(dir=".")


text:
using Test

using Dice
using Dice: CuddMgr, DistBool

code = @dice begin
    function test_rec(a::Int)
        if a <= 0
            flip(0.5)
        else
            flip(0.5) & test_rec(a - 1)
        end
    end

    test_rec(2)
end

@testset "Sugar Tests" begin

    @test infer((@dice if flip(0.5) flip(0.1) else flip(0.3) end), :bdd) ≈ 0.2
    @test infer((@dice if true flip(0.1) else flip(0.2) end), :bdd) ≈ 0.1

    @test infer(code, :bdd) ≈ 0.125

end


text:
using Dispersal, Test, Unitful, Dates, DimensionalData
using Unitful: d, s

struct TestFormulation <: KernelFormulation end
(f::TestFormulation)(d) = 1.0

@testset "dispersal kernel array matches passed in function" begin
    init = [0.0 1.0 0.0;
            1.0 0.0 1.0]
    radius = 2
    dk = DispersalKernel{radius}(formulation=ExponentialKernel(1.0), kernel=init, cellsize=1.0).kernel
    @test size(dk) == (5, 5)
    @test sum(dk) ≈ 1.0
end

@testset "binary dispersal and growth mask" begin
    init = [0 0; 0 1]
    radius = 1
    hood = DispersalKernel{radius}()

    # time sequence of layers
    a = cat([0.1 0.2; 0.3 0.4], [0.5 0.6; 0.7 0.8], dims=3)
    dimz = X(1:2), Y(1:2), Ti(1d:10d:11d)
    suitseq = DimensionalArray(a, dimz)

    # Regular
    ruleset1 = Ruleset(InwardsBinaryDispersal(neighborhood=hood, prob_threshold=0.0),
                       MaskGrowthMap(layerkey=Val(:suitseq), threshold=0.4); timestep=1d)
    # Chained
    ruleset2 = Ruleset(Chain(InwardsBinaryDispersal(neighborhood=hood, prob_threshold=0.0),
                       MaskGrowthMap(layerkey=Val(:suitseq), threshold=0.4)) ; timestep=1d)
    output1 = ArrayOutput(init; tspan=1d:1d:25d)
    output2 = ArrayOutput(init; tspan=1d:1d:25d)

    sim!(output1, ruleset1; aux=(suitseq=suitseq,))
    sim!(output2, ruleset2; aux=(suitseq=suitseq,))

    results = [[0 0; 0 1],
               [0 0; 0 1],
               [0 0; 0 1],
               [0 0; 0 1],
               [0 0; 0 1],
               [1 1; 1 1],
               [1 1; 1 1],
               [0 0; 0 1],
               [0 0; 0 1],]

    @testset "Rounded" begin
        @test output1[[1, 2, 5, 8, 10, 15, 20, 22, 25]] == results
        @test output2[[1, 2, 5, 8, 10, 15, 20, 22, 25]] == results
        output1[[1, 2, 5, 8, 10, 15, 20, 22, 25]]
        output2[[1, 2, 5, 8, 10, 15, 20, 22, 25]]
    end

    # @testset "Interpolated" begin
    #     # All offset by one, because 1 = t0
    #     @test output1[1]  == output2[1]  == [0 0; 0 1]
    #     @test output1[2]  == output2[2]  == [0 0; 1 1]
    #     @test output1[5]  == output2[5]  == [0 0; 0 1]
    #     @test output1[8]  == output2[8]  == [0 0; 1 1]
    #     @test output1[10] == output2[10] == [0 1; 1 1]
    #     @test output1[15] == output2[15] == [1 1; 1 1]
    #     @test output1[20] == output2[20] == [0 1; 1 1]
    #     @test output1[22] == output2[22] == [0 0; 1 1]
    #     @test output1[25] == output2[25] == [0 0; 0 1]
    # end

    @test_throws BoundsError output1[26]
    @test_throws BoundsError output2[26]
end

@testset "binary dispersal simulation with suitability mask" begin

    suit =  [1 0 1 1 0;
             0 0 1 1 1;
             1 1 1 1 0;
             1 1 0 1 1;
             1 0 1 1 1]

    init =  Bool[0 0 0 0 0;
                 0 0 0 0 0;
                 0 0 1 0 0;
                 0 0 0 0 0;
                 0 0 0 0 0]

    test1 = Bool[0 0 0 0 0;
                 0 0 0 0 0;
                 0 0 1 0 0;
                 0 0 0 0 0;
                 0 0 0 0 0]

    test2 = Bool[0 0 0 0 0;
                 0 0 1 1 0;
                 0 1 1 1 0;
                 0 1 0 1 0;
                 0 0 0 0 0]

    test3 = Bool[0 0 1 1 0;
                 0 0 1 1 1;
                 1 1 1 1 0;
                 1 1 0 1 1;
                 1 0 1 1 1]

    # Dispersal in radius 1 neighborhood
    maskgrowth = MaskGrowthMap(layerkey=:suit)
    radius = 1
    hood = DispersalKernel{radius}(; formulation=ExponentialKernel(1.0))

    @testset "inwards binary dispersal fills the grid where reachable and suitable" begin
        inwards = InwardsBinaryDispersal(neighborhood=hood, prob_threshold=0.0)
        rules = Ruleset(inwards, maskgrowth)
        output = ArrayOutput(init; tspan=1:3)
        sim!(output, rules; aux=(suit=suit,))
        @test output[1] == test1
        @test output[2] == test2
        @test output[3] == test3

        # As subrules
        rules = Ruleset(inwards, maskgrowth)
        sim!(output, rules; aux=(suit=suit,))
        @test output[1] == test1
        @test output[2] == test2
        @test output[3] == test3
    end
    @testset "outwards binary dispersal fills the grid where reachable and suitable" begin
        outwards = OutwardsBinaryDispersal(neighborhood=hood, prob_threshold=0.0)
        rules = Ruleset(outwards, maskgrowth)
        output = ArrayOutput(init; tspan=1:3)
        sim!(output, rules, aux=(suit=suit,))
        @test output[1] == test1
        @test output[2] == test2
        @test output[3] == test3
    end

end

@testset "floating point population dispersal simulation with suitability mask" begin

    suit =  [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
             0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
             0.0 0.0 1.0 0.0 1.0 0.0 1.0 1.0 0.0 0.0 0.0;
             0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0 1.0 0.0 0.0;
             0.0 0.0 1.0 0.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0;
             0.0 0.0 1.0 0.0 1.0 1.0 0.0 1.0 1.0 0.0 0.0;
             0.0 0.0 1.0 1.0 1.0 0.0 1.0 1.0 1.0 0.0 0.0;
             0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
             0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0;
             0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
             0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;]

    init =  [0.0 0.0   0.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;
             0.0 0.0   0.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;
             0.0 0.0 100.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;
             0.0 0.0   0.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;
             0.0 0.0   0.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 0.0;
             0.0 0.0   0.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;
             0.0 0.0   0.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;
             0.0 0.0   0.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;
             0.0 0.0   0.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;
             0.0 0.0   0.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;
             0.0 0.0   0.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;]

    test1 = [0.0 0.0   0.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;
             0.0 0.0   0.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;
             0.0 0.0 100.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;
             0.0 0.0   0.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;
             0.0 0.0   0.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 0.0;
             0.0 0.0   0.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;
             0.0 0.0   0.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;
             0.0 0.0   0.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;
             0.0 0.0   0.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;
             0.0 0.0   0.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;
             0.0 0.0   0.0 0.0 0.0 0.0   0.0 0.0 0.0 0.0 0.0;]

    test2 = [0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0;
             0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0;
             0.0  0.0  4.0  0.0  8.0  0.0  4.0  4.0  0.0  0.0  0.0;
             0.0  0.0  4.0  4.0  0.0  0.0  4.0  4.0  4.0  0.0  0.0;
             0.0  0.0  4.0  0.0  8.0  4.0  4.0  4.0  0.0  0.0  0.0;
             0.0  0.0  0.0  0.0  4.0  4.0  0.0  4.0  4.0  0.0  0.0;
             0.0  0.0  0.0  0.0  4.0  0.0  4.0  4.0  4.0  0.0  0.0;
             0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0;
             0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0;
             0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0;
             0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0;]

    test3 = [0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0;
             0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0;
             0.0  0.0  1.28 0.0  1.92 0.0  1.92 1.28 0.0  0.0  0.0;
             0.0  0.0  1.44 1.76 0.0  0.0  2.56 1.76 1.44 0.0  0.0;
             0.0  0.0  1.6  0.0  2.56 2.88 3.2  2.24 0.0  0.0  0.0;
             0.0  0.0  1.12 0.0  1.92 2.24 0.0  1.92 1.6  0.0  0.0;
             0.0  0.0  0.8  1.12 1.44 0.0  2.08 1.44 1.12 0.0  0.0;
             0.0  0.0  0.32 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0;
             0.0  0.0  0.16 0.16 0.32 0.48 0.64 0.48 0.48 0.0  0.0;
             0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0;
             0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0;]

    # Dispersal in radius 1 neighborhood
    mask = MaskGrowthMap(layerkey=:suit)
    radius = 2

    @testset "inwards population dispersal fills the grid where reachable and suitable" begin
        hood = DispersalKernel{radius}(; formulation=TestFormulation(), distancemethod=CentroidToCentroid())
        inwards = InwardsPopulationDispersal(neighborhood=hood)
        rules = Ruleset(inwards, mask; timestep=1d)
        in_output = ArrayOutput(init; tspan=4d:1d:6d, aux=(suit=suit,))
        sim!(in_output, rules)
        @test in_output[1] == test1
        @test in_output[2] == test2
        @test in_output[3] ≈ test3
        # As subrules
        rules = Ruleset(Chain(inwards, mask); timestep=1d)
        sim!(in_output, rules; tspan=4d:1d:6d)
        @test in_output[1] == test1
        @test in_output[2] == test2
        @test in_output[3] ≈ test3
    end

    @testset "outwards population dispersal fills the grid where reachable and suitable" begin
        hood = DispersalKernel{radius}(; formulation=TestFormulation(), distancemethod=CentroidToCentroid())
        outwards = OutwardsPopulationDispersal(neighborhood=hood)
        rules = Ruleset(outwards, mask; timestep=Month(1))
        out_output = ArrayOutput(init; tspan=Date(2001, 1):Month(1):Date(2001, 3))
        sim!(out_output, rules; aux=(suit=suit,))
        @test out_output[1] == test1
        @test out_output[2] == test2
        @test out_output[3] ≈ test3
    end

end


text:
"""
    stp_to_pss(stp_table, det_meta, env, ps_simulator, noise_model, config_name)

Table, PropDict, Env, PSSimulator, NoiseModel, AbstractString -> Table, Table

Simulate waveforms based on stepping info given in <stp_table>,
    LEGEND detector metadata <det_meta>, environment settings given in <env>,
    pulse shape simulation method <ps_simulator>, and <noise_model>

The output is a table with simulated pulses, and a table with simulation truth
    (may be abolished in the future since basically corresponds to stp table)
"""
function stp_to_pss(stp_table::Table, det_meta::PropDict, env::Environment, ps_simulator::PSSimulator, noise_model::NoiseModel,
    cached_name::AbstractString)

    # add fano noise, don't add if data noise is applied later
    # note: current understanding is that Siggen simulation does NOT include fano noise
    @info "//\\//\\//\\ Fano noise"
    stp_table_fano = fano_noise(stp_table, det_meta, env, noise_model)

    @info "_||_||_||_ Simulate detector"
    detector = simulate_detector(det_meta, env, cached_name, ps_simulator)

    @info "~.~.~.~.~ Simulate charge pulses"
    simulate_waveforms(stp_table_fano, detector)
end

function stp_to_pss(stp_filepath::AbstractString, det_meta_fullpath::AbstractString, sim_config_filename::AbstractString)
    stp_table = h5open(stp_filepath, "r") do input
        LegendHDF5IO.readdata(input, "stp")
    end

    stp_to_pss(stp_table, det_meta_fullpath, sim_config_filename)
end

function stp_to_pss(stp_table::Table, det_meta_fullpath::AbstractString, sim_config_filename::AbstractString)
    # construct simulation config based on given inputs and settings
    sim_config = load_config(stp_table, det_meta_fullpath, sim_config_filename)

    stp_to_pss(sim_config)
end

function stp_tp_pss(sim_config_filename::AbstractString)
    # given config file already contains inputs as well
    sim_config = propdict(sim_config_filename)
    # ok this looks kinda stupid i'm in a multiple dispatch loop what do i do
    stp_to_pss(sim_config.input_file, sim_config.detector_metadata, sim_config_filename)
end

function stp_to_pss(sim_config::PropDict)
    @info "---------------------- stp -> pss (pulse shape simulation)"

    det_meta = propdict(sim_config.detector_metadata)
    env = Environment(sim_config)
    ps_simulator = PSSimulator(sim_config)
    noise_model = NoiseModel(sim_config)

    stp_to_pss(sim_config.input_file, det_meta, env, ps_simulator, noise_model, sim_config.pss.cached_name)
end


text:
using IndividualDisplacements, MeshArrays, OrdinaryDiffEq
p=dirname(pathof(MeshArrays))
include(joinpath(p,"../examples/Demos.jl"))

"""
    simple_flow_field(Γ::Dict,np,nz)

Set up an idealized flow field which consists of
[rigid body rotation](https://en.wikipedia.org/wiki/Rigid_body),
plus a convergent term, plus a sinking term.

```
Γ=simple_periodic_domain(12)
u,v,w=simple_flow_field(Γ,12,4)
```
"""
function simple_flow_field(Γ,np,nz)
    γ=Γ["XC"].grid;

    #Solid-body rotation around central location ...
    i=Int(np/2+1)
    u=-(Γ["YG"].-Γ["YG"][1][i,i])
    v=(Γ["XG"].-Γ["XG"][1][i,i])

    #... plus a convergent term to / from central location
    d=-0.01
    u=u+d*(Γ["XG"].-Γ["XG"][1][i,i])
    v=v+d*(Γ["YG"].-Γ["YG"][1][i,i])

    #Replicate u,v in vertical dimension
    uu=MeshArray(γ,γ.ioPrec,nz)
    [uu[k]=u[1] for k=1:nz]
    vv=MeshArray(γ,γ.ioPrec,nz)
    [vv[k]=v[1] for k=1:nz]

    #Vertical velocity component w
    w=fill(-0.01,MeshArray(γ,γ.ioPrec,nz));

    return uu,vv,w
end

"""
    setup_random_flow(;np=12,nq=18)

Set up a random flow field over a gridded domain of size np,nq

```
ϕ,u,v=setup_random_flow()
```
"""
function setup_random_flow(;np=12,nq=18)
    Γ=simple_periodic_domain(np,nq)
    (_,ϕ,_,_)=demo2(Γ)
    ϕ .*= 0.5

    #For the convergent / scalar potential case, ϕ is interpreted as being
    #on center points -- hence the standard gradient function readily gives
    #what we need
    #(u,v)=gradient(ϕ,Γ)
    #return u[1],v[1],ϕ[1]

    #For the rotational / streamfunction case, ϕ is interpreted as being
    #on S/W corner points -- this is ok since the grid is homegeneous,
    #and conveniently yields an adequate substitution u,v <- -v,u; but note
    #that doing the same with gradient() would shift indices inconsistenly
    u=-(circshift(ϕ[1], (0,-1))-ϕ[1])
    v=(circshift(ϕ[1], (-1,0))-ϕ[1])
    return u,v,ϕ[1]
end

function setup_point_cloud(U::Array{T,2},V::Array{T,2};X=[],Y=[]) where T
    np,nq=size(U)
    Γ=simple_periodic_domain(np,nq)
    u=MeshArray(Γ["XC"].grid,[U])
    v=MeshArray(Γ["XC"].grid,[V])
    #vel=dxy_dt
    (u,v)=exchange(u,v,1)
    vel=dxy_dt!

    𝑃 = (u0=u, u1=u, v0=v, v1=v, 𝑇=[0.0,1.0], ioSize=Γ["XC"].grid.ioSize)
    pp=postprocess_xy
    isempty(X) ? X=np*rand(10) : nothing
    isempty(Y) ? Y=nq*rand(10) : nothing

    xy = permutedims([[X[i];Y[i];1.0] for i in eachindex(X)])
    tr = DataFrame(ID=Int[], x=Float64[], y=Float64[], t=Float64[])
    solv(prob) = solve(prob,Tsit5(),reltol=1e-5,abstol=1e-5)

    I=(position=xy,record=tr,velocity=vel,
       integration=solv,postprocessing=pp,parameters=𝑃)

    return Individuals(I)
end

"""
    setup_global_ocean(;k=10,ny=2)

Set up Global Ocean particle simulation in 2D with seasonally varying flow field.

```
𝑃=setup_global_ocean(k=1,ny=2);
```
"""
function setup_global_ocean(;k=1,ny=2)

  #k=10 #choice of vertical level
  #ny=2 #number of simulated years (20 for k>20)
  r_reset = 0.01 #fraction of the particles reset per month (0.05 for k<=10)

  #read grid and set up connections between subdomains
  p=dirname(pathof(IndividualDisplacements))
  γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
  Γ=GridLoad(γ)
  Γ=merge(Γ,IndividualDisplacements.NeighborTileIndices_cs(Γ))

  #initialize u0,u1 etc
  𝑃=set_up_𝑃(k,0.0,Γ,ECCOclim_path);

  #add parameters for use in reset!
  tmp=(frac=r_reset, Γ=Γ)
  𝑃=merge(𝑃,tmp)

  return 𝑃

end

##

"""
    example3_setup(;backward_in_time::Bool=false)

Define gridded variables and return result as Dictionary (`uvetc`).
"""
function example3_setup(;backward_in_time::Bool=false)

   p=dirname(pathof(IndividualDisplacements))
   dirIn=joinpath(p,"../examples/llc90_latlon/")
   γ=gcmgrid(dirIn,"PeriodicChannel",1,
                  [(360,178)], [360 178], Float32, read, write)

   Γ=Dict("XC" => read(γ.path*"XC.latlon.data",MeshArray(γ,Float32)),
   "YC" => read(γ.path*"YC.latlon.data",MeshArray(γ,Float32)),
   "XG" => read(γ.path*"XG.data",MeshArray(γ,Float32)),
   "YG" => read(γ.path*"YG.data",MeshArray(γ,Float32)),
   "DXC" => read(γ.path*"DXC.latlon.data",MeshArray(γ,Float32)),
   "DYC" => read(γ.path*"DYC.latlon.data",MeshArray(γ,Float32)) );

   file = matopen(γ.path*"uv_lonlat.mat")
   u=read(file, "u")
   v=read(file, "v")
   close(file)

   u=dropdims(mean(u,dims=3),dims=3)
   v=dropdims(mean(v,dims=3),dims=3)

   #mask out near edge values to avoid exiting domain
   u[:,1:2].=NaN
   v[:,1:2].=NaN
   u[:,end-2:end].=NaN
   v[:,end-2:end].=NaN

   u=read(u,MeshArray(γ,Float32))
   v=read(v,MeshArray(γ,Float32));

   u[findall(isnan.(u))]=0.0
   v[findall(isnan.(v))]=0.0

   backward_in_time ? s=-1.0 : s=1.0
   u0=s*u; u1=s*u;
   v0=s*v; v1=s*v;

   t0=0.0; t1=86400*366*2.0;

   u0=u0./Γ["DXC"]
   u1=u1./Γ["DXC"]
   v0=v0./Γ["DYC"]
   v1=v1./Γ["DYC"]

   #nr=50; kk=1;

   #mskW=read(γ.path*"hFacW.latlon.data",MeshArray(γ,Float32,nr))
   #mskW=1.0 .+ 0.0 * mask(mskW[:,kk],NaN,0.0)
   #mskS=read(γ.path*"hFacS.latlon.data",MeshArray(γ,Float32,nr))
   #mskS=1.0 .+ 0.0 * mask(mskS[:,kk],NaN,0.0)
   #msk=Dict("mskW" => mskW, "mskS" => mskS)

   𝑃 = (u0=u0, u1=u1, v0=v0, v1=v1, 𝑇=[t0,t1], ioSize=(360,178),
        XC=exchange(Γ["XC"]), YC=exchange(Γ["YC"]))

   return 𝑃,Γ

end

"""
    OCCA_setup(;backward_in_time::Bool=false)

Define gridded variables and return result as Dictionary (`uvetc`).
"""
function OCCA_setup(;backward_in_time::Bool=false)

   γ=GridSpec("PeriodicChannel",MeshArrays.GRID_LL360)
   Γ=GridLoad(γ)
   n=length(Γ["RC"])
   n=10

   delete!.(Ref(Γ), ["hFacC", "hFacW", "hFacS","DXG","DYG","RAC","RAZ","RAS"]);
   backward_in_time ? s=-1.0 : s=1.0

   function rd(filename, varname,n)
   fil = NetCDF.open(filename, varname)
   siz = size(fil)
   tmp = zeros(siz[1:2]...,n)
   [tmp .+= fil[:,:,1:n,t] for t=1:12]
   tmp ./= 12.0
   tmp[findall(tmp.<-1e22)] .= 0.0
   return tmp
   end

   fileIn=OCCAclim_path*"DDuvel.0406clim.nc"
   u=s*read(rd(fileIn,"u",n),MeshArray(γ,Float32,n))

   fileIn=OCCAclim_path*"DDvvel.0406clim.nc"
   v=s*read(rd(fileIn,"v",n),MeshArray(γ,Float32,n))

   fileIn=OCCAclim_path*"DDwvel.0406clim.nc"
   w=s*rd(fileIn,"w",n)
   w=-cat(w,zeros(360, 160),dims=3)
   w[:,:,1] .=0.0
   w=read(w,MeshArray(γ,Float32,n+1))

   fileIn=OCCAclim_path*"DDtheta.0406clim.nc"
   θ=read(rd(fileIn,"theta",n),MeshArray(γ,Float32,n))

#   fileIn=OCCAclim_path*"DDsalt.0406clim.nc"
#   𝑆=read(rd(fileIn,"salt",n),MeshArray(γ,Float64,n))

   for i in eachindex(u)
      u[i]=u[i]./Γ["DXC"][1]
      v[i]=v[i]./Γ["DYC"][1]
   end

   for i in eachindex(u)
      u[i]=circshift(u[i],[-180 0])
      v[i]=circshift(v[i],[-180 0])
      θ[i]=circshift(θ[i],[-180 0])
#      𝑆[i]=circshift(𝑆[i],[-180 0])
   end

   for i in eachindex(w)
      w[i]=w[i]./Γ["DRC"][min(i[2]+1,n)]
      w[i]=circshift(w[i],[-180 0])
   end

   tmpx=circshift(Γ["XC"][1],[-180 0])
   tmpx[1:180,:]=tmpx[1:180,:] .- 360.0
   Γ["XC"][1]=tmpx

   tmpx=circshift(Γ["XG"][1],[-180 0])
   tmpx[1:180,:]=tmpx[1:180,:] .- 360.0
   Γ["XG"][1]=tmpx
   Γ["Depth"][1]=circshift(Γ["Depth"][1],[-180 0])

   t0=0.0; t1=86400*366*2.0;

   for k=1:n
    (tmpu,tmpv)=exchange(u[:,k],v[:,k],1)
    u[:,k]=tmpu
    v[:,k]=tmpv
   end
   for k=1:n+1
    tmpw=exchange(w[:,k],1)
    w[:,k]=tmpw
   end

   𝑃 = (θ0=θ, θ1=θ, u0=u, u1=u, v0=v, v1=v, w0=w, w1=w, 𝑇=[t0,t1],
   XC=exchange(Γ["XC"]), YC=exchange(Γ["YC"]),
   RF=Γ["RF"], RC=Γ["RC"],
   ioSize=(360,160,n))

   return 𝑃,Γ

end

"""
    init_global_range(lons::Tuple = (-160.0, -150.0),lats::Tuple = (35.0, 45.0))

Randomly distribute `np` points over a lon,la region, and
return position in grid index space (`i,j,subdomain`).
"""
function init_global_range(lons::Tuple = (-160.0, -150.0),lats::Tuple = (35.0, 45.0))
    lo0, lo1 = lons #(-160.0, -150.0)
    la0, la1 = lats #(35.0, 45.0)
    np = 100
    lon = lo0 .+ (lo1 - lo0) .* rand(np)
    lat = la0 .+ (la1 - la0) .* rand(np)
    (u0, _) = initialize_lonlat(Γ, lon, lat; msk = Γ["hFacC"][:, k])
    id=collect(1:np)
    return u0
end

"""
    init_global_randn(np ::Int , 𝑃::NamedTuple)

Randomly distribute `np` points over the Earth, within `𝑃.msk`
region, and return position in grid index space (`i,j,subdomain`).
"""
function init_global_randn(np ::Int , 𝑃::NamedTuple)
    (lon, lat) = randn_lonlat(2*np)
    (u0, _) = initialize_lonlat(𝑃.Γ, lon, lat; msk = 𝑃.msk)
    u0[:,1:np]
end

"""
    reset_lonlat!(𝐼::Individuals)

Randomly select a fraction (𝐼.𝑃.frac) of the particles and reset their positions.
"""
function reset_lonlat!(𝐼::Individuals)
    np=length(𝐼.🆔)
    n_reset = Int(round(𝐼.𝑃.frac*np))
    (lon, lat) = randn_lonlat(2*n_reset)
    (v0, _) = initialize_lonlat(𝐼.𝑃.Γ, lon, lat; msk = 𝐼.𝑃.msk)
    n_reset=min(n_reset,size(v0,2))
    k_reset = rand(1:np, n_reset)
    v0 = permutedims([v0[:,i] for i in 1:size(v0,2)])
    𝐼.📌[k_reset].=v0[1:n_reset]
    isempty(𝐼.🔴.ID) ? m=maximum(𝐼.🆔) : m=max(maximum(𝐼.🔴.ID),maximum(𝐼.🆔))
    𝐼.🆔[k_reset]=collect(1:n_reset) .+ m
end

##

"""
    isosurface(θ,T,z)

```
isosurface(𝐼.𝑃.θ0,15,Γ["RC"])
```
"""
function isosurface(θ,T,z)
    d=NaN*similar(θ[:,1])
    nr=size(θ,2)
    for j=1:size(d,1)
        for k=1:nr-1
            i=findall(isnan.(d[j]).&(θ[j,k].>T).&(θ[j,k+1].<=T))
            a=(θ[j,k][i] .- T)./(θ[j,k][i] .- θ[j,k+1][i])
            d[j][i]=(1 .- a).*Γ["RC"][k] + a.*Γ["RC"][k+1]
            i=findall(isnan.(d[j]).&(θ[j,k].<=T).&(θ[j,k+1].>T))
            a=(θ[j,k+1][i] .- T)./(θ[j,k+1][i] .- θ[j,k][i])
            d[j][i]=(1 .- a).*Γ["RC"][k+1] + a.*Γ["RC"][k]
        end
    end
    return d
end


text:
using Distributions
using Pathfinder
using Random
using Test
using Transducers

@testset "ELBO estimation" begin
    @testset "elbo_and_samples" begin
        σ_target = 0.08
        target_dist = Normal(0, σ_target)
        logp(x) = logpdf(target_dist, x[1])
        σs = [1e-3, 0.05, 0.8, 1.0, 1.1, 1.2, 5.0, 10.0]
        rng = Random.seed!(Random.default_rng(), 42)
        @testset for σ in σs
            dist = Normal(0, σ)
            est = @inferred Pathfinder.elbo_and_samples(rng, logp, dist, 1_000_000)
            @test est isa Pathfinder.ELBOEstimate

            # explicit elbo calculation
            r = σ / σ_target
            elbo_exp = (1 - r^2) / 2 + log(r)
            @test est.value ≈ elbo_exp atol = 3 * est.std_err
            @test est.log_densities_target ≈ logp.(eachcol(est.draws))
            @test est.log_densities_fit ≈ logpdf.(dist, first.(eachcol(est.draws)))
            @test est.log_density_ratios == est.log_densities_target - est.log_densities_fit
            @test est.value ≈ mean(est.log_density_ratios)
            @test est.std_err ≈ std(est.log_density_ratios) / sqrt(1_000_000)
        end
    end

    @testset "maximize_elbo" begin
        σ_target = 0.08
        target_dist = Normal(0, σ_target)
        logp(x) = logpdf(target_dist, x[1])
        σs = [1e-3, 0.05, σ_target, 1.0, 1.1, 1.2, 5.0, 10.0]
        dists = Normal.(0, σs)
        if VERSION ≥ v"1.7.0"
            executors = [SequentialEx(), ThreadedEx()]
        else
            executors = [SequentialEx()]
        end
        @testset "$executor" for executor in executors
            rng = Random.seed!(Random.default_rng(), 42)
            lopt, estimates = @inferred Pathfinder.maximize_elbo(
                rng, logp, dists, 100, executor
            )
            @test lopt == 3
            @test estimates[lopt].value ≈ 0
            rng = Random.seed!(Random.default_rng(), 42)
            lopt2, estimates2 = Pathfinder.maximize_elbo(rng, logp, dists, 100, executor)
            @test lopt2 == lopt
            @test getproperty.(estimates2, :value) == getproperty.(estimates, :value)
            @test getproperty.(estimates2, :std_err) == getproperty.(estimates, :std_err)
            lopt3, estimates3 = @inferred Pathfinder.maximize_elbo(
                rng, logp, dists[2:1], 100, executor
            )
            @test lopt3 == 0
            @test isempty(estimates3)
        end
    end
end


text:
# Set up model for testing
using ModelObjectsLH

mutable struct Obj1 <: ModelObject
    objId :: ObjectId
    x :: Float64
    y :: Vector{Float64}
    z :: Array{Float64,2}
end

function init_obj1(objId)
    # objId = ObjectId(:obj1);
    valueX = 7.3;
    # Important to have vector of length 1 as test case
    valueY = fill(1.1, 1);
    valueZ = [3.3 4.4; 5.5 7.6];
    o1 = Obj1(objId, valueX, valueY, valueZ);
    return o1
end

mutable struct Obj3 <: ModelObject
    objId :: ObjectId
    x :: Float64
    y :: Vector{Int}
end

function init_obj3(objId :: ObjectId)
    return Obj3(make_child_id(objId, :obj3), 0.5, [1,2]);
end

mutable struct Obj2 <: ModelObject
    objId :: ObjectId
    a :: Float64
    y :: Float64
    b :: Array{Float64,2}
    obj3 :: Obj3
end

# function Obj2(a, y, b)
#     return Obj2(a, y, b, ParamVector(ObjectId(:pv1)))
# end

function init_obj2(objId)
    # objId = ObjectId(:obj2);
    valueX = 17.3;
    valueY = 9.4;
    valueB = 2.0 .+ [3.3 4.4; 5.5 7.6];
    o2 = Obj2(objId, valueX, valueY, valueB, init_obj3(objId));
    return o2
end

mutable struct TestModel <: ModelObject
    objId :: ObjectId
    o1 :: Obj1
    o2 :: Obj2
    a :: Float64
    y :: Float64
end

function init_test_model()
    objName = ObjectId(:testModel, "Test model");
    o1 = init_obj1(make_child_id(objName, :o1, "Child object 1"));
    o2 = init_obj2(make_child_id(objName, :o2, "Child object 2"));
    return TestModel(objName, o1, o2, 9.87, 87.73)
end

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


text:
using Documenter, Fuji

makedocs(
    format = :html,
    sitename = "Chaplin"
)

deploydocs(
    repo = "github.com/timo-cmd/chaplin.jl.git",
    julia = "0.6",
    osname = "windows",
    deps = nothing,
    make = nothing,
    target = "build"
)


text:
function testNHHill_M()
close all
clear all

% Physical parameters
Param=PhysParameters();

nx=60;
ny=2;
lx=600000;
ly=20000;
x0=-300000;
y0=0;
Boundary.WE='Period';
Boundary.BT='Period';
Param.hS='';
Param.Grid=CartGrid(nx,ny,lx,ly,x0,y0,OrientFaceCart,Boundary,Param);

Param.H=15600;
nz=40;
zP=zeros(nz,1);
z=zeros(nz+1,1);
dz=Param.H/nz;
zP(1)=dz/2;
for i=2:nz
  zP(i)=zP(i-1)+dz;
end
for i=2:nz+1
  z(i)=z(i-1)+dz;
end

Param.Grid.nz=nz;
Param.Grid.z=z;
Param.Grid.zP=zP;
Param.Grid.dz=dz;

% Model
Param.ModelType='Curl';
%Param.ModelType='Div';
Param.Coriolis=false;
Param.Thermo='';
Param.Source=false;
Param.HyperD=0; %1.e15*(nx/30)^3.2;
Param.StrideDamp=6000;
Param.Relax=1.e-1;
Param.Damping=true;
Param.Flat=false;
Param.Coriolis=false;
Param.Buoyancy=true;
Param.hC=400;
Param.x0C=0;
Param.aC=10000;
Param.Th0=300;
Param.uMax=10;
Param.vMax=0;
Param.NBr=1.e-2;
Param.Equation='Compressible';
Param.TopoS='AgnesiCartX';
Param.EndTime=21600;

% Discretization
OrdPoly=4;
OrdPolyZ=1;
[CG,Param]=Discretization(OrdPoly,OrdPolyZ,@JacobiDG3,Param);
LRef=11*1.e5;
Param.HyperVisc=true;
Param.HyperDCurl=2.e4; %1.e14*(dx/LRef)^3.2;
Param.HyperDGrad=2.e4;
Param.HyperDDiv=2.e4;
Param.Upwind=false;

% Initial conditions
Param.NumV=5;
U=zeros(CG.NumG,nz,Param.NumV);
Param.ProfRho='GravityHill';
Param.ProfTheta='GravityHill';
Param.ProfVel='Const';
Param.RhoPos=1;
Param.uPos=2;
Param.vPos=3;
Param.wPos=4;
Param.ThPos=5;
U(:,:,Param.RhoPos)=Project(@fRho,CG,Param);
[U(:,:,Param.uPos),U(:,:,Param.vPos)]=ProjectVec(@fVel,CG,Param);
U(:,:,Param.ThPos)=Project(@fTheta,CG,Param).*U(:,:,Param.RhoPos);

% Output
Slice.Type='XZ';
Slice.y=y0+ly/ny/2;
Slice.iy=1;
Param.SliceXZ=Slice;
Param.SliceXY.Type='XY';
Param.SliceXY.iz=1;
Param.fig=1;
Param.vtk=0;
Param.vtkFileName='Hill_M';
Param.vtk=0;
vtkGrid=vtkCGGrid(CG,@TransCart,@Topo,Param);
Param.cNames(1).s='u';
Param.cNames(2).s='w';
cOut=zeros(CG.NumG,nz,2);

IntMethod='Rosenbrock';
%IntMethod='RungeKutta';
if strcmp(IntMethod,'Rosenbrock')
  dtau=3;
else
  dtau=.4;
end
nIter=6000;
%v = VideoWriter ('Galewsky.avi');
%open (v);
Param.RK=RungeKuttaMethod('RK4');
Param.ROS=RosenbrockMethod('ROSRK3');
CFL=0.125;
time=0;
nIter=Param.EndTime/dtau;
PrintTime=1000;
PrintInt=floor(PrintTime/dtau);
cOut(:,:,1)=U(:,:,Param.uPos);
W=BoundaryWOutput(U,CG,Param);
cOut(:,1,2)=W;
cOut(:,2:nz-1,2)=0.5*(U(:,1:nz-2,Param.wPos)+U(:,2:nz-1,Param.wPos));
cOut(:,nz,2)=0.5*U(:,nz-1,Param.wPos);
Param.vtk=vtkCG(cOut,CG,Param,vtkGrid,Param.vtk);
switch IntMethod
  case 'Rosenbrock'
    tic
    for i=1:nIter
      i
      U=RosenbrockSchur(U,dtau,@FcnNHCurlVec,@JacSchur,CG,Param);
      time=time+dtau;
      if mod(i,PrintInt)==0
        cOut(:,:,1)=U(:,:,Param.uPos);
        W=BoundaryWOutput(U,CG,Param);
        cOut(:,1,2)=W;
        cOut(:,2:nz-1,2)=0.5*(U(:,1:nz-2,Param.wPos)+U(:,2:nz-1,Param.wPos));
        cOut(:,nz,2)=0.5*U(:,nz-1,Param.wPos);
        Param.vtk=vtkCG(cOut,CG,Param,vtkGrid,Param.vtk);
      end
    end
    toc
  case 'RungeKutta'
    for i=1:nIter
      i
      U=RungeKuttaExplicit(U,dtau,@FcnNHCurlVec,CG,Param);
      time=time+dtau;
      if mod(i,PrintInt)==0
        cOut(:,:,1)=U(:,:,Param.uPos);
        W=BoundaryWOutput(U,CG,Param);
        cOut(:,1,2)=0.5*(U(:,1,Param.wPos)+W);
        cOut(:,2:nz-1,2)=0.5*(U(:,1:nz-2,Param.wPos)+U(:,2:nz-1,Param.wPos));
        cOut(:,nz,2)=0.5*U(:,nz-1,Param.wPos);
        Param.vtk=vtkCG(cOut,CG,Param,vtkGrid,Param.vtk);
      end
    end
end
%fig=PlotCG(U(:,1),CG,@JacobiSphere2,Param,fig);
%close(v)
end


text:
export
    global_verbosity,
    set_global_verbosity!

global VERBOSITY = 0

"""
	global_verbosity()

Query the global verbosity setting (VERBOSITY::Int)
"""
global_verbosity() = copy(VERBOSITY)

"""
	set_global_verbosity!(val::Int)

Set the global verbosity to val
"""
set_global_verbosity!(val::Int) = begin global VERBOSITY = val end

function print_styled_header(header,msg...;bold=true,color=136)
	printstyled(header;bold=bold,color=color)
	println(msg...)
end

export @log_info

"""
    @log_info

A helper macro for printing at various verbosity levels.
Usage:
	`@log_info(limit::Int,level::Int,msg...)`
	`@log_info(limit::Int,solver,msg...)`
Args:
* limit::Int - logging threshold
* level::Int - the verbosity level
* msg... - message to print if level > limit
"""
macro log_info(limit,level,msg...)
	filename = split(string(__source__.file),"/")[end]
	ex = :(print_styled_header(log_preamble))
    for x in esc(msg).args[1]
		push!(ex.args,esc(x))
    end
	return quote
		local log_preamble = string("[ logger ](",$filename,"-",$(__source__.line),"): ")
		$(esc(level)) > $(esc(limit)) || global_verbosity() > $(esc(limit)) ?
			$ex : nothing
	end
end


text:
export BrickletSegmentDisplay4x7V2Segments
struct BrickletSegmentDisplay4x7V2Segments
    digit0::Vector{Bool}
    digit1::Vector{Bool}
    digit2::Vector{Bool}
    digit3::Vector{Bool}
    colon::Vector{Bool}
    tick::Bool
end

export BrickletSegmentDisplay4x7V2SPITFPErrorCount
struct BrickletSegmentDisplay4x7V2SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletSegmentDisplay4x7V2Identity
struct BrickletSegmentDisplay4x7V2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletSegmentDisplay4x7V2
"""
Four 7-segment displays with switchable dots
"""
mutable struct BrickletSegmentDisplay4x7V2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletSegmentDisplay4x7V2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_segment_display_4x7_v2")
        deviceInternal = package.BrickletSegmentDisplay4x7V2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_segments
"""
    $(SIGNATURES)

Sets the segments of the Segment Display 4x7 Bricklet 2.0 segment-by-segment.

The data is split into the four digits, two colon dots and the tick mark.

The indices of the segments in the digit and colon parameters are as follows:

.. image:: /Images/Bricklets/bricklet_segment_display_4x7_v2_segment_index.png
   :scale: 100 %
   :alt: Indices of segments
   :align: center
"""
function set_segments(device::BrickletSegmentDisplay4x7V2, digit0, digit1, digit2, digit3, colon, tick)
    device.deviceInternal.set_segments(digit0, digit1, digit2, digit3, colon, tick)
end

export get_segments
"""
    $(SIGNATURES)

Returns the segment data as set by :func:`Set Segments`.
"""
function get_segments(device::BrickletSegmentDisplay4x7V2)
    return device.deviceInternal.get_segments()
end

export set_brightness
"""
    $(SIGNATURES)

The brightness can be set between 0 (dark) and 7 (bright).
"""
function set_brightness(device::BrickletSegmentDisplay4x7V2, brightness)
    device.deviceInternal.set_brightness(brightness)
end

export get_brightness
"""
    $(SIGNATURES)

Returns the brightness as set by :func:`Set Brightness`.
"""
function get_brightness(device::BrickletSegmentDisplay4x7V2)
    return device.deviceInternal.get_brightness()
end

export set_numeric_value
"""
    $(SIGNATURES)

Sets a numeric value for each of the digits. They represent:

* -2: minus sign
* -1: blank
* 0-9: 0-9
* 10: A
* 11: b
* 12: C
* 13: d
* 14: E
* 15: F

Example: A call with [-2, -1, 4, 2] will result in a display of "- 42".
"""
function set_numeric_value(device::BrickletSegmentDisplay4x7V2, value)
    device.deviceInternal.set_numeric_value(value)
end

export set_selected_segment
"""
    $(SIGNATURES)

Turns one specified segment on or off.

The indices of the segments are as follows:

.. image:: /Images/Bricklets/bricklet_segment_display_4x7_v2_selected_segment_index.png
   :scale: 100 %
   :alt: Indices of selected segments
   :align: center
"""
function set_selected_segment(device::BrickletSegmentDisplay4x7V2, segment, value)
    device.deviceInternal.set_selected_segment(segment, value)
end

export get_selected_segment
"""
    $(SIGNATURES)

Returns the value of a single segment.
"""
function get_selected_segment(device::BrickletSegmentDisplay4x7V2, segment)
    return device.deviceInternal.get_selected_segment(segment)
end

export start_counter
"""
    $(SIGNATURES)

Starts a counter with the *from* value that counts to the *to*
value with the each step incremented by *increment*.
*length* is the pause between each increment.

Example: If you set *from* to 0, *to* to 100, *increment* to 1 and
*length* to 1000, a counter that goes from 0 to 100 with one second
pause between each increment will be started.

Using a negative *increment* allows to count backwards.

You can stop the counter at every time by calling :func:`Set Segments`
or :func:`Set Numeric Value`.
"""
function start_counter(device::BrickletSegmentDisplay4x7V2, value_from, value_to, increment, length)
    device.deviceInternal.start_counter(value_from, value_to, increment, length)
end

export get_counter_value
"""
    $(SIGNATURES)

Returns the counter value that is currently shown on the display.

If there is no counter running a 0 will be returned.
"""
function get_counter_value(device::BrickletSegmentDisplay4x7V2)
    return device.deviceInternal.get_counter_value()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletSegmentDisplay4x7V2)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletSegmentDisplay4x7V2, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletSegmentDisplay4x7V2)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletSegmentDisplay4x7V2, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletSegmentDisplay4x7V2, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletSegmentDisplay4x7V2, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletSegmentDisplay4x7V2)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletSegmentDisplay4x7V2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletSegmentDisplay4x7V2)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletSegmentDisplay4x7V2, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletSegmentDisplay4x7V2)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletSegmentDisplay4x7V2)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletSegmentDisplay4x7V2, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
struct mbedtls_mpi
    s::Cint
    n::Csize_t
    p::Ptr{Cuint}
end

struct mbedtls_rsa_context
    ver::Cint
    len::Csize_t
    N::mbedtls_mpi
    E::mbedtls_mpi
    D::mbedtls_mpi
    P::mbedtls_mpi
    Q::mbedtls_mpi
    # More fields follow, but omitted here, since they
    # are not required for this wrapper
end

mutable struct RSA
    data::Ptr{mbedtls_rsa_context}

    function RSA(padding=MBEDTLS_RSA_PKCS_V21, hash_id=MD_MD5)
        ctx = new()
        ctx.data = Libc.malloc(1000)
        ccall((:mbedtls_rsa_init, libmbedcrypto), Cvoid,
            (Ptr{Cvoid}, Cint, Cint),
            ctx.data, padding, hash_id)
        finalizer(ctx->begin
            ccall((:mbedtls_rsa_free, libmbedcrypto), Cvoid, (Ptr{Cvoid},), ctx.data)
            Libc.free(ctx.data)
        end, ctx)
        ctx
    end
end

function mpi_import!(mpi::Ptr{mbedtls_mpi}, b::BigInt)
    # Export from GMP
    size = ndigits(b, base=2)
    nbytes = div(size+8-1,8)
    data = Vector{UInt8}(undef, nbytes)
    count = Ref{Csize_t}(0)
    # TODO Replace `Any` with `Ref{BigInt}` when 0.6 support is dropped.
    ccall((:__gmpz_export,:libgmp), Ptr{Cvoid},
            (Ptr{Cvoid}, Ptr{Csize_t}, Cint, Csize_t, Cint, Csize_t, Any),
            data, count, 1, 1, 1, 0, b)
    @assert count[] == nbytes
    # Import into mbedtls
    @err_check ccall((:mbedtls_mpi_read_binary, libmbedcrypto), Cint,
        (Ptr{mbedtls_mpi}, Ptr{UInt8}, Csize_t),
        mpi, data, nbytes)
end

function mpi_size(mpi::Ptr{mbedtls_mpi})
    ccall((:mbedtls_mpi_size, libmbedcrypto), Csize_t, (Ptr{mbedtls_mpi},), mpi)
end

function pubkey_from_vals!(ctx::RSA, e::BigInt, n::BigInt)
    Nptr = Ptr{mbedtls_mpi}(ctx.data+fieldoffset(mbedtls_rsa_context,3 #= :N =#))
    mpi_import!(Nptr, n)
    mpi_import!(Ptr{mbedtls_mpi}(ctx.data+fieldoffset(mbedtls_rsa_context,4 #= :E =#)), e)
    nptr_size = mpi_size(Nptr)
    unsafe_store!(Ptr{Csize_t}(ctx.data+fieldoffset(mbedtls_rsa_context,2 #=:len =#)),
        nptr_size)
    @err_check ccall((:mbedtls_rsa_check_pubkey, libmbedcrypto), Cint,
        (Ptr{Cvoid},), ctx.data)
    ctx
end

function verify(ctx::RSA, hash_alg::MDKind, hash, signature, rng = nothing; using_public=true)
    (!using_public && rng == nothing) &&
        error("Private key verification requires the rng")
    # All errors, including validation errors throw
    @err_check ccall((:mbedtls_rsa_pkcs1_verify, libmbedcrypto), Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, Any, Cint, Cint, Csize_t, Ptr{UInt8}, Ptr{UInt8}),
        ctx.data,
        rng == nothing ? C_NULL : c_rng[],
        rng == nothing ? Ref{Any}() : rng,
        using_public ? 0 : 1,
        hash_alg, sizeof(hash), hash, signature)
end

function gen_key!(ctx::RSA, f_rng, rng, nbits, exponent)
    @err_check ccall((:mbedtls_rsa_gen_key, libmbedcrypto), Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, Any, Cint, Cint),
        ctx.data, f_rng, rng, nbits, exponent)
    ctx
end

function gen_key(rng::AbstractRNG, nbits=2048, exponent=65537)
    ctx = RSA()
    gen_key!(ctx, c_rng[], rng, nbits, exponent)
    ctx
end

function public(ctx::RSA, input, output)
    @err_check ccall((:mbedtls_rsa_public, libmbedcrypto), Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), ctx.data, input, output)
    output
end

function private(ctx::RSA, f_rng, rng, input, output)
    @err_check ccall((:mbedtls_rsa_private, libmbedcrypto), Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, Any, Ptr{Cvoid}, Ptr{Cvoid}),
        ctx.data, f_rng, rng, input, output)
    output
end

function private(ctx::RSA, rng::AbstractRNG, input, output)
    private(ctx, c_rng[], rng, input, output)
end


text:
function HMMsmooth(v,phghm,ph1,pvgh,alpha=[],beta=[];ReturnPairwiseMarginals=false,ReturnLogLikelihood=false)
    #HMMSMOOTH Smoothing for a Hidden Markov Model (HMM)
    # [phtgV1T,phthtpgV1T]=HMMsmooth(v,phghm,ph1,pvgh)
    # Return the HMM pointwise p(h(t)|v(1:T)) and pairwise posterior p(h(t),h(t+1)|#v(1:T)).
    #
    # Inputs:
    # v : visible (observation) sequence
    # phghm : transition distribution
    # ph1 : initial distribution
    # pvgh : emission distribution p(v|h)
    #
    # Outputs:
    # phtgV1T : smoothed posterior p(h(t)|v(1:T))
    # phthtpgV1T : smoothed pair p(h(t),h(t+1)|v(1:T))
    T=length(v); H=size(phghm,1);

    if isempty(alpha)
        alpha,loglik=HMMforward(v,phghm,ph1,pvgh)
    end
    if isempty(beta)
        beta=HMMbackward(v,phghm,pvgh)
    end

    # smoothed posteriors: pointwise marginals:
    phtgV1T=zeros(H,T)
    for t=1:T
        tmp=alpha[:,t].*beta[:,t]
	phtgV1T[:,t]=tmp./sum(tmp)
    end

    if ReturnPairwiseMarginals
        # smoothed posteriors: pairwise marginals p(h(t),h(t+1)|v(1:T)):
        phthtpgV1T=zeros(H,H,T-1)
        for t=1:T-1
	    atmp=alpha[:,t]
	    btmp=beta[:,t+1]
	    #ctmp = repmat(atmp,1,H).*phghm'.*repmat(pvgh[v[t+1],:].*btmp',H,1)
            ctmp = repmat(atmp,1,H).*phghm'.*repmat(pvgh[v[t+1],:]'.*btmp',H,1)
	    phthtpgV1T[:,:,t]=ctmp./sum(ctmp)
        end
        if ReturnLogLikelihood
            return phtgV1T, phthtpgV1T, loglik
        end
        return phtgV1T, phthtpgV1T
    end

    if ReturnLogLikelihood
        return phtgV1T,loglik
    end
    return phtgV1T
end


text:
function make_household_ptrs!(
    ptrs::AbstractVector{Tuple{Ti,Ti}},
    household_indices::AbstractVector{T} where T<:Integer
    ) where Ti<:Integer

    @assert length(ptrs) == length(household_indices)

    ptrarr = reshape(reinterpret(Ti, ptrs), 2, :) # tuples as 2-by-N array
    headptrs = view(ptrarr, 1, :)
    tailptrs = view(ptrarr, 2, :)

    groupptrs!(headptrs, tailptrs, household_indices)
    ptrs
  end

  make_household_ptrs(household_indices::AbstractVector{T} where T<:Real) =
    make_household_ptrs!(Vector{Tuple{PersonIdx, PersonIdx}}(undef, length(household_indices)), household_indices)


text:
function neighsingles(x, q, atomtype, ilist)

if !isempty(x)
    xi = x[:,ilist];
else
    xi = [];
end
if !isempty(q)
    qi = q[:,ilist];
else
    qi = [];
end

ai = ilist;
ti = atomtype[ilist];

return xi, qi, ai, ti

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using VisualCircuits, Compose

begin
    width, height = 5, 4
    nand1 = Gate(2.5, 1, :NAND)
    nand2 = Gate(2.5, 3, :NAND)

    circuit = begin
        compose(
        # Range of coordinate system
        context(units=UnitBox(0, 0, width, height)),

        # Gates
        render(nand1), render(nand2),

        # Connect to edges (point, x_value)
        connect(nand1[:in1], 1),
        connect(nand2[:in2], 1),
        connect(nand1[:out], width-1),
        connect(nand2[:out], width-1),

        # crossing lines
        connect(
            nand1[:in2],
            (1.5, nand1[:in2][2]),
            (1.5, 1.5), (3.5, 2.5),
            (3.5, nand2[:out][2]),
            nand2[:out]
        ),
        connect(
            nand2[:in1],
            (1.5, nand2[:in1][2]),
            (1.5, 2.5), (3.5, 1.5),
            (3.5, nand1[:out][2]),
            nand1[:out]
        ),

        # connections/dots
        dots(
            (3.5, nand1[:out][2]),
            (3.5, nand2[:out][2])
        ),

        # render text
        render(
            (0.5, nand1[:in1][2], "S"),
            (0.5, nand2[:in2][2], "R"),
            (width-0.5, nand1[:out][2], "Q"),
            (width-0.5, nand2[:out][2], "~Q")
        ),

        # compose options
        stroke("black"), linewidth(2mm)
        )
    end

    # draw(PDF("RS_FlipFlop.pdf", autosize(width, height)...), circuit)
    draw(SVG("RS_FlipFlop.svg", autosize(width, height)...), circuit)
end