text:
#=
A Bayesian Network (BN) represents a probability distribution
over a set of variables, P(x₁, x₂, ..., xₙ)
It leverages relations between variables in order to efficiently encode the joint distribution.
A BN is defined by a directed acyclic graph in which each node is a variable
and contains an associated conditional probability distribution P(xⱼ | parents(xⱼ))
=#

const DAG = DiGraph

function _build_dag_from_cpds(
	cpds::AbstractVector{T},
	name_to_index::Dict{NodeName, Int}
) where {T<:CPD}

	dag = DAG(length(cpds))
	for cpd in cpds
		j = name_to_index[name(cpd)]
		for p in parents(cpd)
			i = name_to_index[p]
			add_edge!(dag, i, j)
		end
	end
	dag
end
function _enforce_topological_order(
	dag::DAG,
	cpds::AbstractVector{T},
	name_to_index::Dict{NodeName, Int},
) where {T<:CPD}

	topo = topological_sort_by_dfs(dag)

	cpds2 = Array{T}(undef, length(cpds))
	name_to_index2 = Dict{NodeName, Int}()
	for (i,j) in enumerate(topo)
		# i is the new index
		# j is the old index

		cpd = cpds[j]
		name_to_index2[name(cpd)] = i
		cpds2[i] = cpd
	end

	dag2 = _build_dag_from_cpds(cpds2, name_to_index2)

	(dag2, cpds2, name_to_index2)
end

mutable struct BayesNet{T<:CPD} <: ProbabilisticGraphicalModel
	dag::DAG # nodes are in topological order
	cpds::Vector{T} # the CPDs associated with each node in the dag
	name_to_index::Dict{NodeName,Int} # NodeName → index in dag and cpds
end
BayesNet() = BayesNet(DAG(0), CPD[], Dict{NodeName, Int}())
BayesNet(::Type{T}) where {T <: CPD} = BayesNet(DAG(0), T[], Dict{NodeName, Int}())
function BayesNet(cpds::AbstractVector{T}) where {T <: CPD}

	name_to_index = Dict{NodeName, Int}()
	for (i, cpd) in enumerate(cpds)
		name_to_index[name(cpd)] = i
	end

	dag = _build_dag_from_cpds(cpds, name_to_index)

	!is_cyclic(dag) || error("BayesNet graph is non-acyclic!")

	(dag2, cpds2, name_to_index2) = _enforce_topological_order(dag, cpds, name_to_index)

	BayesNet(dag2, cpds2, name_to_index2)
end

Base.get(bn::BayesNet, i::Int) = bn.cpds[i]
Base.get(bn::BayesNet, name::NodeName) = bn.cpds[bn.name_to_index[name]]
Base.length(bn::BayesNet) = length(bn.cpds)

"""
Returns the ordered list of NodeNames
"""
function Base.names(bn::BayesNet)
	retval = Array{NodeName}(undef, length(bn))
	for (i,cpd) in enumerate(bn.cpds)
		retval[i] = name(cpd)
	end
	retval
end

"""
Returns the parents as a list of NodeNames
"""
CPDs.parents(bn::BayesNet, name::NodeName) = parents(get(bn, name))

"""
Returns the children as a list of NodeNames
"""
function children(bn::BayesNet, target::NodeName)
	i = bn.name_to_index[target]
	NodeName[name(bn.cpds[j]) for j in outneighbors(bn.dag, i)]
end

"""
Returns all neighbors as a list of NodeNames.
"""
function neighbors(bn::BayesNet, target::NodeName)
	i = bn.name_to_index[target]
	NodeName[name(bn.cpds[j]) for j in append!(inneighbors(bn.dag, i), outneighbors(bn.dag, i))]
end

"""
Returns all descendants as a list of NodeNames.
"""
dst(edge::Pair{Int,Int}) = edge[2] # LightGraphs used to return a Pair, now it returns a SimpleEdge
function descendants(bn::BayesNet, target::NodeName)
	retval = Set{Int}()
	for edge in edges(bfs_tree(bn.dag, bn.name_to_index[target]))
		push!(retval, dst(edge))
	end
	NodeName[name(bn.cpds[i]) for i in sort!(collect(retval))]
end

"""
Return the children, parents, and parents of children (excluding target) as a Set of NodeNames
"""
function markov_blanket(bn::BayesNet, target::NodeName)
	nodeNames = NodeName[]
        for child in children(bn, target)
		append!(nodeNames, parents(bn, child))
                push!(nodeNames, child)
        end
        append!(nodeNames, parents(bn, target))
        return setdiff(Set(nodeNames), Set(NodeName[target]))
end

"""
Whether the BayesNet contains the given edge
"""
function has_edge(bn::BayesNet, parent::NodeName, child::NodeName)::Bool
	u = get(bn.name_to_index, parent, 0)
	v = get(bn.name_to_index, child, 0)
	u != 0 && v != 0 && has_edge(bn.dag, u, v)
end

function enforce_topological_order!(bn::BayesNet)
	dag2, cpds2, name_to_index2 = _enforce_topological_order(bn.dag, bn.cpds, bn.name_to_index)
	bn.dag = dag2
	bn.cpds = cpds2
	bn.name_to_index = name_to_index2
	bn
end
function adding_edge_preserves_acyclicity(parent_list::Vector{Vector{Int}}, u::Int, v::Int)
    n = length(parent_list)
    visited = falses(n)
    propagate = [u]
    while !isempty(propagate)
        x = pop!(propagate)
        if x == v
            return false
        end
        if !visited[x]
            visited[x] = true
            append!(propagate, parent_list[x])
        end
    end
    return true
end

"""
Returns whether the set of node names `x` is d-separated from the set `y` given the set `given`
"""
function is_independent(bn::BayesNet, x::NodeNames, y::NodeNames, given::NodeNames)

	start_node = x[1]
	finish_node = y[1]

	if start_node == finish_node
		return true
	end

	C = Set(given)
	analyzed_nodes = Set()

	# Find all paths from x to y
	paths = []
	path_queue = []
	push!(analyzed_nodes, start_node)
	for next_node in neighbors(bn, start_node)
		if !in(next_node, analyzed_nodes)
			push!(analyzed_nodes, next_node)
			push!(path_queue, [start_node, next_node])
		end
	end
	while (!isempty(path_queue))
		cur_path = pop!(path_queue)
	 	last_node = cur_path[end]
		for next_node in neighbors(bn, last_node)
			if next_node == finish_node
				push!(paths, push!(copy(cur_path), next_node))
			elseif !in(next_node, analyzed_nodes)
				push!(analyzed_nodes, next_node)
				push!(path_queue, push!(copy(cur_path), next_node))
			end
		end
	end

	# Check each path to see if it contains information indicating d-separation
	for path in paths
		is_d_separated = false
		if length(path) == 2
			is_d_separated = true
		else
			# Examine all middle nodes
			for i in 2:(length(path) - 1)
				prev_node = path[i - 1]
				cur_node = path[i]
				next_node = path[i + 1]

				# Check for chain or fork (first or second d-separation criteria)
				if in(cur_node, C)

					# Chain
					if in(cur_node, children(bn, prev_node)) && in(next_node, children(bn, cur_node))
						is_d_separated = true
						break

					# Fork
					elseif in(prev_node, children(bn, cur_node)) && in(next_node, children(bn, cur_node))
						is_d_separated = true
						break
					end

				# Check for v-structure (third d-separation criteria)
				else
					if in(cur_node, children(bn, prev_node)) && in(cur_node, children(bn, next_node))
						descendant_list = descendants(bn, cur_node)
						descendants_in_C = false
						for d in descendant_list
							if in(d, C)
								descendants_in_C = true
								break
							end
						end

						if !descendants_in_C
							is_d_separated = true
							break
						end
					end
				end
			end
		end

		if !is_d_separated
			return false
		end
	end

	# All paths are d-separated, so x and y are conditionally independent.
	return true
end

function Base.push!(bn::BayesNet, cpd::CPD)

	cpdname = name(cpd)
	!haskey(bn.name_to_index, cpdname) || error("A CPD with name $cpdname already exists!")

	add_vertex!(bn.dag)

	push!(bn.cpds, cpd)
	bn.name_to_index[cpdname] = j = length(bn.cpds)

	# add the necessary edges
	for p in parents(cpd)
		i = bn.name_to_index[p]
		add_edge!(bn.dag, i, j)
	end

	!is_cyclic(bn.dag) || error("BayesNet graph is non-acyclic!")
	enforce_topological_order!(bn)
end
function Base.append!(bn::BayesNet, cpds::AbstractVector{C}) where {C<:CPD}
	for cpd in cpds
		push!(bn, cpd)
	end
	bn
end

"""
	delete!(bn::BayesNets, target::NodeName)
Removing cpds will alter the vertex indeces. In particular, removing
the ith cpd will swap i and n and then remove n.
"""
function Base.delete!(bn::BayesNet, target::NodeName)

	if outdegree(bn.dag, bn.name_to_index[target]) > 0
		@warn("Deleting a CPD with children!")
	end

	i = bn.name_to_index[target]
	replacement = name(bn.cpds[end])
	bn.cpds[i] = bn.cpds[end]
	pop!(bn.cpds)
	bn.name_to_index[replacement] = i
	delete!(bn.name_to_index, target)
	rem_vertex!(bn.dag, i)
	bn
end

"""
The pdf of a given assignment after conditioning on the values
"""
function CPDs.pdf(bn::BayesNet, assignment::Assignment)
	retval = 1.0
	for cpd in bn.cpds # NOTE: guaranteed in topological order
		retval *= pdf(cpd, assignment)
	end
 	retval
end

"""
The logpdf of a given assignment after conditioning on the values
"""
function CPDs.logpdf(bn::BayesNet, assignment::Assignment)
	retval = 0.0
	for cpd in bn.cpds # NOTE: guaranteed in topological order
		retval += logpdf(cpd, assignment)
	end
 	retval
end


text:
"""
    average(em, state_0; kwargs)

Compute the average trajectory of an `AbstractEpidemicModel` over multiple stochastic simulations from the initial state `state_0`.

Keyword arguments:
* `nsims=1000`: the number of simulations.
* `nbins:100`: the number of time steps at which the averages should be computed. Time steps are uniformly distributed between `0` and `tmax`.
* `tmax=100.0`: the maximum time of each simulation.
* `nmax=1000`: the maximum number of iterations of each simulation.
* `progressbar=progressbar(em)`: if set to `true`, a progress bar will be displayed while the algorithm is running.
"""
function average(em::AbstractEpidemicModel, state_0; nsims=1000, nbins=100, tmax=100.0, nmax=1000, progressbar=progressbar(em))
    ts = LinRange(0.0, tmax, nbins)
    sums = init_sums(em, nbins)
    if progressbar
        p = Progress(nsims, dt=1.0)
    end
    for n in 1:nsims
        t, states = gillespie(em, state_0, tmax=tmax, nmax=nmax, process_output=false)
        l = 1
        for k in 1:nbins
            while l < length(t) && t[l+1] <= ts[k]
                l += 1
            end
            update_sums!(sums, k, l, states, em)
        end
        if progressbar
            ProgressMeter.next!(p; showvalues = [(:iteration,n)])
        end
    end
    return ts, finalize_sums(sums, nsims, em)
end

"""
    progressbar(::AbstractEpidemicModel)

 Return whether a progress bar should be displayed by default when using `average`.
"""
progressbar(::AbstractEpidemicModel) = false

"""
    init_sums(em, nbins)

Return an array of size `nbins` of `em`'s output type, initialised to 0.
"""
function init_sums(em::AbstractEpidemicModel, nbins)
    dims = output_elem_size(em)
    sums = [zeros(output_type(em), dims) for i in 1:nbins]
    return sums
end

"""
    update_sums!(sums, k, l, states, em)

Add the information from `states[l]` to `sums[k]`.
"""
function update_sums!(sums, k, l, states, em::AbstractEpidemicModel)
    sums[k] .+= states[l]
end

"""
    finalize_sums(sums, nsims, em)

Perform some final operations on the agregated data from `sums`, e.g. dividing by `nsims` to get the mean.
"""
function finalize_sums(sums, nsims, em::AbstractEpidemicModel)
    return float.(finalize_output(em, sums))./nsims
end

ext: jl


text:
module SMRTypes

import Base: show

export cSMRWMrkChannel, SMRWMrkChannel, cSMRContChannel, SMRContChannel,
       cSMREventChannel, SMREventChannel, cSMRMarkerChannel, SMRMarkerChannel,
       cSMRChannelInfo, cSMRChannelInfoArray, SMRChannelInfo, show,
       channel_string

const MARKER_SIZE = UInt8(4)

abstract type SMRCType end
abstract type SMRType end

# =========================================================================== #
struct cSMRWMrkChannel <: SMRCType
    length::UInt64
    npt::UInt64
    timestamps::Ptr{Float64}
    markers::Ptr{UInt8}
    wavemarks::Ptr{Int16}
end

mutable struct SMRWMrkChannel <: SMRType
    timestamps::Vector{Float64}
    markers::Matrix{UInt8}
    wavemarks::Matrix{Int16}

    function SMRWMrkChannel(x::cSMRWMrkChannel)
        self = new()

        self.timestamps = Vector{Float64}(undef, x.length)
        copyto!(self.timestamps, unsafe_wrap(Vector{Float64}, x.timestamps, x.length, own=false))

        self.markers = Matrix{UInt8}(undef, MARKER_SIZE, x.length)
        copyto!(self.markers, unsafe_wrap(Vector{UInt8}, x.markers, x.length * MARKER_SIZE, own=false))

        self.wavemarks = Matrix{Int16}(undef, x.npt, x.length)
        copyto!(self.wavemarks, unsafe_wrap(Vector{Int16}, x.wavemarks, x.length * x.npt, own=false))

        return self
    end
end
# =========================================================================== #
struct cSMRContChannel <: SMRCType
    length::UInt64
    sampling_rate::Float64
    data::Ptr{Int16}
end

mutable struct SMRContChannel <: SMRType
    data::Vector{Int16}
    sampling_rate::Float64

    function SMRContChannel(x::cSMRContChannel)
        self = new()

        self.data = Vector{Int16}(undef, x.length)
        copyto!(self.data, unsafe_wrap(Vector{Int16}, x.data, x.length, own=false))

        self.sampling_rate = x.sampling_rate

        return self
    end
end
# =========================================================================== #
struct cSMREventChannel <: SMRCType
    length::UInt64
    data::Ptr{Float64}
end

mutable struct SMREventChannel <: SMRType
    data::Vector{Float64}

    function SMREventChannel(x::cSMREventChannel)
        self = new()

        self.data = Vector{Float64}(undef, x.length)
        copyto!(self.data, unsafe_wrap(Vector{Float64}, x.data, x.length, own=false))

        return self
    end
end
# =========================================================================== #
struct cSMRMarkerChannel <: SMRCType
    length::UInt64
    npt::UInt64
    timestamps::Ptr{Float64}
    markers::Ptr{UInt8}
    text::Ptr{UInt8}
end

mutable struct SMRMarkerChannel <: SMRType
    timestamps::Vector{Float64}
    markers::Matrix{UInt8}
    text::Vector{String}

    function SMRMarkerChannel(mrk::cSMRMarkerChannel)
        self = new()

        self.timestamps = Vector{Float64}(undef, mrk.length)
        copyto!(self.timestamps, unsafe_wrap(Vector{Float64}, mrk.timestamps, mrk.length, own=false))

        self.markers = Matrix{UInt8}(undef, MARKER_SIZE, mrk.length)
        copyto!(self.markers, unsafe_wrap(Vector{UInt8}, mrk.markers, mrk.length * MARKER_SIZE, own=false))

        self.text = Vector{String}(undef, mrk.length)

        ary = unsafe_wrap(Vector{UInt8}, mrk.text, mrk.length * mrk.npt, own=false)
        for k = 1:mrk.length
            isrt = ((k-1)*mrk.npt)+1
            iend = k*mrk.npt
            if all(x->x=='\0', ary[isrt:iend])
                self.text[k] = ""
            else
                self.text[k] = strip(join(map(Char, ary[isrt:iend])), '\0')
            end
        end

        return self
    end
end
# =========================================================================== #
struct cSMRChannelInfo <: SMRCType
    title::Cstring
    index::Int32
    kind::UInt8
    phy_chan::Int16
end

struct cSMRChannelInfoArray <: SMRCType
    length::UInt32
    ifo::Ptr{Ptr{cSMRChannelInfo}}
end

mutable struct SMRChannelInfo <: SMRType
    title::String
    index::Int
    kind::Int
    phy_chan::Int

    function SMRChannelInfo(x::cSMRChannelInfo)
        self = new()
        self.title = unsafe_string(x.title)
        self.index = Int(x.index)
        self.kind = Int(x.kind)
        self.phy_chan = Int(x.phy_chan)
        return self
    end
end
# =========================================================================== #
function show(io::IO, ifo::SMRChannelInfo)
    str = "\"" * ifo.title *
          "\": " * channel_string(ifo.kind) *
          ", " * string(ifo.index)
    print(io, str)
end
# =========================================================================== #
function channel_string(x::T) where T<:Integer
    sx = string(x)
    if x == 1
        str = "continuous"
    elseif x == 2
        str = "event"
    elseif x == 3
        str = "event"
    elseif x == 4
        str = "event"
    elseif x == 5
        str = "marker"
    elseif x == 6
        str = "adc_marker"
    elseif x == 7
        str = "real_marker"
    elseif x == 8
        str = "text_marker"
    elseif x == 9
        str = "real_wave"
    else
        @warn(string(x) * " is not a valid channel type")
        str = "INVALID"
    end
    return str
end
# =========================================================================== #
end #END MODULE

max_stars_repo_name: scottiealexander/libsmr

max_issues_repo_name: scottiealexander/libsmr

max_forks_repo_name: scottiealexander/libsmr


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
fib_l =
((:lambda, (:u,), (:u, :u)),
    (:lambda, (:u,),
        (:lambda, (:n, :a, :b),
            (:if, (:(==), :n, 0), :a,
                ((:u, :u), (:-, :n, 1), :b, (:+, :a, :b))))))

l2j(s) = s
function l2j(s::Tuple)
    s[1] === :if     && return Expr(l2j.(s)...)
    s[1] === :lambda && return Expr(:(->), Expr(:tuple, s[2]...), l2j(s[3]))
    Expr(:call, l2j.(s)...)
end
macro l(x) l2j(Core.eval(__module__, x)) end

fib= @l fib_l
@show methods(fib).ms[1]
@show fib(21, 0, 1);

# %%
f_l =
((:lambda, (:u,), (:u, :u)),
    (:lambda, (:u,),
        (:lambda, (:n,),
            (:if, (:(≤), :n, 0), (:vcat,),
                (:vcat, ((:u, :u), (:-, :n, 1)), :n)))))

f = @l f_l
f(10.0)

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
module  A2

#this is useful if we have mdoules
#=
weave(Pkg.dir("$(pwd())\\A2.jl"),
  informat="script",
  out_path = "$(pwd())\\A2.html",
  doctype = "md2html")
=#
if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using DataFrames, Distributions, StatsBase, GZip, JLD,
  Gadfly, CTModCopy, JuMP, NLopt, ForwardDiff

const DATA_PATH = pwd() * "\\data" #path where we store the data
const OUTPUT_PATH = pwd() * "\\output"
const DATE_FORMAT_STR = "m/d/y"
const SP500_NAME = "GSPC"
const FOOTER_NAME = "footer.tex" #these are just headers for the tex tables
const HEADER_NAME = "header.tex"

const DEFAULT_ITER = 12000
const DEFAULT_BURN = 2000
const DO_NOT_WRITE = false
const PlotContainer = Union{Plot,Gadfly.Compose.Context}

function preProcessSP500()::Void
  dateFormat::String = DATE_FORMAT_STR
  sp500DF::DataFrame = readtable("$DATA_PATH\\$SP500_NAME.csv")

  #fix the dates
  rename!(sp500DF, :Date, :DATE)
  sp500DF[:,:DATE] = ((s::String)->Date(s, dateFormat)::Date).(sp500DF[:,:DATE])::DataVector{Date}
  sort!(sp500DF, cols = [:DATE])
  sp500DF[:,:returns] = similar(sp500DF[:,:Adj_Close])
  sp500DF[2:end,:returns] = log.(sp500DF[2:end,:Adj_Close] ./  sp500DF[1:(end-1),:Adj_Close])
  #showcols(sp500DF)

  stream::IOStream = open("$DATA_PATH\\$SP500_NAME.jls", "w")
  serialize(stream, sp500DF)
  close(stream)

  return nothing
end

#contains all prior parameters
struct GibbsPriors
  priorNames::Vector{Symbol}
  priorValues::Vector{Float64}
  index::Dict
end

#constructor which forms the index
GibbsPriors(priorNames::Vector{Symbol}, priorValues::Vector{Float64})::GibbsPriors =
            GibbsPriors(priorNames, priorValues,
              Dict(priorNames[i]=>priorValues[i] for i::Int ∈ 1:length(priorNames)))

#make the priors
function GibbsPriors(priorTuples::Vector{Tuple{Symbol,Float64}})::GibbsPriors
  priorNames = Vector{Symbol}(length(priorTuples))
  priorValues = Vector{Float64}(length(priorTuples))

  for i::Int ∈ 1:length(priorTuples)
    priorNames[i], priorValues[i] = priorTuples[i]
  end

  return GibbsPriors(priorNames, priorValues)
end

#this contains parameters and all distributions
mutable struct GibbsProblem
  names::Vector{Symbol}
  values::Vector{Vector{Float64}}
  posteriors!::Vector{Function}
  priors::GibbsPriors
  Y::Vector{Float64} #the data
  paramDistribution::Matrix{Vector{Float64}}
  recordParams::Vector{Bool}

  index::Dict
  targetIterations::Int
  burnIn::Int
  iter::Int
  T::Int
end

#constructor for a Gibbs Problem
#Requires all symbols, posterior distribution functions, initial values, and simulation parameters
function GibbsProblem(names::Vector{Symbol},
            posteriors!::Vector{Function},
            priors::GibbsPriors,
            Y::Vector{Float64};
            initialValues::Vector{Vector{Float64}}=1.0/length(names),
            targetIterations::Int=DEFAULT_ITER,
            burnIn::Int = DEFAULT_BURN,
            recordParams::Vector{Bool} = trues(length(names)))::GibbsProblem

  #this is the space where we will store the parameters
  paramDistribution::Matrix{Vector{Float64}} =
    Matrix{Vector{Float64}}(length(initialValues), targetIterations - burnIn)

  #pre-allocate the space
  for c::Int ∈ 1:size(paramDistribution, 2)
    for r::Int ∈ 1:length(initialValues)
      if recordParams[r]
        paramDistribution[r,c] = similar(initialValues[r])
      end
    end
  end

  index::Dict =   Dict(names[i] => i for i::Int ∈ 1:length(initialValues))

  return GibbsProblem(names, initialValues, posteriors!, priors, Y,
  paramDistribution, #pre-allocated space to store realizations
  recordParams,
  index, #an index of the variable names
  targetIterations,
  burnIn,
  1, #start the iteration counter
  length(Y))
end

#realizes a complete set of parameters
function IterateProblem!(p::GibbsProblem)::Void

  for i::Int ∈ 1:length(p.names)
    p.posteriors![i](p)

    if p.iter > p.burnIn && p.recordParams[i]
      p.paramDistribution[i, p.iter-p.burnIn] .= p.values[i]
    end
  end

  p.iter += 1
  return nothing
end

#performans the Gibbs sampling exercise
function getParameterDistribution!(p::GibbsProblem)::Void
  @fastmath for i::Int ∈ 1:p.targetIterations
    IterateProblem!(p)
  end

  return nothing
end

getSeriesLabels(N::Int, prefix::Symbol = :Y)::Vector{Symbol} = [Symbol(prefix,i) for i::Int ∈ 1:N]

#easy access of values
get(p::GibbsProblem, s::Symbol)::Vector{Float64} = p.values[p.index[s]]
get(priors::GibbsPriors, s::Symbol)::Float64 = priors.index[s]

function getResult(p::GibbsProblem, s::Symbol)::Vector{Float64}

  valIndex::Int = p.index[s]
  result = Vector{Float64}(p.targetIterations - p.burnIn)

  for i::Int ∈ 1:length(result)
    result[i] = p.paramDistribution[valIndex, i][1]
  end

  return result
end

getResults(p::GibbsProblem, names::Vector{Symbol} = p.names)::Dict =
  Dict(s => getResult(p,s) for s::Symbol ∈ names)

#helper function to draw from the posterior mean of a normal distribution
function μPosterior!(p::GibbsProblem)::Void
  ##Priors
  θ::Float64 = get(p.priors, :θ)
  δ²::Float64 = get(p.priors, :δ²)

  #parameters
  σ²::Float64 = get(p,:σ²)[1]
  Z::Vector{Float64} = get(p, :Z)
  ξ::Vector{Float64} = get(p, :ξ)
  pμ::Vector{Float64} = get(p,:μ)

  Δ²::Float64 = 1.0/(p.T/σ²+1.0/δ²)

  #intermediate variables
  YₜMinusξₜZₜOverσ²::Float64 = sum(p.Y .-  (ξ .* Z))/σ²

  #form the distribution
  pμ[1] = rand(Normal(Δ²*(YₜMinusξₜZₜOverσ²+θ/δ²), Δ²^0.5))

  return nothing
end

function μₛPosterior!(p::GibbsProblem)::Void
  ##Priors
  θₛ::Float64 = get(p.priors, :θₛ)
  δ²ₛ::Float64 = get(p.priors, :δ²ₛ)

  #parameters
  σ²ₛ::Float64 = get(p,:σ²ₛ)[1]
  ξ::Vector{Float64} = get(p, :ξ)
  pμₛ::Vector{Float64} = get(p,:μₛ)

  Δ²::Float64 = 1.0/(p.T/σ²ₛ+1.0/δ²ₛ)

  #intermediate variables
  ξOverσ²::Float64 = sum(ξ)/σ²ₛ

  #form the distribution
  pμₛ[1] = rand(Normal(Δ²*(ξOverσ²+θₛ/δ²ₛ), Δ²^0.5))

  return nothing
end

#helper function to draw from the psoterior variance of a normal distribution
function σ²Posterior!(p::GibbsProblem)::Void
  ##Priors
  α::Float64 = get(p.priors, :α)
  β::Float64 = get(p.priors, :β)

  #parameters
  μ::Float64 = get(p, :μ)[1]
  Z::Vector{Float64} = get(p, :Z)
  ξ::Vector{Float64} = get(p, :ξ)
  pσ²::Vector{Float64} = get(p, :σ²)

  #intermediate variables
  YtMinusμ²::Float64 = sum((p.Y .- μ .- (ξ .* Z)).^2.0)

  #form the distribution
  pσ²[1] = rand(InverseGamma(p.T/2.+α,0.5*(YtMinusμ²+2.0*β)))
  #draw
  return nothing
end

#helper function to draw from the psoterior variance of a normal distribution
function σ²ₛPosterior!(p::GibbsProblem)::Void
  ##Priors
  αₛ::Float64 = get(p.priors, :αₛ)
  βₛ::Float64 = get(p.priors, :βₛ)

  #parameters
  μₛ::Float64 = get(p, :μₛ)[1]
  ξ::Vector{Float64} = get(p, :ξ)
  pσ²ₛ::Vector{Float64} = get(p, :σ²ₛ)

  #intermediate variables
  YtMinusμ²ₛ::Float64 = sum((ξ .- μₛ).^2.0)

  #form the distribution
  pσ²ₛ[1] = rand(InverseGamma(p.T/2.+αₛ,0.5*(YtMinusμ²ₛ+2.0*βₛ)))
  #draw
  return nothing
end

#helper function to draw λ
function λPosterior!(p::GibbsProblem)::Void
  ##Priors
  γ::Float64 = get(p.priors, :γ)
  η::Float64 = get(p.priors, :η)

  #parameters
  Z::Vector{Float64} = get(p, :Z)
  pλ::Vector{Float64} = get(p, :λ)

  #intermediate variables
  ΣZ::Float64 = sum(Z)

  #form the distribution and draw
  pλ[1] = rand(Beta(ΣZ+γ,p.T-ΣZ+η))

  return nothing
end

#helper function to draw ξₜ
function ξPosterior!(p::GibbsProblem)::Void
  #parameters
  μ::Float64 = get(p,:μ)[1]
  μₛ::Float64 = get(p,:μₛ)[1]
  σ²::Float64 = get(p,:σ²)[1]
  σ²ₛ::Float64 = get(p,:σ²ₛ)[1]

  Z::Vector{Float64} = get(p,:Z)

  pξ::Vector{Float64} = get(p,:ξ)

  for t::Int ∈ 1:p.T
    Δ²::Float64 = (1.0 / (Z[t] / σ² + 1.0 / σ²ₛ))

    pξ[t] =  rand(Normal(Δ² * ((p.Y[t] - μ) * Z[t] / σ² + μₛ/σ²ₛ), Δ²^0.5))
  end

  return nothing
end

#helper function to draw Zₜ
function ZPosterior!(p::GibbsProblem)::Void

  μ::Float64 = get(p,:μ)[1]
  σ²::Float64 = get(p,:σ²)[1]
  λ::Float64 = get(p,:λ)[1]
  ξ::Vector{Float64} = get(p,:ξ)

  pZ::Vector{Float64} = get(p,:Z)

  for t::Int ∈ 1:p.T
    propZ1::Float64 = λ * exp(-0.5 * (p.Y[t] - μ - ξ[t])^2.0 / σ²)
    propZ0::Float64 = (1.0 - λ) * exp(-0.5 * (p.Y[t] - μ)^2.0 / σ²)
    pZ[t] = rand(Bernoulli(propZ1/(propZ0+propZ1)))
  end

  return nothing
end

function problemA2Script(;refreshData::Bool = true, doNotWrite::Bool = DO_NOT_WRITE)::Void
  if !isfile("$DATA_PATH\\$SP500_NAME.jls") || refreshData
    preProcessSP500()
  end

  #read the binary file
  stream::IOStream = open("$DATA_PATH\\$SP500_NAME.jls")
  sp500DF::DataFrame = deserialize(stream)
  close(stream)
  Y::Vector{Float64} = sp500DF[2:end, :returns]
  T::Int = length(Y)

  #################Priors
  priors = GibbsPriors([
  (:θ, mean(Y)), (:δ², var(Y)), # μ~N(θ,δ²) #use std deviation as upper bound for δ
  (:α, 0.001), (:β, 0.001), #σ²~IG(α,β)
  (:θₛ, 0.0 ),  (:δ²ₛ, 1.0),  #μₛ~N(θₛ,δ²ₛ)
  (:αₛ, 0.001 ), (:βₛ, 0.001), #σ²ₛ~IG(αₛ,βₛ)
  (:γ, 1.0 ), (:η, 1.0)]) #λ~β(γ,η) (use uniform distribution)

  #Form posterior distributions
  names::Vector{Symbol} = Vector{Symbol}()
  posteriors!::Vector{Function} = Vector{Function}()
  initialValues::Vector{Vector{Float64}} =  Vector{Float64}()

  #assign variable labels, posterior functions, and initial values
  push!(names, :μ)
  push!(posteriors!, μPosterior!)
  push!(initialValues, [0.0])

  push!(names, :μₛ)
  push!(posteriors!, μₛPosterior!)
  #push!(posteriors!, μₛPosterior)
  push!(initialValues, [0.0])

  push!(names, :σ²)
  push!(posteriors!, σ²Posterior!)
  push!(initialValues, [0.01])

  push!(names, :σ²ₛ)
  push!(posteriors!, σ²ₛPosterior!)
  push!(initialValues, [0.01])

  push!(names, :λ)
  push!(posteriors!, λPosterior!)
  push!(initialValues, [0.01])

  push!(names, :ξ)
  push!(posteriors!, ξPosterior!)
  push!(initialValues, zeros(T))

  push!(names, :Z)
  push!(posteriors!, ZPosterior!)
  push!(initialValues, zeros(T))

  println("""Beginning Gibbs Problem
    targetIterations=$DEFAULT_ITER,
    burnIn = $DEFAULT_BURN
    length Of Values = $(length(initialValues))
    length Of posteriors! = $(length(posteriors!))
    length Of Names = $(length(names))

  """)

  p::GibbsProblem = GibbsProblem(names, posteriors!, priors, Y,
              initialValues=initialValues,
              targetIterations=DEFAULT_ITER,
              burnIn = DEFAULT_BURN,
              recordParams = [true, true, true, true, true, false, false])

  @time getParameterDistribution!(p)

  results::Dict = getResults(p, [:μ, :σ², :μₛ, :σ²ₛ, :λ])

  μMean::Float64 = mean(results[:μ])
  σ²Mean::Float64 = mean(results[:σ²])
  μₛMean::Float64 = mean(results[:μₛ])
  σ²ₛMean::Float64 = mean(results[:σ²ₛ])
  λMean::Float64 = mean(results[:λ])

  println("""
  Results (Iterations=$(p.targetIterations), Burnin=$(p.burnIn)): (Mean, annualized)
    μ: $(round(exp(μMean)-1.0,6))
    σ: $(round(σ²Mean^0.5,6))
    μₛ: $(round(exp(μₛMean)-1.0,6))
    σₛ: $(round(σ²ₛMean^0.5,6))
    λ: $(round(λMean,6))

    μ (Ann.): $(round(exp(255.0*μMean)-1.0,6))
    σ (Ann.): $(round((255.0 * σ²Mean)^0.5,6))
  """)

  if !doNotWrite
    gstream::GZipStream = gzopen("$DATA_PATH\\A2Out_$(p.targetIterations)_$(p.burnIn).jls.gz", "w")
    serialize(gstream, results)
    close(gstream)
  end

  return nothing
end

function getA2Graphs(;targetIterations::Int = DEFAULT_ITER, burnIn::Int = DEFAULT_BURN)::Void

  gstream::GZipStream = gzopen("$DATA_PATH\\A2Out_$(targetIterations)_$(burnIn).jls.gz")
  results::Dict =  deserialize(gstream)
  close(gstream)

  stream::IOStream = open("$DATA_PATH\\$SP500_NAME.jls")
  sp500DF::DataFrame = deserialize(stream)
  close(stream)

  numPoints::Int = targetIterations-burnIn

  push!(results, :σ=>results[:σ²].^0.5)
  push!(results, :σₛ=>results[:σ²ₛ].^0.5)

  Gadfly.push_theme(:default)
  plots::Vector{PlotContainer} = Vector{PlotContainer}()
  plotNames::Vector{String} = Vector{String}()

  push!(plotNames, "IterationPlots1")
  push!(plots, vstack(
    plot(x=1:numPoints, y=results[:μ], Geom.line, Theme(line_width = 0.25pt),
      Guide.ylabel("μ"), Guide.xlabel("Iteration"),
      Guide.xticks(ticks=collect(0:(numPoints÷5):(numPoints))),
      Guide.title("Draws of μ")),
    plot(x=1:numPoints, y=(results[:σ]), Geom.line, Theme(line_width = 0.2pt),
      Guide.ylabel("σ"), Guide.xlabel("Iteration"),
      Guide.xticks(ticks=collect(0:(numPoints÷5):(numPoints))),
      Guide.title("Draws of σ"))))

  push!(plotNames, "IterationPlots2")
  push!(plots, vstack(
    plot(x=1:numPoints, y=results[:μₛ], Geom.line, Theme(line_width = 0.2pt),
      Guide.ylabel("μₛ"), Guide.xlabel("Iteration"),
      Guide.xticks(ticks=collect(0:(numPoints÷5):(numPoints))),
      Guide.title("Draws of μₛ")),
    plot(x=1:numPoints, y=(results[:σₛ]), Geom.line, Theme(line_width = 0.2pt),
      Guide.ylabel("σₛ"), Guide.xlabel("Iteration"),
      Guide.xticks(ticks=collect(0:(numPoints÷5):(numPoints))),
      Guide.title("Draws of σₛ")),
    plot(x=1:numPoints, y=(results[:λ]), Geom.line, Theme(line_width = 0.2pt),
      Guide.ylabel("λ"), Guide.xlabel("Iteration"),
      Guide.xticks(ticks=collect(0:(numPoints÷5):(numPoints))),
      Guide.title("Draws of λ"))))

  push!(plotNames, "HistogramPlots1")
  push!(plots, vstack(
    plot(x=results[:μ], Geom.histogram, Theme(line_width = 0.2pt),
      Guide.ylabel("μ"), Guide.xlabel("Iteration"),
      Guide.title("Distribution of μ")),
    plot(x=(results[:σ]), Geom.histogram, Theme(line_width = 0.2pt),
      Guide.ylabel("σ"), Guide.xlabel("Iteration"),
      Guide.title("Distribution of :σ"))))

  push!(plotNames, "HistogramPlots2")
  push!(plots, vstack(
    plot(x=results[:μₛ], Geom.histogram, Theme(line_width = 0.2pt),
      Guide.ylabel("μₛ"), Guide.xlabel("Iteration"),
      Guide.title("Distribution of μₛ")),
    plot(x=results[:σₛ], Geom.histogram, Theme(line_width = 0.2pt),
      Guide.ylabel("σₛ"), Guide.xlabel("Iteration"),
      Guide.title("Distribution of σₛ")),
    plot(x=results[:λ], Geom.histogram, Theme(line_width = 0.2pt),
      Guide.ylabel("λ"), Guide.xlabel("Iteration"),
      Guide.title("Distribution of λ"))))

  for i ∈ 1:length(plots)
    draw(SVG("$OUTPUT_PATH\\$(plotNames[i]).svg", 7.5inch, 9inch),plots[i])
  end

  ########################make the summary data table
  summaryDecimals::Int = 6
  summaryCols::Vector{Symbol} = [:μ, :σ, :μₛ,:σₛ, :λ]
  summaryRowNames::Vector{String} =
    ["mean", "median", "Ann. Median", "Std Dev", "5th Percentile",
      "95th Percentile", "N"]
  summaryNumRows::Int = length(summaryRowNames)

  annualizedMedians::Vector{Float64} = [exp(median(results[:μ])*255)-1.0,
    median(results[:σ])*255^0.5,
    exp(median(results[:μₛ])*255)-1.0,
    median(results[:σₛ])*255^0.5,
    round(median(results[:λ])*255,1)]

  #make the content matrix
  summaryContent::Vector{Vector{String}} =
    [Vector{String}() for i∈ 1:summaryNumRows]

  summaryContent[1] = string.([
    round(mean(results[s]),summaryDecimals) for s::Symbol ∈ summaryCols])

  summaryContent[2] = string.([
    round(median(results[s]),summaryDecimals) for s::Symbol ∈ summaryCols])

  summaryContent[3] = string.([
    round(annualizedMedians[i],summaryDecimals) for i::Int ∈ 1:length(summaryCols)])

  summaryContent[4] = string.([
    round(std(results[s]),summaryDecimals) for s::Symbol ∈ summaryCols])

  summaryContent[5] = string.([
    round(quantile(results[s], 0.05),summaryDecimals) for s::Symbol ∈ summaryCols])

  summaryContent[6] = string.([
    round(quantile(results[s], 0.95),summaryDecimals) for s::Symbol ∈ summaryCols])

  summaryContent[7] = ["$numPoints"]

  tableSummary::String = texTable( "Summary of Data",
    """See Tex File""", #caption
    [["\$\\mu\$", "\$\\sigma\$", "\$\\mu_{s}\$","\$\\sigma_{s}\$", "\$\\lambda\$"]], #colNames
    Vector{String}(),#contentRowNames
    Vector{Matrix{String}}(), #content
    summaryRowNames, #descRowNames
    summaryContent, #descContent
    Vector{String}(),
    columnSepPt = -5,
    widthDescContent = [[ones(Int, length(summaryCols)) for i∈ 1:(length(summaryContent)-1)];
      [[length(summaryCols)]]]::Vector{Vector{Int}})

  writeTables2File([tableSummary],
    HEADER_NAME, FOOTER_NAME, path=OUTPUT_PATH,
    outName = "SummaryTable.tex")

  return nothing
end

#enable these to run the program
problemA2Script(refreshData=false)
#getA2Graphs()
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
import PortAudio, FileIO
using GR
using FFTW

function play(c::Channel)
  data, = FileIO.load("Monty_Python.wav")
  stream = PortAudio.PortAudioStream(0, 1)
  spectrum = zeros(Int32, 300, 225)

  offset = 1
  while offset + 1024 < length(data)
    amplitudes = data[offset:offset+1024]
    offset += 1024

    PortAudio.write(stream, amplitudes)

    power = log.(abs.(fft(float(amplitudes))) .+ 1) * 50
    spectrum[1, :] = round.(Int, power[1:225])
    spectrum = circshift(spectrum, [-1, 0])

    put!(c, spectrum')
  end
  put!(c, nothing)
end

function main()
  c = Channel(0)
  task = @task play(c::Channel)
  bind(c, task)
  schedule(task)
  start = time_ns()

  while isopen(c)
    spectrum = take!(c)
    if spectrum === nothing
      break
    end

    if time_ns() - start > 20 * 1000000   # 20ms
      imshow(spectrum, colormap=-13, yflip=true)

      start = time_ns()
    end
  end
end

main()


text:
# ImagePFM.jl
# FileIO.jl integration for PFM image files.
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

"""
    ImagePFM

FileIO.jl integration for PFM image files.
"""
module ImagePFM

using Base: Float32
import Base:
    iterate, read, write

import ColorTypes:
    RGB

using FileIO:
    @format_str,
    DataFormat, File, Stream,
    query, skipmagic, stream

export
    InvalidPfmFileFormat, InvalidRGBEltype

include("exceptions.jl")
include("implementation.jl")

#=
Here we implement the `load` and `save` methods for the PFM data format that will be used by `FileIO`
For more information consult the `FileIO` documentation on adding new formats and implementing load/save behavior
=#

function load(f::File{DataFormat{:PFM}}; kwargs...)
    open(f) do s
        load(s, kwargs...)
    end
end
function load(s::Stream{FMT}; permute_horizontal=false, kwargs...) where {FMT<:DataFormat{:PFM}}
    if permute_horizontal
        return read(stream(s), FMT, kwargs...)
    else
        return permutedims(read(stream(s), FMT, kwargs...), (2, 1))
    end
end

function save(f::File{DataFormat{:PFM}}, image::AbstractMatrix; kwargs...)
    open(f, "w") do s
        save(s, image, kwargs...)
    end
end
function save(s::Stream{FMT}, image::AbstractMatrix; permute_horizontal=false, mapi=identity, kwargs...) where {FMT<:DataFormat{:PFM}}
    imgout = map(mapi, image)
    if permute_horizontal
        return write(stream(s), FMT, imgout, kwargs...)
    else
        perm = ndims(imgout) == 2 ? (2, 1) : error("$(ndims(imgout)) dims array is not supported")
        return write(stream(s), FMT, PermutedDimsArray(imgout, perm), kwargs...)
    end
end

end # module ImagePFM


text:
module BenchDistance
using BenchmarkTools
using ParallelKMeans
using Distances
using Random

suite = BenchmarkGroup()

Random.seed!(2020)
X = rand(3, 100_000)
centroids = rand(3, 2)
d = Vector{Float64}(undef, 100_000)
suite["100kx3"] = @benchmarkable ParallelKMeans.colwise!($d, $X, $centroids)

X = rand(10, 100_000)
centroids = rand(10, 2)
d = Vector{Float64}(undef, 100_000)
suite["100kx10"] = @benchmarkable ParallelKMeans.colwise!($d, $X, $centroids)

end # module

BenchDistance.suite


text:
# Base32 Encoder
# ==============

struct Base32Encoder <: Codec
    table::CodeTable32
    state::State
    buffer::Buffer
end

Base32Encoder(table::CodeTable32) = Base32Encoder(table, State(), Buffer(4))

"""
    Base32Encoder(;hex::Bool=false)

Create a base32 encoding codec.

Arguments
---------
- `hex`: use extended hex alphabet (Table 4 of RFC4648)
"""
function Base32Encoder(;hex::Bool=false)
    if hex
        table = BASE32_HEX
    else
        table = BASE32_STD
    end
    return Base32Encoder(table)
end

const Base32EncoderStream{S} = TranscodingStream{Base32Encoder,S} where S<:IO

"""
    Base32EncoderStream(stream::IO; kwargs...)

Create a base32 encoding stream (see `Base32Encoder` for `kwargs`).
"""
function Base32EncoderStream(stream::IO; kwargs...)
    return TranscodingStream(Base32Encoder(;kwargs...), stream)
end

function TranscodingStreams.startproc(
        codec :: Base32Encoder,
        state :: Symbol,
        error :: Error)
    start!(codec.state)
    empty!(codec.buffer)
    return :ok
end

function TranscodingStreams.process(
        codec  :: Base32Encoder,
        input  :: Memory,
        output :: Memory,
        error  :: Error)
    table = codec.table
    state = codec.state
    buffer = codec.buffer

    # Check if we can encode data.
    if !is_running(state)
        error[] = ArgumentError("encoding is already finished")
        return 0, 0, :error
    elseif output.size < 8
        # Need more output space.
        return 0, 0, :ok
    end

    # Load the first five bytes.
    i = j = 0
    while buffer.size < 4 && i < input.size
        buffer[buffer.size+=1] = input[i+=1]
    end
    b1 = b2 = b3 = b4 = b5 = 0x00
    npad = 0
    status = :ok
    if i < input.size
        b1 = buffer[1]
        b2 = buffer[2]
        b3 = buffer[3]
        b4 = buffer[4]
        b5 = input[i+=1]
    elseif input.size == 0
        # Found the end of the input.
        if buffer.size == 0
            finish!(state)
            return i, j, :end
        elseif buffer.size == 1
            b1 = buffer[1]
            npad = 6
        elseif buffer.size == 2
            b1 = buffer[1]
            b2 = buffer[2]
            npad = 4
        elseif buffer.size == 3
            b1 = buffer[1]
            b2 = buffer[2]
            b3 = buffer[3]
            npad = 3
        elseif buffer.size == 4
            b1 = buffer[1]
            b2 = buffer[2]
            b3 = buffer[3]
            b4 = buffer[4]
            npad = 1
        else
            @unreachable
        end
        status = :end
    else
        # Need more data to encode.
        return i, j, :ok
    end
    empty!(buffer)

    # Encode the body.
    #  01234567 89012345 67890123 45678901 23456789
    # +--------+--------+--------+--------+--------+
    # |< 1 >< 2| >< 3 ><| 4 >< 5 |>< 6 >< |7 >< 8 >|
    # +--------+--------+--------+--------+--------+
    @inbounds while true
        output[j+1] = encode(table, b1 >> 3)
        output[j+2] = encode(table, b1 << 2 | b2 >> 6)
        output[j+3] = encode(table, b2 >> 1)
        output[j+4] = encode(table, b2 << 4 | b3 >> 4)
        output[j+5] = encode(table, b3 << 1 | b4 >> 7)
        output[j+6] = encode(table, b4 >> 2)
        output[j+7] = encode(table, b4 << 3 | b5 >> 5)
        output[j+8] = encode(table, b5)
        j += 8
        if i + 5 ≤ input.size && j + 8 ≤ output.size
            b1 = input[i+1]
            b2 = input[i+2]
            b3 = input[i+3]
            b4 = input[i+4]
            b5 = input[i+5]
            i += 5
        else
            break
        end
    end

    # Epilogue.
    while npad > 0
        output[j-npad+1] = table.padcode
        npad -= 1
    end
    if status == :end || status == :error
        finish!(state)
    end
    return i, j, status
end


text:
using DataStructures: Deque

export AbstractStateSpaceSparse, StateSpaceSparse, expand!, deleteat!, get_state_count, get_sink_count, get_stoich_matrix, get_statedict, get_state_connectivity, get_sink_connectivity

"""
$(TYPEDEF)
"""
abstract type AbstractStateSpaceSparse{NS,NR,IntT<:Integer,SizeT<:Integer} <: AbstractStateSpace end

"""
    $(TYPEDEF)

Basic Sparse FSP State space.

# Fields

$(FIELDS)

# See also
`expand!`,`deleteat!`
"""
Base.@kwdef mutable struct StateSpaceSparse{NS,NR,IntT<:Integer,SizeT<:Integer} <: AbstractStateSpaceSparse{NS,NR,IntT,SizeT}
    "Stoichiometry matrix S = [s₁ ... sₘ] of size N x M where N is the number of species, M the number of the reactions."
    stoich_matrix::Matrix{IntT}

    "Number of sinks"
    sink_count::SizeT

    "Array of CME states included in the subspace"
    states::Vector{MVector{NS,IntT}}

    "Dictionary of states, containing pairs `(x=>i)` for (i,x) in enumerate(states). The implementation must ensure that each state in `states` is a key in `state2idx` and conversely every key in `state2idx` exists in `states`."
    state2idx::Dict{MVector{NS,IntT},SizeT}

    "List of state connectivity information. `state_connectivity[i][k] = j` if xᵢ = xⱼ + sₖ, that is, `states[i] = states[j] + stoich_mat[:, k]`. If there is no existing state that can reach xᵢ via reaction k, the implementation must ensure that `state_connectivity[i][k] = 0`."
    state_connectivity::Vector{MVector{NR,SizeT}}

    "Matrix to store reaction events by which the included states transit to outside of the projected state space"
    sink_connectivity::Vector{MVector{NR,SizeT}}
end

# Getters
"""
$(TYPEDSIGNATURES)

Return the stoichiometry matrix.
"""
get_stoich_matrix(space::StateSpaceSparse) = space.stoich_matrix

"""
$(TYPEDSIGNATURES)

Return number of states.
"""
get_state_count(statespace::StateSpaceSparse) = length(statespace.states)

"""
$(TYPEDSIGNATURES)

Return number of sinks.
"""
get_sink_count(statespace::StateSpaceSparse) = statespace.sink_count

"""
$(TYPEDSIGNATURES)

Return list of states.
"""
get_states(statespace::StateSpaceSparse) = statespace.states

"""
$(TYPEDSIGNATURES)

Return state dictionary.
"""
get_statedict(statespace::StateSpaceSparse) = statespace.state2idx

get_state_connectivity(statespace::StateSpaceSparse) = statespace.state_connectivity

get_sink_connectivity(statespace::StateSpaceSparse) = statespace.sink_connectivity

# Constructors
"""
$(TYPEDEF)

Construct a basic FSP state space with stoichiometry matrix `stoich_mat` and initial list of states `initstates` with their integer entries being stored in type `IntT <: Integer`.

# Examples

```
julia> S = [[1,0] [-1,0] [0,1] [0,-1]]
2×4 Matrix{Int64}:
 1  -1  0   0
 0   0  1  -1
 julia> states = [[0,1], [10, 1], [0, 10]]
 3-element Vector{Vector{Int64}}:
 [0, 1]
 [10, 1]
 [0, 10]
julia> StateSpaceSparse(S, states)
```
"""
function StateSpaceSparse(stoich_matrix::Matrix{IntT}, initstates::Vector; index_type::Type{SizeT} = UInt32) where {IntT<:Integer,SizeT<:Integer}
    species_count = size(stoich_matrix, 1)
    reaction_count = size(stoich_matrix, 2)
    sink_count = size(stoich_matrix, 2)

    states = Vector{MVector{species_count,IntT}}()
    state2idx = Dict{MVector{species_count,IntT},IntT}()
    state_connectivity = Vector{MVector{reaction_count,IntT}}()
    sink_connectivity = Vector{MVector{reaction_count,IntT}}()

    fspstatespace = StateSpaceSparse{species_count,reaction_count,IntT,index_type}(
        stoich_matrix = stoich_matrix,
        states = states,
        state2idx = state2idx,
        sink_count = sink_count,
        state_connectivity = state_connectivity,
        sink_connectivity = sink_connectivity,
    )

    _addstates!(fspstatespace, initstates)
    fspstatespace
end

"""
`StateSpaceSparse(stoich_mat::Matrix{IntT}, init_state::Vector{IntT}; index_type::Type{<:Integer}=UInt32)`

Construct a basic FSP state space with stoichiometry matrix `stoich_mat` and a single state `init_state` with their integer entries being stored in type `IntT <: Integer`. The optional keyword argument `index_type` allows for more customization on internal indexing representations.

# Examples

```jldoctest
julia> S = [[1,0] [-1,0] [0,1] [0,-1]]
2×4 Matrix{Int64}:
 1  -1  0   0
 0   0  1  -1
 julia> x0 = [0,1]
julia> StateSpaceSparse(S, states)
````
"""
function StateSpaceSparse(stoich_matrix::Matrix{IntT}, init_state::Vector{IntT}; index_type::Type{SizeT} = UInt32) where {IntT<:Integer,SizeT<:Integer}
    return StateSpaceSparse(stoich_matrix, [init_state])
end

"""
$(TYPEDSIGNATURES)

Expand the FSP state space to include all states that are reachable from the existing states in `expansionlevel` or fewer reaction events. State exploration direction is restricted to only reaction channels within `onlyreactions` if this argument is non-empty.
"""
function expand!(statespace::StateSpaceSparse{NS,NR,IntT,SizeT}, expansionlevel::Integer; onlyreactions = []) where {NS,NR,IntT<:Integer,SizeT<:Integer}
    if expansionlevel <= 0
        return
    end
    # Reference to mutable fields within stateset
    stoich_matrix = get_stoich_matrix(statespace)
    states = get_states(statespace)
    sink_connectivity = get_sink_connectivity(statespace)

    reaction_count = size(stoich_matrix, 2)
    expandreactions = isempty(onlyreactions) ? (1:reaction_count) : onlyreactions

    explorables = Deque{SizeT}()
    for idx in 1:length(states)
        for ir in expandreactions
            if sink_connectivity[idx][ir] != 0
                push!(explorables, idx)
                break
            end
        end
    end

    candidate_states = Vector{MVector{NS,IntT}}()
    for level in 1:expansionlevel
        resize!(candidate_states, length(explorables) * length(expandreactions))
        i = 1
        while !isempty(explorables)
            idx = pop!(explorables)
            @simd for ir in expandreactions
                candidate_states[i] = states[idx] .+ stoich_matrix[:, ir]
                i += 1
            end
        end
        statecount_old = length(states)
        _addstates!(statespace, candidate_states)
        statecount = length(states)
        for idx in statecount_old+1:statecount
            push!(explorables, idx)
        end
    end
    nothing
end

# Helper function. Returns `true` if all elements of a vector is positive. Otherwise returns `false`.
function _is_all_nonnegative(x::AbstractVector)::Bool
    for xi in x
        xi < 0 && return false
    end
    true
end

#=
Helper function to add a state vector `newstate` to the current state space `statespace`.
If `newstate` is already included in `statespace`, this function will return without modifying any of the input arguments.
=#
function _addstates!(statespace::StateSpaceSparse{NS,NR,IntT,SizeT}, newstates::Vector{VT}) where {NS,NR,IntT<:Integer,VT<:Union{Vector{IntT},MVector{NS,IntT}},SizeT<:Integer}
    states = get_states(statespace)
    state2idx = get_statedict(statespace)
    num_sinks = get_sink_count(statespace)
    stoich_matrix = get_stoich_matrix(statespace)
    state_connectivity = get_state_connectivity(statespace)
    sink_connectivity = get_sink_connectivity(statespace)
    reaction_count = size(stoich_matrix, 2)

    unique_newstates = Vector{MVector{NS,IntT}}()
    sizehint!(unique_newstates, length(newstates))
    newidx = length(states)
    for state in newstates
        if get(state2idx, state, 0) == 0 && _is_all_nonnegative(state)
            newidx += 1
            push!(unique_newstates, state)
            state2idx[state] = newidx
        end
    end

    statecount_old = length(states)
    added_state_count = length(unique_newstates)
    append!(states, unique_newstates)
    resize!(state_connectivity, statecount_old + added_state_count)
    resize!(sink_connectivity, statecount_old + added_state_count)
    statecount = length(states)

    @simd for newidx in statecount_old+1:statecount
        @inbounds state_connectivity[newidx] = zeros(IntT, reaction_count)
        @inbounds sink_connectivity[newidx] = zeros(IntT, reaction_count)
    end

    reachablestate = similar(states[1])
    @simd for newidx in statecount_old+1:statecount
        # Determine connectivity between states
        for ir = 1:reaction_count
            copy!(reachablestate, states[newidx])
            axpy!(-1, stoich_matrix[:, ir], reachablestate)
            # reachablestate = states[newidx] - stoich_matrix[:, ir]
            ridx = get(state2idx, reachablestate, 0)
            if ridx ≠ 0
                @inbounds state_connectivity[newidx][ir] = ridx
                @inbounds sink_connectivity[ridx][ir] = 0
            end
        end

        # Determine states that can transit outside of the current truncated state space
        for ir = 1:reaction_count
            reachablestate = states[newidx] + stoich_matrix[:, ir]
            if _is_all_nonnegative(reachablestate)
                ridx = get(state2idx, reachablestate, 0)
                if ridx == 0
                    sink_connectivity[newidx][ir] = ir
                else
                    state_connectivity[ridx][ir] = newidx
                end
            end
        end
    end
end

import Base: deleteat!

"""
    $(TYPEDSIGNATURES)

Delete states with indices `ids` from the state space.
"""
function deleteat!(statespace::StateSpaceSparse, ids::Vector{T}) where {T<:Integer}
    stoich_matrix = get_stoich_matrix(statespace)
    reaction_count = size(stoich_matrix, 2)
    statecount_old = get_state_count(statespace)
    ids = unique(ids)

    # Map from current indices to new indices when the targeted states are removed
    newidxs = Vector{T}(1:statecount_old)
    newidxs[ids] .= 0
    idx = 1
    for i in 1:statecount_old
        if newidxs[i] ≠ 0
            newidxs[i] = idx
            idx += 1
        end
    end

    # Remove target states from member state list
    states = get_states(statespace)
    state2idx = get_statedict(statespace)
    state_connectivity = get_state_connectivity(statespace)
    sink_connectivity = get_sink_connectivity(statespace)

    for idx in ids
        delete!(state2idx, states[idx])
    end
    deleteat!(states, ids)
    deleteat!(state_connectivity, ids)
    deleteat!(sink_connectivity, ids)

    if length(states) == 0
        return nothing
    end
    # Reindexing
    reachablestate = similar(states[1])
    @simd for i in 1:length(states)
        @inbounds eachstate = states[i]
        state2idx[eachstate] = newidxs[state2idx[eachstate]]
        for ir in 1:reaction_count
            @inbounds state_connectivity[i][ir] = (oldidx = state_connectivity[i][ir]) ≠ 0 ? newidxs[oldidx] : 0
        end
    end
    @simd for i in 1:length(states)
        eachstate = states[i]
        for ir in 1:reaction_count
            if sink_connectivity[i][ir] == 0
                copyto!(reachablestate, eachstate)
                axpy!(1, stoich_matrix[:, ir], reachablestate)
                if _is_all_nonnegative(reachablestate) && (get(state2idx, reachablestate, 0) == 0)
                    @inbounds sink_connectivity[i][ir] = ir
                end
            end
        end
    end
    nothing
end


text:
module CellularPotts

using OffsetArrays, #Allow some arrays to be zero indexed to include medium
      GLMakie, #For plotting and making the GUI (also exports AbstractPlotting)
      Colors, #More color options for the cells (e.g. :Purples)
      StatsBase, Random, #Currently just used for countmap, inverse_rle; shuffle
      SparseArrays, #For creating large adjacency matrices
      LinearAlgebra, #Additional functionality for arrays
      LightGraphs, #Needed for creating graphs that look like graphDimension
      Metis, #lightning fast method for partitioning graphs (only thing in this package that is not Julia)
      TimerOutputs, #looking for allocations
      Catalyst, #Testing reactions DLS
      DiffEqJump, #Discrete Differential Equations
      Plots #Well...for faster plotting

include("BaseFunctions.jl")
include("CellActions.jl")
include("InitializeCells.jl")
include("MarkovStep.jl")
include("Gui.jl")

export

#BaseFunctions.jl
      Hamiltonian,
      AdhesionPenalty,
      VolumePenalty,
      ModelParameters,
      NamedGraph,
      CellAttributes,
      CellPotts,
      UpdateConnections!,

#MarkovStep.jl
      MHStepInfo,
      MHStep!,
      MHStep_naive!,

#CellActions.jl
      CellDivision!,
      CellDeath!,

#Gui.jl
      CellGUI

end


text:
const CPModification = Dict{String, Union{Array{Int}, Array{Bool}, SetModification}}

"""
    merge!(prunedDomains::CPModification, newPrunedDomains::CPModification)

Merge `newPrunedDomains` into `prunedDomains`, concatenating the arrays if concerning the same variable.
"""
function merge!(prunedDomains::CPModification, newPrunedDomains::CPModification)
    for k in keys(newPrunedDomains)
        if haskey(prunedDomains, k)
            prunedDomains[k] = vcat(prunedDomains[k], newPrunedDomains[k])
        else
            prunedDomains[k] = newPrunedDomains[k]
        end
    end
end

function updateChildren!(prunedDomains::CPModification, x::Union{IntVar, BoolVar}, pruned::Union{Array{Int}, Array{Bool}, BitArray})
    for y in x.children
        childrenPruned = childrenValue.([y], pruned)
        if haskey(prunedDomains, y.id)
            prunedDomains[y.id] = vcat(prunedDomains[y.id], Array(childrenPruned))
        else
            prunedDomains[y.id] = Array(childrenPruned)
        end
    end
    return
end

"""
    addToPrunedDomains!(prunedDomains::CPModification, x::IntVar, pruned::Array{Int})

Update the `CPModification` by adding the pruned integers.

# Arguments
- `prunedDomains::CPModification`: the `CPModification` you want to update.
- `x::IntVar`: the variable that had its domain pruned.
- `pruned::Array{Int}`: the pruned integers.
"""
function addToPrunedDomains!(prunedDomains::CPModification, x::Union{IntVar, BoolVar}, pruned::Union{Array{Int}, Array{Bool}, BitArray})
    if isempty(pruned)
        return
    end

    if haskey(prunedDomains, x.id)
        prunedDomains[x.id] = vcat(prunedDomains[x.id], Array(pruned))
    else
        prunedDomains[x.id] = Array(pruned)
    end

    updateChildren!(prunedDomains, x, pruned)
end

"""
    addToPrunedDomains!(prunedDomains::CPModification, x::IntVar, pruned::Array{Int})

Update the `CPModification` by adding the pruned integers.

# Arguments
- `prunedDomains::CPModification`: the `CPModification` you want to update.
- `x::IntVar`: the variable that had its domain pruned.
- `pruned::Array{Int}`: the pruned integers.
"""
function addToPrunedDomains!(prunedDomains::CPModification, x::Union{IntVarView, BoolVarView}, pruned::Union{Array{Int}, Array{Bool}, BitArray})
    if isempty(pruned)
        return
    end

    # Update the ParentVariable only (it will then update all the children variables)
    parentPruned = parentValue.([x], pruned)
    addToPrunedDomains!(prunedDomains, x.x, parentPruned)
end

"""
    addToPrunedDomains!(prunedDomains::CPModification, x::IntVar, pruned::Array{Int})

Update the `CPModification` by adding modified set values.

# Arguments
- `prunedDomains::CPModification`: the `CPModification` you want to update.
- `x::IntSetVar`: the variable that had its values changed.
- `modification::SetModification`: the modified values.
"""
function addToPrunedDomains!(prunedDomains::CPModification, x::IntSetVar, modification::SetModification)
    if isempty(modification.required) && isempty(modification.excluded)
        return
    end

    if haskey(prunedDomains, x.id)
        mergeSetModifications!(prunedDomains[x.id], modification)
    else
        prunedDomains[x.id] = modification
    end
    return
end


text:
using Pkg
println("Activating environment in $(pwd())...")
pkg"activate ."
println("Installing packages..."); flush(stdout);
pkg"instantiate"
pkg"precompile"
println("Done!")


text:
module COSMA

import LinearAlgebra: mul!
using LinearAlgebra: BlasFloat, BlasReal, BlasComplex, Adjoint, Transpose
using DistributedArrays
using Distributed
using MPI
using MPIClusterManagers: MPIManager
using COSMA_jll: cosma

const DMatrix{T} = DArray{T,2}

const mpimanager = Base.RefValue{Union{Nothing,MPIManager}}(nothing)

function use_manager(manager::MPIManager)
    mpimanager[] = manager
end

struct DArrayWithMPIRanks{T,TA,TB}
    array::TA
    ranks::TB
end

function DArrayWithMPIRanks(D::DArray{T}) where T
    ranks = map(pid -> mpimanager[].j2mpi[pid], D.pids)

    return DArrayWithMPIRanks{T,typeof(D),typeof(ranks)}(D, ranks)
end

# Here we could add an underscore for fortran bindings if necessary
macro cosmafunc(x)
    return Expr(:quote, x)
end

struct Block
    data::Ptr{Nothing}
    ld::Cint
    row::Cint
    col::Cint
end

struct JuliaLayout
    rowblocks::Cint
    colblocks::Cint
    rowsplit::Vector{Cint}
    colsplit::Vector{Cint}
    owners::Vector{Cint}
    nlocalblocks::Cint
    localblocks::Vector{Block}
end

function JuliaLayout(A::DArrayWithMPIRanks{T}, rank) where T
    rowblocks = Cint(size(A.ranks, 1))
    colblocks = Cint(size(A.ranks, 2))
    rowsplit = map(i -> Cint(i - 1), A.array.cuts[1])
    colsplit = map(i -> Cint(i - 1), A.array.cuts[2])
    owners = map(Cint, A.ranks)

    # Get the blocks we own
    indices = findall(i -> i == rank, owners)
    @assert length(indices) == 1 # should be enforced by DArray
    block_row, block_col = Cint.(first(indices).I .- 1)

    localblocks = Block(
        Base.unsafe_convert(Ptr{Cvoid}, pointer(localpart(A.array))),
        stride(localpart(A.array), 2),
        block_row,
        block_col
    )

    return JuliaLayout(rowblocks, colblocks, rowsplit, colsplit, vec(owners), 1, [localblocks])
end

struct Layout
    rowblocks::Cint
    colblocks::Cint
    rowsplit::Ptr{Cint}
    colsplit::Ptr{Cint}
    owners::Ptr{Cint}
    nlocalblocks::Cint
    blocks::Ptr{Block}
end

Layout(A::JuliaLayout) = Layout(
    A.rowblocks,
    A.colblocks,
    pointer(A.rowsplit),
    pointer(A.colsplit),
    pointer(A.owners),
    A.nlocalblocks,
    pointer(A.localblocks)
)

for (fname, elty) in ((:dmultiply_using_layout,:Float64),
                      (:smultiply_using_layout,:Float32),
                      (:zmultiply_using_layout,:ComplexF64),
                      (:cmultiply_using_layout,:ComplexF32))
    @eval begin
        function gemm!(C::DArrayWithMPIRanks{$elty}, transa::AbstractChar, transb::AbstractChar, A::DArrayWithMPIRanks{$elty}, B::DArrayWithMPIRanks{$elty}, alpha::$elty, beta::$elty)
            # Get our rank
            comm = MPI.COMM_WORLD
            rank = MPI.Comm_rank(comm)

            A_julia = JuliaLayout(A, rank)
            B_julia = JuliaLayout(B, rank)
            C_julia = JuliaLayout(C, rank)

            GC.@preserve A_julia B_julia C_julia A B C begin
                A_cosma = Layout(A_julia)
                B_cosma = Layout(B_julia)
                C_cosma = Layout(C_julia)

                ccall((@cosmafunc($fname), cosma), Cvoid, (
                    MPI.MPI_Comm #= comm =#,
                    Ref{Cchar} #= transa =#,
                    Ref{Cchar} #= transb =#,
                    Ref{$elty} #= alpha =#,
                    Ref{Layout} #= layout_a =#,
                    Ref{Layout} #= layout_b =#,
                    Ref{$elty} #= beta =#,
                    Ref{Layout} #= layout_c =#,
                ),
                    comm, transa, transb, alpha, A_cosma, B_cosma, beta, C_cosma
                )
            end
        end
    end
end

function gemm_wrapper!(C::DMatrix{T}, transa::AbstractChar, transb::AbstractChar, A::DMatrix{T}, B::DMatrix{T}, alpha::T, beta::T) where {T<:BlasFloat}
    Ad1, Ad2 = (transa == 'N') ? (1,2) : (2,1)
    Bd1, Bd2 = (transb == 'N') ? (1,2) : (2,1)
    mA, nA = (size(A, Ad1), size(A, Ad2))
    mB, nB = (size(B, Bd1), size(B, Bd2))
    if mB != nA
        throw(DimensionMismatch("matrix A has dimensions ($mA, $nA), matrix B has dimensions ($mB, $nB)"))
    end
    if size(C,1) != mA || size(C,2) != nB
        throw(DimensionMismatch("result C has dimensions $(size(C)), needs ($mA, $nB)"))
    end

    # Get unique workers
    pids = unique(vcat(vec(C.pids), vec(A.pids), vec(B.pids)))

    # Julia pids are not mpi ranks :(
    A_mpi = DArrayWithMPIRanks(A)
    B_mpi = DArrayWithMPIRanks(B)
    C_mpi = DArrayWithMPIRanks(C)

    @sync for p in pids
        @spawnat p gemm!(C_mpi, transa, transb, A_mpi, B_mpi, alpha, beta)
    end

    return C
end

function promote_alpha_beta(a, b, ::Type{T}) where {T}
    a_prom, b_prom = promote(a, b, zero(T))
    a_prom, b_prom
end

# Pirate the mul! implementation
mul!(C::DMatrix{T}, A::DMatrix{T}, B::DMatrix{T}, a::Number, b::Number) where T<:BlasFloat =
    gemm_wrapper!(C, 'N', 'N', A, B, promote_alpha_beta(a, b, T)...)
mul!(C::DMatrix{T}, trA::Transpose{T, <:DMatrix{T}}, B::DMatrix{T}, a::Number, b::Number) where T<:BlasFloat =
    gemm_wrapper!(C, 'T', 'N', parent(trA), B, promote_alpha_beta(a, b, T)...)
mul!(C::DMatrix{T}, A::DMatrix{T}, trB::Transpose{T, <:DMatrix{T}}, a::Number, b::Number) where T<:BlasFloat =
    gemm_wrapper!(C, 'N', 'T', A, parent(trB), promote_alpha_beta(a, b, T)...)
mul!(C::DMatrix{T}, trA::Transpose{T, <:DMatrix{T}}, trB::Transpose{T, <:DMatrix{T}}, a::Number, b::Number) where T<:BlasFloat =
    gemm_wrapper!(C, 'T', 'T', parent(trA), parent(trB), promote_alpha_beta(a, b, T)...)
mul!(C::DMatrix{T}, adjA::Adjoint{T, <:DMatrix{T}}, B::DMatrix{T}, a::Real, b::Real) where T<:BlasReal =
    gemm_wrapper!(C, 'T', 'N', parent(adjA), B, promote_alpha_beta(a, b, T)...)
mul!(C::DMatrix{T}, adjA::Adjoint{T, <:DMatrix{T}}, B::DMatrix{T}, a::Number, b::Number) where T<:BlasComplex =
    gemm_wrapper!(C, 'C', 'N', parent(adjA), B, promote_alpha_beta(a, b, T)...)
mul!(C::DMatrix{T}, A::DMatrix{T}, adjB::Adjoint{T, <:DMatrix{T}}, a::Real, b::Real) where T<:BlasReal =
    gemm_wrapper!(C, 'N', 'T', A, parent(adjB), promote_alpha_beta(a, b, T)...)
mul!(C::DMatrix{T}, A::DMatrix{T}, adjB::Adjoint{T, <:DMatrix{T}}, a::Number, b::Number) where T<:BlasComplex =
    gemm_wrapper!(C, 'N', 'C', A, parent(adjB), promote_alpha_beta(a, b, T)...)
mul!(C::DMatrix{T}, adjA::Adjoint{T, <:DMatrix{T}}, adjB::Adjoint{T, <:DMatrix{T}}, a::Real, b::Real) where T<:BlasReal =
    gemm_wrapper!(C, 'T', 'T', parent(adjA), parent(adjB), promote_alpha_beta(a, b, T)...)
mul!(C::DMatrix{T}, adjA::Adjoint{T, <:DMatrix{T}}, adjB::Adjoint{T, <:DMatrix{T}}, a::Number, b::Number) where T<:BlasComplex =
    gemm_wrapper!(C, 'C', 'C', parent(adjA), parent(adjB), promote_alpha_beta(a, b, T)...)
mul!(C::DMatrix{T}, trA::Transpose{T, <:DMatrix{T}}, adjB::Adjoint{T, <:DMatrix{T}}, a::Real, b::Real) where T<:BlasReal =
    gemm_wrapper!(C, 'T', 'T', parent(trA), parent(adjB), promote_alpha_beta(a, b, T)...)
mul!(C::DMatrix{T}, trA::Transpose{T, <:DMatrix{T}}, adjB::Adjoint{T, <:DMatrix{T}}, a::Number, b::Number) where T<:BlasComplex =
    gemm_wrapper!(C, 'T', 'C', parent(trA), parent(adjB), promote_alpha_beta(a, b, T)...)
mul!(C::DMatrix{T}, adjA::Adjoint{T, <:DMatrix{T}}, trB::Transpose{T, <:DMatrix{T}}, a::Real, b::Real) where T<:BlasReal =
    gemm_wrapper!(C, 'T', 'T', parent(adjA), parent(trB), promote_alpha_beta(a, b, T)...)
mul!(C::DMatrix{T}, adjA::Adjoint{T, <:DMatrix{T}}, trB::Transpose{T, <:DMatrix{T}}, a::Number, b::Number) where T <: BlasComplex =
    gemm_wrapper!(C, 'C', 'T', parent(adjA), parent(trB), promote_alpha_beta(a, b, T)...)
end


text:
"""
    Comrade
Radio Observation Sampling Exploration
"""
module Comrade

using AbstractFFTs
using AbstractMCMC
using Accessors: @set
using ArgCheck: @argcheck
using BasicInterpolators
using DocStringExtensions
using ChainRulesCore
using ComradeBase
using ForwardDiff
using FFTW: fft, fftfreq, fftshift, ifft, ifft!, ifftshift, plan_fft
using FLoops
using MappedArrays: mappedarray
using MeasureBase
using NFFT: nfft, plan_nfft
using PaddedViews
using PyCall: pyimport, PyNULL, PyObject
using SpecialFunctions
using Reexport
using Requires: @require
using StructArrays: StructArray
# Write your package code here.

@reexport using ComradeBase
using ComradeBase: visanalytic, imanalytic

export SequentialEx, ThreadedEx

const ehtim = PyNULL()

"""
    $(SIGNATURES)
Loads the [eht-imaging](https://github.com/achael/eht-imaging) library and stores it in the
`ehtim` variable.

# Notes
This will fail if ehtim isn't installed in the python installation that PyCall references.
"""
function load_ehtim()
    #try
    copy!(ehtim, pyimport("ehtim"))
    #catch
    #    @warn "No ehtim installation found in python path. Some data functionality will not work"
    #end
end

export rad2μas, μas2rad, ehtim, load_ehtim

"""
    rad2μas(x)
Converts a number from radians to μas
"""
@inline rad2μas(x) = 180/π*3600*1e6*x

"""
    μas2rad(x)
Converts a number from μas to rad
"""
@inline μas2rad(x) = x/(180/π*3600*1e6)

#include("interface.jl")
#include("images/images.jl")
import ComradeBase: flux, radialextent, intensitymap, intensitymap!
export create_cache
include("observations/observations.jl")
include("models/models.jl")
include("distributions/distributions.jl")
include("distributions/radiolikelihood.jl")
include("visualizations/visualizations.jl")
include("bayes/bayes.jl")
include("inference/inference.jl")

function __init__()
    @require AdvancedHMC="0bf59076-c3b1-5ca4-86bd-e02cd72cde3d" include("inference/advancedhmc.jl")
    @require NestedSamplers="41ceaf6f-1696-4a54-9b49-2e7a9ec3782e" include("inference/nested.jl")
    @require Pathfinder="b1d3bc72-d0e7-4279-b92f-7fa5d6d2d454" include("inference/pathfinder.jl")
    @require GalacticOptim="a75be94c-b780-496d-a8a9-0878b188d577" include("inference/galacticoptim.jl")
end

end


text:
#
#  weight_factors.m  ver 1.0  October 18, 2012
#
function weight_factors_mod(iweight)

fc=zeros(44,1);
wk=zeros(44,1);
wd=zeros(44,1);
wf=zeros(44,1);
wc=zeros(44,1);
we=zeros(44,1);
wj=zeros(44,1);
wb=zeros(44,1);
#
#
fc[1]=0.02;
wk[1]=-90.;
wd[1]=-90.;
wf[1]=-32.33;
wc[1]=-90.;
we[1]=-90.;
wj[1]=-90.;
wb[1]=-90.;
#
fc[2]=0.025;
wk[2]=-90.;
wd[2]=-90.;
wf[2]=-28.48;
wc[2]=-90.;
we[2]=-90.;
wj[2]=-90.;
wb[2]=-90.;
#
fc[3]=0.0315;
wk[3]=-90.;
wd[3]=-90.;
wf[3]=-24.47;
wc[3]=-90.;
we[3]=-90.;
wj[3]= -90.;
wb[3]=-90.;
#
fc[4]=0.04;
wk[4]=-90.;
wd[4]=-90.;
wf[4]=-20.25;
wc[4]=-90.;
we[4]=-90.;
wj[4]=-90.;
wb[4]=-90.;
#
fc[5]=0.05;
wk[5]=-90.;
wd[5]=-90.;
wf[5]=-16.10;
wc[5]=-90.;
we[5]=-90.;
wj[5]=-90.;
wb[5]=-90.;
#
fc[6]=0.063;
wk[6]=-90.;
wd[6]=-90.;
wf[6]=-11.49;
wc[6]=-90.;
we[6]= -90.;
wj[6]=-90.;
wb[6]=-90.;
#
fc[7]=0.08;
wk[7]=-90.;
wd[7]=-90.;
wf[7]=-6.73;
wc[7]= -90.;
we[7]=-90.;
wj[7]=-90.;
wb[7]=-90.;
#
fc[8]=0.1;
wk[8]=-30.11;
wd[8]=-24.09;
wf[8]=-3.16;
wc[8]=-24.11;
we[8]=-24.08;
wj[8]=-30.18;
wb[8]=-32.04;
#
fc[9]=0.125;
wk[9]=-26.26;
wd[9]=-20.24;
wf[9]=-0.96;
wc[9]=-20.25;
we[9]=-20.22;
wj[9]=-26.32;
wb[9]=-28.20;
#
fc[10]=0.16;
wk[10]=-22.05;
wd[10]=-16.01;
wf[10]=0.05;
wc[10]=-16.03;
we[10]=-15.98;
wj[10]=-22.11;
wb[10]=-23.98;
#
fc[11]=0.2;
wk[11]=-18.33;
wd[11]=-12.28;
wf[11]=-0.07;
wc[11]=-12.03;
we[11]=-12.23;
wj[11]=-18.38;
wb[11]=-20.23;
#
fc[12]=0.25;
wk[12]=-14.81;
wd[12]=-8.75;
wf[12]=-1.37;
wc[12]=-8.78;
we[12]=-8.67;
wj[12]=-14.86;
wb[12]=-16.71;
#
fc[13]=0.315;
wk[13]=-11.60;
wd[13]=-5.52;
wf[13]=-4.17;
wc[13]=-5.56;
we[13]=-5.41;
wj[13]=-11.65;
wb[13]=-13.51;
#
fc[14]=0.4;
wk[14]=-9.07;
wd[14]=-2.94;
wf[14]=-8.31;
wc[14]=-3.01;
we[14]=-3.01;
wj[14]=-9.10;
wb[14]=-10.98;
#
fc[15]=0.5;
wk[15]=-7.57;
wd[15]=-1.38;
wf[15]=-13.00;
wc[15]=-1.48;
we[15]=-1.29;
wj[15]=-7.60;
wb[15]=-9.53;
#
fc[16]=0.63;
wk[16]=-6.77;
wd[16]=-0.50;
wf[16]=-18.69;
wc[16]=-0.64;
we[16]=-0.55;
wj[16]=-6.78;
wb[16]=-8.71;
#
fc[17]=0.8;
wk[17]=-6.43;
wd[17]=-0.07;
wf[17]=-25.51;
wc[17]=-0.24;
we[17]=-0.53;
wj[17]=-6.42;
wb[17]=-8.38;
#
fc[18]=1.;
wk[18]=-6.33;
wd[18]=0.10;
wf[18]=-32.57;
wc[18]=-0.08;
we[18]=-1.11;
wj[18]=-6.30;
wb[18]=-8.29;
#
fc[19]=1.25;
wk[19]=-6.29;
wd[19]=0.07;
wf[19]=-40.02;
wc[19]=0.00;
we[19]=-2.25;
wj[19]=-6.28;
wb[19]=-8.27;
#
fc[20]=1.6;
wk[20]=-6.12;
wd[20]=-0.28;
wf[20]=-48.47;
wc[20]=0.06;
we[20]=-3.99;
wj[20]=-6.32;
wb[20]=-8.07;
#
fc[21]=2.;
wk[21]=-5.49;
wd[21]=-1.01;
wf[21]=-56.19;
wc[21]=0.10;
we[21]=-5.82;
wj[21]=-6.34;
wb[21]=-7.60;
#
fc[22]=2.5;
wk[22]=-4.01;
wd[22]=-2.20;
wf[22]=-63.93;
wc[22]=0.15;
we[22]=-7.77;
wj[22]=-6.22;
wb[22]=-6.13;
#
fc[23]=3.15;
wk[23]=-1.90;
wd[23]=-3.85;
wf[23]=-71.96;
wc[23]=0.19;
we[23]=-9.81;
wj[23]=-5.62;
wb[23]=-3.58;
#
fc[24]=4.;
wk[24]=-0.29;
wd[24]=-5.82;
wf[24]=-80.26;
wc[24]=0.20;
we[24]=-11.93;
wj[24]=-4.04;
wb[24]=-1.02;
#
fc[25]=5.;
wk[25]=0.33;
wd[25]=-7.76;
wf[25]=-90.;
wc[25]=0.11;
we[25]=-13.91;
wj[25]=-2.01;
wb[25]=0.21;
#
fc[26]=6.3;
wk[26]=0.46;
wd[26]=-9.81;
wf[26]=-90.;
wc[26]=-0.23;
we[26]=-15.94;
wj[26]=-0.48;
wb[26]=0.46;
#
fc[27]=8.;
wk[27]=0.31;
wd[27]=-11.93;
wf[27]=-90.;
wc[27]=-1.00;
we[27]=-18.03;
wj[27]=0.15;
wb[27]=0.21;
#
fc[28]=10.;
wk[28]=-0.10;
wd[28]=-13.91;
wf[28]=-90.;
wc[28]=-2.20;
we[28]=-19.98;
wj[28]=0.26;
wb[28]=-0.23;
#
fc[29]=12.5;
wk[29]=-0.89;
wd[29]=-15.87;
wf[29]=-90.;
wc[29]=-3.79;
we[29]=-21.93;
wj[29]=0.22;
wb[29]=-0.85;
#
fc[30]=16.;
wk[30]=-2.28;
wd[30]=-18.03;
wf[30]=-90.;
wc[30]=-5.82;
we[30]=-24.08;
wj[30]=0.16;
wb[30]=-1.83;
#
fc[31]=20.;
wk[31]=-3.93;
wd[31]=-19.99;
wf[31]=-90.;
wc[31]=-7.77;
we[31]=-26.02;
wj[31]=0.10;
wb[31]=-3.00;
#
fc[32]=25.;
wk[32]=-5.80;
wd[32]=-21.94;
wf[32]=-90.;
wc[32]=-9.76;
we[32]=-27.97;
wj[32]=0.06;
wb[32]=-4.44;
#
fc[33]=31.5;
wk[33]=-7.86;
wd[33]=-23.98;
wf[33]=-90.;
wc[33]=-11.84;
we[33]=-30.01;
wj[33]=0.00;
wb[33]=-6.16;
#
fc[34]=40.;
wk[34]=-10.05;
wd[34]=-26.13;
wf[34]=-90.;
wc[34]=-14.02;
we[34]=-32.15;
wj[34]=-0.08;
wb[34]=-8.11;
#
fc[35]=50.;
wk[35]=-12.19;
wd[35]=-28.22;
wf[35]=-90.;
wc[35]=-16.13;
we[35]=-34.24;
wj[35]=-0.24;
wb[35]=-10.09;
#
fc[36]=63.;
wk[36]=-14.61;
wd[36]=-30.60;
wf[36]=-90.;
wc[36]=-18.53;
we[36]=-36.62;
wj[36]=-0.62;
wb[36]=-12.43;
#
fc[37]=80.;
wk[37]=-17.56;
wd[37]=-33.53;
wf[37]=-90.;
wc[37]=-21.47;
we[37]=-39.55;
wj[37]=-1.48;
wb[37]=-15.34;
#
fc[38]=100.;
wk[38]=-21.04;
wd[38]=-36.99;
wf[38]=-90.;
wc[38]=-24.94;
we[38]=-43.01;
wj[38]=-3.01;
wb[38]=-18.72;
#
fc[39]=125.;
wk[39]=-25.35;
wd[39]=-41.28;
wf[39]=-90.;
wc[39]=-29.24;
we[39]=-47.31;
wj[39]=-5.36;
wb[39]=-23.00;
#
fc[40]=160.;
wk[40]=-30.91;
wd[40]=-46.84;
wf[40]=-90.;
wc[40]=-34.80;
we[40]=-52.86;
wj[40]=-8.78;
wb[40]=-28.56;
#
fc[41]=200.;
wk[41]=-36.38;
wd[41]=-52.30;
wf[41]=-90.;
wc[41]=-40.26;
we[41]=-58.33;
wj[41]=-12.30;
wb[41]=-34.03;
#
fc[42]=250.;
wk[42]=-42.0;
wd[42]=-57.97;
wf[42]=-90.;
wc[42]=-45.92;
we[42]=-63.99;
wj[42]=-16.03;
wb[42]=-39.69;
#
fc[43]=315.;
wk[43]=-48.00;
wd[43]=-63.92;
wf[43]=-90.;
wc[43]=-51.88;
we[43]=-69.94;
wj[43]=-19.98;
wb[43]=-46.65;
#
fc[44]=400.;
wk[44]=-54.20;
wd[44]=-70.12;
wf[44]=-90.;
wc[44]=-58.08;
we[44]=-76.14;
wj[44]=-24.10;
wb[44]=-51.84;
#
wk=10 .^ (wk/20.);
wd=10 .^ (wd/20.);
wf=10 .^ (wf/20.);
wc=10 .^ (wc/20.);
we=10 .^ (we/20.);
wj=10 .^ (wj/20.);
wb=10 .^ (wb/20.);
#
www=zeros(44,1);
#
for i=1:44
#
    if(iweight==1)
        www[i]=wk[i];
    end
    if(iweight==2)
        www[i]=wd[i];
    end
    if(iweight==3)
        www[i]=wf[i];
    end
    if(iweight==4)
        www[i]=wc[i];
    end
    if(iweight==5)
        www[i]=we[i];
    end
    if(iweight==6)
        www[i]=wj[i];
    end
    if(iweight==7)
        www[i]=wb[i];
    end
#
end
#
imax=length(fc);
#
fw=fc;
#
fwl=fw/(2^(1/6));
#
fwu=zeros(imax,1);
#
for i=1:imax-1
    fwu[i]=fwl[i+1];
end
#
fwu[imax]=fw[imax]*(2^(1/6));
#
return (fwl,fw,fwu,wk,wd,wf,wc,we,wj,wb,www)
end


text:
module JuMPUtils

import JuMP
export get_values, set_solver_attributes!, init_if_hasproperty!, isconverged

get_values(; jump_vars...) = (; map(((k, v),) -> k => JuMP.value.(v), collect(jump_vars))...)

function set_solver_attributes!(opt_model; solver_attributes...)
    foreach(
        ((k, v),) -> JuMP.set_optimizer_attribute(opt_model, string(k), v),
        pairs(solver_attributes),
    )
end

function init_if_hasproperty!(v, init, sym; default = nothing)
    init_value = hasproperty(init, sym) ? getproperty(init, sym) : default
    if !isnothing(init_value)
        JuMP.set_start_value.(v, init_value)
    end
end

function isconverged(opt_model)
    JuMP.termination_status(opt_model) in (JuMP.MOI.LOCALLY_SOLVED, JuMP.MOI.OPTIMAL)
end

end


text:
#import BlockRegistration
using RegisterDeformation
using CoordinateTransformations, Interpolations, ImageCore, ForwardDiff
using StaticArrays, LinearAlgebra, Distributed, Statistics
using AxisArrays: AxisArray
using OffsetArrays
using Test
using RegisterDeformation.RegisterUtilities
using JLD2, HDF5, FileIO
using ImageFiltering

@testset "Nodes" begin
    nodes = (range(1, stop=15, length=5), range(1, stop=11, length=3))
    @test arraysize(nodes) == (15,11)

    nodesall = (nodes[1].*ones(length(nodes[2]))', ones(length(nodes[1])).*nodes[2]')
    for (i, knt) in enumerate(eachnode(nodes))
        @test knt[1] == nodesall[1][i]
        @test knt[2] == nodesall[2][i]
    end
end

@testset "Interpolation and extrapolation" begin
    nodes = (range(1, stop=15, length=5), range(1, stop=11, length=3))
    u = randn(length(nodes), map(length, nodes)...)
    ϕ = GridDeformation(u, nodes)
    ϕi = interpolate(ϕ)
    ϕe = extrapolate(ϕ)
    @test extrapolate(ϕi) == ϕe

    ϕi = interpolate!(copy(ϕ))
    ϕe = extrapolate!(copy(ϕ))
    @test extrapolate(ϕi) == ϕe
end

@testset "Deformation tests" begin
    randrange(xmin, xmax) = (xmax-xmin)*rand() + xmin

    ## Deformations
    for nodes in ((range(1, stop=15, length=5),),
                  (range(1, stop=15, length=5), range(1, stop=8, length=4)),
                  (range(1, stop=15, length=3), range(1, stop=8, length=4), range(1, stop=7, length=2)),
                  (range(1, stop=15, length=5), range(1, stop=8, length=4), range(1, stop=7, length=6)))
        sz = map(length, nodes)
        rng = map(z->convert(Int, maximum(z)), nodes)
        # Zero deformation
        u = zeros(length(nodes), sz...)
        ϕ = interpolate(GridDeformation(u, nodes))
        for i = 1:10
            x = map(z->rand(1:z), rng)
            @test @inferred(ϕ(x...)) ≈ [x...]
            if all(x->x>2, sz)
                y = map(z->rand(2:z-1)+rand()-0.5, rng)
                @test ϕ(y...) ≈ [y...]
            end
        end
        # Constant shift
        dx = randn(length(nodes))
        u = zeros(length(nodes), sz...).+dx
        ϕ = interpolate(GridDeformation(u, nodes))
        for i = 1:10
            x = map(z->rand(1:z), rng)
            @test ϕ(x...) ≈ [x...]+dx
            if all(x->x>2, sz)
                y = map(z->rand(2:z-1)+rand()-0.5, rng)
                @test ϕ(y...) ≈ [y...]+dx
            end
        end
        # "Biquadratic"
        if all(x->x>2, sz)
            N = length(nodes)
            u_is = Vector{Any}(undef, N)
            fs = Vector{Any}(undef, N)
            R = CartesianIndices(sz)
            for i = 1:N
                let cntr = Float64[rand(1:z-1)+rand() for z in rng], q = rand(N)
                    f = x -> dot(q, (x-cntr).^2)
                    ui = zeros(sz)
                    for I in R
                        ui[I] = f(Float64[nodes[d][I[d]] for d = 1:N])
                    end
                    u_is[i] = reshape(ui, 1, sz...)
                    fs[i] = f
                end
            end
            u = cat(u_is..., dims = 1)
            ϕ = interpolate!(GridDeformation(copy(u), nodes), InPlaceQ(OnCell()))
            for i = 1:10
                y = Float64[randrange((node[2]+node[1])/2, (node[end]+node[end-1])/2) for node in nodes]
                dx = Float64[fs[d](y) for d=1:N]
                try
                    @test ϕ(y...) ≈ y + Float64[fs[d](y) for d=1:N]
                catch err
                    @show y nodes
                    @show [((node[2]+node[1])/2, (node[end]+node[end-1])/2) for node in nodes]
                    @show Float64[fs[d](y) for d=1:N]
                    rethrow(err)
                end
            end
        end
    end
end

@testset "tform2deformation" begin
    s = [3.3,-2.6]
    gsize = (3,2)
    A = tformtranslate(s)
    ϕ = tform2deformation(A, map(Base.OneTo, (500,480)), gsize)
    u = reshape(reinterpret(Float64, ϕ.u), tuple(2,gsize...))
    @test all(u[1,:,:] .== s[1])
    @test all(u[2,:,:] .== s[2])
end

@testset "similarϕ" begin
    uref = reshape([1.5, -0.3, 0.8], 1, 3)
    nodes = (0:5:10,)
    ϕref = GridDeformation(uref, nodes)
    unew = [0.7, 0.4, -0.33]
    ϕnew = similarϕ(ϕref, unew)
    ϕi = interpolate!(copy(ϕnew))
    @test ϕi.u(0)  ≈ SVector(0.7)
    @test ϕi.u(5)  ≈ SVector(0.4)
    @test ϕi.u(10) ≈ SVector(-0.33)

    ucoefs = copy(ϕi.u.itp.coefs)
    ϕiref = interpolate!(copy(ϕref))
    ϕi = similarϕ(ϕiref, ucoefs)
    @test ϕi.u(0)  ≈ SVector(0.7)
    @test ϕi.u(5)  ≈ SVector(0.4)
    @test ϕi.u(10) ≈ SVector(-0.33)
end

@testset "warp" begin
    # Ensure there is no conflict between ImageTransformations and RegisterDeformation
    tform = tformrotate(pi/4)
    ramp = 0.0:0.01:1
    nhalf = length(ramp) ÷ 2
    for img in (ramp*ramp', OffsetArray(ramp*ramp', -nhalf:nhalf, -nhalf:nhalf))
        imgw1 = warp(img, tform, axes(img))
        ϕ = tform2deformation(tform, axes(img), (7, 7))
        imgw2 = warp(img, ϕ)
        @test all(x->isnan(x) || abs(x) < 0.02, imgw1 - imgw2)
    end
end

@testset "WarpedArray" begin
    znan(x) = isnan(x) ? zero(x) : x

    p = (1:100).-5

    dest = Vector{Float32}(undef, 10)

    # Simple translations in 1d
    ϕ = GridDeformation([0.0,0.0]', axes(p))
    q = WarpedArray(p, ϕ);
    getindex!(dest, q, 1:10)
    @test maximum(abs.(dest - p[1:10])) < 10*eps(maximum(abs.(dest)))

    ϕ = GridDeformation([1.0,1.0]', axes(p))
    q = WarpedArray(p, ϕ);
    getindex!(dest, q, 1:10)
    @test maximum(abs.(dest - p[2:11])) < 10*eps(maximum(abs.(dest)))

    ϕ = GridDeformation([-2.0,-2.0]', axes(p))
    q = WarpedArray(p, ϕ);
    getindex!(dest, q, 1:10)
    @test maximum(abs.(dest[3:end] - p[1:8])) < 10*eps(maximum(abs.(dest[3:end])))

    getindex!(dest, q, 3:12)
    @test maximum(abs.(dest - p[1:10])) < 10*eps(maximum(abs.(dest)))

    ϕ = GridDeformation([2.0,2.0]', axes(p))
    q = WarpedArray(p, ϕ);
    getindex!(dest, q, 1:10)
    @test maximum(abs.(dest - p[3:12])) < 10*eps(maximum(abs.(dest)))

    ϕ = GridDeformation([5.0,5.0]', axes(p))
    q = WarpedArray(p, ϕ);
    getindex!(dest, q, 1:10)
    @test maximum(abs.(dest - p[6:15])) < 10*eps(maximum(abs.(dest)))

    # SubArray (test whether we can go beyond edges)
    psub = view(collect(p), 3:20)
    ϕ = GridDeformation([0.0,0.0]', axes(p))
    q = WarpedArray(psub, ϕ);
    getindex!(dest, q, -1:8)
    @test maximum(abs.(dest[3:end] - p[3:10])) < 10*eps(maximum(abs.(dest[3:end])))
    any(isnan, dest) && @warn("Some dest are NaN, not yet sure whether this is a problem")

    # Stretches
    u = [0.0,5.0,10.0]
    ϕ = interpolate(GridDeformation(u', axes(p)), Line(OnCell()))
    q = WarpedArray(p, ϕ)
    getindex!(dest, q, 1:10)
    @test abs(dest[1] - p[1]) < sqrt(eps(1.0f0))
    getindex!(dest, q, 86:95)
    @test isnan.(dest) == [falses(5);trues(5)]  # fixme
    dest2 = getindex!(zeros(Float32, 100), q, 1:100)
    @test all(abs.(diff(dest2)[26:74] .- ((u[3]-u[1])/99+1)) .< sqrt(eps(1.0f0)))

    #2d
    p = reshape(1:120, 10, 12)
    u1 = [2.0 2.0; 2.0 2.0]
    u2 = [-1.0 -1.0; -1.0 -1.0]
    ϕ = GridDeformation((u1,u2), axes(p))
    q = WarpedArray(p, ϕ)
    dest = zeros(axes(p))
    rng = (1:size(p,1),1:size(p,2))
    getindex!(dest, q, rng...)
    @test maximum(abs.(znan.(dest[1:7,2:end] - p[3:9,1:end-1]))) < 10*eps(maximum(znan.(dest)))
end

@testset "Composition" begin
    # Test that two rotations approximately compose to another rotation
    gridsize = (49,47)
    imgaxs = centeraxes(map(Base.OneTo, (200,200)))
    ϕ1 = tform2deformation(tformrotate( pi/180), imgaxs, gridsize)
    ϕ2 = tform2deformation(tformrotate(2pi/180), imgaxs, gridsize)
    ϕi = interpolate(ϕ1)
    ϕc = ϕi(ϕi)
    uc = reshape(reinterpret(Float64, ϕc.u), (2, gridsize...))
    u2 = reshape(reinterpret(Float64, ϕ2.u), (2, gridsize...))
    @test ≈(uc, u2, rtol=0.02, norm=x->maximum(abs, x))

    ## Test gradient computation

    function compare_g(g, gj, gridsize)
        for j = 1:gridsize[2], i = 1:gridsize[1]
            indx = (LinearIndices(gridsize))[i,j]
            rng = 2*(indx-1)+1:2indx
            @test g[i,j] ≈ gj[rng, rng]
            for k = 1:2*prod(gridsize)
                if !in(k, rng)
                    @test abs(gj[k,rng[1]]) < 1e-14
                    @test abs(gj[k,rng[2]]) < 1e-14
                end
            end
        end
        nothing
    end

    # A composition function that ForwardDiff will be happy with
    function compose_u(ϕ1, u2, f, shp, imaxs)
        ϕ2 = f(GridDeformation(reshape(u2, shp), imaxs))
        ϕ = ϕ1(ϕ2)
        ret = Vector{eltype(eltype(ϕ.u))}(undef, prod(shp))
        i = 0
        for vu in ϕ.u
            for d = 1:length(vu)
                ret[i+=1] = vu[d]
            end
        end
        ret
    end

    gridsize = (11,11)
    u1 = randn(2, gridsize...)
    u2 = randn(2, gridsize...)
    imaxs = map(Base.OneTo, (100,101))
    ϕ1 = interpolate(GridDeformation(u1, imaxs))
    for f in (interpolate, identity)
        ϕ2 = f(GridDeformation(u2, imaxs))
        ϕ, g = compose(ϕ1, ϕ2)
        u2vec = vec(u2)
        gj = ForwardDiff.jacobian(u2vec -> compose_u(ϕ1, u2vec, f, size(u2), imaxs), u2vec)
        compare_g(g, gj, gridsize)
    end

    # Test identity-composition
    _, g = compose(identity, ϕ2)
    gj = zeros(2*prod(gridsize), 2*prod(gridsize))
    for i = 1:prod(gridsize)
        rng = 2*(i-1)+1:2i
        gj[rng,rng] = Matrix{Float64}(I,2,2)
    end
    compare_g(g, gj, gridsize)

    # Extrapolation (issue #7)
    gridsize = (3,3)
    img = rand(10,10)
    tform = AffineMap(1.3, [0.1, 0.1])
    ϕ = tform2deformation(tform, axes(img), gridsize)
    θ = tform2deformation(tform, axes(img), gridsize)

    ϕ = extrapolate!(ϕ)
    ϕθ = ϕ(θ)
    ϕθi = interpolate(ϕθ)
    θi = interpolate(θ)
    x = SVector(7, 3)
    @test norm(ϕθi(x) - tform(tform(x))) <= 2*norm(ϕ(θi(x)) - tform(tform(x)))
    x = SVector(3, 7)
    @test norm(ϕθi(x) - tform(tform(x))) <= 2*norm(ϕ(θi(x)) - tform(tform(x)))
end

@testset "warpgrid" begin
    nodes = (range(1, stop=100, length=4), range(1, stop=100, length=3))
    gridsize = map(length, nodes)
    ϕ = GridDeformation(5*randn(2,gridsize...), nodes)
    kg = nodegrid(nodes)
    @test eltype(kg) == Bool
    @test size(kg) == (100, 100)
    kg2 = nodegrid(ϕ)
    @test kg2 == kg
    kg = nodegrid(Float64, nodes)
    @test eltype(kg) == Float64
    A = warpgrid(ϕ)
    B = warpgrid(ϕ, scale=1.5)
    @test A != B
    @test isa(A, Matrix{Float32})
    C = warpgrid(ϕ, showidentity=true)
    @test eltype(C) == RGB{Float32}
end

@testset "warp I/O" begin
    # Writing warped data to disk
    o = ones(Float32, 5, 5)
    A = o .* reshape(1:7, (1,1,7))
    img = AxisArray(A, :y, :x, :time)
    fn = tempname()
    # With Vector{GridDeformation}
    ϕs = tighten([GridDeformation(zeros(2,3,3), axes(o)) for i = 1:nimages(img)])
    open(fn, "w") do io
        warp!(Float32, io, img, ϕs)
    end
    warped = open(fn, "r") do io
        read!(io, Array{Float32}(undef, size(img)))
    end
    @test warped == img
    # With Array{SVector}
    uarray = reshape(reinterpret(SVector{2,Float64}, zeros(2,3,3,7)), (3,3,7))
    open(fn, "w") do io
        warp!(Float32, io, img, uarray)
    end
    warped = open(fn, "r") do io
        read!(io, Array{Float32}(undef, size(img)))
    end
    @test warped == img
    # With Array{Real}
    uarray = zeros(2,3,3,7)
    fn = tempname()
    open(fn, "w") do io
        warp!(Float32, io, img, uarray)
    end
    warped = open(fn, "r") do io
        read!(io, Array{Float32}(undef, size(img)))
    end
    @test warped == img
    # Multi-process
    open(fn, "w") do io
        warp!(Float32, io, img, uarray; nworkers=3)
    end
    warped = open(fn, "r") do io
        read!(io, Array{Float32}(undef, size(img)))
    end
    @test warped == img
end

@testset "Saving deformations" begin
    # Saving arrays-of-staticarrays efficiently
    nodes = (range(1, stop=20, length=3), range(1, stop=30, length=5))
    ϕ = GridDeformation(rand(2,map(length,nodes)...), nodes)
    fn = string(tempname(), ".jld2")
    save(fn, "ϕ", ϕ)
    ϕ2 = load(fn, "ϕ")
    @test ϕ2.u == ϕ.u
    @test ϕ2.nodes == ϕ.nodes
    if HDF5.h5_get_libversion() >= v"1.10" && get(ENV, "CI", nothing) != "true"
        str = read(`h5dump $fn`, String)
        @test something(findfirst("SIMPLE { ( 5, 3, 2 ) / ( 5, 3, 2 ) }", str), 0:-1) != 0:-1||
            something(findfirst("SIMPLE { ( 30 ) / ( 30 ) }", str), 0:-1) != 0:-1||
            something(findfirst("SIMPLE { ( 2 ) / ( 2 ) }", str), 0:-1) != 0:-1
    end
    rm(fn)
end

@testset "Temporal interpolation" begin
    u2 = [1.0 1.0]
    u4 = [3.0 3.0]
    nodes = (range(1, stop=20, length=2),)
    ϕ2 = GridDeformation(u2, nodes)
    ϕ4 = GridDeformation(u4, nodes)
    ϕsindex = [ϕ2, ϕ4]
    ϕs = tinterpolate(ϕsindex, [2,4], 5)
    @test length(ϕs) == 5
    @test eltype(ϕs) == eltype(ϕsindex)
    @test all(map(x-> x == @SVector([1.0]), ϕs[1].u))
    @test all(map(x-> x == @SVector([2.0]), ϕs[3].u))
end

@testset "Median filtering" begin
    u = rand(2, 3, 3, 9)
    ϕs = griddeformations(u, (range(1, stop=10, length=3), range(1, stop=11, length=3)))
    ϕsfilt = medfilt(ϕs, 3)
    v = ϕsfilt[3].u[2,2]
    v1 = median(vec(u[1,2,2,2:4]))
    v2 = median(vec(u[2,2,2,2:4]))
    @test v[1] == v1 && v[2] == v2
end

@testset "Regridding" begin
    z = [(x-2.5)^2 for x = 0:5]
    uold = zeros(2,6,6)
    uold[1,:,:] = z .* ones(6)'
    uold[2,:,:] = ones(6) .* z'
    nodes = (range(1, stop=20, length=6), range(1, stop=30, length=6))
    ϕ = GridDeformation(uold, nodes)
    ϕnew = regrid(ϕ, (10,8))
    @test ϕnew.nodes[1] ≈ range(1, stop=20, length=10)
    @test ϕnew.nodes[2] ≈ range(1, stop=30, length=8)
    ϕi = interpolate(ϕ)
    for (iy,y) in enumerate(range(1, stop=30, length=8))
        ys = (y-1)*5/29
        for (ix,x) in enumerate(range(1, stop=20, length=10))
            xs = (x-1)*5/19
            @test ≈(ϕnew.u[ix,iy], [(xs-2.5)^2,(ys-2.5)^2], rtol=0.2)
            @test ϕnew.u[ix,iy] ≈ ϕi.u(x,y)
        end
    end
end

@testset "Equality" begin
    # https://github.com/HolyLab/BlockRegistrationScheduler/pull/42#discussion_r139291426
    nodes = (range(0, stop=1, length=3), range(0, stop=1, length=5))
    u = rand(2, 3, 5)
    ϕ1 = GridDeformation(u, nodes)
    ϕ2 = GridDeformation(u, nodes)
    @test ϕ1 == ϕ2
end

@testset "Transformation of arrays" begin
    using RegisterDeformation: center
    # TODO: get rid of dependency on ImageFiltering for faster CI times
    padwith0(A) = parent(padarray(A, Fill(0, (2,2))))

    for IT in ((BSpline(Constant())),
                    (BSpline(Linear())),
                    (BSpline(Quadratic(Flat(OnCell())))))
        A = padwith0([1 1; 1 2])
        itp = interpolate(A, IT)
        tfm = tformtranslate([1,0])
        tA = TransformedArray(extrapolate(itp, NaN), tfm)
        @test @inferred(transform(tA))[3:4,3:4] ≈ [1 2; 0 0]
        @test tA[3,3] ≈ 1
        @test tA[3,4] ≈ 2
        tfm = tformtranslate([0,1])
        tA = TransformedArray(extrapolate(itp, NaN), tfm)
        @test transform(tA)[3:4,3:4] ≈ [1 0; 2 0]
        tfm = tformtranslate([-1,0])
        tA = TransformedArray(extrapolate(itp, NaN), tfm)
        @test transform(tA)[3:4,3:4] ≈ [0 0; 1 1]
        tfm = tformtranslate([0,-1])
        tA = TransformedArray(extrapolate(itp, NaN), tfm)
        @test transform(tA)[3:4,3:4] ≈ [0 1; 0 1]

        A = padwith0([2 1; 1 1])
        itp = interpolate(A, IT)
        tfm = tformrotate(pi/2)
        tA = TransformedArray(extrapolate(itp, NaN), tfm)
        @test transform(tA)[3:4,3:4] ≈ [1 2; 1 1]
        tfm = tformrotate(-pi/2)
        tA = TransformedArray(extrapolate(itp, NaN), tfm)
        @test transform(tA)[3:4,3:4] ≈ [1 1; 2 1]

        # Check that getindex and transform yield the same results
        A = rand(6,6)
        itp = interpolate(A, IT)
        tfm = tformrotate(pi/7)∘tformtranslate(rand(2))
        tA = TransformedArray(extrapolate(itp, NaN), tfm)
        dest = transform(tA)
        tfm_recentered = AffineMap(tfm.linear, tfm.translation + center(A) - tfm.linear*center(dest))
        tA_recentered = TransformedArray(extrapolate(itp, NaN), tfm_recentered)
        nbad = 0
        for j = 1:size(dest,2), i = 1:size(dest,1)
            val = tA_recentered[i,j]
            nbad += isfinite(dest[i,j]) != isfinite(val)
            if isfinite(val) && isfinite(dest[i,j])
                @test abs(val-dest[i,j]) < 1e-12
            end
        end
        @test nbad < 3
        dest = transform(tA, origin_src=zeros(2), origin_dest=zeros(2))
        nbad = 0
        for j = 1:size(dest,2), i = 1:size(dest,1)
            val = tA[i,j]
            nbad += isfinite(dest[i,j]) != isfinite(val)
            if isfinite(val) && isfinite(dest[i,j])
                @test abs(val-dest[i,j]) < 1e-12
            end
        end
        @test nbad < 3
    end

    A = Float64[1 2 3 4; 5 6 7 8]
    tfm = tformeye(2)
    dest = zeros(2,2)

    itp = interpolate(A, BSpline(Linear()))
    tA = TransformedArray(extrapolate(itp, NaN), tfm)
    @test transform!(dest, tA) ≈ [2 3; 6 7]
    dest3 = zeros(3,3)
    @test all(isapprox.(transform!(dest3, tA), [NaN NaN NaN; 3.5 4.5 5.5; NaN NaN NaN], nans=true))

    itp = interpolate(A, BSpline(Constant()))
    tA = TransformedArray(extrapolate(itp, NaN), tfm)
    @test transform!(dest, tA) == [2 3; 6 7]

    # Transforms with non-real numbers
    using DualNumbers
    a = tformtranslate([0,dual(1,0)])
    A = reshape(1:9, 3, 3)
    At = transform(A, a)
    @test At[:,1:2] == A[:,2:3]
end


text:
struct TransmissionNetwork <: AbstractTN
  external::BitArray{1}
  internal::BitArray{2}

  function TransmissionNetwork(individuals::Integer)
    return new(fill(0, individuals),
               fill(0, (individuals, individuals)))
  end

  function TransmissionNetwork(external::BitArray{1},
                               internal::BitArray{2})
    if !(length(external) == size(internal, 1) == size(internal, 2))
      error("Mismatched BitArray sizes")
    end
    multiple_exposures = findall((sum(internal, dims=1)[:] .+ external) .> 1)
    if length(multiple_exposures) > 0
      error("Multiple exposures detected for individual(s): $multiple_exposures")
    end
    return new(external, internal)
  end

  function TransmissionNetwork(starting_states::Vector{DiseaseState})
    return new(starting_states .!= Ref(State_S), fill(0, (length(starting_states), length(starting_states))))
  end
end


text:
export VaxFloatD

primitive type VaxFloatD <: VaxFloat 64 end

VaxFloatD(x::UInt64) = reinterpret(VaxFloatD, ltoh(x))
function VaxFloatD(x::T) where T <: Real
    y = reinterpret(UInt64, convert(Float64, x))
    part1 = y & bmask32
    part2 = (y >>> 32) & bmask32

    if ENDIAN_BOM === 0x04030201
        vaxpart2 = part1
        ieeepart1 = part2
    else
        vaxpart2 = part2
        ieeepart1 = part1
    end

    e = reinterpret(Int64, ieeepart1 & IEEE_T_EXPONENT_MASK)

    if ieeepart1 & ~SIGN_BIT_64 === zero(UInt64)
        # ±0.0 becomes 0.0
        return zero(VaxFloatD)
    elseif e === IEEE_T_EXPONENT_MASK
        # Vax types don't support ±Inf or NaN
        throw(InexactError(:VaxFloatD, VaxFloatD, x))
    else
        e >>>= IEEE_T_MANTISSA_SIZE
        m = ieeepart1 & IEEE_T_MANTISSA_MASK

        if e === zero(Int64)
            m = (m << 1) | (vaxpart2  >>> 31)
            vaxpart2 <<= 1
            while m & IEEE_T_HIDDEN_BIT === zero(UInt64)
                m = (m << 1) | (vaxpart2  >>> 31)
                vaxpart2 <<= 1
                e -= one(Int64)
            end
            m &= IEEE_T_MANTISSA_MASK
        end

        e += one(Int64) + VAX_D_EXPONENT_BIAS - IEEE_T_EXPONENT_BIAS
        if e <= zero(Int64)
            # Silent underflow
            return zero(VaxFloatD)
        elseif e > (2*VAX_D_EXPONENT_BIAS - 1)
            # Overflow
            throw(InexactError(:VaxFloatD, VaxFloatD, x))
        else
            vaxpart = (ieeepart1 & SIGN_BIT_64) |
                (e << VAX_D_MANTISSA_SIZE) |
                (m <<  (VAX_D_MANTISSA_SIZE - IEEE_T_MANTISSA_SIZE)) |
                (vaxpart2 >>> (32 - (VAX_D_MANTISSA_SIZE - IEEE_T_MANTISSA_SIZE)))
            vaxpart2 <<= (VAX_D_MANTISSA_SIZE - IEEE_T_MANTISSA_SIZE)
        end
    end

    vaxpart_1 = vaxpart & bmask16
    vaxpart_2 = (vaxpart >>> 16) & bmask16

    vaxpart_3 = vaxpart2 & bmask16
    vaxpart_4 = (vaxpart2 >>> 16) & bmask16

    res = htol((vaxpart_3 << 48) |
          (vaxpart_4 << 32) |
          (vaxpart_1 << 16) |
          vaxpart_2)

    return reinterpret(VaxFloatD, res)
end

Base.typemax(::Type{VaxFloatD}) = VaxFloatD(0xffffffffffff7fff)
Base.floatmax(::Type{VaxFloatD}) = VaxFloatD(0xffffffffffff7fff)
Base.typemin(::Type{VaxFloatD}) = VaxFloatD(0x0000000000000080)
Base.floatmin(::Type{VaxFloatD}) = VaxFloatD(0x0000000000000080)

Base.zero(::Type{VaxFloatD}) = VaxFloatD(0x0000000000000000)
Base.one(::Type{VaxFloatD}) = VaxFloatD(0x0000000000004080)

function Base.convert(::Type{Float64}, x::VaxFloatD)
    y = reinterpret(UInt64, ltoh(x))

    vaxpart_1 = y & bmask16
    vaxpart_2 = (y >>> 16) & bmask16
    vaxpart1 = (vaxpart_1 << 16) | vaxpart_2

    vaxpart_3 = (y >>> 32) & bmask16
    vaxpart_4 = (y >>> 48) & bmask16
    vaxpart2 = (vaxpart_3 << 16) | vaxpart_4

    if vaxpart1 & VAX_D_EXPONENT_MASK === zero(UInt64)
        if vaxpart1 & SIGN_BIT_64 === SIGN_BIT_64
            # Reserved floating-point reserved operand
            throw(InexactError(:convert, Float64, x))
        end

        # Dirty zero
        return zero(Float64)
    else
        ieeepart1 = ((vaxpart1 & SIGN_BIT_64) |
                     ((vaxpart1 & ~SIGN_BIT_64) >>>
                      (VAX_D_MANTISSA_SIZE - IEEE_T_MANTISSA_SIZE))) -
                     ((UNO64 + VAX_D_EXPONENT_BIAS - IEEE_T_EXPONENT_BIAS) << IEEE_T_MANTISSA_SIZE)
        ieeepart2 = (vaxpart1 << (32 - (VAX_D_MANTISSA_SIZE - IEEE_T_MANTISSA_SIZE))) |
                    (vaxpart2 >>> (VAX_D_MANTISSA_SIZE - IEEE_T_MANTISSA_SIZE))

        if ENDIAN_BOM === 0x04030201
            out1 = ieeepart2
            out2 = ieeepart1
        else
            out1 = ieeepart1
            out2 = ieeepart2
        end
    end

    res = (out2 << 32) | (out1 & bmask32)

    return reinterpret(Float64,res)
end
Base.convert(::Type{T},x::VaxFloatD) where T <: Union{Float16,Float32, BigFloat, Integer} = convert(T,convert(Float64,x))

Base.promote_rule(::Type{T},x::Type{VaxFloatD}) where T <: Union{AbstractVax, Float16, Float32, Float64, Integer} = Float64
Base.promote_rule(::Type{BigFloat},x::Type{VaxFloatD}) = BigFloat

Base.promote_type(::Type{VaxFloatD}, ::Type{VaxFloatD}) = Float64


text:
struct BanditExperience{T,TA} <: Any where {T,TA}
    action::TA
    blogp::T
    reward::T
end

struct ContextualBanditExperience{T,TX,TA} <: Any where {T,TX,TA}
    obs::TX
    action::TA
    blogp::T
    reward::T
end

mutable struct Trajectory{T,TS,TA} <: Any where {T,TS,TA}
    states::Array{TS,1}
    actions::Array{TA,1}
    blogps::Array{T,1}
    rewards::Array{T,1}
    done::Bool

    function Trajectory(::Type{T}, ::Type{TS}, ::Type{TA}) where {T,TS,TA}
        new{T,TS,TA}(Array{TS,1}(), Array{TA,1}(), Array{T,1}(), Array{T,1}(), false)
    end

    function Trajectory(prob::SequentialProblem)
        TS = get_type(prob.X)
        TA = get_type(prob.A)
        T = eltype(TS)
        Trajectory(T,TS,TA)
    end

    function Trajectory(prob::MDP)
        TS = get_type(prob.S)
        TA = get_type(prob.A)
        T = eltype(TS)
        Trajectory(T,TS,TA)
    end

    function Trajectory(prob::POMDP)
        TS = get_type(prob.X)
        TA = get_type(prob.A)
        T = eltype(TS)
        Trajectory(T,TS,TA)
    end

end

function get_type(t::Tuple)
    return get_type.(t)
end

function get_type(::UnitRange{T}) where T
    return T
end

function get_type(::Array{T,1}) where T
    return T
end

function get_type(::Array{T, 2}) where T
    return Array{T,1}
end

function length(τ::Trajectory)
    return length(τ.rewards)
end

function push!(τ::Trajectory{T,TS,TA}, state::TS, action::TA, blogp::T, reward::T) where {T,TS,TA}
    push!(τ.states, deepcopy(state))
    push!(τ.actions, deepcopy(action))
    push!(τ.blogps, blogp)
    push!(τ.rewards, reward)
end

function finish!(τ::Trajectory)
    τ.done = true
end

function finished(τ::Trajectory)
    return τ.done
end


text:
include("run_cutest.jl")
using Ipopt
TIME_LIMIT = 60.0^2
problem_list = default_list()

if false
  folder_name = "ipopt/very_low_tol"
  if_mkdir("../results/$folder_name")
  run_cutest_problems_on_solver(problem_list, folder_name, print_level=3, tol=1e-2, max_iter=3000, max_cpu_time=TIME_LIMIT, nlp_scaling_method="none", bound_relax_factor = 0.0, acceptable_iter=999999)
end

if true
  folder_name = "ipopt/Jan2"
  if_mkdir("../results/$folder_name")
  run_cutest_problems_on_solver(problem_list, folder_name, 4, 1e-6, 3000, TIME_LIMIT, "none", 0.0, 999999)
end

if false
  folder_name = "ipopt/medium_perturb"
  if_mkdir("../results/$folder_name")
  run_cutest_problems_on_solver(problem_list, folder_name, print_level=3, tol=1e-6, max_iter=3000, max_cpu_time=TIME_LIMIT, bound_relax_factor = 1e-8, nlp_scaling_method="none", acceptable_iter=Inf)
end

if false
  folder_name = "ipopt/medium_high_tol"
  if_mkdir("../results/$folder_name")
  run_cutest_problems_on_solver(problem_list, folder_name, print_level=3, tol=1e-8, max_iter=3000, max_cpu_time=TIME_LIMIT, bound_relax_factor = 0.0, nlp_scaling_method="none", acceptable_iter=Inf)
end

if false
  folder_name = "ipopt/medium_perturb_high_tol"
  if_mkdir("../results/$folder_name")
  run_cutest_problems_on_solver(problem_list, folder_name, print_level=3, tol=1e-8, max_iter=3000, max_cpu_time=TIME_LIMIT, nlp_scaling_method="none", acceptable_iter=Inf)
end

max_stars_repo_name: ohinder/OnePhase

max_issues_repo_name: ohinder/OnePhase

max_forks_repo_name: ohinder/OnePhase